pom_v16 by huanghengdong


									Geophysical Fluid Experiments with the
Princeton Ocean Model

Lie-Yauw Oey
Princeton University
                                                                                            Page 2

                               COVER ILLUSTRATION

A three-dimensional surface of near-inertial energy = 0.03 m2s-2 on Sep/03/12:00,
approximately one week after the passage of the disastrous hurricane Katrina, Aug/25-
30/2005, in the Gulf of Mexico, USA, simulated by the Princeton Ocean Model. This
shows penetration of intense energies to deep layers due to the presence of the warm-core
ring and the Loop Current, both represented by the dark contours in the three cut-away
xy-planes.    The location of an observational mooring where extensive model-
observational analyses have been conducted is shown as vertical dashed line [from: Oey
et al. 2008: “Stalled inertial currents in a cyclone,” in press, Geophysical Res. Lett., with
permission to reproduce].
                                                                              Page 3

PROLOGUE: About This Book

“.. I have no special talents. I am only passionately curious..” (Albert

       This book is about exploring various aspects of fluid motions on a
rotating earth using a popular, relatively simple yet powerful numerical
model - the Princeton Ocean Model (POM). It is a book aimed primarily
for advanced undergraduates and graduates in physical oceanography; but
the book should also be useful for researchers or anyone fascinated by and
intensely curious about oceanic fluid motions. Some knowledge of fluid
mechanics is assumed, equivalent to a two-semester course which usually
covers up to the derivations of conservation laws of viscous fluid motions,
including the boundary layer theory; these prerequisites should not be an
impediment to those interested enough to want to open this book. On the
other hand, we have strived to make the book more or less self-contained by
building each chapter from simple to more advanced concepts. By
conducting geophysical fluid experiments on a computer and analyzing the
results, we hope that the student will gain a solid understanding of
geophysical fluid dynamics (GFD) and physical oceanography; the student
will also learn to critically examine the results (be a skeptic) and to attempt
connecting them to the real-world phenomena. And of course, the student
will learn to use POM, numerical methods and, we hope, to keep exploring
long after he or she surpasses this book.
       Why POM? One reason is our familiarity with the model. More
importantly, however, it is because the model is relatively easy to tweak, say
to suit a different flow problem, or to change the model physics. This
makes it an excellent educational tool because it gives the student a hands-
on experience. It is much like the difference between being a driver and
being a driver as well as a mechanics. Most of us are the former, but only
the privileged few belong to the latter. POM shares many of the same
features of other popular ocean models now available in the scientific
community, and it uses basically the same equations of motions and
conservation of mass etc. Therefore, an accomplished student of this book
should find it easy to transit to other models should he or she so desire.
       The book begins in chapter 1 with the near-surface ocean response to
wind built upon Ekman’s fundamental work. This is not the common
                                                                             Page 4

approach of a book in physical oceanography. On the other hand, most of
us experience the ocean perhaps through a summer swim in a coastal sea or
in a lake, and wind-driven surface motions are the most apparent. For a
majority of flows on the rotating earth, we find that by directly dealing with
viscous boundary-layer (i.e. Ekman’s) flows, one can quickly grasp why
most of the ocean’s interior is nearly inviscid yet why the boundary layers
are so important for the interior flows.

POM and related files for exercises outlined in this book can be downloaded

Download latest public-released POM from:

Download POM User Guide other releases from:
                                                                                            Page 5

CHAPTER 1: Wind-Driven Ocean Currents
       The book begins with the near-surface ocean response to wind.                      We
will for the moment assume that the ocean density  ( 1035 kg m-3) is
constant. Throughout this book, unless explicitly stated, we will use the SI
units. Most of us experience the ocean perhaps through a summer swim in
a coastal sea or in a lake. Wind-driven surface motions in the form of
waves and surfs are familiar. However, these motions can involve large
vertical accelerations and are of very small scales (order of centimeters to
meters, O(cm-m)) which we will not model for now, at least not directly.
Instead, we focus on motions which are of sufficiently large scales in the
horizontal (scale L ~ O(1-1000 km)) compared to vertical (scale D ~ O(10-

1000 m)), D/L << 1, that the pressure field satisfies the hydrostatic equation:

                      ,                                                           (1-1)
where a Cartesian coordinate system xyz is used such that, conventionally, x
and y are the west-east and south-north axes respectively, and z is the
vertical axis with z = 0 placed at the mean sea level (MSL; Figure 1-1); x, y
and z are measured in m. The pressure p is in N m-2 (also called Pascal,
Pa), and g ( 9.806 m s-2) is the acceleration due to gravity.                    In other

words, the motions of interest occur in such thin fluid layers (with D/L <<
1) that as far as the vertical force balance is concerned the fluid may be
treated as being static.1 For such an atmosphere-ocean system, the pressure
at some depth z = z’ in the water is just the weight per unit area of fluid
(water and air) above, or
       p = [mass per area].g = [ (-z’) + <eazea> ].g                           (1-2)

1 The thinness of the oceanic (or atmospheric) layer relative to its horizontal extent is
comparable to that of the thickness of paper of this book relative to its width.
                                                                                          Page 6

Figure 1-1. A sketch of water of mean depth H(x,y) and free surface (x,y,t).

Here, z = (x,y,t) is the free surface;  and/or z’ can be positive or negative,

and z’ ( 0) represents the distance from the sea-surface to the point z = z’

in the water.       The term <eazea> represents the mass per area of the
atmosphere, it is given by
       <eazea> =     
                                                                             (1-3)

where a is the air density. Thus ea and zea may be thought of as the
effective density and height of the “air-containing” atmosphere. The main
bulk (about 80%) of the atmosphere’s mass is within the lowest layer (the
troposphere) of about 10 km thick such that the atmospheric pressure at the
sea surface, pa = <eazea>.g  105 N m-2 = 1 bar.2 Equation (1-2) can be
written as

2 An easy way to remember this is to use ea  1 kg m-3, zea  10 km and g  10 m s-2. Note
also the familiar weather reporting of “1000 millibar” etc, which is 1 bar, the approximate sea-
level pressure.
                                                                                Page 7

      p = g(-z’) + pa(x,y,t)                                          (1-4)

which is also obtained by integrating (1-1) from z = z’ to z = ; in general,
the pa is a function of the horizontal position and time.
      If both a and pa do not vary with (x, y) and either are steady or at
most only vary so slowly with time that the fluid remains in hydrostatic
equilibrium, i.e. (1-1) is satisfied, then the fluid remains motionless provided
that initially it is also at rest. A simple system to imagine is water in an
annulus channel (Figure 1-2). If the axis of the annulus is far from the
center, so that r/Ranu << 1, where r is the width of the channel, and Ranu is the
radius of curvature of the annulus, then one can approximate the
(motionless) water to be in a straight channel with x directed along the
channel axis, y across the channel (Figure 1-2) and a vertical slice along the
axis gives Figure 1-1. Modelers refer to this kind of channel a “periodic
channel” since, starting from a yz-section at any x-location and going around
along the axis, the field variables return to the same values. When
conducting geophysical fluid experiments, a periodic channel is a convenient
configuration to use because it often allows the extraction of essential flow
physics while at the same time alleviates the modeler from having to
formulate more complicated boundary conditions.
                                                                               Page 8

Figure 1-2. A circular annulus channel of radius Ranu much larger than the
width of the annulus r. This is used to illustrate along-channel periodic
fluid motion within the annulus.

  1-1: A Simple Shearing Flow by the Wind
       Consider therefore the straight channel (Figure 1-1) in which we will
further assume that all variables are independent of the cross-channel axis
“y,” and that the earth’s rotation is nil. The flow is described by the three
components of velocity (u, v, w) and the pressure p. Since there can be no
flow across the channel wall, i.e.
      v = 0,            at y = 0 and y = r,    t.                     (1-5)

the cross-channel velocity v must be zero everywhere; the symbol t means
“at all time.”   The channel is further idealized by letting its depth H =
constant. Initially, at t = 0, the water is at rest. A wind stress  (N m-2) is

then applied uniformly at the surface, so that  is a function of time only:

       = (t)                                                         (1-6)
Similarly, we could stipulate that pa is also a function of time only; but for
simplicity we will set
        pa = 0.                                                      (1-7)
It follows that, the along-channel velocity u cannot vary with x. Therefore,
at any point, there can be no accumulation (convergence or divergence) of
mass, and since w = 0 at z = H, the vertical velocity is nil everywhere, and
the sea-surface remains flat. Thus,
      u = u(z, t),   w = 0 = /x = p/x.                            (1-8)
       Under these very specialized conditions, a parcel of fluid experiences
acceleration only due to the vertical shear stress. The momentum balance
is then:
                                                                      (1-9)
                                                                               Page 9

where zx denotes shear stress (force per unit area) in the x-direction acting
on the fluid elemental face that is perpendicular to the z-axis, and D(.)/Dt is
the material derivative which for a fluid property S is given by
                                          .                           (1-10)
For the present specialized case, the last three terms are nil, and with S = u,
we have      =     . A loose analogy is a stack of poker cards which are
glued with a special adhesive that is never dry and thus remains sticky.
The top card is then “pulled” parallel to the card’s surface a small distance
and it drags upon the card below it. In our fluid system, the wind is doing
the “pulling” by transferring air momentum onto the water’s surface (and we
assume that no ripples or wind waves are produced!). Momentum is
vertically transferred from the surface to fluid layers below – upper fluid
drags upon the lower fluid which in turn drags upon the layer further below,
and so on (Figure 1-3).

Figure 1-3. The ocean is treated as a (vertical) stack of thin layers each of
thickness z : shown are layers “1” and “2;” layer “0” is air. The Fmn is
(shearing) force per unit (horizontal) area due to layer “m” on layer “n.”
The formulae show that the net force per unit area in layer “1” is F1net = F01
+ F21 =z.zx/z. Therefore, the net shearing force per unit mass is F1net x
                                                                                     Page 10

y/(xyz) = (zx/z)/. This, in the idealized case described in the text,

is = u/t, the acceleration in the x-direction.

For the so called Newtonian fluid such as water, the shear stress zx is
proportional to the vertical rate of change of fluid velocity (i.e. the vertical
“rate of strain”):
        zx =                                                               (1-11)

where      is the viscosity which to a good approximation may be taken as a

constant and moreover is rather small; it is  10-3 kg m-1 for water and
1.7×10-5 kg m-1 for air. Equation (1-11) is valid for laminar or slowly-
moving viscous flows void of turbulent movements we usually see in say, a
mountain stream or in swirling air vortices around a house. However, as
will be seen below, one often model oceanic and atmospheric flows using a
similar formula, except that an “eddy viscosity”         e   is used to represent the
aggregated effects of small eddies on the large-scale flows. The eddy
viscosity is much larger than the molecular viscosity and is moreover not a
constant – it depends on the flow.         The “ e” will be ascribed below but for
now we will proceed with our model formulation using equation (1-11).
Equation (1-9) then becomes:

where     e   =   e/,   the kinematic eddy viscosity. This is a simple “heat-
diffusion” equation that can be easily solved subject to some initial
conditions and boundary conditions at z = 0 and z =H; for example:
        u = 0,                    t=0                                       (1-13)
        u = 0,                    z = H                                    (1-14a)
                                                                                         Page 11

                        ,      z = 0,                                           (1-14b)

where o is a constant wind stress.
Exercise 1-1A: Use POM to solve (1-12) subject to (1-13) and (1-14) in a
periodic straight channel. Compare your solution with exact (analytical)
solution. Experiment with different grid sizes as well as with more
complicated initial and/or boundary conditions. Analyze and discuss your
results (e.g. use Matlab or IDL etc).
Exercise 1-1B: Repeat Exercise 1-1A using a                  e-profile   that decays with

depth, e.g.   e   =   eo   exp(z), where   eo   and  are constants.
Exercise 1-1C: Repeat Exercise 1-1A using an annulus channel (Figure 1-
2) instead of the straight (periodic) channel. Use a curvilinear grid (in
POM) for this exercise. Experiment with various values of Ranu and r.
Compare and discuss your results. Are the results substantially different
from the straight-channel case when Ranu is “not too large” (define what this
means), why? Is equation (1-12) still valid then? Why (or why not)?

1-2: Effects of Earth’s Rotation
       In the above straight-channel flow example, we saw that the cross-
channel flow is nil. This is no longer true if the channel is placed on a
turntable and rotates about a vertical axis. We now expand the study of
wind-driven shearing flow to the case when the earth’s rotation cannot be
ignored. We can imagine that the earth is the “turntable.” If our channel
is placed at the (celestial) pole, then it rotates about the earth’s rotation axis
once every day (0.997258 day to be more exact) anticlockwise at the north
pole and clockwise at the south pole (Figure 1-4). Thus the rotational
frequency  = 2/(86163.09 s) = 7.292  10-5 s-1. If the channel is placed
on the (celestial) equator, then there is no rotation (though the channel is
being carried west to east at a great speed = Rearth  1673 km/hour; Rearth 
6371 km is the radius of a sphere having the same volume as the earth).               At
a given latitude, the rotational frequency is:
      Rot() = .sin(), /2    +/2                                        (1-15)
                                                                             Page 12

with the convention that positive (negative) is anticlockwise (clockwise).

Figure 1-4.    The Earth's axial tilt (or obliquity) and its relation to the
rotation axis and plane of orbit (from Wikipedia).

       We are anchored (by gravity) to the earth and therefore also rotate as
the channel does. It is most convenient to describe phenomenon referenced
to this rotating coordinate frame rather than to a fixed, inertial frame. For
scalar quantities such as the temperature or salinity of a (generally moving)
fluid parcel, a change in the frame of reference clearly will not change their
value. The parcel’s momentum is a different story. While the phenomena
are still independent of the frame of reference, our perception of the
phenomena as well as their descriptions are altered. An example is offered
in Figure 1-5.
                                                                               Page 13

Figure 1-5.    A turntable rotating at a constant  radians/s with two
observers “A” and “B” on it. At time t o “A” throws the ball (red) to “B.”
A time-interval “t” later, the relative position of both observers are
unchanged (and neither perceive any change). Let the ball’s speed be such
that it arrives at the original B’s position in the same time-interval. To
“A,” the ball has been forced to the right of “B.” This apparent force is
called the Coriolis force. It is only perceived by “A” (and “B”) in the
rotating frame but not felt by an outside observer in an inertial frame.

      We now derive the momentum equation on a rotating frame. We
examine first the rate of change of an arbitrary vector (e.g. the momentum)
B = B1i1 + B2i2 + B3i3 in a rotating Cartesian frame of reference with unit
vectors (along x, y and z) i1, i2 and i3 [Pedlosky, 1979]. By chain rule:
where repeated index k means summation over all k = 1, 2 and 3, and the
subscript I denotes rates of change as seen from an observer in the inertial
frame, respectively.     Clearly, the first term on the RHS of (1-16) is
i.e. it is the rate of change of Bk when ik is fixed, i.e. when the observer is
himself rotating, hence is not able to sense the rotation of the coordinate
axes. The second term on the RHS of (1-16) involves              , the rate of
change of a vector of fixed length (i.e. a unit vector), hence its change is
completely due to its changing direction as the coordinate axes rotate.
From Figure 1-6, we see that
                                                                     (1-18)

where      is the angular velocity at which the constant-length vector A
rotates. Setting A = ik and using (1-17) and (1-18), equation (1-16) gives
                                                                     (1-19)
                                                                              Page 14

Figure 1-6. Rate of change of a vector A of fixed length caused by it being
rotated by . This rate of change is seen from the sketch to be a vector of

length |||A|sin(), where  is the angle between A and , directed

perpendicular to the plane formed by A and , i.e. it is simply the vector

product A.

Thus the rate of change of the vector B in the inertial (non-rotating) system
consists of two parts. The first part,        , is the rate of change as seen by
an observer (that is us) “glued” to the rotating frame of reference (that is the
earth). To this must be added the second part, B, in order that the

correct rate of change,       , can be used in Newton’s law to describe fluid

motion. This equates the rate of change            of momentum following a
fluid element in an inertial frame of reference to the sum of forces acting on
the element:
                                                                                   Page 15

                                                                        (1-20)
where the subscript I again reminds us that the equation is valid only in an
inertial frame of reference. The RHS is the sum of the pressure gradient
force  p, the body force            , and the viscous force F(uI) which for

Newtonian fluids with molecular viscosity         is

These equations are derived, for example, in Batchelor [1967]. Applying
(1-19) to the position vector r of the fluid element, and also to uI (i.e. setting
B = r and B = uI in turn), we obtain,
                                                                         (1-22a)

                                                                         (1-22b)

where in (1-22a) we have set                  and                . As seen in the
inertial frame (from a “fixed star” outside the planet earth), the fluid velocity
uI consists of the velocity uR measured by an earth-bound observer, plus an
amount equal to the rate r (in m/s) at which the fluid element is being
carried along by the rotating earth. Since uR is the velocity we directly
observe, we substitute (1-22a) into (1-22b), (1-21) and (1-20) to express the
equation of motion entirely in term of uR. Thus, (1-22b) becomes:
                                                          ,            (1-22b)
where we have assumed that            is small compared to other terms.    Figure

1-7 explains that (r) is a centrifugal acceleration that can be expressed

as the gradient of a centrifugal potential  c; it therefore can be grouped

with the body force on the RHS of (1-20), i.e. ( + c) = g, the acceleration
due to gravity. The direction of g is conveniently taken to align with the z-
                                                                                  Page 16

axis, positive upwards. Since |(r)| < ||2 Rearth  3.410-2 m s-2, the
centrifugal acceleration makes but a small correction to the radial direction
that is “vertical” were the earth not rotating.       Note that (r) is zero at
the poles and is maximum at the equator; it accounts for the slight bulge of
the earth around the equator. Substituting (1-22a) into (1-21), one can
show that F(uI) = F(uR). Thus substituting (1-22b) into (1-20) gives:
                                                                     (1-23)
which is entirely in terms of quantities in the rotating frame of reference.
As written, equation (1-23) is valid also for non-constant . Henceforth the
subscript R will be dropped.

Figure 1-7.     (a) Sketch showing the centrifugal acceleration vector
(r) = ||2|r|sin()nc       of a fluid material volume element with
                                                                                         Page 17

position r on the earth’s surface (at latitude = /2 - ), directed along the unit
vector nc perpendicular to and away from the axis of rotation of angular
momentum . For convenience, the origin of r is at the earth center, but

this is irrelevant. (b) Sketch showing the vector sum  ( + c), of (i)                   ,
the gravitational acceleration vector directed towards the center of the earth,
where       is the gravitational potential, and (ii)         c   =(r) =||2r [see

sketch (a) for meaning of r], where           c   is the centrifugal potential. Let r =

xk ik , then     c   = (||2x2k )/2 as direct differentiation using          = ik (/xk)

shows. Note also that           c   = (||2|r|2)/2 = |r|2/2. In geophysics, the

summed vector           ( + c) = g is called the acceleration due to gravity and
the direction of g (i.e. “upwards”) is conveniently aligned with the z-axis.

The Equation of Motion Applied to Large-Scale Oceanic Flows:
     Written in component form, (1-23) is:
                                                                            (1-24a)

                                                                            (1-24b)

                                                                           (1-24c)
In arriving at these equations, we have used the incompressibility condition
to eliminate the second term on the RHS of (1-21). Equation (1-25) is also
often called the continuity equation, a term that we will also use, though this
is correct only if effects of density change on mass balance are negligible.
Instead of           , we have also used the more conventional notation                  to
denote the substantive derivative:
                                                                                   Page 18

i.e. the rate of change of a field variable G following the fluid’s material
volume element. This is permitted since the vector r (Fig.1-7) denotes the
position of the fluid element. A localized coordinate has been chosen such
that the direction of 23 points in the z-direction and its magnitude is (c.f.

Figure 1-7) 2sin() = f, say, where  is the latitude of the fluid element.
The “f” is called the Coriolis parameter; it is zero at the equator, and is a
maximum (=2) at the pole as noted previously when discussing equation

(1-15). For large-scale, thin-fluid flows (i.e. D/L << 1) that we will mostly
focus on, it can be shown [e.g. Pedlosky, 1979; Gill, 1982] through a scaling
analysis that the         terms in (1-24) are generally quite small compared
to other terms in the same equation; these terms will therefore be dropped.
The molecular viscous terms       (    u) are also very, very small; though their
importance increases as the scale of motion becomes very small and
eventually they provide the necessary energy dissipation for the fluid
system. In models, effects of turbulence are usually parameterized by
invoking an eddy or effective viscosity which in the vertical (z) is denoted by
M (or diffusivity KH for heat and salt, chapter 3; two excellent texts are
Hinze, 1961, and Schlichting, 1963); these have much enhanced values than
their molecular counterparts: KM >> , a result of “eddy mixing” on the
small scales.3 Unlike the molecular viscosity, the eddy viscosity depends
on the flow itself, hence is a function of time and space, and in general takes
different forms in the vertical and horizontal directions (i.e. non-isotropic).
Various models of turbulence are available for KM (and KH). For example,
the Mellor and Yamada’s [1982] parameterization scheme is very popular
and is used in many oceanic and atmospheric circulation models including

3 A proper treatment involves Reynolds averaging – see Hinze, 1961. Equation (1-24)
is then called the Reynolds equation.
                                                                               Page 19

the Princeton Ocean Model. While parameterizations for the vertical eddy
viscosity and diffusivity are relatively well-known, those for enhanced
horizontal eddy viscosity and diffusivity are in comparison still a subject of
much debate. Unless otherwise stated, we will assume that horizontal
viscosity and diffusivity are small, and include them only for the purpose of
numerical computations: stabilization of the numerical schemes and/or
elimination of grid-point (“2x” oscillatory) noise [e.g. Roach, 1972].
Excellent reviews of turbulence mixing in geophysical flows can be found in
the various chapters of the book edited by Baumert et al. [2005].
      With the above simplifications, equation (1-24c) becomes
                                                                     (1-27)
which is the hydrostatic equation (1-1) except that it is seen to be valid also
for non-uniform .      Indeed, (1-24) admits a “static” solution

      u = 0,  = o(z) and p = po(z)                                  (1-28)
which however is not very interesting. More interesting flows are
produced as a result of deviations of the pressure field from this hydrostatic
       = o(z) + ’,          p = po(z) + p’                         (1-29)

where the ’ and p’ are perturbation density and pressure respectively. In

the ocean   o  oc = 1025 kg m-3 say, where oc is a constant density, and

’ is generally a small fraction of this: ’/o  0.01 and smaller.

Subtracting (1-28) from (1-24), letting   o in the acceleration terms, and

dividing through by o, we obtain:


                                                                             Page 20

                                                                     (1-30c)

with an error in (1-30a,b) that is at most O(’/o).

1-3: Wind-Driven Homogeneous Flows with Rotation
      Equations (1-25) and (1-30) are insufficient to solve for u, v, w, ’ and
p’. Additional equations involving the heat and salt equations as well as
the equation of state are required. These will be discussed in chapter 3.
Here, we study the case of a homogeneous fluid of constant density (i.e. ’

is known, ’ = 0), for which equations (1-25) and (1-30) are then complete
and may be solved for u, v, w, and p’, given appropriate initial and boundary
      For this homogeneous fluid system, while ’ = 0, p’ is not and it may

be written in terms of the free surface . Assume then a free-surface z =

(x, t) in an idealized ocean of depth z = H but with no lateral boundaries.

There can be no flow across the surface and bottom, i.e. z =  and z = H are

material surfaces, and D(z)/Dt and D(z+H)/Dt are both = 0:
                             
                                                                     (1-31a)

These equations constitute the upper and lower boundary conditions for our
problem. From (1-30c) with ’ = 0, we have p’ = p’(x, y, t); also the first of

(1-29) gives  = o, a constant. Integrating (1-27) and using the second of
(1-29) yield:
      p = ogz + function(x, y, t) = po(z) + p’(x, y, t),            (1-32a)
                                                                              Page 21

from which we may let
      po = ogz                                                     (1-32b)

Across the air-sea interface, z = , the total pressure p must be continuous
and equal to pa, so that (1-32) gives:
      p’ = og + pa.                                                (1-33)

The incompressibility condition (1-25) can now be used to relate  (or p’) to

the velocity. Integrating (1-25) from z = H to z =  and applying (1-31):
Equation (1-34) can also be re-written in term of a vertically-averaged
velocity                                 , where                ,   so   that
                        . The form used in (1-34) shows explicitly however
that the equation is valid for general (u, v) that may also be a function of z
(as well as of x, y and t).
       Substituting (1-33) into (1-30a,b), the resulting equations and (1-34)
are grouped together here for convenience:
                                             .                       (1-35c)

These constitute three equations to be solved for (u, v, ) given appropriate
initial and boundary conditions. These are in fact the equations solved
numerically in POM. After (u, v) are solved, the incompressibility
condition (1-25) is then used to solve for w. Before we describe the
numerical solution however, it is instructive to examine more closely how
one can approach the problem analytically.
                                                                             Page 22

1.3.1 Boundary-Layer and Interior Flows:
       For analytic treatment, we closely follow Pedlosky [1979]. We
study a problem in which the ocean is confined in the vertical by a free
surface z = (x,y,t) and a bottom z = H(x,y), but is unbounded in the
horizontal. It is more direct to use (1-25) and (1-30), which are repeated



Under each term, we indicate its order of magnitude by using the following
scale for each variable:
      (x,y)  L;    z  D;      (u,v)  U;     KM  K;                (1-37a)

      t  L/U;     w  W ~ UD/L;                                      (1-37b)

      p’  P ~ oLfoU; f  fo                                         (1-37c)

where “” means “have/has the scale of.”         Expression (1-37a) sets the
basic scales for the horizontal (L) and vertical (D) distances, and also the
horizontal velocity (U) and the eddy viscosity (K). The time scale (L/U) in
(1-37b) is chosen to be the “advective time scale,” i.e. the time taken for a
parcel of fluid to cover a distance of about L. If f is constant, then:
       f/fo = ±1                                                        (1-37d)
where the plus (minus) sign is for the northern (southern) hemisphere. The
vertical velocity scale W is chosen as follows. In order to satisfy the
continuity equation (1-36a), horizontal divergences or convergences
                                                                                Page 23

           must be balanced by vertical motion         such that their net sum
is zero; thus we set terms (1) and (2) in (1-36a) to be of the same order, and
W ~ UD/L. It follows then that the scale of term (2) of equations (1-36b,c),
WU/D ~ U2/L, which is also the scale of term (1) of these equations;
although |w| is in general much smaller than |u| or |v|, it multiplies the
vertical gradient       or      which is generally much larger than the

horizontal velocity gradients such as        . The scale for pressure in (1-
37c), P, is such that the horizontal pressure gradients balance the Coriolis
acceleration (see below). It is easy to see now that when the variables in
(1-36) are non-dimensionalized by the above scales, the equations become
[after dividing (1-36b,c) through by the scale (foU)]:



where now all variables are dimensionless, f = ±1 for constant Coriolis and
are the two parameters of the problem (the factor “2” in Ev is for
convenience, below).   The non-dimensionalized kinematic boundary
conditions at z =  and z = H have the same form as equations (1-31a,b).
      Equations (1-38b,c) show that the Coriolis and pressure gradient
terms are of the same order, both of these terms in (1-38b,c) are multiplied
by “1.” Thus by choosing the scales as those in (1-37c), we implicitly
assume that the Coriolis acceleration terms are important, and are in fact of
the same order of magnitude as the pressure-gradient terms. Clearly this
choice is only valid if “f” is not zero; the scaling (and balance of terms) will
need to be changed close to the equator. Observations of most oceanic
(and atmospheric) flows do indicate that away from the equator, pressure
gradients approximately balance Coriolis accelerations. From (1-38b,c) we
                                                                                 Page 24

see that the conditions for this to be true are that the Rossby and Ekman
numbers are small. The former requires that either the flow is not strong
and/or the horizontal scale L is sufficiently large – for example, U  0.1 m/s,

L > 10 km for fo  10-4 s-1. The latter requires that the Ekman scale

(2K/fo)1/2 is much smaller than the ocean depth D; if K  10-4~10-1 m2 s-1, this
requires that D > 10~50 m.
      However, close to a boundary, no matter how small the eddy viscosity
KM is, at some point the shear terms                          become important.
This is most apparent at the bottom boundary where the fluid parcel cannot
have a relative motion with respect to that boundary, i.e. the “no-slip”
condition must be satisfied:
        nt u = 0, or approximately u = 0,    at z = H(x,y)             (1-40)
where nt is the unit vector tangential to the bottom. Therefore, unless the
flow is trivially zero, the interior (i.e. away from the bottom) velocity of a
real-fluid (i.e. fluid with viscosity) flow must decrease in value until it
becomes zero at the bottom.4 Since Ev is small, the only way the viscous
terms                    can become large enough to balance the Coriolis
acceleration term is if the z-derivatives of (u, v) become very large, i.e. if the
decrease of velocity from the interior to the bottom occurs within a thin
layer in the immediate neighborhood of the bottom. This thin layer is
called a boundary layer, or a bottom boundary layer in the present case. A
similar thin layer also exists near the surface where wind stresses are

Geostrophic Nearly-Inviscid Interior Flows:
       Far away from boundaries, then, one expects that the flow is nearly
inviscid and one may approximate equation (1-38b,c) by dropping the term

4 An inviscid flow needs only to satisfy the no-normal flow condition
equation (1-31b), so that a fluid parcel can ‘slip’ past the bottom.
                                                                                   Page 25

involving Ev. Just exactly how far is “far” will be determined below. If
is also small, equations (1-38b,c) then give (for f = ±1):

Oceanographers call the horizontal velocity (ug, vg) that satisfy (1-41)
geostrophic velocity, and the equation geostrophic relation. In the northern
hemisphere, a higher pressure in the south (east) than in the north (west),
p’/y < 0 (p’/x > 0) would produce an eastward (northward) geostrophic
velocity ug > 0 (vg > 0). The geostrophic flow around a center of high
(low) pressure is therefore clockwise (anticlockwise) or anticyclonic
(cyclonic). The direction of the flow is reversed in the southern
hemisphere (f < 0). For homogeneous fluid flow, since p’ is independent of
z (see 1-38d), (ug, vg) is therefore also independent of z. From (1-41), a
geostrophic flow also has zero horizontal divergence:
                      ,                                                   (1-42)
so that the corresponding vertical velocity wg satisfies:
All three components of the geostrophic velocity (ug, vg, wg) are therefore
independent of z, the rotation axis. This is the Taylor-Proudman theorem
which is valid for inviscid homogeneous slow (small ) flow. If the upper
or lower boundary is flat, so that wg = 0 there, then wg = 0 through the water
column. The motion is then two dimensional, (ug, vg)  0, but wg = 0.
Figure 1-8 shows an example of this type of flow.      (See Exercises).

Strictly speaking, the model solution in Fig.1-8 is still evolving in time, so it
is not steady and does not exactly satisfy equations (1-41), (1-42) and (1-43).
This is an important point – purely geostrophic flows satisfying exactly (1-
41) through (1-43) are indeterminate. Any p’(x,y) satisfies the hydrostatic
relation (1-38d), and (ug, vg) can be constructed to satisfy the two
(approximate) momentum equations (1-41a,b); the unknown (ug, vg) in turn
                                                                                      Page 26

is horizontally non-divergent (equation 1-42) which simply requires that the
vertical velocity wg is independent of z. We know all these information
(hence the Taylor-Proudman theorem) but cannot pin down what the flow is!
To arrive at the unique model solution shown in Fig.1-8, we actually
included the small terms on the left hand side of momentum equations (1-
38b,c), the time-dependent terms in particular.5 Another way to pin down
a unique solution is to include the shear terms                              which of
course have the added benefit that we can then also satisfy the physically
relevant no-slip bottom condition (1-40) and the conditions at the sea surface
where windstress may be applied.

5 There were additionally horizontal viscous terms on the RHS of (1-38b,c) included for
numerical stability, but these were small and unimportant to the arguments presented
                                                                          Page 27

Figure 1-8. Nearly-steady homogeneous (constant-density) flow in an x-
periodic channel, 600km by 300km, of depth 200m except at the channel’s
center where a cylinder rises 50m above the bottom. Color is sea-surface
height in meters and vectors are velocity at (A) z=0m (i.e. surface), (B)
z=90m and (C) z=180m. This model calculation was carried out for 100
days when the flow has reached a nearly steady state. The experiment is
meant to illustrate Taylor-Proudman theorem: flow below the cylinder’s
height goes around the cylinder while above it the flow also tends to go
around as if the cylinder extends to the surface. The velocity does not vary
with “z” and the vertical velocity (which is not shown) is nearly zero.

Viscous Boundary-Layer Flows:
                                                                             Page 28

      It was mentioned before that no matter how small the Ekman number
Ev is, at some point close to the boundary the shear terms
become large enough that they balance the Coriolis and pressure gradient
terms which by our scaling (see equation 1-38) are of O(1); purely
geostrophic flows then break down. Thus in a thin layer close to the
boundary, the small Ev multiplies                 to make their product O(1).

We express this mathematically by changing the coordinate z to , where

                                                                   (1-44)
so that the shear terms become
                                          ,                         (1-45)
                                     

which is of O(1). The heuristic interpretation is, since                     is
large, of O(Ev-1), inside the boundary layer because a small change in z
results in a large change in (u, v), the change is made more gradual in the
stretched coordinate  (i.e. ‘one moves more slowly’ in  than in z) and one
can resolve or “see” the rapid velocity change so that the boundary condition
(1-40) (and a corresponding one at the surface, see below) can be satisfied.
The “z” measures distance changes outside the boundary layer, i.e. “z” is an
outer coordinate (see equation 1-37a, in which the dimensional z* (where
the asterisk denotes dimensional) is non-dimensionalized by the water depth
D >> boundary layer depth), and the geostrophic relation (1-41 through 1-
43) is an outer solution, though at the moment an incomplete one. Inside
the boundary layer, vertical distances are measured by    scaled by a much
smaller (localized) length scale l* which from (1-44) is
                       .                                            (1-46)

The  is called an inner coordinate, and the solution to be derived below is
called an inner solution. Indeed, equation (1-38) belongs to a class of
problems that are amendable to solution by the singular perturbation
                                                                                      Page 29

methods, the early development and usage of which is in aerodynamics [van
Dyke, 1964].
      With the new vertical scale l*, we have                       , which would
appear to be much larger than (and therefore cannot be balanced by) the
horizontal divergence term                in the continuity equation (1-38a).
This apparent dilemma is resolved by recognizing that in actuality the
vertical velocity w is very small. Indeed, inside the boundary layer, the
continuity equation requires that W*/l* ~ U/L, where W* is the scale of w
inside the boundary layer; thus
where we have used (1-46) and (1-37b; for W). Therefore, because of the
thinness of the boundary layer, hence a much smaller aspect ratio l*/L than
the D/L previously assumed in (1-37), the appropriate vertical velocity scale
W* is correspondingly much smaller, by O(Ev1/2), than the scale W that we
have assumed. In the absence of other processes, the reader should satisfy
himself or herself to show that w* is the scale of the vertical velocity that is
valid both inside and outside the boundary layer.
      If we now write the newly rescaled vertical velocity              as

            = w*/W* = w.W/W* = w/Ev1/2                                       (1-48)

using (1-47), the            terms on the RHS of the momentum equations (1-
38b,c) become:
and equation (1-38) becomes:

                                                                 

                                                                     

                                                                              Page 30

where (1-45) has been used in the viscous shear terms in (1-50b,c), f = ±1
for constant Coriolis and the hydrostatic relation is unchanged but for a
trivial change of the vertical coordinate. Because of (1-50d) (or 1-38d), the
thinness of the boundary layer means that the pressure-gradient force (per
unit mass), (p’/x, p’/y)/o, may be represented by their values just
outside the boundary layer in the geostrophic interior of the water column,
i.e., by f(vg, ug) from (1-41) (f = ±1).
     Now that equation (1-50) is properly scaled for flows inside the
boundary layer, we may simplify it (as we did with equation (1-38) when
deriving the geostrophic relation) by dropping the terms involving :


      0                                                              (1-51b)
                                       

      0                                                              (1-51c)
                                           
where (1-41) has been used in (1-51b,c) to replace the pressure-gradient
terms by the geostrophic velocity. Also, for constant Coriolis, f = ±1,
where the first (second) or top (bottom) sign of          or   applies to the
northern (southern) hemisphere with positive (negative) Coriolis.   Note that
since (ug, vg) is independent of z (or ), the (u, v) appearing in the -
derivatives in equations (1-51b,c) may be replaced by
      (uE, vE) = (u-ug, v-vg),                                     (1-52)
the Ekman velocity; this is in fact what we will do below.
      A wind stress is applied at the surface, and the boundary condition in
dimensional form is:
                                 , at z* = 0                         (1-54)

where the kinematic wind stress (Xs*, Ys*) = (o*/o) (Xs, Ys), and o* is the
scale of the wind stress in N m-2. In non-dimensional form, (1-54) is:
                                        ,       at   =0              (1-55a)
                                                                               Page 31

where s = (Xs, Ys) is the non-dimensionalized kinematic wind stress, and
                           .                                          (1-55b)
The boundary condition at the bottom is the no-slip condition:
      (u, v) = (0, 0),         at z = H                              (1-56)
Surface Boundary Layer with Constant Eddy Viscosity:
       We assume constant eddy viscosity KM (= 1, non-dimensional), and
also that the depth of the channel is sufficiently deep that the surface and
bottom boundary layers do not interact. Equations (1-51b,c) are solved by
first re-writing them in terms of A = uE + ivE, where i = (1). Thus {(1-

51c)  i (1-51b)} gives (for f = ±1):

                       .                                              (1-57)
We first consider the surface boundary layer, for which the solution that
satisfies (1-55a) and that vanishes as  ~  (so that u ~ ug and v ~ vg as one
moves into the geoestrophic region outside the boundary layer) is:
                                                                     (1-58)

In dimensional form, equation (1-59a,b) are:


where                                     . Equation (1-59) shows that over a
distance of about
      E* = l* =                                                      (1-61)
                                                                                 Page 32

 the velocity profiles (u, v) ~ (ug, vg), the geostrophic velocity in the ocean’s
interior outside the boundary layer.    The E* is called the Ekman depth; it
measures the thickness of the boundary layer which is therefore thicker for
larger K and/or smaller fo. At the surface,  = 0, the vector product (uE,

vE) (Xs, Ys) gives the angle w between the velocity (uE, vE) and wind-stress
vector (Xs, Ys) as:
       w = sin-1[cos(/4)] = /4;                                      (1-62)
the positive (negative) sign is for the northern (southern) hemisphere and it
means that the Ekman velocity is directed to the right (left) of the wind
stress, at an angle of 45o.
       To obtain w, we use the continuity equation (1-51a), and also (1-59):

which gives upon integration           and setting                 :

In dimensional form, using equation (1-48), we have:
                                                                      (1-64b)

where  = z*/(2K/fo)1/2.

       Outside the surface boundary layer,  ~ , the vertical velocity is:

                                                                       (1-65)
In terms of the original variables “w” and “z,” we note from equation (1-44)
that the  ~  means a finite but small z-distance below the surface, i.e.

        ~  is equivalent to z ~ 0 and Ev1/2 << 1,                   (1-66)
so that:
                                                                              Page 33

In dimensional form, this is:
                                                                    (1-67b)
Equation (1-67) is the main result from the analysis of surface boundary
layer: the effect of the wind is to induce a small vertical velocity which is
directly proportional to the wind stress curl:
For positive wind stress curl, there is upwelling (downwelling) in the
northern (southern) hemisphere. The vertical velocity is small, but as we
shall see in a simple example, it has a profound effect on the interior flow
field. Before we give this example, it is necessary to conduct a similar
analysis on the bottom boundary layer near z =H.

Bottom Boundary Layer with Constant Eddy Viscosity:
      The analysis for the bottom boundary layer follows much the same
way as for the surface layer. Instead of the wind-stress condition (1-55),
we now apply the no-slip condition (1-56). Let
      zb = z + H                                                     (1-69)
where the subscript “b” (here and below) is used to denote the bottom
boundary layer, to distinguish the solution from that obtained previously for
the surface boundary layer. The bottom boundary layer is near zb > 0.
Equation (1-56) then becomes:
      (u, v) = (0, 0),   at zb = 0                                   (1-70)
The boundary-layer coordinate (1-44) becomes:
                                                                    (1-71)

Equation (1-57) is the same, but a solution that decays as b ~ + is chosen:

      uEb + i vEb = (ug + i vg) exp[(1 i)b ]                      (1-72)
which also satisfies the no-slip condition (1-70). Thus,

                                                                                   Page 34

and the corresponding dimensional form simply has the u’s and v’s replaced
by corresponding variables with asterisks:

The vertical velocity is again obtained from the continuity equation (1-51a),
and also (1-73):

which gives upon integration            and setting                   :

In dimensional form, using equation (1-48), we have:


where b = (z*+H*)/(2K/fo)1/2.
       In terms of the original variables “w” and “z,” we note from equation
(1-71) that the b ~ + means a finite but small z-distance above the bottom,
       b ~ + is equivalent to zb ~ 0+, or z = H+ and Ev1/2 << 1,       (1-77)
so that (1-76a) and (1-48) give:
In dimensional form, this is:


Equation (1-78) is the main result from the analysis of bottom boundary
layer. Unlike the surface boundary layer, it is the interior geostrophic flow
that is driving the bottom boundary layer. In the northern (southern)
hemisphere, a positive curl in the geostrophic velocity produces upwelling.
                                                                                Page 35

A Simple Problem: Quasi-Geostrophic Interior Solution:
       Equipped with the above boundary-layer formulae, we are now ready
to revisit the interior flow and attempt to close the otherwise indeterminate
(geostrophic) solution (equations 1-41 through 1-43). The correctly scaled
equation that governs the interior flow is (1-38). We will find an
approximate solution which as will become clear is called quasi-geostrophic.
We already saw that a first approximation is the geostrophic relation, and
indicated that the solution in Figure 1-8 was obtained through the inclusion
of small terms on the left-hand-side of the momentum equations. These
small terms represent small imbalances to the purely geostrophic relation.
Equations (1-38b,c) suggest then that the “degree” of this imbalance, or
ageostrophy, is of O( ). It is reasonable to assume then that the solution
consists of predominantly the geostrophic part, plus a small correction that is
of O( ):

      u = ug + u1 + O( 2);         v = vg + v1 + O( 2); etc.           (1-79)

and similarly for other variables also.   Higher order correction terms , O( 2)
and higher are also indicated in equation (1-79), since there is no reason to
expect that the geostrophic relation plus the O( ) correction (i.e. u1 etc) will
provide the complete solution. It is important to understand that just as the
(ug, vg, …) are properly scaled of O(1), the (u1, v1, …), (u2, v2, …), etc are
also all scaled to be of O(1). Only then it makes sense to assume the form
of (1-79) for (u, v, …). Equation (1-79) is called a perturbation expansion,
which by the way may not be necessarily convergent [van Dyke, 1964].
From the discussion prior to Figure 1-8 (Taylor-Proudman Theorem) the
geostrophic velocity components are all independent of z: (ug,vg,wg)/z =
(0, 0, 0). In particular, for “w,” we already know from the boundary-layer
solutions at either the surface or the bottom boundary (see 1-67a and/or 1-
78a) that w ~ O(Ev1/2).         Therefore, since Ev1/2 << 1 ~ O( ) say, the

corresponding expansion (i.e. 1-79) for w = wg + w1 + … shows that
      wg = 0,      for all z;                                          (1-80)
                                                                               Page 36

otherwise, the solution for w will be inconsistent with its small values near
the surface and bottom boundaries. In quasi-geostrophic approximations,
the vertical velocity component is very small, O( ; Ev1/2), and represents
only the ageostrophic portion of the flow. However, this ageostrophy is
very important.
      Substituting (1-79) into (1-38b,c), and collecting and equating terms
that are of O(1) and O( ) separately we obtain for the O(1) terms the

geostrophic relation as before, while the O( ) terms give:





Note that the viscous terms                          and the vertical advection

terms                are missing in (1-81) because                   of course,

but in the latter case also because wg = 0 from (1-80). Thus the O( )
dynamics are also hydrostatic and incompressible, but are not in perfect
geostrophic balance. The ageostrophy is due to the unsteadiness and
advective terms by the geostrophic flow (i.e. the left-hand-side of 1-81b,c).
We eliminate the pressure gradients from (1-81) by subtracting the y-
derivative of (1-81b) from the x-derivative of (1-81c), i.e. by [(1-

81c)/x(1-81b)/y], and using (1-42) and also (1-81a):
                          
                                                                               Page 37

where g = vg/xug/y is the relative vorticity of the geostrophic flow (not

to be confused with the scaled vertical coordinate  of the boundary layer,

e.g. equations (1-44) or (1-71)). Since g is also independent of z, equation

(1-82) can be integrated          to yield, upon using (1-67a) and (1-78a):
                          
                                                                     (1-83)

Given the wind stress curl, s, this equation now allows the determination

of the geostrophic vorticity field g, hence also the geostrophic velocity (ug,
vg). It is the consideration of the viscous boundary layers near the surface
(providing the wind stress curl) and also near the bottom (providing the
‘damping’, g) that allows the existence of a small vertical velocity and the
determination of the O(1) geostrophic field.      A simple example further
illustrates these ideas.

A simple problem:
      Consider a homogeneous fluid in a channel of constant depth = H that
extends to     in the x-direction; all fields (and forcing) are assumed

independent of the y-direction. A y-directed kinematic wind stress s = (Xs,
Ys) is applied at the surface:
        (Xs, Ys) = (0, x)                                             (1-84a)
The solution is seeked only in the region L*  x*  +L*, where L* = 500 km,

or in terms of the non-dimensional “x,” in 1  x  +1. Taking the scale of

the kinematic wind stress o*/o = 10-4 m2 s-2:

                                             ;                        (1-84b)
                                                                             Page 38

we then have wind stresses of about    0.1 N m-2 at x* =   500 km. Since the
flow is independent of y, the geostrophic relation (1-41) requires that ug = 0.
Assuming also steady state, the LHS of equation (1-83) is then zero, yielding
the simple solution:
                                                                     (1-85a)
or in dimensional form:


The solution (equation 1-85) states a simple balance between the applied
wind stress at the surface and the frictional damping dues to the viscous
effects at the bottom. Knowing the geostrophic velocity in the interior, the
complete solution is supplemented near the surface by equations (1-59) and
(1-64a) (or in dimensional forms, equations 1-60 and 1-64b) and near the
bottom by equations (1-73) and (1-76a) (or in dimensional forms, equations
1-74 and 1-76b). Some details of this solution are given in Figures 1-9 and
1-10 (see exercises).
Exercise 1-3A: Use POM to solve the simple problem (above) of wind-
driven flow in a channel. (See Figures 1-9 and 1-10). Modify the code to
solve also the case with no-slip conditions at the bottom, and replot the
figures for this latter case.
Exercise 1-3B: Use POM to solve the Taylor-Proudman column (Figure
1-8). Note that in this case POM is modified without vertical viscosity.
Solve also the case of a seamount instead of the cylinder.
Exercise 1-3C: Use POM to solve a case that has wind as well as the
seamount: i.e. a combination of exercises 1-3A and 1-3B above. Derive
the corresponding analytical solution.
                                                                          Page 39

Figure 1-9. Nearly-steady homogeneous (constant-density) wind-driven
rotational flow in a channel of depth 200m and theoretically unbounded in x,
but in practice an x-length = 1000 km is used in the solution. The flow is
assumed independent of y. The wind stress is specified in y-direction only:
(0, 10-4).(x/500km  1) m2 s-2. Panel (A) is for the y-directed velocity and
panel (B) the x-directed velocity. In (B), profiles of u are also plotted at
the four indicated x-locations with scales 0.1 m/s given along the bottom
for the first and fourth locations. Detailed plots comparing the profile at
x=205 km (i.e. the first x-location) with analytical solution are given in
Figure 1-10. The Coriolis parameter is constant, fo = 610-5 s-1, and the

vertical eddy viscosity is also constant, K = 510-3 m2 s-1. This numerical
solution is at time = 100 days, and differs slightly from the analytical one
discussed in the text; instead of the no-slip condition at the bottom, a
                                                                        Page 40

matching of the velocity to the law-of-the-wall log layer is used [see the
POM manual by Mellor, 2004].
                                                                         Page 41

Figure 1-10. Profiles of (A) u, (B) v and (C) w at the x = 205 km location
described in Figure 1-9b. Solid line indicates analytical solution assuming
no-slip condition at the bottom, while dashed line the corresponding
numerical (i.e. Figure 1-9) solution assuming log-layer at the bottom.
Black line indicates the Ekman velocity while red (in panels (A) and (B))
indicates the total (Ekman + geostrophic) velocity. Note that in panel (B),
total v/4 is plotted (so the profile fits in the same panel).
                                                                                                Page 42

CHAPTER 2: Wind-Driven Homogeneous Ocean: Boundary Currents

        In this chapter we study the effects of wind on the general circulation of the
ocean. We will examine how the large-scale oceanic gyre may be explained by the
theories of Sverdrup [1947], Stommel [1948] and Munk [1950]. The circulation
consists of a slow equatorward drift in the mid-ocean and over the eastern portion of the
ocean basin, and a swift and narrow poleward current along the western boundary. We
then numerically simulate these large-scale features with simple model experiments using
POM. We will compare the numerical results against theories.

        We still assume that the fluid is homogeneous ( = constant), and that the motions
have horizontal scales much larger than the vertical scales; in other words, we still
assume that the hydrostatic equation (1-30c) is satisfied. In fact the scales are so large
(~1000 km’s) that the Coriolis parameter “f” is no longer a constant, i.e. equation (1-37d)
is not valid, and

        f = 2.sin(), /2    +/2                                                  (2-1)

can be a function of the latitude 

2-1: The -plane & the QG-Potential Vorticity Equation

        For the moment, unprimed symbols are dimensional variables; primed symbols
(such as u’ etc) denote dimensionless variables. Suppose we focus on a certain latitude,
 = o where f = fo, and let “y” be the northward (distance) coordinate measured from this
latitude (so that y = yo at  = o), then for “small” distances away from this latitude, we

        f  fo + [df/dy][yyo] = fo + [(df/d)/ro][yyo] = fo + [2.cos(o)/ro][yyo]

        f  fo + o(yyo)                                                               (2-3a)


        o = 2.cos(o)/ro                                                              (2-3b)
                                                                                               Page 43

Equation (2-3) constitutes the so called “-plane,” and is valid for |yyo| not too large (i.e.
a few thousands of kilometers, for o  2×10-11 m-1s-1 near mid-latitude.) We now
derive a vorticity equation similar to equation (1-82) but taking into account that f is a
function of y, i.e. is given by (2-3). In terms of dimensionless variables, denoted by
primes, (2-3) is:

       f’ f’o + ’(y’y’o),    ’ = (oL/|fo|) << 1.                                  (2-4)

where f’o = fo/|fo| = ±1 denotes northern or southern hemisphere. Note that ’ is small;
we will in fact let ’  , and from (1-39):

        = ’/ = oL2/U  O(1).                                                        (2-5)

With (2-4), the Coriolis terms in equation (1-38) (which is dimensionless) becomes:

          f'v’ = [f’o + ’(y’y’o)][v’g + v’1 + O( 2)]
               = f’ov’g + f’ov’1 + ’(y’y’o) v’g + O( ’,       2
                                                                     )                 (2-6a)

         f'u’ = f’ou’g  f’ou’1  ’(y’y’o) u’g + O( ’,              2
                                                                             )]        (2-6b)

With these expressions, substituting (1-79) into (1-38b,c), and collecting and equating
terms that are of O(1) and O( , ) separately we obtain for the O(1) terms the geostrophic
relation as before, while the O( , ) terms give (c.f. equation 1-81):


                                                                                   (2-7b)

                                                                                   (2-7c)


where fo = ±1,  = ’/ , and dimensionless variables are now unprimed (except for the
perturbation pressure p’1). Taking the curl, i.e. [(2-7c)/x(2-7b)/y] now gives
(instead of 1-82):
                                                                                            Page 44

                                                                                 (2-8a)



Since (ug, vg) hence g are independent of z, equation (2-8) can be integrated              to

yield, upon using (1-67a) and (1-78a) (c.f. 1-83):

                                                                                (2-9)

In dimensional form [with asterisks *; using (1-39) and (1-55b)], this is:

                                                                                (2-10)

where           is the z-component of the curl of the kinematic wind stress               (see
equation 1-55a) and

                     ;                                                          (2-11)

is the bottom friction coefficient; l* is the Ekman depth (from 1-46); note that an asterisk
has been placed on D, i.e. D* to avoid confusion with non-dimensional variables.

2-2: A Model of Large-Scale Ocean Gyre & Western Boundary Current

         In this section, unless otherwise stated, all variables will be dimensional.
Consider a “box” ocean with sides Lx and Ly on a -plane (Figure 2-1). A zonal wind
stress is applied:

                                                                                 (2-12)

so that
                                                                                         Page 45

                                                                             (2-13)

where wo (which will be set = 10-4 m2s-2 below) is the magnitude of the (kinematic) wind
stress. Initially, the fluid is at rest:

        (u, v, w) = (0, 0, 0), and also, free surface  = 0                    (2-14)

At the four side walls, both the normal and tangential components of the horizontal
velocities are zero:

        (u, v) = (0, 0) at x = 0, x = Lx and also at y = 0, y = Ly.            (2-15)

At the ocean bottom and at the free surface, equations (1-31) are satisfied:

                              
                                                                              (2-16a)


The continuity equation and the x and y momentum equations are then (c.f. 1-35):



                                              .                                (2-17c)



is the depth-integrated transport vector (per unit length). These equations differ from
(1-35) in that the nonlinear advective terms in momentum equations (2-17b,c) have been
                                                                                         Page 46

Figure 2-1. A “box” ocean circulation driven by eastward (northern half of basin) and
westward (southern half) wind stress distribution.

       Depth-integrate (2-17b,c), and again drop nonlinear terms such as u etc:



where w and b are wind and bottom stress vectors respectively. These equations
together with equation (2-17a) constitute there equation for three unknowns (, U, V).
For our purpose, it is more convenient to write them in terms of the depth-averaged

                                                                              (2-20)

so that (2-17a) and (2-19) become:
                                                                                         Page 47




where r is the bottom friction coefficient with dimension m/s, we have modeled the
bottom friction as:

       (                     ,                                                 (2-22)

and we have also added a source term Q to the right hand side of (2-21a). This we will
use later to model the net effect of precipitation ( ) minus evaporation ( ), i.e.

       Q=(         –   )/o (unit is m/s).                                     (2-23)

       Let H = constant, a vorticity equation similar to (2-10) may now be derived from
(2-21). Taking the curl of (2-21b,c), we obtain:

                        
                                                                             (2-24a)

Apart from the terms involving Q and /t, this equation is equivalent to (2-10). The
term /t does not appear in (2-10) because it is neglected when we derived the vertical
velocity at the surface (see equation 1-67). We will explain in the next section when
this term may be deleted; assuming that it can be, then (2-24a) becomes:

                                                                             (2-24b)

The appearance of -Qf together with curlz(w) means that as far as the large-scale ocean
circulation is concerned, effects of wind stress curl and net precipitation minus
evaporation are equivalent. The reason is that in both cases, the ocean layer beneath the
surface boundary layer (where wind and precipitation etc act) ‘sees’ only vertical mass
flux coming in (or out) of the surface layer. In the case of wind stress curl, this is
“Ekman pumping” (Chapter 1) caused by convergence or divergence in the surface layer,
                                                                                            Page 48

while in the case of precipitation etc, it is more direct – due to the flux of water from the
atmosphere. We will model these effects and give a physical interpretation below.
Also, while H is the ocean depth, it may be (and for our purpose will be) interpreted as
the depth of the main thermocline in which the predominant large-scale currents reside.
In that case, the friction      is then considered as the friction at the base of the main

Sverdrup Relation:

        In the open ocean away from side boundaries, the relative vorticity  is assumed
weak, in fact   , and the main vorticity balance in steady state is, from (2-24b)
(omitting Q for simplicity):

                                                                                 (2-25)

Figure 2-2. A schematic to explain the equatorward Sverdrup transport V caused by a
negative wind stress curl, curlz(w) < 0.
                                                                                             Page 49

Equation (2-25) is the Sverdrup relation; it shows that in steady motion, a negative
curlz(w) produces southward transport. Figure 2-2 explains the physics involved. The
negative curlz(w) produces convergent flow in the thin surface (Ekman) layer (Chapter 1)
and hence a downward pumping that squashes the vortex tube in the deeper interior; this
is equivalent to a divergent subsurface flow. Squashing increases the cross-sectional
area of the tube, and if it stays at the same latitude a negative relative vorticity (i.e. a
reduced spin) would be produced,           . To keep the same   0, and assuming that
the tube remains approximately vertical (and it does), the tube must move south along the
curved earth surface where it can stretch and retain its initial zero spin. The curved
surface in this case is the physical equivalence of the -effect, i.e. the Coriolis parameter f
varies with latitude. Gill [1982; p.465] gives another, equivalent explanation. The
divergent subsurface flow beneath the negative curlz(w) increases the cross-sectional
area of the vortex tube, so the absolute vorticity (i.e. f + ) decreases. But since f >> ,
the only way that this decrease can be accomplished is that the tube moves southward
where f is smaller.

         Figure 2-2 illustrates two other features. One is that the effect of (positive) Q [=
( – )/o] is similar to (negative) curlz(w): this is because Q ‘pumps’ water into the
subsurface layer and again squashes vortex tube. The other one is that while the surface
is slightly elevated ( > 0) it is not important in the argument given above; in steady
state, it is in fact = 0.

       If as shown in Figure 2-1 the Sverdrup transport occurs over most portion of the
ocean basin, then the total southward transport is:

       SV                             (m3 s-1).                                  (2-26)
                                                                                         Page 50

Figure 2-3. Annual mean distribution of Ekman pumping, positive = upwelling, in unit
of 10-6 m s-1 (or  0.1 m/day). From Tomczak and Godfrey.

Figure 2-3 shows a global map of Ekman pumping – note large areas of downwelling in
the subtropical regions, suggesting equatorward Sverdrup transports there. From (2-26),
we have roughly SV  curlz(w).Lx/o  2      Lx/(Lyo)  (2×10-4/10-11)(Lx/Ly) m3 s-1 
20 Sv (1 Sv = 106 m3 s-1) for Lx  Ly,          10-4 m2 s-2 and o  10-11 m-1s-1. The
corresponding Ekman pumping (downwelling) is from (1-64b) = curlz(w)/fo 
2     /(Lyfo)  1~2×10-6 m s-1 for Ly  2000 km and fo  6×10-5 s-1, roughly equal to the
values shown in Figure 2-3. This Ekman velocity (or Figure 2-3) may be compared
with contours of Q [= ( – )/o] shown in Figure 2-4, which gives local maxima of
roughly 10-7 m s-1, about 10 times weaker than the Ekman pumping.
                                                                                        Page 51

Figure 2-4. Annual mean distribution of precipitation minus evaporation, Q = ( – )/o,
in m yr-1 ( 3×10-8 m s-1); shaded regions are positive.

Stommel’s Model of the Western Boundary Current (WBC):

        In our idealized closed basin (Figure 2-1), the southward Sverdrup transport must
be returned northward. Since this return cannot be in the mid-ocean (i.e. open ocean)
because of the Sverdrup constraint (as seen above), it must occur near the western and/or
eastern boundaries where the Sverdrup balance is broken. In other words, other terms in
the potential vorticity equation (2-24b), such as friction and/or nonlinearity (which is
neglected in 2-24b), become important near the boundary. Stommel assumes that
(bottom) friction is important so that (2-24b) becomes, in steady-state:

                        .                                                      (2-27)

Including curlz(w) adds unnecessary complexity without adding physical insight so it is
omitted. Also, near the western or eastern boundary,         . The solution is:

                                                                              (2-28)

where C(y) is an arbitrary function of y to be determined. Equation (2-28) shows a
northward jet which is maximum at the boundary (for an eastern boundary, the “x” is
                                                                                             Page 52

replaced by “xLx”) and either decays (if western boundary) or grow (if eastern
boundary) exponentially away from the boundary. Therefore, the northward return flow
must be along the western boundary and crucially depends on a non-zero o; the width of
this return jet may be estimated as the e-1 x-scale:

       xwbc = r/(Ho)                                                               (2-29)

This is  50 km for r  10-4 m s-1, H  200 m and o  10-11 m-1 s-1. The C(y) is
determined by requiring that the total northward transport carried by the WBC is equal to
the southward Sverdrup transport at the same latitude (i.e. same “y”). Thus,

                              , i.e.

                                                                                  (2-30)

With the aforementioned values of SV  20 Sv, o  10-11 m-1 s-1 and r  10-4 m s-1, C 
2 m s-1.

        The Stommel WBC solution (2-28) therefore consists of a jet concentrated against
the western boundary. Since this jet has to return northward all the Sverdrup transport
across a narrow width xwbc, the jet’s velocity is large. The Sverdrup velocity scale is,
from (2-25),  curlz(w)/(Ho)  2         /(HLyo)  (2×10-4/10-11)/(HLy) m s-1  0.1 m s-1
for H = 200 m and Ly = 2000 km; while the Stommel velocity scale is C  2 m s-1.
Stommel’s jet is maximum at the coast x = 0 and decays exponentially offshore to
become very weak in the ocean’s interior where the southward Sverdrup transport
prevails. Note that the Stommel’s jet, being wholly positive, cannot join smoothly to
the Sverdrup flow, nor is it expected to since it is based on a frictional equation that
ignores the Ekman pumping by the wind; but this is not important. Nor is it important
that the Stommel’s jet cannot satisfy the no-slip condition at the coast. The important
thing is that: (1) beta (o) is crucial for the existence of the western-intensified boundary
current, and (2) the WBC jet provides a mechanism that closes the ocean-basin mass

Munk’s Model of the Western Boundary Current (WBC):
                                                                                            Page 53

       A slightly more complicated model can be constructed to satisfy the no-lip
condition at x = 0. Instead of the bottom friction we include horizontal viscosity terms
on the RHS of the momentum equations:




where AH is the horizontal (eddy) viscosity coefficient (unit is m2 s-1), assumed a constant.
If we now take the curl etc as in the derivation of (2-24b) we obtain:

                                                                                (2-32)

Define a stream function         (unit m2 s-1) such that:


and therefore

                                                                                 (2-34)

In the WBC, we now let the beta term balances the horizontal viscosity term, equation (2-
32) in steady state is then (c.f. 2-27):

                                                                                 (2-35a)


            = (o/AH)1/3x.                                                        (2-35b)
                                                                                                Page 54

Equation (2-35a) is a fourth-order ordinary differential equation (with “y” as the ‘floating’
parameter not directly involved), and a solution of the form exp(m ) maybe assumed.
Four roots are:

m1 = 0, m2 = 1, m3 = exp(i/3)=(1+i3)/2, m4 = exp(i2/3) =(1i3)/2               (2-36)

The solution then is of the form

           = C1 + C2exp(m2 ) + C3exp(m3 ) + C4exp(m4 )                                (2-37)

where the C’s are arbitrary function of “y” yet to be determined. The m2 root must be
deleted - by choosing C2 = 0, because it leads to a solution which exponentially grows
with (or x). The root m1 leads to C1 – which is chosen to be C1 = I(y), where the
subscript “I” denotes “interior” (of the ocean away from the boundary) and from the
interior Sverdrup solution (2-25) with             ,

                                                                                   (2-38)

is the stream function of the interior’s Sverdrup flow; o(y) is an arbitrary function.
This choice is so that ~ I for ~ , since m3 and m4 give exponentially decaying
solutions for large .

Substituting m3 and m4 from (2-36) to (2-37), we obtain then:

           = I +      [C3           + C4            ]                                (2-39)

the C3 and C4 here differ from those in (2-37) but the detail is not important; they are still
arbitrary functions of “y.” The corresponding velocity components from (2-33) are:



The conditions that there can be no normal flow (         ) and no slip (      ) at    = 0 give:
                                                                                              Page 55

       C3 = I(0,y) + K, where K = constant, and                                   (2-41a)

                    .                                                               (2-41b)

Therefore, (2-39) gives:

         = I(x, y)  I(0,y)           [              +                ]

                    +K          [              +               ]                    (2-42)

From equation (2-35b), we note that (AH/o)1/3 has the unit of length (m). Its value is 
100 km (for AH = 103 m2 s-1 and o  10 -11 m-1 s-1) which is thin compared to the width of
the basin  2000 km (say); it is the e-folding decay scale of the western boundary current
near x  0. We are therefore making only a small error if I(0,y) in (2-42) is replaced by
I(x,y). Another way of saying the same thing is that, since I varies over a large scale
>> (AH/o)1/3, I(x, y)  I(0,y) near the western boundary. To determine K, we note
that, since our basin is closed, the stream function around the boundary is zero.
Therefore, for to be also = 0 at x = 0, (2-42) shows that K must also = 0. Therefore,
(2-42) becomes:

         = I(x, y) {1             [              +               ]}               (2-43)

The northward velocity is  /x:

                                                         ,       = (o/AH)1/3x,   (2-44)

in which we have neglected the interior velocity contribution I/x, which is much
smaller. Note that the jet form is again apparent by the exponential decay term. But in
contrast to the Stommel’s solution equation (2-28), because of the sine function, the
northward velocity is zero at x = 0. Also because of the sine function the velocity
becomes slightly negative (i.e. reversed flow) at the offshore edge of the jet.

Class: estimate the magnitude of            from (2-44).
                                                                                            Page 56

Dear Jane, Eda and Alu, here (section 2-3) is a derivation of the quasi-geostrophic
potential vorticity equation used in Sheremet’s paper (his eqn.1). He used a so called
reduced-gravity form of the equation. Pls. read and work out the steps as much as you
can, and hopefully I can help you on Thursday. Have a good weekend. leo

2-3: The Equivalent Reduced-Gravity Model

        In this section non-dimensional variables will be denoted by subscripts “n.” The
above model has one layer and involves the free surface . We now show that with
slight reinterpretations of the variables, the equations of that model are exactly the same
as those governing the motion of a two-layer ocean in which one of the layers is infinitely
thick and is quiescent. Before we do this, we find out first the condition when the term
involving /t in equation (2-24a) may be neglected. This is necessary in order to
understand the restrictions on the spatial and temporal scales under which the resulting
equation(s) apply. Instead of (2-24a), we will non-dimensionalize the vorticity equation
that uses the horizontal viscosity instead of the bottom friction, i.e. we use (2-32) with
f(/t)/H retained:

                        
                                                                                 (2-45)

Define the following non-dimensional variables:

                                                    

          w = wno,   Q = Qn Qo,   AH = AHn AHo.                                 (2-46)

Note that the scale for , i.e. fUL/g, is determined from balancing the pressure gradient
and Coriolis terms in equation (2-31b,c). Substitute (2-46) into (2-45):

                               
                                                                              


where R = (gH)1/2/f is the Rossby radius of deformation. We see that dropping the
term involving n/tn amounts to saying that R >> L, i.e. the spatial scale of motion is
much smaller than the Rossby radius. For H  1000 m and f  6×10-5 s-1, R  1500 km,
                                                                                         Page 57

so that motions with scales smaller than about 1000 km will not be much affected by the
production of vorticity due to stretching associated with the motion of the free-surface.
This kind of motion is called a barotropic motion. On the other hand, as we will now
show, R can be quite small for the reduced-gravity model, in which case the stretching
term becomes important.

The Reduced-Gravity Model from the Two-Layer Model:

       We now derive the reduced-gravity (or 1.5-layer) equations from a two-layer
model as sketched schematically in Figure 2-5:

Figure 2-5. An illustration in the xz-plane of an ocean model with two (inmiscible fluid)
layers with subscripts “1” and “2” denoting the upper and lower layers respectively.
Symbols are: 1 and 2 = constant densities of the upper and lower layers respectively, 2
> 1; H1 = constant undisturbed upper-layer depth; H2(x,y) = undisturbed lower-layer
depth which may be a function of (x,y) for variable bottom topography; h(x,y,t) = layer
depth; u(x,y,t) = the x-component velocity; v(x,y,t) = the y-component velocity; p(x,y,z,t)
= pressure; 1(x,y,t) = disturbed free-surface height measured from z = 0; 2(x,y,t) =
disturbed inter-layer height measured from z = -H1. Note that 1 and 2 can be positive
(upward) or negative (downward), and that h1 = H1 + 1  2.
                                                                                            Page 58

The caption of Figure 2-5 explains all the variables to be used. Since in each layer the
density is constant, the continuity (conservation of volume) equation takes exactly the
same form as that for a homogeneous fluid (see equation 2-31a):

                              , i = 1, 2.                                         (2-48)

where i = 1 is for layer-1 and i = 2 is for layer-2, and for simplicity the source/sink term
Q is omitted. (The subscripts ‘i’ in equation 2-48 are not to be confused with x and y
components, and repeated subscripts do not mean summation!). Except for the pressure
gradient terms, the momentum equations (2-31b,c) will also remain essentially the same.
To evaluate the pressure gradients, we will assume as before that the horizontal scales are
much larger than the vertical scales, so that the pressure is hydrostatic (see equations 1-1
and 1-2). Thus the pressure at some depth z = z’ in the water is just the weight per unit
area of fluid above. For a point in layer 1, the pressure p1(x,y,z,t) is:

       p1(x,y,z,t) = pa(x,y,t) + g1(1  z), for 1  z  H1+2,                (2-49a)

where pa = atmospheric pressure (i.e. = weight per unit area of air above z = ).
Similarly, for a point in layer 2, the pressure p2(x,y,z,t) is:

       p2(x,y,z,t) = p1(at z =H1+2) + g2(H1 + 2  z),
                                           for H1+2  z  (H1+H2),             (2-49b)

Then the horizontal pressure gradients are (for example /x):

       p1/x = pa/x + g11/x, for 1  z  H1+2,                          (2-50a)

       p2/x = p1/x + g(1h1)/x, for H1+2  z  (H1+H2),               (2-50b)

where  = (2  1), and we have used (2-49a) to evaluate p1(at z =H1+2), and also 2
= H1 + 1  h1 in obtaining the last term in (2-50b) (see Figure 2-5).

The momentum equations (2-31b,c) applied to layers 1 and 2 separately lead to (2-51b,c)
and (2-52b,c) below:

                                                                                            Page 59




                                               

                                               

Note that we have generalized /t to D/Dt in the momentum equations. These are the
equations for our two-layer ocean. More complicated forms with additional terms exist
in the literature; for our purpose, the above suffices. Note that the lower layer
experiences both the layer-1 pressure gradient (p1/x)/o as well as the additional force
(per unit mass) (g/o)(1h1)/x; and similarly for the y-component.

In the reduced-gravity model, we let one of the layers to be infinitely thick. Supposing
we let layer-2 to be infinitely thick, then from the continuity equation (2-52a) the deep-
layer currents (u2, v2) must approach zero, and h2/t ~ 0 also. Subtracting (2-52) from
(2-51) then gives:




in which we have put back the “Q” on the RHS of the continuity equation, and have
made a further approximation that 1x << h1/x and similarly 1y << h1/y,
which means that the variations of the free surface is negligible in comparison to the
variations of the upper-layer depth. Equations (2-53) are identical to the equations that
govern the motion of a single-layer (homogeneous) fluid of constant depth H1 (i.e. the
topography) with a free surface except that the acceleration due to gravity g is reduced to:

       g’ = g/o                                                                (2-54)
                                                                                             Page 60

which appropriately is called the “reduced gravity.”      To see the equivalence, set

       h1 = H1 + ,                                                                 (2-55)

where now the “” is interpreted as the departure of the upper-layer depth from its
‘mean’ of H1. For typical /o  10-4~10-2, the corresponding shallow-water
baroclinic phase speed Ci = (g’H1)1/2, hence also the baroclinic Rossby radius of
deformation Ri = (g’H)1/2/f, is some 10 to 100 times smaller. Therefore, for reduced-
gravity (baroclinic) motions, the time-dependent stretching term (L/Ri)2 h1n/tn that
produces vorticity in equation (2-47) can be significant. In fact, this term overwhelms
the other time-dependent term 1n/t in the limit of

       (L/Ri)2 >> 1.                                                                (2-56)

The reduced-gravity equivalence of (2-45) is:

                                                                                 (2-57)

in which the nonlinear advection term u1. 1 (from the D/Dt terms in 2-53b,c) is
included. But, consistent with the quesi-geostrophic approximations of Chapter 1, we
may replace the u1 in this nonlinear term by its geostrophic values:

       H1v1 = +(Ci2/f)h1/x,     H1u1 = (Ci2/f)h1/y                             (2-58)

Equation (2-58) suggests a stream function 1 such that:

       1/x = H1v1 and 1/y = H1u1                                             (2-59)

By comparing (2-58) and (2-59), we then have:

       1 = (Ci2/f)h1 = fRi2 h1                                                     (2-60)

Equation (2-57) then becomes:
                                                                                         Page 61

                                                   
                                                                             


where J(A,B) = A/x.B/yB/x.A/y is the Jacobian (of the functions A and B).

Eda, Jane and Alu, this equation (2-61) is identical to Sheremet’s equation (1), except
that his  is actually my 1/H1, and he sets the wind stress curl and Q to zero.

2-4: An Application of the Reduced-Gravity Model: Jet & Eddy Interaction

       Figure 2.6a shows an eddy, either cyclonic or anticyclonic, drifting westward and
approaching a straight meridional jet. The jet’s velocity is a maximum along its center,
and exponentially decays symmetrically on either side. Such a velocity profile is
obtained by specifying a constant potential vorticity (PV) jump across the jet. The PV is
defined as an anomaly from a rest state. For a reduced gravity model  an upper layer
of depth h = H +  (where H = undisturbed layer depth and  = perturbed layer height)
and velocity (u, v) overlying an infinitely deep resting lower layer  the PV is:

       PV = (f + )/h  fo/H = (  fo/H + y)/(H + ),      for a -plane,   (2-62a)


       PV = (  fo/H)/(H + ),            for an f-plane.                    (2-62b)

Then the PV-distribution of the jet is:

       PVJ     = +QJ/2,       x > 0,
               = QJ/2,       x < 0.                                          (2-63)

Note PV-jump|Jet = QJ. Similarly, for the eddy, assumed to be of a circular shape with
radius REddy, we have,

       PVEddy    = QE,        r < REddy,
                 = 0,          r > REddy,                                      (2-64)
                                                                                       Page 62

and the PV-jump|Eddy = QE.    Here, r denotes the radial distance from the center of the

Inverting the PVJ and PVEddy to Obtain the Corresponding Velocity Fields:

Figure 2-6a. A sketch of the jet with PV-jump|Jet = QJ, and the approaching eddy with
PV-jump|Eddy = QE.

Figure 2-6b. A sketch of a weak (cyclonic) eddy and/or a strong jet, QJ > QE. The
eddy does not cross the jet and is being advected along with the jet.

       For a weak eddy and/or a strong jet, such that QJ > QE (Figure 2.6b), the eddy
introduces only a small perturbation on the jet. The eddy does not cross the jet and is
being advected along with the jet [Stern and Flierl, 1987; Bell, 1990].
                                                                                              Page 63

Figure 2-6c. A sketch of the evolution at four indicated times for the combined system of
a strong (cyclonic) eddy and/or a weak jet, QE > QJ.

If the eddy’s strength is increased and/or the jet’s strength is weakened, QE > QJ,
where  > 1, the eddy interacts nonlinearly with the jet. In other words, the eddy’s
velocity field can significantly influence the jet, and the jet’s velocity field in turn affects
the eddy. The evolution of the combined system of a cyclone and jet is illustrated in
Figure 2.6c for four time instances, t1 < t2 < t3 < t4. Time t1 is when the jet begins to
“feel” the presence of the eddy “E,” when the latter comes within 2Ro of the jet’s axis,
where Ro = (g’H)1/2/fo is the Rossby radius. The sketch at time = t2 shows that the
eddy’s velocity field downstream (upstream) pushes the jet westward (eastward) to
develop the downstream (upstream) meander “DM” (“UM”). Vorticity inside “DM”
(“UM”) is anticyclonic (cyclonic). The “dipole” formed by DM and E moves
southwestward due to their mutually-induced velocity shown by the green vector, while
UM and E also have a mutually-induced velocity shown by the blue vector. The
resultant vector is westward indicated by the red vector in a direction that tends to cross
the jet. Figure 2-6c shows how this jet-crossing may occur at later time = t3, during
which the E and UM merge, but the DM-E dipole continues to “peal” itself from the jet.
This “pealing” process is complete at time = t4 when the original eddy E has completely
crossed the jet.
                                                                                       Page 64

         If eddy “E” is an anticyclone, then the same argument as Figure 2-6c will show
that the eddy cannot cross the jet because the resultant (red) vector would be eastward.
However, the jet may become baroclinically unstable (this is not modeled by the reduced-
gravity model because baroclinicity requires an active lower layer), and jet-crossing may
still occur.


Bell, G.I., 1990: Interaction between vortices and waves in a simple model of
geophysical flows. Phys. Fluids A 2 (4), 575-586.

Stern, M.E. and G.R. Flierl, 1987: On the interaction of a vortex with a shear flow. J.
Geophys. Res, 92 (C10), 10733-10744.

Eda, start here please. Thank you.
                                                                                                           Page 65

Pressure Compensation

In Smith and O’Brien [1983; JPO, 13, 1681-1697], you will hear many times the word
“compensation.” All this means is that the lower-layer has very small velocity. In my
note “pom_v14.pdf” Fig.2-5, this just means that the 1 and 2 are of opposite signs. To
see this, we use lower-layer equation (2-52b,c) and set the pressure gradient terms on the
RHS to zero. After integration, and note that there is an arbitrary constant “C” (has to
be a constant – why?), we get:

         1 + C = (  or

         1 = (

choosing C = (. When the free surface convex upwards, 1 > 0, the isopycnal
separating layers 1 and 2 concave downward, 2 < 0.

2-5: An Application of the Two-Layer Model

Isobe [2004; JPO, 34, 1839-1855; equations will be prefixed with an”I”] provides a nice
application of the 2-layer model. The 2-layer system was previously shown in Fig.2-5
(of the “GFD book”), copied here:

Figure 2-5. An illustration in the xz-plane of an ocean model with two (inmiscible fluid) layers with
subscripts “1” and “2” denoting the upper and lower layers respectively.        Symbols are: 1 and 2 =
constant densities of the upper and lower layers respectively, 2 > 1; H1 = constant undisturbed upper-layer
depth; H2(x,y) = undisturbed lower-layer depth which may be a function of (x,y) for variable bottom
                                                                                                          Page 66

topography; h(x,y,t) = layer depth; u(x,y,t) = the x-component velocity; v(x,y,t) = the y-component
velocity; p(x,y,z,t) = pressure; 1(x,y,t) = disturbed free-surface height measured from z = 0; 2(x,y,t) =
disturbed inter-layer height measured from z = -H1.     Note that 1 and 2 can be positive (upward) or
negative (downward), and that h1 = H1 + 1  2.

The equations of a two-layer system were derived in section 2-3, equations (2-51) and (2-
52). The x and y-momentum equations (2-51b,c) for layer-1, ignoring wind stresses,
atmospheric pressure and horizontal viscosity, are:



Taking the curl: (2-5.1c)/x - (2-5.1b)/y] then gives:

         1/t + ( 1u1)/x + ( 1v1)/y + f .u1 = 0                                           (2-5.2)

For layer-2 (see equation (2-52)), the x and y-momentum equations are:

                                         

                                         

where g’ = g/o is the reduced gravity. Taking the curl, we obtain exactly the same
equation (2-5.2) with subscript “1” replaced by “2.” Thus, we write:

         k/t + .( kuk) + f .uk = 0                                                       (I2.1) (2-5.4)

The continuity equations are from section 2-3, equations (2-51a) and (2-52a):

        hk/t + .(hk uk) = 0                                                              (I2.2) (2-5.5)

Split each field variable into time-mean and fluctuation:
                                                                                                        Page 67

         u = <u> + u’, h = <h> + h’, = < > + ’                                                (2-5.6)

where <.> denotes a mean and primes denote fluctuations such that <A’> = 0 for any
variable A. Substitute (2-5.6) into the continuity equation (2-5.5), then take time-

           .(<hk> <uk>) =  .(<hk‘uk’>)                                                  (I2.4) (2-5.7)

The LHS represents divergence of the mean mass flux field, and this is balanced by the
divergence of the mass flux due to eddies – the RHS.

Analogous to heat diffusion by turbulent mixing, say -<w’T’> which we model as:

         -<w’T’> = (KHT/z)/z                                                              (2-5.8)

see Fig.2-5.1,

Fig.2-5.1 On average, turbulent velocities bring warmer water    downward and cooler water upward, as
shown in the sketch.   Therefore,   -<w‘T’>   > 0. This is often called down-gradient diffusion (down the
gradient of T – i.e. heat goes from warm to cold – and is commonly often modeled as <w‘T’> = KH dT/dz,
where KH is the turbulent eddy diffusivity.

we also model the RHS of (2-5.7) as:

         -<hk‘uk’> = k <hk >                                                            (I2.5) (2-5.9)

where k is the diffusivity for layer-thickness (which you can think of also as “heat” since
thicker h means warmer water). Equation (2-5.7) then becomes:
                                                                                           Page 68

         .(<hk > <uk>) = .(k <hk >)                                          (I2.6) (2-5.10)

By using the idea of down-gradient diffusion, i.e. using “diffusivity,” we avoid the very
difficult task of computing turbulence. In the case of large-scale geophysical fluid
flows (e.g. Kuroshio meanders and frontal eddies), we then avoid directly analyzing
‘eddies’ – which are called “geostrophic turbulence.” Recall in GFD course that eddies
in geostrophic turbulence are rigidly constrained by rotation because of the Taylor-
Proudman theorem, so that their motions are “planar” in the horizontal plane
perpendicular to the rotation axis. In the case of Kuroshio frontal eddies, equation (2-
5.9) models in a crude way how eddies transport heat from the Kuroshio to the (cooler)
shelf, and how they also transport cooler shelf water to the Kuroshio – i.e. mixing. The
development of these eddies are say, via baroclinic instability. In Isobe (2004), these
eddies (or some idealized form of eddies) are actually numerically computed.

Similar to time-averaging the continuity, the time-average of the vorticity equation (2-
5.4) gives:

         .(< k> <uk>) + .(< k‘uk’>) + f .<uk> = 0                             (I2.7) (2-5.11)

(note that the time-average of the time-dependent term  k/t is 0). The second term
  .(< k‘uk’>) can again be modeled using some kind of eddy viscosity times the gradient
of the mean vorticity. But it is usually small, and is neglected. Eddy-mixing of
relative vorticity in (2-5.11) is not as important as eddy-mixing of heat in equation (2-
5.10) because | k’| is generally smaller than |f| (and also smaller than |< k>|). Physically,
what is important is the mixing of heat across the Kuroshio, and less so for the relative
vorticity. Mixing of heat (h’) causes layer to change which intuitively is also important
because the change directly affects potential vorticity: e.g. (f + 1)/h1  f/h1.

The f .<uk> in (2-5.11) can be obtained from (2-5.10):

       f .<uk> = f<uk> <hk>/<hk> + f .(k <hk >)/<hk>

                         = +<uk><hk>. (f /<hk>) + f .(k <hk>)/<hk>;(assume f –plane)

                         = +<Uk>. (f /<hk>) + f .(k <hk>)/<hk>;(let           <Uk>        =
                                                                                                Page 69

Substituting this into (2-5.11), we obtain:

           .(< k><uk>) + <Uk>. (f /<hk>) = f .(k <hk>)/<hk>           (I2.8) (2-5.12)
Now the down-gradient portion on the RHS of the previous equation, “ <hk>”, can be
approximated by “thermal wind” for this 2-layer system which can be obtained from the
geostrophic form (time-averaged) of (2-5.1) and (2-5.3). We first do (2-5.1) minus (2-
5.3) (i.e. upper-layer minus lower-layer), in vector form:

          z×<ug1-2>=(g’/f) (<h1>-<1>)(g’/f) <h1>             (since |1| << h1 see fig.2-5)

so that

           <h1> = z×<ug1-2>(f /g’)                                              (2.9) (2-5.13)

Substituting this into (2-5.12) for the upper layer (k = 1), and noting that:

           .( z×<ug1-2>)=<ug1-2>. zz. ×<ug1-2>             (Use .(a×b)=b.( ×a)-a.( ×b))

                                       =      0                Curlz(<ug1-2>)(Since z = 0)

we obtain from (2-5.13) and (2-5.12):

 .(< 1><u1>) + <U1>. (f /<h1>) = f21g’-1 Curlz(<ug1-2>)/<h1>                  (I2.10)(2-5.14)

To derive the lower-layer vorticity equation, we similarly get the “thermal-wind”
equation by doing (2-5.3) minus (2-5.1) (i.e. lower-layer minus upper-layer), and noting
from Fig.2-5 that (1-h1) = 2 = (h2-H2), where H2(x,y) is the bottom topography:

          z×<ug2-1>=(g’/f) ( <1>-<h1>)=(g’/f) (<h2>H2)                           (2-5.15)

from which we obtain:

           <h2> = z×<ug2-1>(f /g’) + H2                                             (2-5.16)

Substituting this into (2-5.12) for the upper layer (k = 2), we obtain:

 .(< 2><u2>)+<U2>. (f/<h2>)=f22g’-1Curlz(<ug2-1>)/<h2>f .(2 H2)/<h2>
                                                                                             Page 70


However, if we assume in (2-5.15) that <h2> is >>        H2, i.e. that the topographic slope
varies slower than the interfacial slope, then:

 .(< 2><u2>) + <U2>. (f /<h2>) = f22g’-1 Curlz(<ug2-1>)/<h2>              (I2.11)(2-5.18)

2-6: Heuristic derivation of the quasigeostrophic (QG) potential vorticity (PV)
equation for a stratified fluid

The inviscid x, y & z-momentum and continuity equations are:

       u/t + uu/x + vu/y + wu/z -fv = -r-1(p/x)                        (2-6.1a)
       v/t + uv/x + vv/y + wv/z +fu = -r-1(p/y)                        (2-6.1b)
       p/z = -g                                                                (2-6.1c)
       u/x + v/y + w/z = 0                                                  (2-6.1d)

Also, a beta-plane is assumed such that:

       f = fo + oy                                                               (2-6.2)

where fo is the Coriolis parameter at y=0 (around where later we will develop the
instability analysis), o = df/dy, and it is also assumed that |oy| << |fo|. For o  10-11
m-1 s-1, and |fo|  5×10-5 s-1 or larger (poleward of 20oN/S), we restrict |y|  1000 km or

The symbols have the usual meanings, and r = reference density which is a function of z
only. The basic QG-approximation idea is that the motions are close to being
geostrophic, so we will use the geostrophic velocity assuming constant f  fo (which is
‘zeroth-order’) to evaluate the “difficult” terms which in (2-6.1a,b) are the non-linear and
“o” terms.

To the zeroth-order approximation, the motions are nearly geostrophic, so that:

       -fovo = -r-1(po/x)
       +fouo = -r-1(po/y)                                                      (2-6.3a,b)
                                                                                                                     Page 71

          so that

                    uo/x + vo/y = -wo/z  0                                                         (2-6.4)

          Therefore, wo  0 since it is  0 at the surface. In any case, the zeroth-order or
          geostrophic approximation of w = wo is small compared to |uo| or |vo|.6

          Substitute the (uo, vo, wo0) into the non-linear and beta parts of (2-6.1), we then obtain
          the following equations for the next-order term with subscripts “1” (u1, v1, p1):

                    do/dt{uo} – fov1  oyvo = -r-1(p1/x)                                              (2-6.5a)
                    do/dt{vo} + fou1 + oyuo = -r-1(p1/y)                                              (2-6.5b)


                    do/dt = /t + uo/x + vo/y.                                                       (2-6.6)

          Taking the (z-component of the) curl of (2-6.5) to eliminate the pressure, and using the
          first of (2-6.4) that uo/x + vo/y  0, we obtain:

                    do/dt{o} + ovo + fo     H.u1   = 0, where    H=i/x+j/y     and u1=(u1,v1);

          or,       do/dt{o + oy} = fow1/z                                                            (2-6.7)

          after using (2-6.6) and (2-6.1d), i.e. u1/x + v1/y + w1/z = 0.

          For  = constant, by vertically integrating (from z=-H to z=0) equation (2-6.7) we
          recover the QGPV equation for a homogeneous fluid, which we studied previously in
          section 2-1 (see equation (2.10)). For stratified fluid, we need to evaluate w1/z using
          the density (or buoyancy b = -g/r) equation [see section 11.1.1 of Marshall and Plumb,
          2008] which in the absence of heat/salt sources and sinks on the RHS is:

                    /t + u/x + v/y + w/z = 0                                                  (2-6.8)

6   Typically, w = Ekman pumping near the surface  ×o/(fo)  0.1(N/m2)/[104(m).10-4(s-1).103(kg/m3)]  10-4 m/s  10
m/day, for a strong wind stress curl of 0.1 N/m2 over a distance of 10 km.   Thus comparing to typical |uo|  5×10-2 m/s, the
“w” is indeed small.    In the open ocean, the “w” is typically smaller than 10 m/day.
                                                                                         Page 72

The density is given by:

         = r(z) + ’(x,y,z,t)                                                (2-6.9)

and the hydrostatic equation is then assumed also for the perturbation density ’ (and
pressure po):

        po/z = -g’   (also of course pr/z = -gr)                         (2-6.10a)

or      (po/z)/r = -g’/r = b’                                             (2-6.10b)

Using the QG-approximation on (2-6.8), we get:

        do/dt{’} + w1(dr/dz) = 0,                                            (2-6.11a)

or multiplying by –g/r, we get do/dt{-g’/r} + w1(-gdr/dz)/r = 0, i.e.

        do/dt{b’} + w1N2(z) = 0                                                (2-6.11b)

where N2(z) = (-gdr/dz)/r is the squared buoyancy (or Brunt-Vaisala) frequency. We
see that (in the absence of sources and sinks) the vertical velocity w1 is related to the
isopycnal movement due to the geostrophic flow:

        w1 = do/dt{b’/N2(z)}, or also = do/dt{’}/{-dr/dz}.                  (2-6.11c)

The last form shows, since -dr/dz > 0, that deepening isopycnal (i.e. do{’}/dt < 0)
corresponds to downwelling, and shallowing isopycnal corresponds to upwelling.

We now use w1 in (2-6.7):

        do/dt{o + oy} = fow1/z = fodo/dt{[b’/N2(z)]/z}

or upon using the hydrostatic equation (2-6.10b):

        do/dt{o + oy} =  fodo/dt{[(po/z)/(rN2)]/z}

i.e.,   do/dt{o + oy + fo[(po/z)/(rN2)]/z} = 0                          (2-6.12a)
                                                                                                   Page 73

or,    do/dt{(       2
                         po)/(for) + oy + fo[(po/z)/(rN2)]/z} = 0                (2-6.12b)

since from (2-6.3a,b), we have for the geostrophic vorticity:

       o = vo/xuo/y =            2
                                           po/(for)                                    (2-6.13)

For the ocean, r can be assumed to be  constant in (2-6.12b), so that equation becomes:

       do/dt{       2
                        (po/for) + oy + [(po/for)/z)(fo /N2)]/z} = 0            (2-6.14)

From the geostrophic relation (2-6.3), the po/(for) is equivalent to geostrophic stream

        = po/(fo r)                                                                   (2-6.15)

so that (2-6.14) is usually expressed in terms of :

       do/dt{       2
                         + oy + [(/z)(fo /N2)]/z} = 0                            (2-6.16a)

The quantity inside the {..} is the QGPV for a stratified fluid:

       Q=       2
                     + oy + [(/z)(fo /N2)]/z                                     (2-6.16b)


We use the following scales:

Table 2-6.1
Scales                  Variables               So that… (primes denote nondimensional variables)
    L                     (x,y)                                   (x,y) = L (x’,y’)
   D                        z                                          z = Dz’
   U                      (u,v)                                   (u,v) = U(u’,v’)
  L/U                       t                                        t = (L/U)t’
fo rUL                    po                      po = fo rULpo’, using geostrophic eqn.(2-6.3)
   UL                                                       = UL’, using eqn. (2-6.15)
   Ns                       N                        N = NsN’, for some typical Ns value of N

Then the non-dimensionalized form of (2-6.16) is (dropping all primes):
                                                                                              Page 74

       do/dt{       2
                         + y + [S-1(/z)]/z} = 0                             (2-6.17a)

and    Q=       2
                     + y + [S-1(/z)]/z                                      (2-6.17b)

where  = oL2/U, S = (LD/L)2, and LD = NsD/fo is the baroclinic Rossby radius based on
ocean’s depth D. Note that Ns2D ~ g/r ~ g’, the reduced gravity, so that Ns2D ~ g’D
~ baroclinic phase speed squared. The nondimensional  = o/[(U/L)/L] = o/(o/L) is a
measure of the ratio of gradient of planetary vorticity “o” to the gradient of flow
vorticity “o/L”. By assuming that  ~ O(1) in (2-6.17), we are examining an important
GFD case in which the planetary vorticity gradient contributes equally with the relative
vorticity gradient to the overall vorticity balance.

If the scale of motion is small compared to the Rossby radius, L << LD so that S >> 1, the
effect of vortex-tube stretching to the PV-balance (i.e. the [S-1(/z)]/z term) is small,
and the flow’s relative vorticity (i.e. 2) becomes important. If on the other hand the
scale of motion is large, L >> LD, then S << 1, and the flow’s vorticity becomes relatively
unimportant, i.e. the flow looks horizontally more uniform.

2-7: Baroclinic instability

Fig.2-7.1 suggests that baroclinic instability may occur along special paths “(iii)” across
slanting isopycnals which of course imply the existence of vertical velocity shear by the
thermal-wind equations (from geostrophic eqn.(2-6.3a,b) and hydrostatic eqn.(2-6.10b)):

       k × fo uo/z =  b’                                                        (2-7.1a)

or     fo uo/z = k × b’                                                          (2-7.1b)

We will now find the condition – what type of vertical shears can produce baroclinic
instability? Or, in general, what type of vertical and horizontal shears can produce
baroclinic and barotropic instabilities? We will use the non-dimensionalized QGPV
equation (2-6.17) so that all variables are from here on non-dimensional.
Dimensional variables will have subscript “*”. Consider an initial or background
state of purely zonal flow Uo(y,z) with stream function (y,z):

       Uo(y,z) = /y (note that Vo = /x = 0)                                 (2-7.2)
                                                                                             Page 75

Introduce perturbation “ ”, so that the total, time-dependent stream function  is:

       (x,y,z,t) = (y,z) + (x,y,z,t)                                             (2-7.3)

Fig.2-7.1. Schematized diagrams showing various scenarios when exchange of fluid
parcels “A” and “B” along the blue path line either leads to increased potential energy so
that the mass center of the system moves up (stable), is unchanged (neutral) or to
decreased potential energy so that the mass center of the system moves down (maybe
unstable). It is easy to see that only for the special path within the “wedge” formed by the
isopycnal (red line) and the horizon (dashed line) in “(iii)” can the exchange possibly
lead to instability (such that after the exchange the two parcels move away from each
other). Such an instability is called baroclinic instability.

The function is perturbation to the initial state ; it represents the structure of the
evolving perturbation field. Substitute (2-7.3) into (2-6.17):

       (/t + Uo/x    y/x   +   x/y){q   + 2/y2 + y + [S-1/z]/z} = 0
                                                                                              Page 76

or     (/t + Uo/x)q + J( ,q) +     x/y   =0                                (2-7.4)

where q(x,y,z,t) is the perturbation PV:

       q=    2
                 + [S-1( /z)]/z                                                (2-7.5a)

and    o = 2/y2 + y + [S-1/z]/z                                          (2-7.5b)

is the background or initial state’s PV, and its meridional gradient is (using eqn(2-7.2)):

       o/y = 2Uo/y2 +   [S-1Uo/z]/z                                    (2-7.5c)

Notice how the thermal-wind relation leads to the vertical shear Uo/z that then appears
in this last equation for o/y.

How does the structure of Uo(y,z) determine the evolution of the perturbation field ?
That is, given a particular background or initial state Uo(y,z), will the perturbation
“injected” on the flow grows or decays? If grow, then the initial state is unstable with
respect to the perturbation . To show that Uo is stable, we must check all possible ’s.
On the other hand, to show instability, we only need to find one perturbation to which the
initial state Uo is unstable.

Linear Stability Analysis:

We assume that | | << 1 so that the J( ,q) in (2-7.4) is dropped:

       (/t + Uo/x)q +    x/y   =0                                          (2-7.6)

The boundary conditions are that vertical velocity w* = wo* + w1* + O( 2) is zero at z=0
(surface) and z=-1 (ocean’s bottom). Here = U/(foL) is the Rossby number which is
assumed to be small, and we already noted previously (see eqn.(2-6.4) and discussion)
that wo* = 0, so that w*  w1* The non-dimensional w1 can then be expressed in terms of
 using equations (2-6.11c), the hydrostatic relation (2-6.10b), and (2-6.15) which gives:

       b*’ = fo*/z*, hence w1* = do*/dt{fo*/z*/N2}                          (2-7.7a,b)

and    w1 = do/dt{ S-1/z} = 0 at z=0 & -1.                                     (2-7.8)
                                                                                                Page 77

A perturbation energy equation can be derived and it can be shown that:

       s (o/y)(<2>/t) dydz = 0                                              (2-7.9)

where  is the meridional displacement of fluid elements defined by:

       /t + Uo/x =  /x                                                       (2-7.10)

and <.> is a zonal-averaging. From (2-7.9), we see that if there is to be a growth in the
displacement of fluid elements in time, i.e. if <2>/t > 0, then o/y must be
somewhere positive and somewhere else negative in the yz-plane, or o/y can also be
identically zero everywhere. In other words, o/y must vanish on a line in the yz-
plane. Clearly, this condition (i.e. that o/y must vanish somewhere) is not sufficient.

In Tulloch et al.’s paper (JPO, 2011), they use a  that is a function of x,y & z:

       (x,y,z) = V(z)x  U(z)y

Substituting this into (2-6.17) then gives their eqn.2 (instead of eqn.2-7.6 above). Their
second equation is just equation (2-7.8). Their equation (1) is just our equation (2-7.5c)
above with 2Uo/y2 = 0 and their last term (in their eqn.1) is a consequence of applying
the boundary condition (2-7.8) at the surface z=0.
                                                                                      Page 78

Simple Ideas of Frictional Spindown and Buoyancy Shutdown
Ocean is spun down by friction  frictional spindown. But with stratification and
bottom slope, it is possible that bottom mixing (w/stratification+slope) produces
horizontal density gradients, hence thermal-wind shear in the mixed layer so that the
interior flow “sees” a slippery bottom and Ekman pumping & friction are shutdown 
buoyancy shutdown.

Rough derivations:
To compute /y due to bottom slope h/y and bottom mixing (downwelling case), we
consider an isopycnal =constant; then  = 0 = /z h + /y y              so that
/y = -/z.h/y so that (g/o)/y = N .hy, where N = -g(/z)/o. Or,
                                                   2             2

with b = -g/o we have: b/y = -N2hy where N2 = b/z.
By geostrophy within the mixed layer, we have fu = -p/y/o, so that fu/z = +(g/o)
/y = -b/y = N2hy from previous equation. So that integrating across ML, from z=-
h to z=-h+, where  = ML thickness:
ui – u = (-h+-z) N2hy assuming that N2 is constant. Therefore: u=ui-(-h+-z) N2hy.
(a) buoyancy shutdown arises only if bottom has slope h/y0;
                                                                                         Page 79

(b) if hy is =0, then u = ui extending all the way to bottom, and the interior flow is then
    subject to frictional spindown (because then b = -rui operates);
(c) basically, at buoyancy shutdown, the ui is canceled by thermal-wind part and b=0.
Garrett, MacCready & Rhines, Ann. Rev. Fluid Mech, 1993, 291-323.
Chapman, JPO, 2002, 336-352.
                                                                                Page 80

Chapter 3. Barotropic and Baroclinic Instabilities

Barotropic Instability: Rayleigh-Kuo's Theorem

If               , then the jet is stable (For eddy, if          )

Consider potential vorticity equation:


Consider parallel basic flow in x direction, U=U(y)
                       uT  U  u
                       vT  0  v
                       T  Z 

                       Z 

After linearization and neglect quadratic and high order terms, (I1) becomes:

                  Z
             U    v     v  0                                                (I4)
          t    x    y


 2            U       2
  (        )              2v                                                (I5)
t   Z y      Z y x

Integral over a wave length, 2L=wave length
                                                                                                 Page 81

        2            U           L

   L   Z y )dx    Z y            2  vdx                                                   (I6)
t                               L

The 2nd term of LHS is zero, because of the periodicity, the RHS will be:

       v u  1 v 2    uv   v   1 v 2    uv 1 u 2
v  v(  )         (     u )         (           )
       x y  2 x       y   y   2 x      y 2 x
           1  2
                                            L
Lvdx L 2 x (v  u )  y (uv) dx   L y (uv)dx

                                  

Thus, (I6) becomes:

       2                 L
t    Z y )dxdy  2L y (uv)dydx
     (                                                                                               (I8)

since uv=zero at y=±b are walls, so RHS of (I8)=0
       2
     (    )
     Zy
             dxdy  0                                                                                (I9)

                                   2
For there to be an instability,       >0, then (I9) only vanish if β+Zy changes sign somewhere in "y"

Suppose that

                                                               (  Z y ) 1    -1   2 y5
1 y  2                                             (I10)

0  y 1
                                                                                                    Page 82

Then, clearly α changes sign at some "y". Therefore, if ∂/∂t(ζ2)>0, flow is unstable, then
∫α∂/∂t(ζ2)dy=∫αdy =0, so that (I9) is satisfied.

In summary, we say that if the flow is unstable, i.e. if ∂/∂t(ζ2)>0, it is necessary for α to change sign; or
we can also say that "a change in sign in α" is necessary for flow to be unstable.

But, "a change in sign in α" is not sufficient for the flow to be unstable. For example, the integral
∫α∂/∂t(ζ2)dy=-∫αdy =0, then, ∂/∂t(ζ2)<0, it is stable for the flow!

Again, we have (I4):

         Z
   U     v(   )  0                                                                               (I4)
t    x     y

u v               
    0u      , v                                                                              (I11)
x y         y      x

where ψ=perturbation stream function (note: ζ= 2ψ)
so that, (I4) becomes:

                              o
        U             v(
                                    )0       (Rayleigh-Kuo Equation2)                                (I12)
t             x               y

                           o Z
where (Pedlosky)                     U yy                                                      (I13)
                           y   y

Let   real{  y e ik ( x ct) }   (real part of it)                                                (I14)
                                                                                                                                                         Page 83

be a disturbance of wave number k in x (periodic), and

σ=-ikc=-ik(cr+ici)=ikcr+kci                                                                                                                               (I15)

which is -i*kc(frequency), such that

kcr=real frequency                                                                                                                                        (I16a)
kci=growth rate                                                                                                                                           (I16b)

Thus, flow is unstable if the imaginary part of c, i.e. ci>0
From (I14), only take the real part:

u  y   y e ik ( xct)
v  x  ik e ik ( xct)                                                                                                                         (I17a, b, c)
           2
                   k       2
                                            yy   e   ik ( x ct )

(I12) then becomes:

 ikc  k 2               yy    Uik k                2
                                                                         yy    ik   
                                                                                                               oy
(U  c)  k          2
                                  yy      
                                                  0, or                           k   2
                                                                                                 yy    (U  c)     0

Multiply (I18) by                                *, which is the conjugate of                                 , so that    *   =abs( )2=real, and also

                                                     2
    *                  ( *             )                , then integrate                             , we have
        yy                         y              y

                                                                        dy    (U c) 
Y                                  Y                                              Y

Y y ( *
                      y ) dy                          k2                                                   dy  0
                                                                      2                      oy          2
                                             y                                                                                                            (I19)
                                  Y                                             Y

The 1st term on LHS is zero, since we will assume either
                                                                                                                             Page 84

(i) Y and -Y are walls so that                                  = *=0 there, or

(ii) Y and -Y are infinite and                                     decays → 0 as y →±∞

so (I19) becomes:

                                 dy    U c
Y                                        Y                                    
                                                              (U  cr  ici ) 
                  k2                                                                       dy  0
                               2                       oy                               2
          y                                                 2
Y                                       Y                                  

where U  c  (U  cr  ici )(U  cr  ici )  real , so

                                  U c
Y                                                                             Y
                                                                                    oy ci
                                                                        dy  i 
                   k2                                (U  cr )                                           dy  0
                               2                                    2                                 2
                                                                              Y U  c
                                                  2                                           2

It means real=0, and imaginary=0. We focus on the imaginary part since we want to find out more
about "ci" (which determines instability etc.), so

      U c                        dy  0
ci                     2

Since "c" is complex, it occurs in pair in conjugate, so that it has both positive ci (growing) and
negative ci (decaying) parts. So all we need to show for instability is that ci≠0

                                                                                    U c                      dy  0
(I21) can only be satisfied if either ci=0 or                                                     2

If there is an instability, i.e. ci≠0, then

          oy
 U c                         dy  0

But since both U  c                                               are nonzero,  oy must change sign somewhere in y: -Y≤ y≤ Y.
                                              2                2

This is the necessary condition that we obtained in (I9).
                                                                                                      Page 85

                                                                              oy
On the other hand, if  oy does not change sign anywhere in "y", then    U c
                                                                                            dy cannot be 0, so

that ci must be 0, and the flow is stable. So the condition that

"  oy does not change sign in y" is sufficient to ensure stability.                                   (I24)

This proves the Rayleigh-Kuo Theorem.
                                                                                                    Page 86

Chapter 4 (Gill's Chapter 4; Gill’s Figures and Equations are denoted by G*):

1. What is a material element? [p.63]. [material element = infinitesimally small fluid sample that
retains its identity.]

2. Why is it that if molecular diffusion is =0, the material element will consist of the same
particles? What are "particles" in this case? [p.64]. [If diffusion=0, then molecules (=particles)
cannot spread, so material element consists of the same particles.]

3. Why is (G4.2.1) valid (i.e. why can the LHS be broken up into the 3 terms on the RHS)?
[Use d(ab)/dk = adb/dk+bda/dk).]

4. Derive (G4.2.3), (G4.2.4) & (G4.2.5).

5. What is the difference between advection and convection? [p.66]

6. Derive (G4.1.8) from (G4.3.2) and (G4.2.4) (== (G4.2.3)).

7. Derive (G4.2.5) [Use Fig.G4.2].

8. Explain the flux vector F in (G4.3.3) - as advective + diffusive fluxes.

9. Note how (G4.3.4) [or G4.3.7] is similar to the Salinity equation used in POM. Is it
identical? Why (yes or not)?

10. The entropy of a material element is fixed if the element is moved
      (a) without exchanging heat with its surrounding and if
      (b) it retains its fixed composition.
    The motion is then called isentropic.

11. So, the salinity or 'salt' equation is fairly straight-forward, since the Qv on p.67 = rho*s = mass of
salt per unit volume; and the "mass of salt" is a 'measurable' scalar quantity - i.e. a "mass."

12. Temperature is somewhat different however. It is still a scalar of course, but it is closely tied to the
agitation of molecules and therefore also to the entropy or "heat content." Temperature is derived from
the "heat equation." To understand this, we need to explore a little thermodynamics.
                                                                                                Page 87

13. The "heat" equation is obtained from the 2nd law of thermodynamics:

         dE = Td  pdvs                                                                   (G3.2.1)

       E = internal energy per unit mass
       Vs = specific volume (volume per unit mass)
       T and p = temperature and pressure respectively
        = specific entropy (entropy per unit mass)

Energy increases because of (i) agitation (molecules) d, and (ii) compression dvs < 0. We assume no
phase change and salinity is fixed – i.e. the fluid is of fixed composition; so  = (p, T).

       Td  change in heat content per unit mass
        pdvs = work done by p in compressing fluid’s volume per unit mass.

14. Eqn.(G3.2.6) is used often later – so we will derive it. From (G3.2.1):

       Td = dE + pdvs                                                                     (N4.1)

So, increase in heat content for a change in T while keeping p constant is then, from (N4.1):

       Cp  T(/T)p = (E/T)p + p(vs/T)p                                              (N4.2)

The Cp is called the specific heat at constant pressure.

Similarly, we obtain the following for increase in heat content for a change in p while keeping T

       T(/p)T = (E/p)T + p(vs/p)T                                                   (N4.3)

Now, recall that  = (p, T), the Td (i.e. LHS of (G3.2.6)) is:

       Td = T(/T)p dT + T(/p)T dp
                = CpdT + T(/p)T dp                                                      (N4.4)

where the definition of Cp from (N4.2) has been used for the first term. For the second term,
T(/p)T, we may use (N4.3), but we can get rid of the “E” as follows. We can
[T(N4.3) p(N4.2)]:
                                                                                                   Page 88

(/p)T +T( /pT)T( /Tp) = ( E/pT)( E/Tp)+p( vs/pT)p( vs/Tp)(vs/T)p

so that,

           (/p)T = (vs/T)p                                                            (G3.2.5)

which can be used in (N4.4) to give (G3.2.6):

           Td = CpdT  T(vs/T)pdp.                                                       (G3.2.6)

15. A process is reversible if it changes very, very slowly due to infinitesimal changes in some property
of the system, so that no dissipation occurs. For example, if a gas in an insulated cylinder is
compressed by pushing a piston (see Fig.N4.1) and there is no friction between the piston and the
cylinder wall, then if the piston is let go the compressed gas would expand back so the piston is
reversed to its original position. If there is friction, then heat is permanently loss and the process is

In practice, no process is reversible. But if the applied changes are slow and/or if the system’s
response is fast, then the process may be considered as being reversible. Gill assumes this [p.42,
approx. line 6].
                                      gas     cylinder


16. Since Td = increase in heat content (of the system), then if the process is isentropic and reversible,
then it is also called adiabatic. In our case, we assume “reversible processes,” so isentropic =

Thermodynamic Definition: An adiabatic process a process in which no heat is transferred to or from
the working fluid.

17a. When a parcel of fluid sinks adiabatically, say into the deep ocean, it experiences increased
pressure and therefore warms up. Similarly, when air parcel rises (to mountain top) it experiences less
pressure and cools down. In this case, “adiabatic” means that there is no turbulent or molecular
diffusion and no radiation across the parcel’s boundary, and also no change in composition or phase.
                                                                                                    Page 89

Thus, in the case of adiabatic sinking or rising fluid parcel, because of pressure compression
(expansion), DT/Dt  0; it would be nice to define a quantity such that the D(.)/Dt = 0.

Also, in both cases (sinking or rising parcel), the (vertical) column of fluid (air) is apparently unstable
(i.e. “warm” is under “cool”). But this is not true – i.e. the fluid is actually stable. To avoid this
apparent confusion, and to make D(.)/Dt = 0, the quantity potential temperature  is introduced.

The potential temperature “” is the temperature of the parcel at pressure p would have if the parcel
were brought adiabatically (i.e. isentropically; i.e. without exchange of heat with its surrounding) to a
reference (i.e. common) pressure pr.

By using , we eliminate the apparent warming or cooling due to changes in pressure. Since now the
various temperatures at different depths have been brought to a common pressure, we can compare
them without worrying that the ’s represent false temperatures. The variation of  with height can
therefore be used to judge if the vertical column is statically stable or not.

The situation is a bit similar to trying to judge if person (A) on Jupiter is heavier than person (B) on
Pluto (Fig.N4.2). It would appear that (A) is heavier than (B). But to more fairly judge them, we
can move them (without loosing even a single hair; constant entropy) all to a reference place, say the
earth, then re-weigh them. The weight measured on earth may then be called “potential weight.”

                  (A) I’m 2 tons                        (B) I’m

           Jupiter                              Pluto


17b. An example (for air) and some exercises will make the concept of potential temperature clearer.
It is convenient to do this with air than with water because air behaves like an “ideal gas” and explicit
formulae may then be obtained.

Suppose we follow an air parcel from some height with (p,T) which sinks adiabatically (i.e.
isentropically) from a high mountain to the ground. For isentropic process, we have from (G3.2.6)
with d = 0:
                                                                                                  Page 90

       dT = T(vs/T)p,s dp/Cp = T(/T)p,s dp/(2Cp).                                    (N4.5)

Assume that air obeys the ideal-gas law:

        = p/(RT)                                                                           (N4.6)

[note that high p gives rise to high T and also smaller volume (per unit mass) vs = -1], we have then,

       (/T)p,s = p/(RT2)                                                                (N4.7)

and therefore (using (N4.6) again):

       (/T)p,s /(2Cp) = +(R/Cp)/p = (2/7)/p   =  /p,                                  (N4.8)

where  = (R/Cp) = 2/7 for dry air from Gill’s (G3.3.2). Equation (N4.6) then becomes:

       dT/T =  dp/p                                                                        (N4.9)

Integrating (N4.9) from a height where pressure is p and temperature is T to the reference level where
the pressure is pr and temperature is, by definition, the potential temperature , we have

       ln()ln(T) =  [ln(pr)ln(p)], or

       ln(/T)      =      ln(pr/p)2/7, or

       /T          =     (pr/p)2/7.                                                        (N4.10)

It is clear that potential temperature is only meaningful when we have prescribed a reference level; the
most common one is the surface, where the pressure is  1013 mb. It is convenient to simply let pr =
1000 mb.

Exercise 17b1: Air parcel “A” at a height where p = pA = 100 mb has an in situ temperature TA = 175
K. Another parcel “B” closer to the ground where p = pB = 400 mb has a TB = 220 K. Is the air column
stable? [Hint: calculate the corresponding .]

Please do not read beyond this line before you attempt exercise 17b1.
                                                                                                   Page 91

More Thoughts on Exercise 17B1: In the above exercise, TA < TB and according to (N4.6), the density
A seems to be greater than B, i.e. A > B, i.e. air column appears to be unstable. But (N4.6) cannot
be used this way to determine the densities because the pressures pA and pB are different. In fact, since
pA < pB, it may be possible for A to be less than B. This is why we introduce the potential
temperature  by “moving” both “A” and “B” adiabatically to the common level where p = pr. By
doing this “moving” adiabatically, we ensure that both parcels’ heat-gain (and rise in temperature) as
they are moved closer to the ground is “kept inside” the parcels. Therefore, from (N4.10), A  338 K
> B  286 K, so that the air column is in fact stable.

Exercise 17b2: If we increase TB, the air column would then tend to the state of “less stability.” (we
will learn the precise definition later). Calculate the value of TB such that the air column is neutrally

Exercise 17b3:    In Exercise 17b1, is the air column stable if TB = 300 K?

Exercise 17b4: Figure N4.3 plots T as a function of p for three values of ’s. Think about this in
view of the exercises you have just done. By choosing a few numbers, convince yourself that you
have understood the concept of .

Fig.N4.3. The drop in local temperature (T) with height (or pressure) for various potential temperatures
 = 240, 280 and 320 K, computed from (N4.10), rewritten as T =  (p/pr)2/7, pr = 1000 mb.
                                                                                                     Page 92

Answers to Exercises 17B2&3: TB_neutral = 260 K; unstable.

18. Static Stability (section G3.6).

Suppose the fluid is at rest; to judge a parcel’s static stability, we move the parcel up or down
isentropically, calculate its density change, d, then compare this change with the change in the
surrounding (i.e. ambient) density, da. If the parcel is moved up, then for stability, d must > da so
that the parcel can sink back to its original depth when let go.

The general formula for the change in density d(p,T,s) of the fluid (ambient or parcel fluid) as a result
of a change in height “dz” is, by chain-rule:

         d(p,T,s)  (d/dz).dz = [ (/p)T,s (dp/dz) + (/T)p,s (dT/dz) + (/s)p,T (ds/dz) ].dz

where s = salinity. For a parcel of fluid that is moved isentropically, with  = (p, T), and fixed s; see
“13”, (G3.6.8) becomes:

         d = (/p)T,s dp + (/T)p,s dT                                                    (N4.11)

We can replace “dp” and “dT” with “dz” as follows.           From (G3.2.6), for isentropic process (d = 0):

         dT = T(vs/T)p,s dp/Cp = T(/T)p,s dp/(2Cp) = Tdp/(Cp) = dz                  (G3.6.1)

where  = thermal expansion coefficient and  = adiabatic lapse rate:

                    = -1 (/T)p,s > 0                                             (G3.6.2)

                    = gT/Cp,                                                                 (G3.6.5)

and the last equality in (G3.6.1) (involving ) is obtained by using the hydrostatic equation:

                   dp/dz = g                                                                 (N4.12)

Substituting “dT” from (G3.6.1) and “dp” from (N4.12) into (N4.11), which becomes:

         d =  [g(/p)T,s + ] dz                                                         (G3.6.7)
                                                                                                                           Page 93

Reminder: equation (G3.6.7) is the change in density, dpar, of the fluid parcel as it is moved. On the
other hand, equation (G3.6.8) is general, and gives (for example) the change, damb, of the ambient (i.e.
surrounding) fluid density with z. Suppose that the parcel is moved upward (i.e. dz > 0), we must
have, for stability, that dpar > damb, so that the parcel can move back (i.e. downward) to its original
position when let go:

           [ + dT/dz]  s ds/dz > 0                                                                                (G3.6.9)


          Cp-1 2gT +  dT/dz  s ds/dz              > 0                                                             (G3.6.10)


          s = -1 (/s)p,T > 0                                                                                     (G3.6.11)

is the expansion coefficient for salinity (I use s to avoid possible confusion with the planetary beta ).

Exercise 18.1: Show that the same stability condition is obtained if similar arguments as above are
applied but the parcel is moved downward instead of upward.

Answers: If parcel is moved down, then for stability we must have dpar < damb.   Then from (G3.6.7) and (G3.6.8):
          [g(/p)T,s + g(/p)T,s +  + dT/dz  sds/dz] dz < 0.
But since dz < 0, the […] must be > 0, and (G3.6.9) is obtained.

19. The buoyancy frequency (or BruntVaisala frequency; section G3.7) squared is:

          N2 = g [ ( + dT/dz)  s ds/dz] = g[Cp-1 2gT +  dT/dz  s ds/dz].                                      (G3.7.1)

When N2 > 0 (< 0), then the medium is stable (unstable). In a stable medium, N is real and represents
the frequency of (vertical) oscillation of a fluid parcel – closely linked to internal waves. In the ocean
N  10-4 to about 10-2 s-1, or periods of about 10 minutes to a few hours.

20. Since  is the temperature that a parcel has if it is moved isentropically to a reference pressure (see
“17a”), then a particular  has a particular entropy . Therefore  = constant when  = constant.
So if s = fixed,  = (), i.e.  is a function of  only. The function can be obtained by fixing p = pr
(reference pressure) in (G3.2.6) which then becomes, since T is then = :
                                                                                                      Page 94

        d/d = Cp(pr, )/                                                                    (G3.7.6)

21. Assume that s = constant (i.e. fixed composition).      Instead of p and T, we can use p and , so that

         = (p, ).

Then for a parcel that is moved isentropically (i.e. d = 0), the change in its density is (c.f. N4.11)

        dpar = (/p),s dp + (/)p,s d                                                  (N4.13)

The change of the surrounding density is from (G3.6.8) with T replaced by :

        damb(p, ,s)  (d/dz).dz = [ (/p),s (dp/dz) + (/)p,s (d/dz) + (/s)p,T (ds/dz) ].dz


For stability, dpar  damb > 0 (if dz > 0) as in “18”, then:

        g-1N2 = ’(d/dz)  ’s (ds/dz) > 0                                                    (G3.7.9)

where       ’ = -1(/)p,s     and ’s = -1(/s)p, . In the form written as (G3.7.9), if s =
fixed, then the stability condition becomes very simple, i.e. the fluid medium is stable if the potential
temperature increases upwards. This fact was used in Exercise 17b.1-4.

22. The Heat Equation (isentropic motion):

If fluid elements do not exchange heat with their surroundings, and further more if they retain a fixed
composition, then their motions are adiabatic or isentropic, and we have:

        D/Dt = 0 = Cp(pr, )-1 D/Dt                                                         (G4.4.1)

the last equality is from (G3.7.6). Thus for adiabatic motion, the heat equation takes on this simple

        D/Dt = 0,                                                                             (N4.15)
                                                                                                   Page 95

i.e. the potential temperature is constant following the parcel’s movement. The heat equation looks a
little more complicated if written in terms of the in situ temperature. From (G3.2.6), we have

          TD/Dt  CpDT/Dt  (T/)Dp/Dt = 0.                                              (N4.16)

The heat equation can also be written in terms of the internal energy E (per unit mass; why?    See the
pressure work term in eqn.G3.2.1). From (G3.2.1):

          TD/Dt  DE/Dt + pDvs/Dt = 0.

This last equation can be re-written by making use of the continuity equation (G4.2.4) which is:

          /t + .(u) = 0,                                                               (G4.2.4)

so that

          DE/Dt + E*LHSof(G4.2.4) = t + u. E + E*LHSof(G4.2.4) = (E)/t + .(uE); or:

          DE/Dt = (E)/t + .(uE).                                                      (G4.3.6)

Then,      times the second equation of (N4.17) gives:

          DE/Dt  (E)/t + .(uE) = pvs-1Dvs/Dt = p .u,                               (G4.4.3)

in which (G4.3.6) and (G4.2.2) have been used.     Equation (G4.4.3),

          (E)/t + .(uE) = p .u,                                                       (N4.18)

says that the rate of increase of energy per unit volume following a fluid parcel, D(E)/Dt > 0, is equal
to the rate at which work is being done due to the parcel’s compression (i.e. convergence) .u < 0.
Similarly, there is a decrease in E if the parcel expands (i.e. divergence) .u > 0.

By the way, one can use Figure G4.2 and its derivation of the time rate of change, /t, of mass per unit
volume (= ), to apply to energy per unit volume = E.

23. The Heat Equation (non-isentropic motion):
                                                                                                Page 96

Then additional terms are added to account for various heat transfers and sources: Frad, k T and QH.
Similar to the “advective heat flux density” Eu, we can identify Frad and k T as the “radiative and
conductive heat flux densities” respectively. In other words, they are “fluxes” that ‘flow’ in and out
of the side faces of the rectangular volume element of Figure G4.2 (see the comment of the last
sentence of “22” above). By contrast, the “QH” (= rate of heating per unit volume) is a heat source (or
sink) inside the volume element and is therefore grouped together with the pressure work p .u. The
result is therefore equation (G4.4.4):

       (E)/t + .(Eu + Frad  k T) = QH p .u.                                         (G4.4.4)

An alternate form using the in situ temperature T is obtained by equating  times (N4.16) to  times
(N4.17), and then using (G4.3.6) and (G4.4.4), to obtain:

       CpDT/Dt  T Dp/Dt = .( k T  Frad) + QH.                                         (G4.4.6)

Yet another alternative form using the potential temperature  is simpler, using (G4.4.1), (N4.16) and

       TCp(pr, )-1 D/Dt = .(k T  Frad) + QH.                                         (G4.4.7)

The Princeton Ocean Model (POM) uses  (so do most other ocean models), hence it uses equation
(G4.4.7). In fact an approximate form is used, in which QH = 0, the k is replaced by the
Smagorinsky’s horizontal diffusivity coefficient Asmag (why??), and Frad represents the radiative
penetration of light from the ocean’s surface and is significant only near the surface z > 30 m (where
  T):

       D/Dt = .(Asmag )  (Frad)/z.(Cp)-1,                                           (N4.19)

where .(Asmag ) is supposed to represent mixing effects by small-scale eddies; we replace T by ,
and also lump the (Cp)-1 into Asmag (so that it has a unit of m2s-1) thus acknowledging that we don’t
know much about this small-scale eddy mixing process, and therefore the errors due to these
“replacement” and “lumping” of terms are inconsequential. In actual numerical simulation, we set
Asmag to be very small, in fact often zero.

24. Momentum Balance (or the Equation of Motion):

This is (G4.5.20) [or pom_v13.pdf (or *.doc) equations (1-24a,b,c)]:
                                                                                                   Page 97

        (Du/Dt    + 2 × u) =  p’  ’g + .(        u),                                    (G4.5.20)


        u = (u, v, w);     = (1, 2, 3);       g = (0, 0, g).                             (N4.20)

25. Mechanical Energy (Kinetic Energy) Equation:

This is (G4.6.1) [with (G4.6.2)]; it is derived by taking the scalar product of “u” with eqn.(G4.5.20) 
i.e. u.(G4.5.20). Term by term, we have:

        u.Du/Dt = u.(u/t + uu/x + vu/y + wu/z) = D(u.u/2)/Dt        [decompose into components]

        u.(2 × u) = 0 [the 2 vectors are perpendicular, so their dot product = 0]

        u. p’ =  .(u p’) + p’ .u [by chain rule]

        u. (’g) = wg’

        u.[ .(    u)]    u.[ .( u)] =     .(u. u)  ( u . u)
                                         =    .( (u.u/2))   [(u/x)2+(u/y)2+(u/z)2]

Putting all these together, we have then for u.(G4.5.20):

D(u2/2)/Dt = wg’  .(u p’) +       .( (u2/2))   [(u/x)2+(u/y)2+(u/z)2] + p’ .u


        D(u2/2)/Dt = wg’ + .(u p’ +        (u2/2))   + p’ .u                           (G4.6.1)


         =    [(u/x)2+(u/y)2+(u/z)2] = dissipation rate                                (G4.6.2)

26. It is also useful to think of the rate of change of K.E. for a fixed volume – Gill’s equations (G4.6.3-
G4.6.4). Also understand (G4.6.5) and interpretations based on Fig.G4.6.
                                                                                                    Page 98

27. For a large volume, we have eqn.(G4.6.7):

       dK/dt + Fn’dS =   dxdydz,                                               (G4.6.7)

where Fn’ = n.F’ for F’ = energy flux density given by (G4.6.4), and n is the unit vector normal to the
bounding surface dS.

28. Importance of Viscosity!

We normally say (or hear the statement) that viscosity is not important on the large scales (say a few
kilometers or larger) typical of the energy input (say by wind). But let’s think about this more
carefully. Let’s take the simple case of a homogeneous ocean. Input of energy is due to the pressure
and viscous stress (see equation G4.6.4) at the free-surface only, since no energy can pass through the
solid bottom and side walls. Yet despite this constant energy input, we know that the ocean’s energy
is not continually increasing, i.e. dK/dt  0 in (G4.6.7). Therefore, in that equation, the large-scale
energy flux input Fn’dS can only be balanced by the small-scale energy dissipation  dxdydz by
molecular viscosity! How small is small? Let’s call this small scale = l. The l can be assumed to
depend only on and :

       l = ( 3/ )1/4  O(10-3 m)                                                              (N4.22)

Therefore, because of the very disparate scales, the only way that the large-scale energy can be
dissipated by the small-scale (l ) viscous dissipation is that there will be a transfer of energy from large
to small. This transfer is accomplished by the nonlinear terms (e.g. uu/x etc; if u ~ eikx, then uu/x
~ ei2kx etc) – “breaking” down large eddies into smaller ones, which are in turn broken down again into
even smaller ones etc until the smallest eddies have scales of O(l ) when molecular viscosity can act to
dissipate the (kinetic) energy into heat.

The fact that we have disparate scales creates big problem for numerical models which typically must
use some kind of eddy viscosity >> molecular viscosity in order to dissipate input of large-scale (say
wind) energy. However, this “eddy viscosity” parameterization does not guarantee that we are
removing the energy in a realistic way – thus we have (for example) Mellor-Yamada scheme or the
Smagorinsky viscosity to try to model this “energy-removing” process. If we have a super-super
computer that allow a grid resolution of 10-3 m, then we will not need to use Mellor-Yamada or
                                                                                             Page 99

Chapter 5. Centrifugal and Slantwise Instabilities: Applications to Mixed Layer

Recall that when warm fluid “A” is below cold fluid “B,” the system is unstable if the
potential temperature of “A” is higher than that of “B.” In this case, the system is
unstable because of gravity – i.e. lighter fluid tends to rise and heavier fluid tends to sink.
This is called convective instability. This topic was discussed when we studied Gill’s
Chapter 4. An interesting analogous instability occurs in rotating flows. Recall that on
a rotating earth, centrifugal acceleration may be “lumped” together with gravity such that
a gravitational potential may be defined (see Fig.1-7). This suggests that in rotating
flow, centrifugal force acts like gravity and there may exist a similar kind of instability.
In section 5-1 we study the nature of this instability for a homogeneous fluid (constant
density) – centrifugal instability. We then extend the ideas to the case with density
stratification in section 5-2  slantwise or symmetric instability (Emanuel, 1994). In
section 5-3, we then utilize the same parcel theory used in section 5-1 to re-derive the
growth rates of both instabilities, and also to outline same for baroclinic instability.
Finally, in section 5-4, we discuss the applications of the theories to instabilities that can
happen in a mixed layer where horizontal stratification can exist and the vertical
stratification is weak.

5-1 Centrifugal Instability (CI)

We examine axis-symmetric flow – a circular eddy for example – in which there exists an
azimuthal velocity v(r,t), positive cyclonic and a radial velocity u(r,t), positive outward
from the center. Here r is the radius from the axis of rotation and t is the time. All
variables are independent of the azimuthal angle  (Fig.5.1).
                                                                                            Page 100

Figure 5-1. Axis-symmetrically rotating homogeneous fluid with angular velocity  =
v/r. The radial displacement of a circular ring is also sketched.

The motion is characterized by  and M:

        = angular velocity = v/r                                                 (5-1.1)

       M = angular momentum = rv = r2                                            (5-1.2)

The  here should not be confused with the same symbol used for earth’s rotation rate in
the previous chapters. Assuming inviscid fluid, the azimuthal momentum equation is:

       Dv/Dt + uv/r = o (p/)/r                                               (5-1.3)

and the radial momentum equation is:

       Du/Dt  v2/r = o (p/r)                                                 (5-1.4)

where o = 1/o (not to be confused with the thermal expansion coefficient used in
Chapter 4,  = -1(/T)p,s). [If you are familiar with the y and x-momentum
equations involving the Coriolis parameter “f,” then it is easy to remember (5-1.3,4), just
replace the “f” by “” and use equation (5-1.1); also replace “y” by “r” and “x” by
“r.”] We now show that, following a fluid parcel, the angular momentum “M” is
conserved, i.e. DM/Dt = 0, for axis-symmetric flow. To show this, we multiply (5-1.3)
by r and expand the D/Dt written in (r,) coordinate [i.e. D/Dt  /t + u/r + (v/r)/]:

       rv/t + ruv/r + vv/ + uv = o (p/)                               (5-1.5)

We then note that DM/Dt = D(vr)/Dt is:

       DM/Dt = M/t + uM/r + (v/r)M/
            = rv/t + uv + ruv/r + vv/                                      (5-1.6)

which is the same as the LHS of (5-1.5). Therefore, (5-1.5) gives :

       DM/Dt = o (p/) = 0,                                                   (5-1.7)

for axis-symmetric flow, since then “p” (also all other variables) is independent of .
                                                                                            Page 101

When we analyzed convective instability (Chapter 4), we examined a parcel’s movement
in the vertical (z) within a background of (ambient) density stratification. We then
moved the parcel vertically and adiabatically to determine if at the new position the
parcel is lighter or heavier than the ambient fluid. A similar analysis is now carried out
for a parcel undergoing axially symmetric displacements. Instead of “z” we now have
“r,” and instead of potential temperature, we now have that the angular momentum (M) is
conserved. The parcel in the present axially symmetric situation is a “ring.”

So, what is the ambient, equivalent “hydrostatic” state? This is obtained from (5-1.4) by
setting Du/Dt = 0 (similar to setting “Dw/Dt = 0” in hydrostatic case):

       Mo2/r3 = o (dpo/dr)                                                     (5-1.8)

where Mo = rvo, and po is the pressure distribution that maintains the balanced azimuthal
flow having an angular momentum = Mo. Suppose now that the flow is initially at this
balanced state, and consider the ring-parcel radially displaced by r from r = ro. For
small displacement, we assume also that “p” remains at its initial pressure distribution po.
Equations (5-1.4) and (5-1.8) then give, for the ring-parcel:

       Du/Dt = [M2  Mo2(ro+r)]/r3                                               (5-1.9)

Here, u and M represent the fluid state after the displacement at r = ro+r.    Since M is
conserved (i.e. equation 5-1.7), we have:

       M2 = Mo2(ro)                                                               (5-1.10)

since the ring-parcel initially (i.e. before being displaced) has M = Mo.      By Taylor’s
expansion about ro assuming that r is small, we also have,

       Mo2(ro)  Mo2(ro+r)  [(dMo2/dr)|ro].r                                   (5-1.11)

where the vertical-bar symbol “|ro” means “at ro.” Equations (5-1.10) and (5-1.11) then
give M2Mo2(ro+r)  [(dMo2/dr)|ro], and equation (5-1.9) then becomes:

       [Du/Dt]|(ro+r)  [(dMo2/dr)|ro].r/ro3                                   (5-1.12a)

or simply (with the “at” symbol “|” being understood):
                                                                                              Page 102

        Du/Dt  (dMo2/dr).r/ro3                                                       (5-1.12b)

Equation (5-1.12) is the result we want. It shows that as the ring-parcel is displaced
outward, r > 0, it is further accelerated outward (and therefore moves further away from
its initial position ro), Du/Dt > 0, if the squared ambient angular momentum Mo2
decreases with r, i.e. if (dMo2/dr) < 0. Thus,

        if (dMo2/dr) < 0, then the flow is unstable                                     (5-1.13a)

This kind of instability is called centrifugal instability (CI).   On the other hand,

        if (dMo2/dr) > 0, then the flow is stable                                       (5-1.13b)

Suppose the ambient aximuthal velocity varies with “r” as:

        vo = rn                                                                         (5-1.14)

then Mo2 = r2(n+1), and

        dMo2/dr = 2(n+1).r2n+1 < 0 if n < 1, (for vo = rn)                             (5-1.15)

Therefore, an (axially symmetric) eddy which “does not spin fast enough,” i.e. whose
velocity distribution decays faster than 1/r, is unstable to axially symmetric perturbations.


Go to a playground with your friend(s) or someone you may trust. Stand on a “merry-
go-round (MGR; see fig.5-2).” Ask the trusted friend to spin the MGR, then jump off
while the MRG has spun to an almost steady rotational speed. Describe and explain
what you experience the moment you touch the firm (almost stationary) ground. By the
way, why is the ground “almost stationary”? (Warning: this exercise may not be
suitable for teachers at or over 55).
                                                                                          Page 103

Figure 5-2.   A merry-go-round.

Localized Centrifugal Instability (LCI):

The above theory requires axially symmetric perturbations. But as the Merry-Go-
Round exercise suggests, localized perturbation (the exerciser jumping off) may also
induce localized centrifugal instability (LCI). This is in fact true for ocean (and
atmosphere). To develop the relevant formulae, we now define localized conditions for
CI based on the inviscid equations of motion on an f-plane. We will develop the theory
carefully, step by step.

In equation (1-24), we

(a) ignore the horizontal components of the earth’s rotation: 1 and 2;

(b) align the z-axis so that it is perpendicular to the plane tangent to the local earth’s
    surface where the latitude is o, i.e. so that z is locally “upward” and 23 = 2sin(o)
    = f, where  = magnitude of the earth’s angular velocity;

(c) assume motions of scales small enough that “f” does not vary significantly, i.e. the
    motions are confined to the neighborhood of latitude “o”; and,

(d) ignore the viscosity .
                                                                                         Page 104

We then divide (1-24) through by the density , and set  = 1/; so that equation (1-24)

                                                                                 (5-1.16a)

                                                                                 (5-1.16b)

                                                                                 (5-1.16c)

Let us assume that in our flow field, there exists a direction yg such that p/yg = fug =
constant (could be zero!), we can then align the “y” in (5-1.16) in this direction, i.e.:

       p/y = fug = constant                                                   (5-1.17)

where ug is therefore the (constant) geostrophic flow in the (new) x-direction. Equation
(5-1.16) then becomes:

                                                                                 (5-1.18a)

                         )                                                        (5-1.18b)

                                                                                 (5-1.18c)

We can now translate the coordinate system in the x-direction at the constant speed ug, so
that with respect to this translating coordinate system the new “u” and its material
derivative are given by:

       unew = u  ug, and      Dunew/Dt = Du/Dt                                   (5-1.19)

The other velocity components, v and w, in the directions perpendicular to the x-axis, are
clearly not affected by the uniform x-translation (assuming negligible Einstein’s
relativistic mechanics!). Also xnew = x  ug.t, so that /xnew = /x. Equation (5-1.18)
is then (dropping the “new” etc):

                                                                                 (5-1.20a)
                                                                                             Page 105


                                                                                    (5-1.20c)

In this new aligned and translating coordinate system, we note that the y-momentum
equation (5-1.20b) has no pressure forces, so it becomes:

        D(v + f x)/Dt = DM/Dt = 0                                                    (5-1.21)


        M=v+fx                                                                       (5-1.22)

is the absolute momentum. Therefore the absolute momentum M is conserved for
inviscid motions on an f-plane on which one component of the geostrophic velocity (ug in
our case) is constant. While ug is constant, the geostrophic velocity in the y-direciton,
i.e. vg is quite free, and in general is a function of x and z, and is defined as (see fig.5-3):

        vg =  (p/x)/f                                                             (5-1.23)
                                                                                          Page 106

Figure 5-3. A sketch showing ug and vg, and also the tube-parcel elongated in the y-
direction and given a perturbed displacement x in the x-direction.

We consider then the flow situation in Fig.5-3. Suppose a perturbation is introduced.
In order that M (= v + f x) remains conserved, the perturbation must be such that the ug
remains constant. From (5-1.17), it follows then that

        p’/y = 0                                                                  (5-1.24a)

where p’ is the perturbation pressure.    It then follows that

       (perturbation)/y = 0                                                       (5-1.24b)

and the perturbation must be highly (in fact infinitely) elongated in the y-direction as the
tube-parcel of fluid sketched in Fig.5-3 shows. When the tube-parcel is displaced in the
x-direction, the governing acceleration is obtained from (5-1.20a) which upon using (5-
1.23) becomes:

                                                                                   (5-1.25)

where the approximation is from neglecting the perturbation pressure from p/x (i.e. the
pressure remains at its initial pressure distribution that produces vg) and

       Mg = vg + f x                                                                (5-1.26)

is called the geostrophic absolute momentum (in analogous to the absolute momentum in
equation 5-1.22). In (5-1.25), the u, M and Mg are the values at the tube-parcel’s
displaced position x + x. Then comparing (5-1.25) with (5-1.9) and following the
analysis that leads to (5-1.12b), we now have:

       Du/Dt  (Mg/x).x.f                                                       (5-1.27)

which is the result we want. It shows that as the tube-parcel is displaced to the right, x
> 0, it is further accelerated to the right, i.e. Du/Dt > 0, if (Mg/x) < 0. Thus,

       if (Mg/x) < 0,     i.e. if   vg/x + f < 0,   then the flow is unstable   (5-1.28)
                                                                                          Page 107

The growth rate may be estimated by noting that u  x/t, so that (5-1.27) then gives:

       cf2 = (t)-2  f.( vg/x + f)                                          (5-1.29)

Equation (5-1.28) shows that only anticyclones can have centrifugal instability,
especially at lower latitudes.
                                                                                         Page 108

Gill's Chapter 5: Adjutment under Gravity in a Nonrotating System

Section 5.1:

Many fluid mechanical phenomena (including GFD) involve studying how the fluid
system (e.g. the atmosphere-ocean system) tends to adjust to equilibrium; i.e. how it
approaches a steady state once it is initially perturbed. What are the characteristics of
the equilibrium (e.g. a bucket of water or a small lake at rest; the Kuroshio or the
atmospheric jet-stream in perfect geostrophic balance; etc)? How long does the
adjustment take? How may the adjustment process be understood?

In this chapter, we study the adjustment without the complications due to the earth's
rotation and curvature, and furthermore we only consider small departures from the
hydrostatic state - i.e. we will linearize the equations.

Fig.5.1 -- left (side X) has density 2, right (side Z) has density 1 < 2. The response is
similar to the familiar "estuarine circulation;"

Fig.5.2a:   pA  pB = c gh;                                                    (N5.1)

Fig.5.2b: Let "F" be the z-position of the interface that is on the left side of the
partition. Then:

       pA = pF + 2 gh, and
       pB = pF + 1 gh; so that:
       pA  pB = (2  1) gh = 2 g’h                                          (N5.2)

where g' = (2  1) g/2 is the reduced gravity.

Exercise N5.1. Why is the adjustment slower in the case of Fig.5.2b when compared with
that of Fig.5.2a?

Section 5.2:

Consider an ocean of constant depth H and constant density c with an atmosphere of
zero density above, so that the in-situ density of the system is:

        = c for H < z  0
                                                                                             Page 109

        = 0 for z > 0.

A Cartesian coordinate (x, y, z) is chosen such that z = 0 is at the free surface at rest, x
points eastward and y northward. The z-momentum equation applied to this equilibrium
(i.e. at rest) fluid is then:

       dpo/dz = g, or po(z) = gz.                                              (5.2.1)

Section 5.3 -- Surface Waves:

The governing equation is just the Laplace Equation (5.2.7). Note that this is NOT a
wave equation – it is not even time-dependent! It is in fact just a statement of
incompressibility condition (5.2.4).

The “wave” part comes from the surface boundary condition equation (5.2.10) which
gives the time-dependent motion to the water beneath to generate waves.

The bottom boundary condition (5.2.8) restricts wave motion at the ocean’s floor. But
actually, this condition is less important physically especially when the ocean is very
deep – why?

Also, why are there no boundary conditions on the sides?         What implicit assumption is
Gill making by ignoring side boundaries?

To see this, derive the dispersion relation (5.3.8) in section 5.3.

       2 = g.tanh(H).                                                           (5.3.8)

Note that the phase speed is c = / = [(g/)(tanhH)]1/2, so that long waves with small 
travels faster than short waves with large .

For H << 1, we have tanh(H) ~ H, so that:

       c = /  (gH)1/2                                                         (N5.3)

For H >> 1, we have tanh(H) ~ 1, so that:

       c = /  (g/)1/2.                                                       (N5.4)
                                                                                       Page 110

What is the meaning of H << 1 and H >> 1?

Derive the expressions for u and v:

Equation (5.2.5) and (5.3.6) give:

        u/t = gkosin(kx+lyt)cosh[(z+H)]/cosh(H),

and similarly for v; so that:

        (u,v) = (k, l).(go/)cos(kx+lyt)cosh[(z+H)]/cosh(H)             (N5.5)

Let’s consider only waves propagating in the xz-plane (i.e. the y-wavenumber l = 0).
Also consider infinitely deep ocean, i.e. H ~ . Then (N5.5) becomes:

        u  (kgo/).cos(kxt).exp(kz)                                      (N5.6a)

and (5.3.7) gives:

        w  (kgo/).sin(kxt).exp(kz)                                      (N5.6b)

Exercise: Plot on the xz-plane paths experienced by fluid parcels due to the wave motion.
Hint: x ~ u.dt and z ~ w.dt.

Some Useful Things to Remember for Deep-Water Waves:

Condition for “deep water” is: H >> 1, or

        H >> /(2)                                                          (N5.7)

i.e. water depth deeper than about one wavelength . For example, at z = /2, equation
(N5.6) gives the amplitudes for u and w to be approximately exp(2)  0.04, or 4% of
the amplitude at the surface, and at z = , exp(2)  0.002, which is only 0.2% of the
amplitude at the surface.

From (N5.4), we have for deep-water wave phase speed in m/s:
                                                                                     Page 111

       c = /  (g/)1/2 = [g/(2)]1/2 1/2  1.25 1/2                   (N5.8a)

and also for period tp in seconds:

       tp = /c  0.80 1/2                                                (N5.8b)

The Table below summarizes the above results.

                      Table 5.1 Typical values of deep-water waves
   (m)               c (m/s)            tp (s)         Water deeper than H (m) ~ /2
     1                 1.25               0.8                        0.5
    10                    4               2.5                         5
   100                 12.5                 8                        50
  1000                   40                25                       500
                                                                                       Page 112

Gill's Chapter 6:

It is perhaps best to approach the concept of normal modes in a continuous system, then
we will go back to the two layer model. I will rely on your knowledge about the Fourier
theory. Mr. Fourier basically asserted that a continuous (maybe piecewise) function can
be represented by a simple sum of infinite number of terms comprising of the sines and


where the an and bn are the Fourier coefficients. The sines and cosines are called the
basis functions; together they form a complete basis – which loosely means that together
they are sufficient to represent F(x). The sines and cosines are solutions of many
boundary/initial-value differential equations – simple harmonic oscillator, seiches in
channel etc. We will find that normal modes for motions in a stratified sea of constant
depth play the same role as the sines and cosines in these other problems.

Consider the following linearized equations:

       u/t  fv = P/x + x/z                                           (6.2a)
       v/t + fu = P/y + y/z                                           (6.2b)
       (p/z) /o = g/o = +b                                              (6.2c)
       /t + wdo/z = 0, or b/t + N2w = 0                                (6.2d)
       u/x + v/y + w/z = 0                                              (6.2e)

where P = p/o is the “kinematic” pressure, x and y are also kinematic stresses, i.e.
stresses divided by o, N2 = (g/o)(do/dz) is the buoyancy frequency, o(z) =
background density, and (u, v, w, p, , b) are all small perturbations representing
departures from a static fluid at rest.

We consider that the motion is initiated from rest, so that

       (u, v, w, p, , b) = 0,   at t=0                                       (6.3)

Equation (6.2d) then suggests that we define

                                                                                             Page 113

so that b = N2 , and (6.2c) becomes

       (p/z) /o = N2                                                            (6.5)

Now, the total pressure PT is the sum of a hydrostatic part and the perturbation:

       PT = gz + P                                                                 (6.6)

At the free surface z = (x,y,0,t) = (x,y,t), PT is the atmospherics pressure, taken to be
zero. So that (6.6) gives,

       P = g, and P/t = g/t,     at z  0                                     (6.7)

Integrating (6.5) w.r.t. z from z to 0 then gives:

                                                                                   (6.8)

In addition to the boundary condition at z=0, we also have at the bottom:

       w=0        at   z = H (H=constant)                                          (6.9)

Normal Modes:

We seek separable solutions of the homogeneous part of (6.2) (i.e. without forcing ):

       (u, v) = (Uk(x,y,t), Vk(x,y,t)).dQk/dz,       k=1,2,3,…                      (6.10a,b)

where “Q(z)” is a function of z only. From (6.2e), we have:

       (Uk/x + Vk/y).dQk/dz = w/z                                            (6.11)

which suggests that the “z” part of w be “Q”, or since w is related to through (6.4), that

         = Sk(x,y,t).Qk(z)                                                          (6.10c)
                                                                                           Page 114

Substitute (6.10) into /z of (6.2a) (for example; with  = 0) gives, using also P/z
from (6.8):

       LHS of (6.2a)  [Uk/t  fVk].d2Qk/dz2                                    (6.12a)

       RHS of (6.2a)  [Sk/x].N2Qk(z)                                           (6.12b)

For both sides to be consistent, the “z” part must match to an arbitrary constant factor,
say “k”, so that:

       d2Qk/dz2 + N2k.Qk(z) = 0                                                  (6.13)

This is a simple ODE, and is supplied with 2 BC’s one at z=0 (i.e. 6.7) and the other one
at z=H (i.e. 6.9). A particularly simple case is, instead of (6.7), we also take w = 0 at z
= 0:

       Qk = 0 at z = 0 and at z = H                                              (6.14)

The solution of (6.13) consists of sines and cosines, which are the normal modes or
eigenvectors of the linearized, shallow-water equations (6.2) in an ocean of constant

Substituting (6.10) into (6.2) [/z the momentum eqns. first], and making use of (6.13):

       Uk/t  fVk = ck2 Sk/x                                                 (6.15a)
       Vk/t + fUk = ck2 Sk/y                                                 (6.15b)
       Uk/x + Vk/x = Sk/t                                                  (6.15c)

where we have written ck2 = 1/k. These equations are exactly the same as the
equations the depth-integrated equations of a homogeneous fluid in an ocean basin of
constant depth. In fact, we can identify Sk to be the “” and

       ck2 = gHk                                                                  (6.16)

which can be named the square of the equivalent shallow-water wave phase speed for the
mode “k” in an ocean of equivalent depth Hk. Thus the solution to this stratified fluid
problem consists of just solving a set of shallow-water equations, one for each “k”, then
summing them all up to get the complete solution.
                                                                                            Page 115

In Fourier series, we calculate the coefficients an and bn by using the orthogonality
property of of sine and cosine:

                                      
                                       
                                                                                 ; (6.17a)

                                                  ;                                 (6.17b)

                           
                            
                                                            .                      (6.17c)

The first equality on the LHS of (6.17a) is derived by first noting that     

then using ei(m+n)x = [cos(mx).cos(nx)sin(mx).sin(nx)] + i[sin(mx).cos(nx)+
cos(mx).sin(nx)] to equate the real part. The second equality, i.e. that the integral is = 0,
is shown by integration by parts, say of                 
                                                                             , which gives

                           
                            
                                                  , so that in order to be consistent with the

first equality we must have, if m  n, each of the two integrals = 0.
Equating the imaginary part then gives (6.17b). Equation (6.17c) is easily derived.
                                                                     
We first have                                                         

integration by parts. Also, since cos(m+m)x = cos2mx  sin2mx = 1  2sin2mx, we have
                                                                     
                            
                                                         ; and also    

Using equations (6.17), the Fourier coefficients an and bn can be determined for the
function F(x) in (6.1). The requirement is that the basis functions (sines and cosines) be
complete and orthogonal in the interval [-, ].

An exactly parallel approach can be used for the normal modes. We now show that they
are orthogonal. Multiply (6.13) by Qm, take the integral from H to 0, integrate by parts:

       Qmd(dQk/dz) + N2kQmQkdz = 0, or,

       (dQm/dz).(dQk/dz)dz + N2kQmQkdz = 0.
                                                                                          Page 116

Interchange “m” and “k” in the above equation and subtracting the two resulting

       (km). N2 QmQkdz = 0, or

       N2 QmQkdz = 0, if k  m.                                               (6.18a)

so that from the line before “Interchange “m” and “k”..”, we also have

       (dQm/dz).(dQk/dz)dz = 0, if k  m.                                     (6.18b)

Comparing this with (6.17), we see that the normal-mode eigenfunctions Qm (or dQm/dz)
are orthogonal in the interval z = [H, 0]. Thus any (piecewise continuous) function in
“z” may be expressed as a sum of Qm’s with coefficients determined similar to the way
the Fourier coefficients are determined. In particular, the frictional force, x/z and
y/z may be expressed in terms of the dQk/dz’s:

       x/z = g k        k
                                . dQk/dz,

       y/z = g k        k
                                . dQk/dz                                         (6.19)

where the kx,y can be functions of (x,y,t). Given the functions x/z and y/z, the
 k     are obtained just like we obtained the Fourier coefficients an and bn:

        k         = x,y/z. dQk/dz dz/(g( dQk/dz )2dz).                      (6.20)

Thus instead of (6.15), we now have:

       Uk/t  fVk = ck2 Sk/x + g         k
       Vk/t + fUk = ck2 Sk/y + g         k
       Uk/x + Vk/x = Sk/t                                                 (6.21c)

This completes the formulation of the linearized shallow-water equations (6.2) in water
of constant depth H in terms the normal modes. Instead of the velocity (u,v) and density,
we have transformed the problem into one that requires the solution of the “transports”
(Uk, Vk) and “elevation” Sk for each of the mode “k.” We will now show that in a two-
layer system that we studied in Gill’s chapter 6, we use similar “transport” and “elevation”
                                                                                           Page 117

Chapter 7: Notes on Island Rule

Island rule (Godfrey, 1989) is derived by an application of the Kelvin’s circulation
theorem. We will first learn this theorem, following Pedlosky’s book (1979), then we
will learn island rule by choosing a very special circuit.

Kelvin’s Theorem

Take a curve C enclosing a surface A in the fluid. Then the line integral of the velocity
u around C is:

                                                                                 (7.1)

Note that dr is on C. When C is a material circuit, meaning that it is composed of the
same fluid elements at all time and therefore it moves with the fluid. Then we can d/dt
eqn.(7.1) and bring the d/dt inside the integral:

                                                   =                              (7.2)

Since u.du = d|u2|/2, the integral                      .   Now we substitute the momentum

       du/dt + 2  u =  p/ + g + F                                             (7.3)

into (7.2) which then gives:

                                                                                (7.4)

Meaning of Each Term:

1.                       --- consider a u that is perpendicular to .           (7.5a)

2.                                                                                (7.5b)
                                                           

3.              --- this is the frictional term.                                  (7.5c)
                                                                                            Page 118

In application to the Island Rule, the fluid is single-layer and homogeneous, so the
friction F is given by:

       F = (o  b)/H,                                                            (7.6)

where H is water depth and o and b are wind and bottom stresses respectively. The
velocity u is then a depth-averaged velocity (u, v). Since the motion is barotropic, the p
and  surfaces coincide, then the      p  = 0, and the contribution from                  in

(7.5b) vanishes. The x and y momentum equations can be written as:

       u/t + (+f) k  u =  p/ + |u2|/2) + (o  b)/H.                       (7.7)

where k is the unit upward vector and  = v/x  u/y is the z-component of the relative
vorticity. This equation looks somewhat different from the usual form. But by
expanding  (k  u) and (|u2|/2):

        (k  u) =  (iv + ju) = i v(v/xu/y) + j u(v/xu/y)             (7.8a)


           (|u2|/2) = i [(u2+v2)/2]/x + j [(u2+v2)/2]/y
                   = i [uu/x + vv/x] + j [uu/y + vv/y]                     (7.8b)

and adding them (and including the f also), we have:

 (+f) (k  u) + (|u2|/2) = i [uu/x + vu/y  fv] + j [uv/x + vv/y + fu]    (7.9)

which (in (7.7)) gives the usual form.

Integrate (7.7) now around a circuit CI that surrounds an island (as in 7.2 and 7.4):

                                      , on circuit CI.                          (7.10)
                                                                                            Page 119

The circuit-integrals of other terms, (+f) k  u and p/ + |u2|/2), vanish, first because
the normal component of the velocity to the island is zero, and secondly because the
circuit integral of (…) = 0.

Consider that the island is Taiwan at longitude xI and latitudes ys to yn; the island is
located between China at xw and America at xe. Then consider a new circuit CR in the
Pacific Ocean between Taiwan and America: (xe, ys)  (xe, yn)  (xI, yn)  (xI, ys)  (xe,
ys). The Bernoulli’s terms again integrate to zero, but the vortex terms are not zero
along the two latitude lines (xe, yn)  (xI, yn) and (xI, ys)  (xe, ys):

                                                                      , on CR,   (7.11)

where n = unit normal to CR, and the circuit integral is positive cyclonic (hence the minus
sign in front of fn). Now,


is just the Sverdrup transport + Kuroshio transport and is in fact the stream function on
the island, and is equal to the transport between China and Taiwan = Taiwan Strait
transport. To see this, consider the definition of /x = v, and integrate to calculate for
 from China to Taiwan then to America – try it.

Adding the circuit integrals on CI and CR (7.10) and (7.11), we then have:

                                                , on circuit CI + CR.          (7.13)

which in steady state and assuming that b  0:

                                   , on circuit CI + CR.                         (7.14)

which is the Island Rule.
                                                                                                Page 120

                            3. Homogenization of 2 Fluid Layers by Wind

                                    L.-Y. Oey (lyo@princeton.edu)

We wish to derive Eqn.2 of Oey et al. 2006 (Loop Current warming by hurricane Wilma, GRL 33,
L08613, doi:10.1029/2006GL025873, 2006). The followings have variously been derived by several
people: J. Turner (original?), John Simpson (also original?), Chris Garrett, Peter Rhines and perhaps

Consider two fluid layers that are mixed into one layer by turbulence due to, e.g. wind, as shown
schematically by the following figure:

Figure 3-1 Left: two layers of different densities  and  -  ( > 0) before the onset of mixing
(by wind). Right: wind mixes the fluids and homogenizes the two layers to one layer of some
intermediate density = am.

The physical situation may be that in summer, warm fluid overlies cold fluid (Figure 3-1 left).
A tropical cyclone (typhoon or hurricane) comes and mixes the two fluids into one layer (Figure
3-1 right). In this idealization, we assume that there is little agitation of the fluid below the
lower layer of Figure 3-1 (left).
                                                                                                   Page 121

What is am, the final density of the fluid after mixing, in terms of the original densities of the
two fluids before mixing? If h1 = h2 then clearly am will be a simple average of  and ,
i.e. am =   /2. In general, am is a weighted-average of  and . In other words, by
mass conservation (i.e. mass after mixing = mass before mixing), we have per unit horizontal
(xy) area:

               am (h1 + h2) = .h2 + ().h1,                                          (3-1a)

Which shows that am is a weighted-average of  and :

               am = [.h2 + ().h1]/[h1 + h2]                                         (3-1b)

Wind pumps energy into the fluid by doing work through turbulent mixing, and thereby
increases the potential energy of the fluid. The rate (per unit area) at which the wind does work
is simply equal to the frictional stress (the drag force per unit area acted by the ocean surface on
the air above) = aCd|W|2 times the air parcel’s speed (assume that the wind speed >> ocean
current speed), thus:

               (Wind Work)/Area = aCd|W|3 in J/(m2.s) or Watts/m2                        (3-2)

where J = Joules (= Newton  meters). If the wind blows for a time period = , then the energy
(per unit area) supplied by the wind to the water column, or the power dissipation by the wind,

                                         , in J/m2.                                     (3-3)

Here, the factor  accounts for the fact that not all the wind work is used to produce mixing in
the water column – some of it is used to create waves, sprays etc. In fact, only a relatively
small fraction ~ 10% (or less) is used to produce mixing in the fluid below. However, for our
purpose (see below), we need not worry about this.

To see that the wind work raises the potential energy (PE; again per unit area) of the fluid, we
need to compare the PE after mixing, PEAM, with the PE before mixing, PEBM. The PE is
simply the work done in raising the fluid through the water column, i.e. ~ the integral with
respect to z of the weight of the fluid parcel as it is raised through the water column:

                                                                                     (3-4)
                                                                                                    Page 122

where the primes denote general or dummy variables. Evaluating the last two integrals in (3-4),
we obtain for the potential energy of the fluid, per unit area, before mixing:

               PEBM = (g/2)[h22 + 2h1h2 + (  )h12                                    (3-5)


                                                                    
                                                                                           (3-6)

where equation (3-1a) has been used to substitute for “am (h1 + h2)” in the final expression.

The change in PE, PE = PEAM  PEBM, is then given by:

                                                                         
                                                                                   


                                     , in J/m2,                                          (3-7)

which is positive, confirming that the PE of the mixed fluid is indeed increased.

The ratio,

                                                                                       (3-8)

which is equation (2) of Oey et al. (2006), is then a ratio between the amount of energy required
to thoroughly mix two stratified fluid layers and the amount of energy from wind that is actually
available to do the mixing. If is “large,” then either the wind power is not sufficiently strong
(or the wind has not acted long enough) to thoroughly mix the fluid layers, or that the
stratification between the two layers is too strong, or the layer thicknesses are too large for the
given wind to mix the layers. The opposite is true if          is “small,” wind is very strong, or
stratification is weak, or layers are thin.
                                                                                                 Page 123

It is difficult to exactly say just what value of is to qualify it as “large” or “small.” (One of
the difficulties is of course our very poor knowledge of the value of .) However, one can
estimate , , h1 and h2, and calculate (and plot) the for particular situation, say as a time
series, and attempt to correlate (visually is the simplest!) it with, say, the time-series of sea-
surface temperature or better still the temperature a few meters below the surface. If you
suspect a correlation, you can then try to estimate what the ‘critical”       is for your particular

Exercise: Carry out the integrations in (3-4) to obtain (3-5).
                                                                                                  Page 124

Gill’s Chapter 7: Effects of Rotation
The principle of conservation of angular momentum
Why do global winds display bands of easterlies and westerlies (see Gill’s fig.2.2)?

To answer, we use the principle of conservation of angular momentum which is valid when there is no
friction. Let τx(φ) = eastward wind stress = force/area acting on the earth’s surface. This stress
transfers a tourque = rate of transfer of angular momentum per unit area of the earth’s surface:

       torque = rE.cos(φ).τx(φ)

Multiply this by an elemental area of zonal strip = 2πrE.cos(φ).rE.dφ, inegrate from south pole to north
pole to cover the entire earth’s surface, then set the integral to zero because there can be no net torque
(otherwise the earth’s spin may vary with time and it will be hard for us to plan a daily schedule):

or,                                   .


Therefore, since cos2(φ) ≥ 0, τx(φ) must change sign, i.e. must alternate between westerly and easterly
(figure). The principle of conservation of angular momentum is used over and over again in GFD.
                                                                                                Page 125

7.2 The Rossby Adjustment Problem

Previously in chapter 5 we studied the outward propagation of an initial discontinuity in free surface
when there was no rotation. By method of characteristics we were able to calculate that the
discontinuity spreads at a propagation speed = (gH)1/2. We now examine what happens to the
discontinuity if there is rotation. In that case, we have eqn.(7.2.4).

Exercise 1a: Using hydrostatic equation ∂p/∂z = −ρg, show that p = patm + po(z) + ρgη, where po(z) =
−ρgz is the equilibrium (static) pressure and patm is the atmospheric pressure. (Often we set patm = 0,
or at least patm = 0. In numerical modeling of storm-generated currents, the patm may be specified
from weather-chart and can be important).

Exercise 1b: State all conditions and assumptions under which equations (7.2.1), (7.2.2) and (7.2.3) are

Exercise 2: Derive eqn. (7.2.4):

       ∂2η/∂t2 – c2(∂2η/∂x2 + ∂2η/∂y2) + fHζ = 0                                           (7.2.4)

For f = 0, this equation can be solved given initial and boundary conditions for η. For f ≠ 0, we need
another equation that relates ζ and η. (Alternatively, of course, one can solve for (u,v,η) from the 3
                                                                                                   Page 126

equations (7.2.1)−(7.2.3); but it is often easier to solve just 1 equation (7.2.4) for η when f=0, or derive
another one when f≠0 to solve for η and ζ).

Exercise 3: Derive eqn. (7.2.8), the conservation of (linearized) potential vorticity:

       ∂(ζ/f – η/H)/∂t = 0                                                                    (7.2.8)

Thus for a homogeneous rotating fluid, the quantity Q’ is invariant with time:

       Q’ = (ζ/H – fη/H2)                                                                     (7.2.9)

Note that Q’ may be and in general is a function of the spatial coordinates (x,y).    But at each point
(x,y), it retains its initial value for ALL time:

       Q’(x,y,t) = Q’(x,y,0).

Exercise 4: The full PV = (f + ζ)/(H + η). Show that in linearized form (i.e. ignoring products such as
(ζη) terms and higher), this becomes f/H + perturbation, where the perturbation is Q’.

Eqns.(7.2.9) and (7.2.10) give:

       fHζ = fH2Q’(x,y,0) + f2η

or after dividing by f2H:

       ζ/f – η/H = HQ’(x,y,0)/f = ηosgn(x)/H

Substituting this into (7.2.4) then gives:

       ∂2η/∂t2 – c2(∂2η/∂x2 + ∂2η/∂y2) + f2η = − fH2Q’(x,y,0)                                 (N7.2)

Consider now our initial condition of a discontinuous free surface:

       η(x,y,0) = −ηo.sgn(x),      and       u = v = 0,                                       (7.2.11)
                                                                                                  Page 127

so that from (7.2.9) we have Q’(x,y,0) = fηosgn(x)/H2 which after substituting into (N7.2) gives:

       ∂2η/∂t2 – c2(∂2η/∂x2 + ∂2η/∂y2) + f2η = − fH2Q’(x,y,0)
                                          = − f2ηosgn(x)                             (7.2.13)

Exercise 5: Carefully follow the above steps, then derive the appropriate equations on the β-plane, i.e.
for the case when f = fo + βy. Then look back notes discussed previously on Sverdrup transport etc,
compare and discuss.

7.2.2 The Steady Solution: Geostrophic Flow

       We wish to know the steady state that the initial discontinuity in elevation (eqn. 7.2.11) may
evolve into after a long time. It seems reasonable that we find the steady solution from the steady
equations – i.e. from equations (7.2.1)−(7.2.3) with all the ∂/∂t terms set to zero − so let’s try that:

       fv = g∂η/∂x, fu = −g∂η/∂y

Exercise 6: Describe what you find after you try it, and are you then able to find the steady-state
solution for (u,v,η)? Why yes (if you could find the solution) or why not (if you could not)?

Exercise 7: Based on the conclusions of Exercise 6, look back pom-book’s chapter 2 and find out what
we did when we encountered a similar dilemma with the Ekman problem.

         In Exercise 6, you should find that the geostrophic relation (7.2.14) gives a velocity field (u,v)
that identically satisfies the (steady-state) continuity equation (7.2.3). Thus any η(x,y) would produce
(u,v) (from 7.2.14) that satisfies all the steady-state equations. In other words, we cannot find a
unique solution.

        In Exercise 7, you should find that for the steady Ekman problem of pom-book’s chapter 2, we
relied on viscosity (i.e. Ekman layer) to overcome the difficulty and found the unique Ekman-layer
solution. In the present case, we do not have viscosity. But recall that in the Taylor-Proudman
column problem we were able to find the solution numerically, and we mentioned then that that was
accomplished by solving the pom-model by a time-dependent integration. So now we can do the
same, but we do it analytically. Actually, we have done so when we integrated (7.2.8) in time to
obtain (7.2.10), and finally we obtained (N7.1b).
                                                                                                   Page 128

Exercise 8: In (N7.1b), approximate the ζ using the geostrophic velocity from (7.2.14), then show that:
(g/f2)(∂2η/∂x2 + ∂2η/∂y2) – η/H = ηosgn(x)/H, so that:

        −c2(∂2η/∂x2 + ∂2η/∂y2) + f2η = −f2ηosgn(x)

which is the steady-state form of (7.2.13).

        Equation (N7.3) governs the form of η, given some conditions on η for large |x|. Since the
RHS (which is the ‘forcing’) is not a function of ‘y’ we may look for a solution η(x), function of ‘x’
only, so that:

       −c2d2η/dx2 + f2η = +f2ηo, for x < 0, and

       −c2d2η/dx2 + f2η = −f2ηo, for x > 0

The solution that is finite at large |x| is then:

        η/ηo =   1.exp[+x/(c/f)]   + 1, for x < 0, and                                           (N.7.5a)

        η/ηo =   2.exp[−x/(c/f)]   − 1, for x > 0.                                         (N.7.5b)

Matching the two solutions at x =0 where the η(0) needs to be continuous, =0, we choose      1   = −1 =
− 2, so that finally:

        η/ηo = −exp(+x/a) + 1, for x < 0
                = +exp(−x/a) − 1, for x > 0,                                               (7.2.22)


     a = c/f = (gH)1/2/f                                                                          (7.2.23)

is the Rossby radius of deformation.       From the geostrophic relation, we obtain:

        u = 0, and v = −(gηo/fa).exp(−|x|/a)                                            (7.2.24)
                                                                                                   Page 129

The flow is not in the direction of the pressure gradient, but is perpendicular to it along line of constant
η (i.e. contours of η), as shown in Gill’s fig.7.1. That geostrophic flow is along η-contours can be
shown by noting that

       u. η = (g/f).[−∂η/∂y.∂η/∂x + ∂η/∂x.∂η/∂y] = 0,

so that the velocity vector u and the direction normal to the contours η are perpendicular; i.e. u and η-
contours are parallel.

An alternative derivation:

I mentioned in the class that the reason we could obtain a unique solution is that we retain the small
∂/∂t terms, which then leads to (7.2.10). We will now do this more systematically by first scaling the
equations. Scaling reveals which terms are small and which are larger. Let U be the velocity scale,
H the scale for η, L the scale for (x,y), and εL/U the scale for time, where ε = U/fL is the Rossby
number. The reason we scale the time “t” in this way will become clear shortly. We then define
nondimensional variables (denoted by subscript n) as follows:

       u = un.U, v = vn.U, x = xn.L, y = yn.L, t = tn/(εf), η = ηn.H                          (N7.6)

By scaling, we anticipate that all non-dimensionalized variables are of O(1). We see that, since ε is
expected to be small, tn ~ O(1) ~ ε.ft, and we therefore require that the time “t” varies very slowly over
many inertial periods ~ 1/(εf). Then equation (7.2.1) becomes (dropping the subscript “n”):

       fUε (∂u/∂t) – fU.v = −(gH/L).∂η/∂x, or

       ε ∂u/∂t – v = −(R2/εL2).∂η/∂x,

       ε ∂u/∂t – v = −(1/ε) ∂η/∂x,

where R = c/f is the Rossby radius (Gill uses “a” for this), and we have chosen L = R. Similarly for
eqn.(7.2.2), we have:

       ε ∂v/∂t + u = −(1/ε) ∂η/∂y.

The continuity equation (7.2.3) becomes:
                                                                                               Page 130

       ∂η/∂t + (∂u/∂x + ∂v/∂y) = 0.

We try to find the solution to (N7.7a-c) by using power series expansion in ε:

       u = uo + ε.u1 + ε2.u2 + …; v = vo + ε.v1 + ε2.v2 + …; where un & vn ~ O(1). (N7.8a,b)

For η, since the flow is expected to be nearly geostrophic, (N7.7a,b) suggests that:

       η = ε.η1 + ε2.η2 + …; where ηn ~ O(1).

Substituting (N7.8) into (N7.7), we obtain for the zeroth-order equations:

       vo = ∂η1/∂x, uo = −∂η1/∂y, and ∂uo/∂x + ∂vo/∂y = 0 .                               (N7.9a,b,c)

The leading (i.e. zeroth) order then satisfies the degenerate geostrophic balance and the corresponding
velocity (uo, vo) cannot be uniquely determined because it satisfies the (leading order) continuity
equation (N7.9c); i.e. we only have two equations for three unknowns (uo, vo, η1). We therefore
examine the O(1) equations:

      ∂uo/∂t − v1 = −∂η2/∂x,
    ∂vo/∂t + u1 = −∂η2/∂y,
      ∂η1/∂t + ∂u1/∂x + ∂v1/∂y = 0.

Use (N7.10a,b) to eliminate u1 and v1 from (N7.10c):

       ∂[η1 – ζo]/∂t = 0,

where ζo = ∂vo/∂x − ∂uo/∂y = 2η1 from (N7.9a,b). Equation (N7.11) is the PV-conservation equation
(7.2.8) in non-dimensionalized form. Integrating (N7.11):

       ζo – η1 = Q(x,y,t) = Q(x,y,0), or
                                                                                                    Page 131

       − 2η1 + η1 = −Q(x,y,0)

which is the steady-state form of (7.2.13) (i.e. eqn.7.2.20 when ∂/∂y = 0) in non-dimensionalized form
(check it!). We see clearly that through our systematic scaling and series expansion, we obtain the
same steady-state equation (N7.12) which is not degenerate, and can be solved. What have we learnt
from this systematic approach to solve large-scale, slowly-varying flows with rotation?

1. The scale for length is the Rossby radius, L ~ R = (gH)1/2/f;
2. The scale for time is 1/(εf) ~ many inertial periods if ε << 1, which is the case;
3. η*/H is small ~ O(ε); but our chosen velocity scale ‘U’ is such that u*/U ~ O(1); here the asterisks
   *’s denote dimensional quantities.
4. To a good approximation (i.e. to the leading order), the velocity (uo, vo) is geostrophic;
5. The vorticity ζo can be evaluated assuming geostrophy; but the equation that relates it’s temporal
   variation (i.e. ∂ζo/∂t) to the slow time-dependent variation in η*/H (i.e. ∂η1/∂t) can only be obtained
   by considering the O(ε) balance in the series expansion. This leads to the important PV-
   conservation equation (N7.12) that may be used to solve for η1, after which the (uo, vo) may be
   computed through geostrophy – equations (N7.9a,b).

The above approach is general, and the resulting approximation is called quasi-geostrophic
approximation (or QG-approx). “Quasi” means “almost like but not exactly so” – in our case, we
used geostrophy for (u,v) but we did not use it when we wanted to relate η with ζ. Another important
thing to learn is that we can in fact include nonlinear advection terms in the momentum equations and
we will obtain a very similar equation as (N7.12) and geostrophy for (uo, vo) (see exercise 9). The
thing is, we never had to directly deal with ∂u/∂t or ∂v/∂t, just ∂ζ/∂t.

Exercise 9: Show that if nonlinear advection terms are included, we obtain instead of (N7.7a,b):

      ε Du/Dt – v = −(1/ε) ∂η/∂x,
      ε Dv/Dt + u = −(1/ε) ∂η/∂y.

and if H is still kept constant (it is not hard to also include a non-constant but slowly- varying H, i.e.
gentle slope ~ O(ε)) the continuity equation (7.2.3) becomes:

       Dη/Dt + (1 + η)(∂u/∂x + ∂v/∂y) = 0.                                                       (N7.13c)
                                                                                               Page 132

Now show that by substituting the power series expansion (N7.8a,b,c) as before, we obtain (N7.9a,b,c)
as the leading order as before, and then the O(ε) equations are still (N7.10a,b,c) as before. Therefore,
the PV-equation (N7.12) remains unchanged.
                                                                                                 Page 133

Lecture 8. Boundary Waves

Kelvin Wave:

Consider a homogeneous fluid (or the equivalent normal-mode equation, or reduced gravity), the x and
y momentum equations are (Gill’s (7.2.1) and (7.2.2)):

        u/t  fv = g /x,                                                              (8.1a)
        v/t + fu = g /y.                                                                       (8.1b)

Also, the mass conservation equation (H = constant):

        /t + H u/x = 0.                                                                         (8.1c)

Let y = 0 (or y = constant) be the straight coast. Kelvin discovered a type of rotational wave that can
propagate along the coast. He reasoned that, near the coast, v  0, so that (8.1a,b) becomes:

        u/t = g /x,                                                                            (8.2a)
        fu = g /y,                                                                               (8.2b)

Thus the y-momentum equation (8.2b) is geostrophic, i.e. the along-shore velocity “u” is always in
geostrophic balance with the cross-shore sea-level gradient. But the alongshore momentum (8.2a) is
highly ageostrophic. In fact, the x-momentum “u” is driven by pressure gradient in the same
direction (as the momentum) – as if the flow has no rotation. Thus alongshore momentum (8.2a) and
mass conservation (8.1c) in every y = constant plane (i.e. parallel to the coast) are exactly the same as
in non-rotating flow, and may be combined to give:

         t2 = c2  x2                                                                         (8.3)

Noting that “y” is kept constant in (8.3), so its solution is:

         = F(x+ct,y) + G(x-ct,y),                                                                   (8.4a)

which can be substituted in (8.2a) to give u:

        u = (g/H)1/2 [F(x+ct,y)  G(x-ct,y)].                                                       (8.4b)
                                                                                                 Page 134

Substitute  and u from (8.4) into the y-momentum equation (8.2b) now gives, since “F” and “G” are
independent solutions:

       F/y = F/R, and        G/y = G/R,                                        (8.5a,b)

where R = (gH)1/2/f is the Rossby radius of deformation. For f > 0 and y points seaward, only (8.5b)
with the exponentially decaying solution “G” is allowed, so that:

       G = A(x-ct).exp(y/R)                                                                (8.6a)

which represents a wave propagating in the positive x-direction – see fig.8.1.   If f < 0 and y again
points seaward, then (8.5a) with the “F(x+ct,y)” solution must be chosen:

       F = A(x+ct).exp(+y/R)                                                                (8.6b)

Figure 8.1   Kelvin wave.

Substitute (8.6a) into (8.4) then gives for  and u:

        = A(x-ct).exp(y/R)                                                                (8.7a)
       u = (g/H)1/2 A(x-ct).exp(y/R) = (g/c).                                             (8.7b)

where c = fR. If we think R as being always > 0, then c > 0 in northern hemisphere, and c < 0 in
southern hemisphere. Then (8.7a) and the second form of (8.7b) are valid for both f > 0 and f < 0.

 In summary:

   1. Kelvin wave travels at the speed “c”, exactly the same as for non-rotating shallow-water waves;
                                                                                                                       Page 135

      2. The wave form “A” is the same on every y=constant plane (parallel to the coast), except that its
         amplitude decreases exponentially with distance from the coast;
      3. Therefore the wave has significant amplitude only within a distance R  c/f from the coast;
      4. The velocity normal to the coast is zero: v = 0;
      5. The velocity parallel to the coast, u, is in geostrophic balance with the sea level  (show this
         from (8.7)). Therefore, from the second form of (8,7b), u is > 0 (< 0) in the northern
         (southern) hemisphere where coastal sea-level |y=0 is > 0. And the signs of u are reversed
         where the |y=0 is < 0.

Things to discuss on Friday Apr/30/2010:

      1.   Kelvin wave energy flux..
      2.   UE/t .. total atmos+ocean mass transport = 0;
      3.   Why T/z << u/z .. mixed layer for temperature?
      4.   Why barotropic geostrophic flows must be parabathic?
      5.   Read Gill 10.9 (Storm Surge) for next week.

Coastal Setup (or Setdown)by Wind:

Consider an infinite coast bounding an ocean of uniform depth H (Fig.8.1) on a rotating earth with
constant f. At t = 0, a uniform kinematic wind stress ox blows in x-direction. Since coast is infinite
and wind is uniform, we consider the case /x = 0. Then instead of (8.1), we have:7

           u/t  fv = (oxbx)/H                                                                                      (8.8a)
           v/t + fu = g/y  by/H                                                                        (8.8b)
           /t + Hv/y = 0                                                                                            (8.8c)

where (bx,by) is bottom frictional stress assumed = 0. Wind blows in x-direction, so there is Ekman
transport, hence v  0. We can (8.8b)/t, then use (8.8a) for u/t, and use (8.8c)/y for 2/yt:

           2v/t2 + f2v c22v/y2 = fox/H                                                                            (8.9)

After a time longer than inertial period, t >> 1/f, transient inertial oscillations die away, the 2/t2 term
is small, i.e. steady state.8 The solution consists of exponentials, and we get rid of the growing, e+ part
and also make the solution satisfy v = 0 at y = 0:

7   Derive this by integrating the linearized equations with vertical shear stress terms, and assuming  = constant.
                                                                                                                     Page 136

           v = ox/(fH) [1  ey/R].                                                                                     (8.10)

Divergence (u/x + v/y) of water near the coast is (since /x = 0), from (8.10):

           v/y = ox/(RfH)ey/R                                                                                        (8.11)

If ox > 0, we then have convergence i.e. v/y < 0 and sea-level rises from (8.8c).                    If ox < 0, we have
divergence i.e. v/y > 0 and sea-level falls. From (8.8c), the rise and fall is:

            = ox/(Rf)ey/R.t                                                                                            (8.12)

Therefore, although v is steady,  is not, and linearly increases with time. Since from (8.8b) the “u” is
in geostrophic balance with the sea-level, it too linearly increases with time:

           u = +(ox/H).ey/R.t                                                                                           (8.13)

The PV-equation is:

           (  f/H) = 0                                                                                                (8.14)

so that the (relative) vorticity also increases indefinitely.

In practice, as “u” increases, friction eventually becomes large on the RHS of (8.8a).9

Kelvin wave energy flux:

For a fixed y, Kelvin wave behaves as in flow without the earth’s rotation – see eqns.(8.3) & (8.4).                            So
we derive first the energetic for f = 0, and then apply it to Kelvin wave. The equations are:

           /t + (Du)/x + (Dv)/y = Q                                                                     (8.15a)
           u/t  fv = g/x  ru                                                                           (8.15b)

8    To see this, define non-dimensional variables: tn = ft, yn = y/R, vn = v/U, and onx = ox/(fUH), where << 1, then
     vn/tn2 + vn  2vn/yn2 = onx, and the first term can be neglected.
    2 2

9 Effects of friction are seen by considering very near-coast region where “v” is small.    Then (8.8a) gives ox  bx, etc.
                                                                                                    Page 137

        v/t + fu = g/y  rv                                                              (8.15c)

Here, D = H+, Q = source (in m s-1) and r = constant bottom friction coefficient.       Multiply (8.15a) by

        (g2/2)/t + g[(Du)/x + (Dv)/y] = gQ                                         (8.16)

Note that

        g2/2 = PEA                                                                           (8.17a)

is potential energy per unit (horizontal) area, since it is = gz dz’, where the integral is from z’ = 0 to
. For our case,  = constant since the fluid is assumed to be homogeneous. Equation (8.16) is
therefore the potential energy equation.

Multiply (8.15b) by Du and (8.15c) by Dv, and then sum the resulting 2 equations, we obtain, with
KE = (u2+v2)/2:

        D.(KE)/t     = g(uD/x+vD/y) rD(u2+v2)
                       = g{ [(uD)/x+(vD)/y]  [(Du)/x+(Dv)/y] }  2rDKE

which after using (8.16) and (8.17a) and also approximating D.(KE)/t = (D.KE)/t  KE./t 
(D.KE)/t, gives:

        (KEA + PEA)/t = g[(uD)/x+(vD)/y] + gQ  2rDKE                              (8.18)


        KEA = D.KE = D. (u2+v2)/2                                                             (8.17b)

can be seen to be the kinetic energy per unit (horizontal) area of the column of water (of depth “D”).

For Kelvin wave, the solution is (8.7a,b).     Let’s consider the special case when the amplitude is of a
simple cosine form:

        A(x-ct) = cos(x-ct) = cos( ),        = x-ct.                                           (8.19)

So (8.7a,b) become:
                                                                                                  Page 138

        = o cos( ) exp(y/R)                                                                (8.20a)
       u = o (g/H)1/2 cos( ).exp(y/R) = (g/c).                                             (8.20b)

Let <.> denotes the mean over all “ ”, so that since cos2 .d /(2    where the integration  is from
0 to 2, we have from (8.17a,b) that the mean PE or KE per unit x-length (PEx or KEx) is:

       PEx = g<2>/2.dy = go2R/8 = H<u2>/2.dy = KEx                                     (8.21)

where the integral over ‘y’ …dy is taken from y= 0 (coast) to y = . Therefore, the mean kinetic and
potential energies per unit x-length are equal. This is the same conclusion reached for f = 0 step-
problem in Gill’s Chapter 5 that the energy of motion is partitioned equally between KE and PE.

The mean energy flux at any “x” along the coast is given by (see the 1st term on the RHS of 8.18):

       Efluxx = gH<u>dy = g2Ho2/(4|f|)                                                   (8.22)

For Kelvin wave propagating along the coast, the 2nd term on the RHS of (8.18), (vD)/y, is small
compared to the 1st term, so that the energy flux is Efluxx given by (8.22). Now, if friction can be
neglected, we have by integrating the total energy equation (8.18) over an x-section of the coast where
Kelvin wave passes, say from x1 to x2, that:

       (<KEA> +< PEA>)dxdy}/t = 0 = gH[<u>|x1<u>|x2].dy                    (8.23)

The LHS is = 0 since the total energy of Kelvin wave over (x1, x2) is constant. Therefore, <u>|x1 =
<u>|x2 and the mean energy flux Efluxx is constant. Therefore, if at any “x2” the water depth H2 is
shallower than the section previously x = x1 < x2 and H2 < H1, then the amplitude at x2 must increase:

       [o]|2/[o]|1 = (|f2|/|f1|).(H1/H2)1/2 > 1, if |f| does not change much.               (8.24)

Therefore, flooding and drying can occur when a Kelvin wave (generated say by a distant storm)
propagates into a shallow coastal region; made worse if propagation is poleward (so that |f2| > |f1|).


   1. Step-function alongshore wind produces inertial oscillations that radiate away;
   2. The wind then drives steady Ekman transport towards or away from the coast;
   3. The convergence or divergence causes coastal sea-level  to rise or fall linearly with time;
                                                                                                Page 139

   4. The alongshore velocity “u” also increases linearly with time, in geostrohpic balance with ;
   5. Both “u” and “” decays exponentially away from the coast over a distance ~ R.
   6. Eventually bottom friction is important, and balances the wind.

Problem 8.1. Modify the plume-on-shelf problem in http://pom-training.pbworks.com/ (File:
runpom08_01_assignment5b_horcon050_nadv1 in
ftp://aden.princeton.edu/pub/lyo/pom_gfdex/wmo09training/anIntroCourseNumOceanExpsUsingPOM ).
Domain is a square or rectangular basin, coastal walls on all sides. Then ramp (in time) an alongshore
wind in a finite coastal strip 50km alongshore and 10km cross-shore say at the center portion of the
western boundary. Examine the resulting Kelvin wave.

   a.   No rivers, no tides, and ocean is initially homogeneous, bottom is flat;
   b.   Same as ‘a’ but ocean has shelf/slope all around;
   c.   Same as ‘a” but ocean is vertically stratified;
   d.   Same as ‘b” but ocean is vertically stratified;
   e.   Same as ‘a’ but a circular cyclonic wind system (typhoon) propagating from open ocean (east;
        i=im) to land (west; i=1) along the center latitude (j=jm/2) of the domain at a speed of 5 m/s.

Problem 8.2. Design an experiment to verify (8.24).
                                                                                                    Page 140

9. Dynamics of Long Rossby Wave (Propagation) Forced by Wind Stress Curl

Consider the linearized reduced-gravity equations:

       h/t + (hu)/x + (hv)/y = Q  h                                                  (9.1a)
       u/t  fv = g’h/x + ox/H                                               (9.1b)
       v/t + fu = g’h/y + oy/H                                               (9.1c)

where (x, y, t) = (west-east, south-north, time), h = upper-layer depth, H = mean upper-layer depth, Q =
a source term, = Newtonian cooling coefficient, (u,v) = horizontal velocity, g’ = g/o is reduced
gravity, f = Coriolis parameter, and (ox, oy) = wind stress vector. For slow motions with time scales
much larger than 2/f ~ days, the /t terms underlined in (9.1b,c) are small compared to the Coriolis
terms, so may be neglected. Then the motion is nearly geostrophic:

        vg  (g’/f) h/x                                                                   (9.2a)
        ug  (g’/f) h/y                                                                  (9.2b)

where the subscript “g” denotes “geostrophic.”   From (9.2), we have

       ugh/x + vgh/y = 0.                                                               (9.3)

Equation (9.1a) is:

       h/t + (uh/x + vh/y) + h(u/x + v/y) = Q  h                                 (9.4)

The term “(uh/x + vh/y)” is nearly zero by (9.3), while the “h(u/x + v/y)” term is obtained by
taking the curl of the momentum equations (9.1b,c), i.e. by “(9.1c)/x  (9.1b)/y”, which gives:

       h(u/x + v/y)  H(u/x + v/y) = Hv/f + ×o/f
                                            Hvg/f + ×o/f
                                           = [(Hg’)/f2]h/x + ×o/f
                                           = R2 h/x + ×o/f                             (9.5)

where  = df/dy and R = (Hg’)1/2/f is the baroclinic Rossby radius.   Using (9.5), equation (9.4) then
becomes, with CR = R2, the long Rossby wave speed:

       h/t  CR h/x =  ×o/f + Q  h.                                                  (9.6)
                                                                                                  Page 141

Eqn.(9.6) describes the westward propagation of (long) Rossby wave which without forcing (RHS=0)
has the solution of the form

       h = G(x+CRt)                                                                       (9.7)

But the wave is also being continuously forced by some distribution of sources and sinks on the RHS
given by

       (i)      ×o/f = Ekman pumping (note water is added to upper-layer if this is <0), and
       (ii)    Q representing some waters entering (>0 i.e. entrainment) or leaving (<0 i.e.
               detrainment) the upper layer,

Note that mathematically the source/sink terms ×o/f and Q are indistinguishable from each other.
At the same time the wave is also damped by

       (iii)    h = Newtonian cooling or heating,

‘cooling’ or ‘heating’ because “h” can physically be imagined as representing the upper-layer warm
water. It is “damping” because of the minus sign, so that when h, this term tends to produce h/t,
and vice versa.

Note that also the wind stress acts on a thin surface Ekman layer (~E), and its effects on the upper-
layer motion is mainly through Ekman pumping as a result of differential Ekman transports inside E
(as we will shortly see).
                                                                                                 Page 142

Lecture 10. El Niño, La Niña and the Southern Oscillation

0. Historical Introduction:

El Niño is the name given to the irregular development of warm ocean surface waters along the coast of
Ecuador and Peru. It develops when there is change in the circulation of the atmosphere across the
tropical pacific. Globally, the development of El Niño is also associated with a number of other
climatic changes in other parts of the world.

To explain El Niño, it is necessary to explain how the ocean adjusts to the changes in the surface winds.
Responses (wind and ocean currents) near the equator are rapid. Recall Gill’s adjustment with step
problem – non-rotational gravity waves (gH)1/2, and also Rossby waves R2 (for example) are faster
because R = Ci/f is larger. So ocean responds fast to winds, and air-sea coupling is particularly
efficient: slight change in one leads to immediate changes in the other and vice versa.

Gilbert Walker (1923~1937): knew that interannual pressure fluctuations over the Indian Ocean and
eastern tropical Pacific are out of phase: “When pressure is high in the (eastern) Pacific Ocean it tends
to be low in the Indian Ocean from Africa to Australia.” He called this irregular fluctuation the
Southern Oscillation (SO). The SO is correlated with major changes in the rainfall patterns and wind
fields over the tropical Pacific and Indian Oceans, and with temperature fluctuations in southeastern
Africa, southwestern Canada, and southeastern USA.

El Niño is that phase of the SO when the trade winds are weak and when pressure is low over the
eastern and high over the western tropical Pacific.
                                                                                                Page 143

La Niña is apposite when sea surface temperatures in the central and eastern tropical Pacific are
unusually low and when the trade winds are very intense.

                                                                             Page 144

CHAPTER 2???: Internal Waves and Rays in the Ocean
      In chapter 1 we studied how wind generates ocean currents and sea-
level variations assuming homogeneous fluids ( = constant). We also
focused on motions having horizontal scales that were much larger than the
vertical scales; in other words, we assumed also that the hydrostatic equation
was satisfied. These simplifications allowed us to examine rudimentary
boundary-layer dynamics and to understand their profound effects on the
large-scale (interior) oceanic flows. In this chapter, we remove these
restrictions. We examine internal fluid motions in a stratified fluid ( is
not constant, and generally a function of the vertical coordinate z). We do
not assume hydrostasy, but we restrict ourselves to motions which are of
small amplitudes that the equations are linear, and which are also of small
horizontal scales and short time scales that the Coriolis force can be
neglected. We will study rays – which are paths (generally in three
dimensions) along which internal wave energy propagate. We then attempt
to numerically model these waves (and rays) in two-dimensional vertical-
section plane (xz-plane, say). We conclude the chapter by showing how
the study of internal-wave rays may be applied to an entirely different
physical system in the ocean: that of the study of ray propagation of
hydrostatic and rotational topographic Rossby waves. We show that under
certain conditions, the two systems are mathematically equivalent, and
demonstrate a surprising physical consequence of rotation.

2-1: The Spring-Mass System - Generalized Inertia & Stiffness

      An oscillatory system such as the wave motion can be viewed as a
spring-mass system:

                                                                    (2-1)
                                                                             Page 145

where m is the mass,  is the spring’s stiffness, and z is distance by which
the spring has been pulled or compressed (z = 0 is where the spring is in its
natural position with zero tension). The solution to (2-1) consists of
cosines and sines and the frequency squared is:

      2 = /m                                                       (2-2)

Equation (2-1) can be rewritten as


Now, the potential energy (PE) is = z2/2 while the kinetic energy (KE) is =

       , so that (2-3) is simply a statement that the total energy is constant
in a system without dissipation. Equations (2-2) and (2-3) suggest that:


We will now use this idea to derive the frequency of oscillations for surface-
interfacial waves, then internal waves also.

2-2: Surface & Interfacial Waves

      Consider then a thin (horizontal) interface that separates a heavier
fluid of uniform density 2 below from a lighter fluid of uniform density 1
above. We assume that that the interface is far from the free surface above
and also from the ocean bottom below, so that the free surface z = zT and
ocean bottom z = h do not interfere too much with the interfacial motion;
                                                                                  Page 146

here z = 0 at the interface. Let the interface be disturbed by a (small)
displacement . The PE of the lower layer per unit area is:

                               
                                                                       (2-5a)

Thus the lower-layer PE is increased by 2g2/2 from its undisturbed state

(i.e. when  = 0). Similarly, the PE of the upper layer per unit area is:

                                                                      (2-5b)

and in this case the upper layer has its PE decreased, 1g2/2 from its

undisturbed state. Therefore the total PE per unit area due to  is:

                                                                      (2-6)

and the generalized stiffness, the coefficient of 2/2, is seen to be:

      generalized stiffness = g(2  1)                                 (2-7)

      The disturbed interface displacement  produces currents and
therefore kinetic energy. To derive an expression for it, we need to
examine surface (or interfacial) waves.

Free-Surface Waves on Water with Constant Density:
                                                                              Page 147

      Consider then a homogeneous fluid with a free surface  (say we

think of the lower layer only with uniform density 2, and in the followings

think of “o” as “2”). Equations (1-24a,b,c) written in vector form and
dropping the nonlinear, Coriolis and viscous terms become:

      o u/t =  p’                                                 (2-8)

where p’ is the pressure in which a hydrostatic (undisturbed) pressure
distribution po(z) (integrating equation 1-27 from z to 0; see equation (1-

      po(z) = patmos  ogz                                           (2-9)

has been removed (from the total pressure p):

      p’ = p  patmos + ogz                                          (2-10)

and the  that multiplies the acceleration on the LHS of (2-8) has been

approximated by o.      Taking the curl of (2-8) shows then that the vorticity

field u is independent of time:

      ( u)/t = 0                                                   (2-11)

Therefore, the rotational part of the velocity field induced by this stationary
vorticity field is steady, hence its corresponding pressure is zero by (2-8),
and it cannot disturb . The remaining part of u is irrotational (i.e. it has
                                                                             Page 148

zero vorticity) and therefore can be written as (for simplicity we still call it
u, and ignore the rotational part):

      u=                                                               (2-12)

the gradient of a scalar function, so called the velocity potential. Only this
part exhibits the fluctuations associated with wave propagation. From (2-
12), we see that, since water is incompressible (equation 1-25), we have:

            =0                                                         (2-13)

This Laplace equation cannot by itself support wave motions if the
boundaries are steady. But it can if the free surface is undulating in time.
From (2-10), we have, setting p = patmos at z = :

      p’ = og,     at z =   0                                      (2-14)

the last approximation is because          has been assumed to be small

compared to the total water depth. Using (2-12) and (2-8) (giving o /t

= p’), equation (2-14) gives:

       /t = g,     at z = 0                                        (2-15)

      At the free surface, equation (1-31a) also holds (i.e. the free surface is
a material surface), thus:

      /t = w =  /z,     at z = 0                                  (2-16)

Equations (2-15) and (2-16) then give:
                                                                          Page 149

      2 /t2 = g /z,             at z = 0                       (2-17)

       Equation (2-13) describes the motion in the fluid interior away from
boundaries, and with equation (2-17) the system is almost complete. We
assume that side boundaries are sufficiently far (in fact they are infinitely
far) that waves propagate horizontally unimpeded. Besides the surface
condition (2-17) we then also have the condition at the ocean bottom (1-31b)
which for a flat bottom is simply that the vertical velocity is zero:

       /z = 0,             at z = h                              (2-18)

We will, however, relax the matter even further: we will assume that the
bottom is also sufficiently deep that a solution (satisfying 2-13 and 2-17)
that decays sufficiently rapid with depth will be adequate.

      Consider horizontal propagation in the x-direction.      The Laplace
equation admits a separable solution of the form:

        =     (z) exp[i(tkx)]                                     (2-19)

where  is frequency and k wave number. Substituting into (2-13) and

solving for     (z):

        (z) =     o    ekz                                          (2-20)

This solution decays with depth z ~ . Substituting (2-20) into (2-17)

now gives the dispersion relation connecting  and k:
                                                                          Page 150

      2 = gk                                                       (2-21)

      The solution (2-20) allows us now to compute the total kinetic energy:

                                    

where equation (2-13) and the divergence theorem are used to obtain the
second integral over the bounding surface S consisting of the surface and the
bottom only. The bottom contribution vanishes, either because of (2-18) or
because the solution decays exponentially with depth according to (2-20).
At the surface,  /n =  /z = k-1(/t)2 after using (2-19), (2-20) and
then (2-16). Equation (2-22) then gives for the kinetic energy per unit area:

      KE = (o/2) k-1(/t)2                                       (2-23)

Comparing with equation (2-4), then, for this free-surface system, the
generalized inertia is:

      generalized inertia = o k-1                                  (2-24)

We also have, from (2-6), PE = go2/2, so that

      generalized stiffness = og                                   (2-25)

Equation (2-4) then gives 2 = gk in agreement with (2-21).

Dispersion Relation for Interfacial Waves:
                                                                            Page 151

      We saw from (2-6) that the disturbance  gives rise to the PE
available for wave motion.         The same disturbance gives rise to kinetic
energy given by (2-23) with the o replaced by the lower-layer density “2,”

and similarly in the upper layer with o replaced by “1.” The total

      KE = [(1 + 2)/2] k-1(/t)2                                   (2-26)

Equation (2-4) then gives,

      2 = gk(2  1)/ (2 + 1)                                      (2-27)

Comparing with the frequency of the wave on a free surface, equation (2-
21), we see that the interfacial wave frequency (and hence also the
corresponding phase speed = /k) is smaller, by a factor proportional to the
square root of the density difference between the two layers. Note that (2-
27) leads to (2-21) if 2 >> 1.

2-3: Internal Waves

      Extensions of the above ideas to the case when the density o(z) is
continuous require an exposition to the idea of a vertically stable stratified
fluid when the density decreases with height. If a fluid parcel at z is moved
up by a small amount to z + , the density of the surrounding fluid is now

      o + (o/z), with o/z < 0,  > 0.                          (2-28)
                                                                              Page 152

Through the hydrostatic equation (1-27), po(z) = o(z)g, pressure of the
surrounding fluid is also reduced to:

        po(z)  p = po(z)  o(z)g.                                   (2-29)

This pressure drop pressed upon the fluid parcel results in its density being
reduced. Assuming isentropic process, the parcel’s density is reduced to:

        o(z)  o(z)g/co(z)2,                                         (2-30)

where co(z)2 = (p/)s is the square of the speed of sound, and the subscript
“s” denotes constant entropy. The parcel’s excess weight (per unit volume)
over the surrounding buoyancy is then g[(2-30)  (2-28)], or,

        g[o(z)g/co(z)2 + (o/z)] = oN2(z)                        (2-31)


        N2 = [g2/co(z)2 + g(o/z)/o]                                (2-32)

In order that there is a restoring gravitational force on the fluid parcel, this
excess weight must be positive, which means that N is real and also that the
vertical density gradient must be (from 2-31):

        o/z < o(z)g/co(z)2                                         (2-33)

This shows that for stable stratification, it is not sufficient to merely require
that o/z < 0.
                                                                              Page 153

      For N2 > 0, the restoring force (2-31) is like the force exerted by the
spring in the spring-mass system described above. The  is also similar to
the perturbed interfacial displacement described previously, but since the
fluid is continuously stratified, the  may be viewed as a generalized vertical
coordinate for each fluid parcel. Then the potential energy per unit volume
that results from  is:

                                                                      (2-34)


      generalized stiffness = o(z)N2(z)                                (2-35)

What is the corresponding kinetic energy per unit volume? Equation (2-23)
gives KE per unit area for a surface wave having a discontinuity in density,
and a penetration depth scale ~ k-1. With continuous stratification, we do
not know a priori what the penetration depth scale is (in fact there is none(!)
 as we shall see, there exists wave-propagation in the vertical also). It is
more useful, therefore, to consider KE per unit volume which from (2-23) is:

      KE/volume = (o/2)(/t)2                                        (2-36)

With surface or interfacial waves, the “” can only move vertically, and “the
vertical” therefore represents the only generalized coordinate.             The
generalized inertia of these purely vertical oscillations is from (2-36) = o(z).
But this must be the minimum inertia in the case of internal waves in
continuously stratified fluid. This is because in a continuously stratified
fluid, it is possible for a fluid parcel to oscillate along a line that makes an
                                                                            Page 154

angle T to the horizontal (or at  = /2  T to the vertical); the

corresponding displacement is       , say (Figure 2-1).      These types of

oscillations clearly increase the kinetic energy which becomes ( o/2)( /t)2,

hence also the generalized inertia (the coefficient of (/t)2/2), which is

now increased to a value = o(z)/sin2(T). The purely vertical oscillations,

for which T = /2, therefore represent a state of minimum inertia, hence
maximum possible frequency.

Figure 2-1. Parcel fluid oscillations (displacement = , say) inclined at an

angle T to the horizontal, or at  = /2  T to the vertical. Thus           =


In general then, (2-36) is replaced by:

      KE/volume = (o/2)(/t)2/sin2(T)                             (2-37)
                                                                            Page 155

and equation (2-4) gives, with (2-34) and (2-37):

      2 = (generalized stiffness)/(generalized inertia) = N2 sin2   (2-38)


       = N.sin(T) = N.cos(),  = /2  T                          (2-39)

Figure 2-2 illustrates schematically how an IW excited in a thermocline may
be trapped within it. It makes use of the idea expressed in (2-39) that the
maximum possible frequency of oscillations cannot exceed the local N.

Figure 2-2. Left: a typical N-profile in the ocean’s thermocline.     Right: an
internal wave (IW) with  = N1 excited within the thermocline cannot
escape to depths (beyond upper and lower boundaries) where the local N <


2-3A: Find out what the typical values of N are in the ocean. Then estimate
typical (minimum) IW periods and compare them with the Coriolis periods.
                                                                           Page 156

2-3B: Re-derive (2-33) by letting  < 0 instead of positive.

2-3C: Study the following alternative derivation of  = N sin(T) from the
equations of motion. Then close the book, and re-derive yourself.

Derivation of (2-39) Using the Equations of Motion:

In equation (1-24), we drop rotation, viscosity and the nonlinear advective
terms, and replace the  on the LHS (i.e. the acceleration terms) by o (this
is the first part of the so called Boussinesq approximation), to obtain in
vector form the following momentum equation:

      o u/t =  p + g                                            (2-40)

where g = (0, 0, g).         Subtract from this the hydrostatic relation (see
equations 1-27, 1-28 and 1-29):

        po = og                                                     (2-41)

we obtain:

      o u/t =  p’ + ’g                                          (2-42)

The conservation of mass (i.e. continuity equation) is (see e.g. “Fluid
Dynamics” by G.K. Batchelor, 1966):

      /t +      (u) = 0                                          (2-43)
                                                                             Page 157

Now write /t +       (u) = ’/t +   [(o+’)u] and use the second part of

the Boussinesq approximation that density fluctuations ’ (which of course
affects motions, from equation 2-42) in the conservation of mass are
neglected. In other words, the fluid motions of concern evolve sufficiently
slowly that the rate of change of ’ as well as its spatial gradients do not
affect the mass balance very much. Therefore,

         (ou) = 0                                                      (2-44)

Also, since ’ is the excess density of a fluid parcel over its surrounding

o(z) (see equation 1-29), we have from the derivation of (2-31) that

      ’ = oN2/g                                                      (2-45)

for the (excess) density when the fluid parcel is displaced a small distance 

from its initial position; note that (2-45) is valid for positive or negative 
(exercise 2-3B).

We now manipulate (2-42), (2-44) and (2-45) to derive an equation solely in
terms of one variable which we will choose to be ow = upward component
of mass flux.

Take the divergence of (2-42) and use (2-44),

            p’ = g ’/z                                              (2-46)

then eliminate p’ from the z-component of (2-42), which is:
                                                                                    Page 158

       (ow)/t + p’/z = ’g                                               (2-47)

we obtain,

            (ow)/t = g2’/z2  g       2
                                                ’ = g(2’/x2 + 2’/y2)   (2-48)

The ’ on the RHS is now eliminated by first taking the time derivative, (2-

48)/t, then noting from (2-45) that

       g’/t = o(/t)N2 = owN2,                                          (2-49)

 to finally obtain,

            (2q/t2) = N2   h
                                  q                                            (2-50)

where q = ow, and            h
                                      = (2/x2+2/y2) is the horizontal Laplacian
operator. Incidentally, equation (2-49) written as

       b/t = wN2, b =  g’/o = buoyancy,                                  (2-51)

may be recognized as being derivable from the density equation

       D/Dt = 0                                                               (2-52)

after linearization about the base (rest) state o(z) etc. (e.g. Gill, 1982).

The Dispersion Relation:
                                                                        Page 159

      For constant N, we may look for a (separable) wave-like solution of
equation (2-50):

      q = ow = qo exp[i(tkxlymz)],                           (2-53)

where k = (k, l, m) is the wavenumber vector. Substituting this into (2-50)

      2 = N2(k2+l2)/(k2+l2+m2)                                   (2-54a)


       = N(k2+l2)1/2/(k2+l2+m2)1/2                               (2-54b)

which shows, in agreement with (2-38) or (2-39), that the maximum possible
frequency of oscillations is N.
                                                                              Page 160

Figure 2-3. Incompressibility means that the fluid parcel oscillation with
displacement vector       and the internal wave vector k are perpendicular;
assuming a plane-wave solution. (See also Figure 2-1).

       In fact, equations (2-54b) and (2-39) are equivalent (apart from the
fact that we did not explicitly assume constant N in 2-39; see however,
later). This may be proven as follows. Since the fluid is incompressible,
we have from (2-44) that for a plane-wave solution such as (2-53),

        ku=0                                                          (2-55)

meaning that the wavenumber vector k and the fluid velocity u are
perpendicular to each other. Since in Figure 2-1, the u is in the same
direction (and plane) as the oscillation displacement vector , we have then

k and     are also perpendicular as shown in Figure 2-3. Equation (2-39)

then follows by noting that the angle T made to the horizontal by the fluid

parcel’s displacement vector       during the oscillations is also the angle
spanned by the wavenumber vector k with the vertical axis (Figure 2-3).

Expressions for p’ and u in terms of qo:

    For completeness, we give here the corresponding formulae for p’ and
u. For plane wave, we have (c.f. equation 2-53):

        (p’, u) = (po, uo).exp[i(tkxlymz)].                       (2-56)

Substitute the p’ into /t[(2-46)] and using (2-49) and (2-53), we obtain,

        p’ = m(k2+l2)-1qoexp[i(tkxlymz)].                       (2-57)
                                                                               Page 161

Substitute (u, v) into the (x, y)-component of (2-42) and using (2-57) gives

      ou = [km(k2+l2)-1, lm(k2+l2)-1, 1]qoexp[i(tkxlymz)]      (2-58)

where the ow part just repeats (2-53).

A Summary of the Governing Equations for IW’s:

      Before deriving the corresponding dispersion relation for TRW’s in
the next section, it is useful to summarize here the governing equations for
IW’s so that one may better see the differences with the TRW equations.
In IW’s the momentum equations are (2-42):

      o u/t = p’/x                                              (2-59a)

      o v/t = p’/y                                              (2-59b)

      o w/t = p’/z  ’g                                        (2-59c)

The conservation of mass (i.e. the “continuity equation”) is

      /t +    (u) = 0                                             (2-60a)


      D/Dt +  u = 0                                                 (2-60b)

which we used for IW’s in the form of equation (2-44) keeping the o inside

the divergence operator ( ) mainly for conciseness so we obtained an
                                                                             Page 162

expression for the vertical mass flux ow. However, for incompressible
fluid, it may be shown from (2-60b) that the first term is of O(1/co2) so it is
very small and therefore

         u=0                                                           (2-61)

which is the impressibility condition stated previously in equation (1-25).
The important point is that in deriving the (oceanic) IW equations we have
assumed that the fluid is incompressible. Finally, the density equation is
(2-51), repeated here:

      b/t = wN2, b =  g’/o = buoyancy.                           (2-62)

Equations (2-59a,b,c), (2-61) and (2-62) give five equations for the five
unknowns (u, v, w, ’, p’).

2-4: Topographic Rossby Waves

      In this section, we show the curious equivalence for the two-
dimensional propagation of internal waves studied in the previous sections
with the propagation of hydrostatic bottom-trapped topographic Rossby
waves in the horizontal xy-plane. In other words, we show that one can
study the propagation of non-hydrostatic internal waves using a hydrostatic
model, or vice versa!

      We derive the corresponding dispersion relation for TRW’s. We
consider fluid motions over a sloping bottom. We will assume that the
scales are large that the earth’s rotation and hence also the Coriolis
parameter f (Chapter 1) are important, but not too large so that we may
ignore the variations of f with latitudes, i.e. f/y = 0. The scaling for these
large-scale and relatively slow (compared to the earth’s rotation) flows have
been previously detailed in chapter 1 (see, e.g. equation 1-38). It is more
instructive for the present purpose of pointing out the similarities and
                                                                         Page 163

differences between TRW’s and IW’s to formulate by modifying equations
(2-59), (2-61) and (2-62). In TRW’s we also assume that the fluid is
incompressible so that equation (2-61) holds. The density equation (2-62)
also remains the same (of course!). The x and y momentum equations,
however, are modified by adding the Coriolis terms, and the z-momentum
equation is reduced to the hydrostatic equation (even though the fluid is in
motion!); thus equations (2-59a,b,c) become:

      u/t  fv = P’/x                                         (2-63a)

      v/t + fu = P’/y                                         (2-63b)

      P’/z = b                                                   (2-63c)

      u/x + v/y + w/z = 0                                    (2-63d)

      b/t = wN2                                                 (2-63e)

where P’ = p’/o, the buoyancy b is used in place of (g’/o) in the
hydrostatic approximation (2-63c), and the incompressibility and density
equations are repeated from (2-61) and (2-62) respectively. We assume
fluid bounded only at top and bottom but unbounded in the horizontal x and
y directions. The boundary conditions are:

      w = 0 at z = 0                                               (2-64a)
      w = uh/x  vh/y                                         (2-64b)

Take the curl of (2-63a,b) gives:

      /t + f (u/x + v/y) = 0                                (2-65)

Taking the divergence of (2-63a,b) gives:
                                                                          Page 164

      /t (u/x + v/y)  f =                h
                                                          P’        (2-66)

and this becomes after eliminating the horizontal divergence term from (2-

      2/t2 + f2  f              2
                                    h P’   =0                       (2-67a)

or, assuming that the time scales of motions are much larger than the earth’s
rotation period,

      f    h
                     P’ = 0                                         (2-67b)

We now (/z) then (/t) equation (2-63c), and then use (2-63e) to

substitute for (2b/zt) to obtain:

      (/t)(2P’/z2) = N2w/z                                   (2-68)

Use now the incompressibility equation (2-63d) to replace w/z by the
horizontal divergence, and then use (2-65) to obtain:

      (/t)(2P’/z2) = (N2/f) /t                              (2-69)

Eliminate  by using (2-67b) to obtain, finally,

      (/t)[        h
                             P’ + (f2/N2) 2P’/z2] = 0             (2-70)

The boundary conditions are written in terms of P’ by using (2-63c) and (2-
                                                                                         Page 165

      2P’/zt = 0 at z = 0                                                        (2-71a)

      2P’/zt = (N2/f) (h/y P’/xh/x P’/y) at z = h                      (2-71b)

where also fv = P’/x and fu = P’/y from (2-63a,b) have been used,
again assuming that the time scales of motions are much longer than the
earth’s rotation period.

Substituting as before a separable plane-wave solution:

      P’ =       (z).exp[i(tkxly)]                                               (2-72)

(2-70) gives:

        zz      2
                     = 0,         2
                                      = (N2/f2)(k2 + l2)                            (2-73)

The solution is:

         = Ae z + Be        z

where A and B are constants. Equation (2-71) requires                  zt,   which is:

        zt   = i    z   = i [Ae z  Be z]                                        (2-75)

Applying (2-71a) at z = 0 gives A = B =                    o,   say.   Equation (2-71b)

        z    =  [N2/(f)][hyk  hxl] =  [N2/(f)][Kh  h]|z, at z = h
                                                                        Page 166


Equaitons (2-75) and (2-76) then give:

       (/N) tanh( h) = sign(f) (nk  h)|z                        (2-77)

Now, tanh( h) = tanh(h/htrap), where

       htrap =        = |f|/(N|Kh|), where |Kh| = (k2 + l2)1/2    (2-78)

so that if

       h/htrap > 1                                               (2-79)

then tanh( h)  1 and (2-77) becomes:

        = sign(f) N | h| sin(T)                                 (2-80)

where T is the clockwise angle that the wave vector makes with   h (Figure
2-4). Equation (2-80) is in the same form as the corresponding IW
dispersion relation, equation (2-39).
                                                                                   Page 167

Figure 2-4. A sketch that describes the relation between topography h, its gradient h,
TRW wave number vector K and group velocity Cg for the northern hemisphere. The
T is the clockwise angle the K-vector makes with h.
                                                                       Page 168


Batchelor, G.K., 1967: Fluid Dynamics.
Baumert et al., 2005:
Gill, A., 1982: Atmospheric & Oceanic Dynamics
Mellor, G.L., 2004: POM User Guide.
Mellor G.L. and Yamada, T., 1982: Review of Turbulence Models.
Pedlosky, J., 1979: Geophysical Fluid Dynamics.
Roach, P., 1972: Computational Fluid Dynamics.
Van Dyke, M., 1964: Perturbation Methods in Fluid Mechanics.
Stewart, Robert H., 2000: Introduction to Physical oceanography. Columbia
University Press, http://www.earthscape.org/r3/stewart/index.html.

To top