A Coriolis Tutorial: Part 2: A Rotating Shallow Water Model, Potential Vorticity and Geostrophic Adjustment

a Coriolis tutorial
Part 2: a rotating shallow water model, potential vorticity and

geostrophic adjustment
James F. Price
Woods Hole Oceanographic Institution
Woods Hole, Massachusetts, 02543
www.whoi.edu/science/PO/people/jprice jprice@whoi.edu
time = 1.5 days

1
ηo = 100 m
0.5
η/ηo
600
400
200
0 500
−200
0
−400
north, km −600 −500
east, km
Version 7 May 30, 2013
Abstract: This is the second of a three-part introduction to the effects of Earth’s rotation on the fluid
dynamics of the atmosphere and ocean. Part 1 derived the Coriolis force, ∝ − f k×V, where f is the
Coriolis parameter, k is the vertical unit vector, and V is horizontal wind or current, and went on to
analyze some of its basic properties in the context of a very simple, single parcel model. The goal and
plan of this Part 2 is to develop further insight into the consequences of the Coriolis force by analyzing
1
a sequence of geostrophic adjustment experiments in which a mass anomaly of horizontal scale L is
released from rest and allowed to evolve under the influence of gravity and the Coriolis force. These
experiments are posed in a single layer fluid model, often called the shallow water model.
The initial gravitational slumping produces gravity waves (see the cover graphic). If there is no
rotation, i.e., if f = 0 in a model experiment, gravity waves will disperse the mass anomaly in a time
L/C, where C is the gravity wave speed. When rotation is present, and if L exceeds several times the
radius of deformation, Rd = C/ f , and if the eddy has a potential vorticity anomaly, then the Coriolis
force will arrest the slumping and yield a geostrophically balanced eddy (modified by curvature), an
anti-cyclone if there is a mass excess and high pressure. If f is presumed constant and if there is no
friction, then a geostrophically balanced eddy could be exactly steady.
These and other low frequency phenomenon are often best interpreted in the context of potential
vorticity conservation, the geophysical fluid equivalent of angular momentum conservation. Earth’s
rotation contributes planetary vorticity = f , that is generally considerably larger than the relative
vorticity of winds and currents. Small changes in the thickness of a fluid column will cause significant
changes in the relative vorticity of winds and currents.
Cover graphic: A snapshot taken 1.5 days after the start of a numerical experiment in which a
thickness anomaly of a dense fluid layer was released from rest onto an f −plane centered on 30o N. The
thin red circle expands at the rate of a pure, baroclinic gravity wave, roughly 300 km day−1 in this case.
This and other geostrophic adjustment experiments are described in Secs. 3 and 4. If you are viewing
online, click this link www.whoi.edu/jpweb/ga2d-h100-f.mp4 to download an animation.
Contents
1 Large-scale flows of the atmosphere and ocean 3
2 An ideal shallow water model 7

2.1 Mass and momentum balance of an ideal fluid layer . . . . . . . . . . . . . . . . . . . . 8
2.2 Solving and diagnosing the shallow water system . . . . . . . . . . . . . . . . . . . . . 14
2.2.1 Wave velocities and fluid velocities . . . . . . . . . . . . . . . . . . . . . . . . 14
2.2.2 Energy balance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
2.2.3 Potential vorticity balance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
2
1 LARGE-SCALE FLOWS OF THE ATMOSPHERE AND OCEAN 3
2.2.4 Linear and nonlinear; finite amplitude effects . . . . . . . . . . . . . . . . . . . 22

2.3 Boundary and initial conditions define the problem . . . . . . . . . . . . . . . . . . . . 23
2.4 Appendix to Sec. 2: Normal modes of a two layer model ocean . . . . . . . . . . . . . 25
2.4.1 Barotropic normal mode . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
2.4.2 Baroclinic normal mode . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
2.4.3 Reduced gravity approximation of the baroclinic normal mode . . . . . . . . . . 32
3 Gravitational adjustment, f = 0 34
3.1 Just gravity waves . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
3.2 An exact solution of the linear, one-dimensional wave equation . . . . . . . . . . . . . . 35
3.3 The choice of scales . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
3.4 Energy and potential vorticity balances . . . . . . . . . . . . . . . . . . . . . . . . . . 38
3.5 Finer details of the solution; finite amplitude effects . . . . . . . . . . . . . . . . . . . . 40
4 Geostrophic adjustment on a rotating, flat Earth, f = f0 42

4.1 Gravity waves plus inertial and geostrophic motions . . . . . . . . . . . . . . . . . . . . 42
4.2 Dispersion relation for waves on the f -plane . . . . . . . . . . . . . . . . . . . . . . . . 44
4.3 Intrinsic scales of the f -plane . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
4.4 The steady state inferred from potential vorticity conservation . . . . . . . . . . . . . . 48
4.5 If there is no potential vorticity anomaly . . . . . . . . . . . . . . . . . . . . . . . . . . 50
5 Geostrophic adjustment in two dimensions 53
6 Summary 56
7 Supplementary material 58
1 Large-scale flows of the atmosphere and ocean
This essay is the second of a three-part introduction to fluid dynamics on a rotating Earth. It is written
for students who have some preparation in classical dynamics and applied mathematics and who are
beginning a quantitative study of geophysical fluid dynamics. Part 1 examined the origin and
fundamental properties of the Coriolis force, and went on to consider a few of its consequences for the
motion of a single, dense parcel released onto a slope. The Coriolis force, often called simply ’rotation’,
admits two new modes of motion, a free oscillation usually called an inertial oscillation, and if there is a
steady or time-mean forcing, a steady or time-mean motion that is an analogue of geostrophic motion in
the atmosphere and ocean.
The broad goal of this Part 2 is to develop further insight into the consequences of rotation in the
context of a fluid model, called the shallow water model, derived in Sec. 2. The specific goal is to
understand: What circumstances lead to quasi-steady, nearly geostrophic balance? The plan is to
solve and analyze a sequence of geostrophic adjustment experiments in which a mass anomaly is
released into a rotating environment and allowed to evolve freely under the influence of gravity and
rotation.1 As we will see first in Sec. 3.2, a nearly steady, nearly geostrophic state will follow
spontaneously when a mass anomaly is released from rest (and so has a potential vorticity anomaly)
provided that the mass anomaly has a sufficiently large horizontal scale, about 100 kilometers or greater
at mid-latitudes. On the other hand, mass anomalies that do not have a potential vorticity anomaly
(because of their initial velocity) or mass anomalies having a smaller horizontal scale will be more or
less dispersed by gravity waves before reaching a geostrophic balance.
The discussion of the numerical experiments touches on four noteworthy topics or methods that
you will almost certainly encounter over and over again in your study of the atmosphere and ocean:
Wave Dynamics As previewed in the cover graphic, most of the phenomena that arise in these
experiments are wave-like. There are gravity waves that are more or less modified by rotation
depending upon their wavelength and f , and nearly geostrophic eddies2 that may also exhibit wave-like
propagation (when the β -effect is included in Part 3). The language and concepts appropriate to waves
— dispersion relations, phase speed, group speed — are thus essential for the description of the
experiments.
Linear and Nonlinear, Finite Amplitude Waves arise in elementary form implied above in linear
systems, or in the limit that amplitude goes to zero in the solution of a nonlinear system. On the other
hand, many geophysical phenomenon have an amplitude that is large enough that linear dynamics
would not appear to be strictly valid a priori. Nevertheless, linear theory will often be a very useful
1 Much of the pioneering research on the topics discussed in this essay appeared in a series of classic papers by Carl G.
Rossby and colleagues published in the late 1930s. A collection of Rossby’s highly readable papers is available online at
http://www.aos.princeton.edu/WWWPUBLIC/gkv/history/general.html
2 The term ’eddy’ is widely used in fluid dynamics, and with a wide range of meanings. Here eddy means a flow feature
having a more or less circular planform, but with no particular dynamics or time-dependence implied. To say that a movement
is ’propagation’ is suggestive of wave propagation, which is intended.
Aviso, SSH, 2007, week 40

o
60 N
2
o
45 N
1.5
o
30 N
1
o
15 N
0.5
0o
100oW 80oW 60oW 40oW 20oW 0o
Figure 1: A snapshot of SSH over the North Atlantic ocean from week 40 of 2007 (thanks to the Aviso
project). Compared with the year-long mean of Fig. (2), Part 1, this field shows considerable variability
on scales of several hundred kilometers, often termed the oceanic mesoscale (meso is Greek for middle)
especially in subtropical and higher latitudes. If you are viewing online you can animate the field by
clicking on this link: www.whoi.edu/jpweb/Aviso-NA2007.flv You will see that the mesoscale eddies
persist as identifiable features for many weeks. The amplitude of mesoscale SSH perturbations is about
± 0.1 m, and the corresponding thermocline displacements are about ± 50 m (not evident in this figure).
The specific goal of this essay is to understand how Earth’s rotation leads to long-lived eddies and gyres
that are in approximate geostrophic balance.
adjunct to more comprehensive (numerical) solution methods, even though linear theory may not be
adequate for every aspect of an analysis.
Potential Vorticity Balance The dynamics of low frequency, quasi-geostrophic atmospheric and
oceanic phenomena is revealed most effectively by analyzing the balance of potential vorticity,
essentially the angular momentum of a rotating fluid, Sec. 2.3.2, rather than by the balance of linear
momentum, which was all that was possible for the point-like parcel of Part 1. This use of potential
vorticity balance was introduced by Rossby and colleagues1 in the late 1930s, and remains a (and is
perhaps the) central concept applied to the analysis and prediction of large scale and low frequency
geophysical flows today. Potential vorticity balance is an essential tool for understanding the
phenomena discussed here and in Part 3.
Numerical Solutions The shallow water system used here is kept as simple as possible but is
nevertheless nonlinear, as are all more or less complete fluid models. Solving this nonlinear system,
even for an idealized problem, usually requires more than just pencil and paper, viz., a numerical model
that serves to time-step a finite difference representation of the shallow water equations. In this and in
many other ways, numerical models are an essential part of atmospheric and oceanic science, and for
example, nearly all of the PhD students in the MIT/WHOI Joint Program (Physical Oceanography) use
numerical models during their thesis research. For many of these students, that is also their first,
hands-on experience. An implicit goal of this essay is to introduce the practice of numerical modelling,
especially hypothesis testing and solution interpretation, at an introductory level. To that end, some of
the homework assignments suggested here will require the generation of new numerical solutions.
Source codes that are directly applicable to these problems are available on an anonymous ftp site
linked in Sec. 7. By introducing some of the methods (and limitations) of numerical research, this essay
is intended to supplement the excellent GFD texts by Cushman-Roisin, Gill or Pedlosky that treat many
of the same gesotrophic adjustment problems via mainly analytic solution methods.3
3 An introduction to geophysical fluid dynamics at about the level of this essay is by B. Cushman-Roisin, Introduction
to Geophysical Fluid Dynamics (Prentice Hall, Engelwood Cliffs, New Jersey, 1994). Somewhat more advanced and highly
recommended for the topic of geostrophic adjustment is A. E. Gill, Atmosphere-Ocean Dynamics (Academic Press, NY, 1982)
and for waves generally, J. Pedlosky, Waves in the Ocean and Atmosphere (Springer, 2003).
2 AN IDEAL SHALLOW WATER MODEL 7
Figure 2: Schematic of a single fluid

layer (light blue) sitting on a variable-
depth lower, solid boundary, e.g., the
sea floor. The horizontal dashed lines
are level surfaces (perpendicular to
gravity). This layer could be a homo-
geneous ocean, in which case the fluid
above is the atmosphere and δ ρ ρo .
In that case η is the displacement of
the sea surface (as in Figs. 1 and 2).
Or, it could be a relatively dense layer
within a stratified ocean, in which
case the fluid above is just slightly
less dense, δ ρ ρo , and η is the dis-
placement of a density surface. In ei-
ther case, it is presumed that the fluid
above the active layer does not im-
pose a horizontally-varying pressure
upon the active layer. Up until Sec.
4.4.3 b = 0, and hence the layer thick-
ness is just h = H + η .
2 An ideal shallow water model
The horizontal variability of winds and currents is of particular interest here, and hence the fluid model
is taken to be a single fluid layer that varies in the horizontal, (x, y), but not in the vertical, z.4 This layer
is presumed to have a nominal thickness, H, and rests upon a lower, solid boundary, e.g., the sea floor,
that is at a depth z = −b(x, y) (Fig. 2). The vertical displacement of the upper surface (relative to a level
surface) is η (x, y, t) and will vary with horizontal position and with time. The thickness of the layer is
h(x, y, t) = H + b(x, y) + η (x, y, t) in general, but from here until Sec. 4.4.3 only the case b = 0 (flat
bottom) will be considered. The fluid above is presumed to be at rest, and to have a uniform density ρo ,
and pressure, Pa , that is uniform horizontally. The fluid within the active layer has a greater density than
the fluid above, ρo + δ ρ where δ ρ is a specified constant noted below.
Compared with the single parcel model of Part 1, this single layer model is much closer to being a
realistic model of the atmosphere or ocean. There are, however, two important idealizations. The first is
that the horizontal velocity V is implicitly depth-independent since it is represented by a two
component vector field i.e., V(x, y, t) = [u(x, y, t), v(x, y, t)], vs. V(x, y, z, t), a bigger help than it may first
4 Readers who already know the shallow water model may, of course, skip this rather lengthy section, but should take a
look at Secs. 2.2 and 2.3.
appear. For this to be a valid approximation, the nominal layer thickness, H, must be somewhat less
than the horizontal scale of the motion, e.g., a wavelength, and hence a single layer model of this kind is
often and appropriately termed a shallow water model, even when H is taken to be the full depth of the
ocean. The most convincing justification for the shallow water model comes from the analysis of a two
layer model that includes both external (barotropic) and internal (baroclinic) normal modes (see the
appendix to this section). The two layer model makes clear that a shallow water model is appropriate to
a study of either barotropic or baroclinic normal modes, but not both at once. This is not a significant
shortcoming given the broad goal of this study and, moreover, it wouldn’t matter which normal mode
was chosen. But given that a specific goal is to model the westward propagation of oceanic mesoscale
eddies seen in Fig. 2, then as nominal values, H = 500 m and δ ρ = 2 kg m−3, are appropriate to the
baroclinic motions of the ocean’s main thermocline. Ocean mesoscale eddies have a radius L ≈ 100 km,
and hence this layer thickness is indeed ’shallow’ in the sense that H/L 1. It is notable that oceanic
eddies are mainly upper ocean phenomenon, i.e., strongest above the main thermocline. They certainly
are not resting on the sea floor, as does the single layer model derived here. The reduced gravity
approximation discussed in Sec. 2.4.3 shows that these two cases will be equivalent, provided that the
gravity wave speed of the single layer model is appropriate to the upper ocean phenomenon of interest.
The second important idealization is that the only physical processes allowed in are those of an
ideal fluid, viz., pressure and transport by the fluid flow, although an external stress F, is also allowed.
The wide range of physical processes associated with the thermodynamic properties of a real fluid, e.g.,
compressibility, diffusion and viscosity, are all omitted. This may seem a bit high-handed, but is
consistent with building the simplest model that will help us understand some of the consequences of
Earth’s rotation.
2.1 Mass and momentum balance of an ideal fluid layer
The start is a prescription for
I) Mass conservation in an ideal fluid: the mass of the layer at a point can change only by
virtue of mass fluxes associated with the (horizontal) fluid velocity within the layer.
This is the simplest statement of mass conservation that could be said to be fluid-dynamical in that the
consequences of fluid flow are accounted, though nothing else.5 To find the corresponding
5 Ifthe intent was to utilize the shallow water model to make a realistic simulation of a given, observed flow phenomenon,
then three approximations or simplifications would have to be valid: 1) that the initial state was free of vertical shear, ∂ V/∂ z =
Figure 3: A control volume defined

by thin dashed lines. The average
thickness is h̄, the length is d, and
the width in the y-direction is c (y is
normal to the page and is not shown).
The statement of mass conservation
appropriate to the problems consid-
ered here is that the mass of the ac-
tive layer within the control volume,
M, can change only because of mass
fluxes due to fluid flow normal to the
sides of the control volume, −ρ AV·n,
where A is the area. The layer thick-
ness variation shown here is exagger-
ated; changes of layer thickness of
even a few percent can be significant.
The F at top is an external stress, wind
stress or bottom frictional drag.
mathematical form, the mass balance will be evaluated over a control volume that spans the full
thickness of the active layer, Fig. (3). The mass of the layer contained within the control volume is
M = ρ h̄dc, (1)
where ρ is the constant density within the layer, and h̄ is the average thickness over the interval x to
x + d and c is the width (the dimension into the page). The mass flux due to fluid flow through a given
surface of area A is just
mass flux = −ρ AV·n, (2)
where n is the outward unit normal of the surface. For the left-facing side of the control volume in Fig.
(3), n = (−1, 0, 0) and the mass flux = ρ hcu, where u is the x-component of the fluid velocity and hc is
the area. These are evaluated at the x of the left-facing side. Noting that the unit normal on the
right-facing side at x + d is n = (+1, 0, 0), then the sum of the mass fluxes through both x-facing sides is
∂M
= ρ (h(x)u(x) − h(x + d)u(x + d))c. (3)
∂t
0, 2) that the fluid of interest was outside of frictional boundary layers, and 3) that the horizontal density variation within the
layer was effectively zero. Most real geophysical flows violate all three of these conditions to some degree, and so a realistic
simulation is likely to require a good deal more resolution in the vertical, e.g., a number of shallow layers stacked one on top
of another to represent vertical shear, boundary layers, stratification, etc. One layer is enough for now.
By our construction of mass conservation, this net mass flux is equal to the time rate of change of mass
within the control volume. The mass M(x, y, t) varies with position and time, but notice that the time
derivative of (3) is a partial derivative, since M is the mass of a control volume that is fixed in space.
Inserting the definition of M, and dividing through by ρ , c and d gives
∂ h̄ 1
= (h(x)u(x) − h(x + d)u(x + d)).
∂t d
Now let d go to zero, and the right-hand side also becomes a partial derivative
∂h ∂ (hu)
=− , (4)
∂t ∂x
the one-dimensional, differential statement of mass conservation (though since the constant density and
area are divided out, the units are length per time). Thus, the mass inside the control volume will
change if and only if there is a divergence of the mass flux associated with the fluid flow; the mass flux
alone is not relevant insofar as the M at a given position is concerned. Had variations in the y-direction
(normal to the page) been considered, then there would arise an identical term, −∂ (hv)/∂ y, and the
two-dimensional form of mass conservation

∂h ∂ (hu) ∂ (hv)
=− + = −∇·(hV). (5)
∂t ∂x ∂y
This way of writing mass conservation is sometime said to be the flux form or conservative form, and is
generally preferred for implementation in a numerical model.
Momentum conservation involves a little more, but follows along much the same lines.
II) Momentum conservation for an ideal fluid observed in a rotating reference frame: the
momentum of the layer at a given point, ρ hV, can change only by momentum fluxes
associated with the fluid velocity, by pressure variations, and by the Coriolis force, an
inertial force that acts throughout the body of the fluid. There may also be an external force,
F, (not to be confused with the lower case f ) that acts as a stress, a tangential force per area
on the upper or lower boundaries.
The x component of momentum (per unit area) is ρ hu, and the momentum flux due to the
x-component of the velocity is just ρ hu2 , which, notice, is proportional to kinetic energy. The pressure
flux can be computed given the known density and layer thickness, and presuming that the vertical
2
accelerations of the motion are very, very small compared to the acceleration of gravity, ∂∂ t 2h g. In
that case the vertical component of momentum balance is the hydrostatic pressure relation,
∂p
= −gρ ,
∂z
and the hydrostatic pressure within the layer (z < h) is
P(x, y, z) = Pa (z) − gρ0 z + gδ ρ (h(x, y, t) − z). (6)
Pa (z) is the z-dependent but horizontally uniform pressure in the fluid above the active layer and the
second term on the right is the hydrostatic pressure due to the constant density ρo . Neither of these
z-dependent terms contribute to horizontally-varying momentum fluxes and may be ignored in what
follows. The third term represents the pressure anomaly associated with the horizontally-varying layer
thickness, h and notice that a relatively large layer thickness indicates a high pressure anomaly within
the layer. The momentum flux due to the anomaly of hydrostatic pressure on an x-facing side of the
control volume having width c is then
h2
Z h
gδ ρ c (h − z)dz = gδ ρ c .
0 2
The x-component momentum flux on an x-facing side thus has two terms:
h2
momentum flux = ρ chu2 + gδ ρ c . (7)
2
The first term on the right is the flux of u-momentum due to the u-component of the velocity and so is
proportional to u2, and the second term is the layer-averaged anomaly of the hydrostatic pressure ∝ h/2
acting over the face of the control volume, ch. When the body force due to the Coriolis force is included,
−f × V = ( f v, − f u), and when the operations noted above to go from Eqn. (3) to Eqn. (5) are repeated
here, the conservative (or flux) form of the shallow water momentum (per mass × area) equations are
∂ (hu) ∂ (hu2 + 12 g0 h2 ) (huv)

=− − + f hv + Fx , (8)
∂t ∂x ∂y
∂ (hv) ∂ (huv) ∂ (hv2 + 21 g0 h2)

=− − − f hu + Fy, (9)
∂t ∂x ∂y
where g0 = gδ ρ /ρ is the reduced gravity (also buoyancy per unit volume). If the upper surface of the
layer is the sea surface, then δ ρ ≈ ρ , and g0 ≈ g; if the upper surface is within the comparatively weakly
stratified ocean, then δ ρ ρ and so g0 g. Unless assumed otherwise for modelling purposes, the f
of these equations varies with latitude. The consequences of this are the principal topic of this essay.
The external force shown here as F = (Fx , Fy ) could represent wind stress, in which case we would
probably say that F = F(x, y) is a specified function of space only. Or, it could represent a bottom drag,
in which case F ∝ −rV as in the single parcel model of Part 1. Notice that F must have dimensions of
force/mass, though it will be referred to as ’force’. Our main interest here is what the fluid itself does,
rather than F, and so unless it is noted otherwise, presume that F = 0.
For purposes other than numerical integration it may be useful to expand the divergence operator
and combine the partial time derivative and the advection term into the material derivative,
D( ) ∂ ( ) ∂( ) ∂( )
≡ +u +v , (10)
Dt ∂t ∂x ∂y
about which more below. The mass (or layer thickness, since ρ = constant) and momentum balances
written using the material derivative are

Dh ∂ h ∂h ∂h ∂u ∂v
= +u +v = −h + , (11)
Dt ∂t ∂x ∂y ∂x ∂y
Du ∂ u ∂u ∂u ∂h Fx
= +u +v = −g0 + fv+ . (12)
Dt ∂t ∂x ∂y ∂x h
Dv ∂ v ∂v ∂v ∂h Fy
= + u + v = −g0 − fu+ (13)
Dt ∂t ∂x ∂y ∂y h
The vector equivalents are
Dh ∂ h
= + V·∇h = −h∇·V (14)
Dt ∂t
DV ∂ V F
= + V·∇V = −g0∇h − f k×V + (15)
Dt ∂t h
It may be helpful to pause for a moment and consider how the momentum equation (15)
appropriate to a fluid layer compares with the equation of motion appropriate to the single parcel of Part
1, Sec. 5 (repeated here but omitting the bottom friction term from the latter and ignoring F of the
former),
dV
= g0 ∇b − f k × V.
dt
1) The Coriolis terms are identical in the two systems. The Coriolis force is an inertial force
that depends solely upon the fluid velocity relative to the rotation vector, and makes no
distinction of the physical properties of the material. And, of course, absent that
fundamental property, it would not have been appropriate to use a single parcel model as
the starting point for this study.
2) The pressure gradient term of the shallow water model, −g0 ∇h, has the form of the
buoyancy force on a dense parcel sitting on a slope, g0 ∇b (sign aside). A crucial difference
is that the bottom slope of the single parcel model is prescribed and fixed, while the
gradient of layer thickness in the shallow water model is a dependent variable, dependent
upon space and time via the continuity (mass balance), Eqn. (14). The continuity equation
has no counterpart in the single parcel model. In an ideal fluid model such as this one, the
variable pressure is the only way that a fluid parcel interacts with the rest of the fluid
domain. For example, if fluid begins to pile up (converge) at a given location, then the layer
thickness and the hydrostatic pressure will go up until the pressure gradient is sufficient to
push fluid away from that location. The rate at which the fluid responds to a convergence
determines the speed of gravity wave propagation (Sec. 3.1).
3) The time derivatives have quite different meanings implicit in the use of different
symbols: the ordinary time derivative d/dt of the single parcel model and the material
derivative D/Dt of the shallow water model. The former is the time rate change of a
specific, material parcel, which in a fluid dynamics context would often be termed a
Lagrangian coordinate system. The material derivative is equal to the time rate of change
following a moving parcel at the instant the parcel is coincident with the spatial position
where D/Dt is evaluated. The independent coordinates of the shallow water model are fixed
spatial coordinates and time, which in fluid dynamics is often referred to as an Eulerian
coordinate system. The fluid present at a fixed location changes due to the fluid flow itself,
the advection process represented by the term, V·∇( ), and which is the nonlinear product
of two dependent variables, velocity times thickness gradient in Eqn. (14) or velocity times
the velocity gradient in Eqn. (15). Advection and the associated nonlinearity are right at the
heart of fluid dynamics and are a significant reason why numerical methods and numerical
models are central to atmospheric and oceanic sciences.6
6
What is the interpretation of d/dt = 0 in a Lagrangian system? How about ∂ /∂ t = 0 and D/Dt = 0 in an Eulerian
system? The material derivative and some aspects of advection are discussed in greater detail in ’Lagrangian and Eulerian
representations...’. Price, James F., 12.808 Supplemental Material, Topics in Fluid Dynamics: Dimensional Analysis, the
Coriolis Force, and Lagrangian and Eulerian Representations, http://ocw.mit.edu/ans7870/resources/price/index.htm
2.2 Solving and diagnosing the shallow water system
Equations (14) - (35) are a coupled set of nonlinear, partial differential equations in three dependent
variables, the two components of the horizontal velocity, V = (u, v), and the thickness, h. The
parameters g and δ ρ are constants for any given problem, and the Coriolis parameter f is a specified
function of the north-south coordinate, y, including in one experiment f = 0, about which much more
below.
The shallow water model is nonlinear because of the advection terms of the material derivative
noted above, e.g., the east-west advection of north-south momentum is u ∂ v/∂ x, which is the product of
two dependent variables. Solutions of a shallow water system including these nonlinear terms are
necessarily generated by numerical methods: the one-dimensional experiments shown here were solved
by the Matlab script geoadj 1d.m, and the two-dimensional experiments by the Fortran code
geoadj 2d.f, both of which are linked in Sec. 7.7
In common with any numerical model, these models produce (putative) solutions in the form of
very large data files, u(xi , yi , t j ), v(xi , yi , t j )..., where xi , ti are discretized position and time. To verify that
this mass of data is faithful to the model equations and to the boundary/initial conditions, and then to go
on and learn something useful requires as much thought and effort as does generating the model and
solution in the first place. To wit, we will seek to 1) construct useful visualizations of the solutions (this
can be the artful side of numerical modelling), 2) diagnose the balances of energy and potential
vorticity, 3) interpret the wave-like properties of the solution using the dispersion relation of the
corresponding, linearized system, and then 4) form and test hypotheses regarding parameter
dependence by conducting further numerical experiments (that will be the homework, actually).
2.2.1 Wave velocities and fluid velocities
A qualitative difference between a shallow water model and the single parcel model of Part 1 is that a
shallow water model supports wave motions — relatively fast-moving gravity waves will arise in all of
our experiments, and much lower frequency Rossby waves will arise when there is a north-south
variation of f (Sec. 4). If the layer thickness anomaly is a simple harmonic motion varying with x and t
only, η (x, t) ∝ ηo sin(kx − ω t), then a constant phase propagates at the rate of the
ω
wave phase speed, C p = , (16)
k
7 The details of the numerical methods are all-important in setting the efficiency and the accuracy of the solutions, but not
something that will be discussed here. An excellent, concise reference on numerical methods suitable for the shallow water
model is http://www.mathworks.com/moler/exm/chapters/water.pdf
where wave frequency is ω = 2π /(wave period), or waves per time, and the wavenumber is
k = 2π /(wave length), or waves per space interval. This wave propagation velocity is often readily
apparent (as in the cover graphic). In many cases the thickness anomaly will be the result of a
superposition of waves, i.e., waves of different wave numbers and frequencies, and in that case the
envelope of the superposition will move at the
∂ω
wave group speed, Cg = , (17)
∂k
This makes clear that the dispersion relation, ω (k), is going to be very important. The dispersion
relation is determined by the physics of the wave medium, specifically the relationship of restoring
force to wavelength. This differs greatly between pure gravity waves (Sec. 3.1), gravity waves in the
presence of rotation (Sec. 3.2), and Rossby waves (Part 3).8
The motion of the fluid itself (i.e., the material) is, in general, qualitatively different from the phase
speed of waves that may be propagating through the fluid; the fluid velocity is usually much less than
the gravity wave speed and the fluid velocity is often much greater than the Rossby wave speed. To see
the fluid motion one can simply plot the field of the instantaneous velocity vectors, which shows the
Eulerian fluid velocity, V(xi , yi , t j ) = (u(xi , yi , t j ), v(xi , yi , t j )), (18)
at the fixed positions (xi , yi ) and the times, t j . This Eulerian velocity is the direct output of the numerical
model, useful in itself, and the starting point for much else. The fluid velocity is generally proportional
to the amplitude of the motion; in the problems discussed here, V ∝ ηo , where ηo is the initial interface
displacement. Wave speeds, on the other hand, are independent of amplitude in linear approximation,
and generally only weakly dependent upon amplitude in nonlinear models (more on this below).
To see the transport of fluid over a long term we can compute the evolution of a tracer, say s, that is
carried along with the flow without in any way altering the flow, i.e., a passive tracer. This tracer may be
embedded in the fluid at the starting time, so = s(x, y, to ) in any way that will serve to highlight the
features of interest. This tracer is then presumed to be conserved following the flow,
Ds
= 0. (19)
Dt
The tracer value at a fixed location will thus change in time due solely to advection, expanding the
material derivative,
∂s ∂s ∂s
= −(u + v ) = −V · ∇s. (20)
∂t ∂x ∂y
8 If phase speed and group speed are not already familiar to you, then this very brief discussion will not suffice. An
excellent primary reference is Chapters 1 and 2 of Pedlosky’s ’Waves in the Ocean and Atmosphere’ (footnote 3). For a
quick refresher, you might take a look at the script twowaves.m (Sec. 7), which allows you to define an arbitrary dispersion
relationship between two waves that are then superimposed. The envelope of the superposition, and thus a wave form or pulse,
propagates at the group speed, Cg = δ ω /δ k, where δ is the difference between the two waves.
Diffusion of tracer is omitted in this ideal fluid model, though some inadvertent ’numerical diffusion’
will always be present in numerical solutions. Thus Eqn. (19) would hold exactly only in a perfect
numerical solution (and you may never see one).
The trajectories of discrete fluid parcels, or ’floats’, make a useful complement to the continuous
tracer. The trajectory of the ith float, (xi , yi), is found by integrating the (Eulerian) fluid velocity at the
moving location of the float,
Z t Z t
i
x (t) = xio + i i
u(x , y , t)dt and i
y (t) = yio + v(xi , yi , t)dt. (21)
to to
Which specific fluid parcel or float is being tracked in this way is given by the initial position, i.e., the
ith float starts at xio = xi (t = to). The float will likely be found between the discrete grid points of the
numerical model’s Eulerian velocity field and so in practice this will require some interpolation. The
trajectory is the fundamental dependent variable of a Lagrangian description (rather than the velocity as
in the Eulerian description). But if needed, the Lagrangian velocity may be computed from the
trajectories via
∂ xi (t) ∂ yi (t)
Lagrangian fluid velocity, Ui(t) = ( , ). (22)
∂t ∂t
Tracer evolution by advection Eqn. (19) and float trajectories by Eqn. (21) are very closely related, and,
for example, the tracer concentration at the moving location of a float should be exactly conserved
(numerical diffusion and tracking errors aside).
The relationship between the Eulerian and the Lagrangian velocities depends entirely upon the
length of time that a float is tracked, the t = t0 + ∆t of Eqn. (21). If ∆t is so brief that the float moves
only an infinitesimal distance compared to the scale over which the Eulerian velocity varies in space,
then the Lagrangian velocity will converge to the Eulerian velocity at the starting point. Much more
interesting is the case that ∆t is very long, in Sec. 4.1 a year, in which case the float may wander into
regions having very different Eulerian velocity, which is also likely changing in time. In that event, and
aside from special cases, the Lagrangian velocity will likely bear no simple relationship to the Eulerian
velocity.
2.2.2 Energy balance
Energy conservation is a fundamental physical law that can sometimes be of use in analyzing a fluid
flow, especially if, as here, dissipation and thermal (internal) energy may be neglected. In that case the
mechanical energy, E = KE + PE, is conserved, where E is the sum of kinetic energy (per unit mass),
KE = 12 hV2 , and potential energy (per unit mass), PE = g0 0η zdz = 12 g0 η 2 + const. To find the rate of
R
change of KE, take the dot product of the momentum equation with the velocity times thickness, and
the rate of change of PE is found by multiplying the continuity equation by the thickness anomaly times
reduced gravity. The mechanical energy balance is
DE
= −g0 ∇·(η hV) + F·V (23)
Dt
The flux term on the right-hand side is the product of pressure anomaly, thickness and velocity, i.e.,
pressure work (actually, pressure work rate, or power). For example, if the pressure anomaly is positive
on the (imaginary) boundary of a control volume where the velocity is directed outward, say, then the
fluid inside the control volume will do pressure work on the fluid outside of the control volume. The
energy within the control volume will thus decline, while the energy outside the control volume will
increase. The pressure-work term thus accounts for the outward energy transport associated with wave
radiation through an open boundary, for example. The fluid flow can also transport energy at the rate V,
a process accounted by the advection term, V·∇E, of the total derivative. Finally, note that an external
force F can either increase or decrease energy depending upon whether it has a component parallel or
anti-parallel to the fluid velocity.
2.2.3 Potential vorticity balance
Throughout Part 1, our perspective on dynamics was through the linear momentum balance (linear here
in the geometric sense), and rotation appeared by way of the Coriolis force. That was all that was
possible given a single parcel model. Momentum balance is always relevant, but may not always be
highly revealing. Given the present fluid model, there is another and complementary point of view,
angular momentum balance, that has proven immensely fruitful for understanding some of the most
important low frequency (frequencies less than f ) phenomena of geophysical fluid dynamics. Two main
reasons, to get a little ahead in this short story, are that 1) Earth’s rotation provides a very large
background angular momentum that is made visible by small changes in the thickness or latitude of a
fluid column, and, 2) the angular momentum balance amounts to a kind of filter that eliminates high
frequency gravity wave motions and so serves to highlight the processes that cause departures from
geostrophic balance.
Back a step or two..... to analyze the motion of a rotating, solid object, say a gyroscope, you might
begin by computing the linear momentum balance of the component pieces. This would require an
accounting of radial accelerations and internal stresses on each piece and would likely be a fairly
arduous task. Assuming that the gyroscope is not at risk of breaking up, then at some point you might
decide to take the structure for granted, and focus your effort toward analyzing the angular (azimuthal)
momentum balance of the gyroscope as a whole. As a first step you would define a coordinate system
that gave the most compact and least complex accounting of the moment of inertia of the gyroscope and
then consider the processes that cause the angular momentum to change with time. The physical
content of your angular momentum analysis would not be fundamentally different from the linear
momentum description, but it would likely be a great deal simpler in the same way that choosing
appropriate variables and an appropriate coordinate system can facilitate any mathematical analysis.
The same considerations apply to an analysis of a fluid flow: as we will see in examples here, an
angular momentum analysis and description will often (but not always!) be a good deal simpler and so
provide a great deal more insight than does the otherwise equivalent linear momentum analysis.
The fluid flow equivalent of angular velocity is the curl of the fluid velocity,
ξ ≡ ∇×V,
called the vorticity. If the fluid velocity is a three-dimensional vector, then the vorticity is also a
three-dimensional vector. In the special case of the shallow water model, the velocity varies only in the
two horizontal dimensions, and so the shallow water vorticity,
∂v ∂u
ξ = ∇×V = ( − )z , (24)
∂x ∂y
has a vertical component only and is effectively a scalar. You can visualize vorticity as the rotation of
small (but not quite point-like) two-dimensional parcels, e.g., cylinders, that make up the fluid.9 As we
will see, however, vorticity in a plane wave will arise from horizontal shear in the direction normal to
the velocity, e.g., ξ = ∂ v/∂ x or ξ = −∂ u/∂ y are quite possible, and so spinning cylinders are not
always apropos. If the direction of rotation is the same as Earth’s rotation, the vorticity is said to be
cyclonic (from the Greek kyklon, for circular motion). Cyclonic rotation is thus counterclockwise in the
northern hemisphere and clockwise in the southern hemisphere. Anticyclonic is the reverse. Notice that
vorticity has units of inverse time, or frequency, the same as the Coriolis parameter, f . There is an
important sense in which ξ may be compared to and even added to f , as discussed below.
The vorticity of a fluid is unlike the angular momentum associated with a solid object in that it is
defined at every point in a fluid, i.e., in principle there is a vorticity field that accompanies every fluid
flow (although it could be zero), just as there is a velocity field and a thickness field in the shallow water
model. The governing equation for this vorticity field may be found by taking the curl of the
momentum equation, ∂ /∂ x of the y-component minus ∂ /∂ y of the x-component. A consequence of
9 An essential resource for all students of fluid mechanics is the collection of educational films available online
at http://web.mit.edu/hml/ncfmf.html These are admittedly a bit old-fashioned, but nevertheless provide excellent vi-
sualizations of many key concepts, including vorticity and Eulerian and Lagrangian coordinate systems; these films
are very highly recommended. A more modern film collection, many of which emphasize rotational effects, is
http://planets.ucla.edu/featured/spinlab-geoscience-educational-film-project/
Figure 4: Some illustrations of q-conservation,

Dq/Dt = 0, by fluid columns that have a con-
stant volume but may have a variable thick-
ness, h and that are free of external torques,
∇×F = 0. (upper) If h is variable, then the
absolute vorticity, f + ξ , will change in pro-
portion to h. This is true whether the abso-
lute vorticity ≈ ξ , because f ξ as in small
scale or engineering flows, or if ξ f and
ξ is effectively zero, as in a gyre-scale flow,
illustrated in the lower panel. (middle) The
columns perched at upper and right on the
Earth are intended to show that the planetary
vorticity due solely to Earth’s rotation is given
by ∇×VΩ = 2Ω Ω ·n = 2Ωsin(latitude) = f . The
planetary velocity VΩ and thus the planetary
vorticity f are apparent to an inertial observer,
implicit in this view from space. We Earth-
bound observers will see only the relative ve-
locity V, i.e., winds and ocean currents, and the
relative vorticity, ξ . The planetary vorticity is
not directly observable (though the stars turn-
ing overhead are a mirror image). The three
columns in the North Atlantic are shown end-
on to illustrate the sense of relative vorticity
that would be acquired by a q-conserving col-
umn that was displaced north or south away
from a reference site where ξ = 0 while hold-
ing h constant. This and the next mechanism
involve f (latitude) that will be taken up in Part
3. (lower) If the relative vorticity ξ is much,
much less than f and remains effectively zero,
as in a gyre-scale flow, then a change in the
latitude of a column will be accompanied by
a change in the thickness of a column in the
sense shown here, i.e., smaller h in the direc-
tion of smaller | f |. The solid lines are depths
of the 27.0 and 27.4 potential density surfaces
along a plausible trajectory through the main
thermocline of the eastern subtropical North
Atlantic. These data are after Fig. 4 of Luyten
J. R. et al., 1983: ’The ventilated thermocline’,
J. Phys. Oceanogr. 13, 292 - 309.
applying the curl operator is that all of the forces that are derivable from a potential are eliminated, most
notably the pressure gradient in the shallow water model, i.e.,
∂ 2η ∂ 2η
∇ × ∇(η (x, y)) = − = 0.
∂ x∂ y ∂ y∂ x
A divergence term ∂ u/∂ x + ∂ v/∂ y will arise and may be eliminated using the thickness (continuity)
equation, (14). After a little further rearrangement, the result is a balance equation for a scalar, q, called
the potential vorticity,10
Dq ∂ q 1 F
= + V·∇q = ∇× (25)
Dt ∂t h h
where
f +ξ
q= (26)
h
Up until Part 3 the interest will be the case F = 0, in which case we have a q-conservation law
Dq
= 0. (27)
Dt
The terms in potential vorticity are:
ξ = ∇×V is often called the relative vorticity in this context, since it is the vorticity of the
relative velocity, i.e., the winds and currents observed from an Earth-attached reference
frame. In Part 1, Sec. 4 the relative velocity was denoted by V0 , though by now the prime
superscript has been dropped.
f = ∇×VΩ is the planetary vorticity of the planetary velocity, VΩ = Ω × X, due to Earth’s

rotation, Fig. (4), where X is the position vector (Sec. 4.3, Part 1). f may be written in
several forms, f = ∇×Ω Ω × X = 2Ω Ω · n = 2Ωsin(latitude),
f + ξ = ∇×(VΩ + V) is appropriately termed the absolute vorticity, since it is indeed the

10 You will feel much more at ease with the important energy and potential vorticity conservation laws if you will
take the time to derive them yourself. There are no special techniques or approximations required, though the full
shallow water equations entail a fair degree of manipulation. It will be helpful to start with the linear shallow water
equations (see the next section) and then add just one of the advection terms, say u∂ h/∂ x. And if you need to peek:
www.eos.ubc.ca/courses/eosc477/content/pv deriv.pdf Performing this derivation does not aid applications, but it does help
answer two important questions. 1) Does the linear shallow water system yield a self-consistent energy and vorticity balance?
2) Does the q conservation law include all of the phenomena of the shallow water system? How about energy conservation?
(This second question will be answered in Sec. 3.1.)
vorticity of the absolute velocity, VΩ + V, i.e., the velocity that would be observed from an
inertial reference frame, and finally,
h is the layer (or column) thickness.
The conservation equation (27) states that q is conserved following an ideal fluid parcel (or fluid
column since this is a single layer model) aside from external torques due to wind stress or bottom
friction. What is most important is to notice what is missing: there is no process or term comparable to
the pressure work term found in the energy balance, Eqn. (23), that transmits energy at the speed of
gravity waves; q is transported with the fluid velocity, exactly as a passive tracer, Eqn. (19), and not at
the gravity wave velocity as is energy. Thus the potential vorticity of an ideal fluid has the conservation
property of a tracer (or dye) that, once put into a fluid, remains with the fluid no matter how complex
the flow may be. This includes during the process of geostrophic adjustment. But unlike a passive
tracer, the potential vorticity includes an important part of the velocity field, the part having relative
vorticity. In some cases it may be sufficient to calculate the evolution of a flow from q-conservation
only, the steady state of an adjustment process is an important example in Sec 4.4, which is a marked
simplification over solving the full shallow water system. However, the pure gravity wave experiment in
Sec. 3.1 will show that q conservation can be irrelevant if the dynamics is such that there is no vorticity.
Potential vorticity is a generalized angular momentum insofar as it accounts for a variable moment
of inertia (variable thickness) as well as the planetary vorticity f due to Earth’s rotation. The planetary
vorticity is extremely important in this regard, because f is considerably larger than is the ξ of large
scale flows, usually by a factor of 10 and often much more. A fluid column having potential vorticity
will exhibit a kind of gyroscopic rigidity in the sense that it will respond to changes in the parcel
configuration (thickness) or to changes in latitude and thus f . For example, suppose that a fluid column
has absolute vorticity ( f1 + ξ1) and thickness h1 in an initial state, and is then stretched from h1 to h2 .
Assuming that the stretching occurs without frictional or external torques, then the conservation law
(27) applies and the absolute vorticity of the column will change to ( f2 + ξ2) = ( f1 + ξ1)h1 /h2 (Fig. 4).
Whether the change in absolute vorticity is due mainly to a change in the relative vorticity or due to a
change in the planetary vorticity (latitude) can not be told without some additional information;
sometimes knowing just the horizontal scale of the motion will suffice (Sec. 2, Part 3).11
11 Some q-conservation problems for you: 1) A right cylinder having a radius r and height h has a moment of inertia
I = Mr2 /2 where M is the mass of the cylinder. The angular momentum due to rotation at a rate ω about the central axis is
L = I ω . Show that conservation of angular momentum under changes of h and r that are mass (volume) conserving can be
summarized with a vorticity conservation law like Eqn. (27). 2) Because f ξ , generally, a change in thickness of 10% is
often highly significant for the relative vorticity, as is a change in latitude of only a few degrees. Assume that a fluid column
having a radius of 50 km and a thickness of 700 m is moved from 40o N to 25o N and that it conserves q. If all of the change
You are probably wondering why we should make such a fuss over potential vorticity when the
numerical model that solves the more general shallow water equations. It is true that the winds or
currents might be inferred from a q-conservation argument can always be computed by the shallow
water model and so q-conservation is indeed redundant insofar as computation alone is concerned. The
value of potential vorticity becomes evident when it is time to describe and interpret (understand) a
numerical model solution or a set of field observations: a potential vorticity-based description will often
be far simpler and yield far more insight than would the corresponding momentum plus continuity
description. Of course, this presumes that potential vorticity concepts are a part of our working, fluid
dynamics vocabulary. We can take a useful step in that direction by using potential vorticity balance to
help describe and interpret the geostrophic adjustment experiments that follow in Sec. 4 and in Part 3.12
2.2.4 Linear and nonlinear; finite amplitude effects
The phrases ’nonlinear’ and ’finite amplitude’ are often used interchangeably. Here, linear and
nonlinear will be reserved for model equations, which are either one or the other with no gradation. The
shallow water system is nonlinear because of the six terms that are the products of dependent variables,
e.g., in the thickness equation (11) the terms u∂ h/∂ x and h∂ u/∂ x are nonlinear and so the model
system is nonlinear when these are retained, while ∂ h/∂ t is linear.13
in f is accounted by changes in ξ (no thickness change), estimate the magnitude of the resulting current. Now suppose that
the change in f is accommodated entirely by a change in thickness.....how much?
12 Before ending this discussion it should be noted that the shallow water potential vorticity is not the most general form of
potential vorticity, just as the shallow water model is not the most general fluid model. If the fluid velocity is a three component
vector, which it generally is in the ocean and atmosphere, then so too is the vorticity. The three-dimensional vorticity equation
includes a term that represents a change in the direction of the vorticity vector, e.g., from the horizontal into the vertical, often
called ’tipping’. If the fluid flow is baroclinic, which the atmosphere and ocean generally are, then there is an additional term
that arises from the cross-product of the density and pressure gradients, called the ’solenoidal’ term. The corresponding three-
dimensional and baroclinic potential vorticity that takes account of these additional processes is often called the Ertel potential
vorticity; it reduces to the shallow water potential vorticity when the flow is appropriately two-dimensional and barotropic. By
and large, the simpler shallow water potential vorticity will be adequate for an analysis of most large (horizontal) scale phe-
nomenon, e.g., gyre-scale flows and even mesoscale eddies, once they are formed. However, it will be entirely inadequate for
analysis of small scale phenomenon, e.g., boundary layer roll vortices and three-dimensional turbulence. Knowing where this
transition may occur requires a thorough understanding of the more general, Ertel potential vorticity. A valuable discussion
of some of these more advanced concepts is available online http://www.atm.damtp.cam.ac.uk/mcintyre/papers/ENCYC/epv-
times.pdf
13 Terms like u∂ h/∂ x that are the product of two dependent variables are sometimes said to be ’bilinear’ or ’semi-linear’.
Model equations having such terms are nonlinear in the most important sense that a superposition of valid solutions does
not yield a new valid solution. However, some solution methods are applicable to semi-linear systems, e.g., the method of
characteristics, that are not applicable to more comprehensively nonlinear systems.
Finite amplitude effects are the departures from linear dependence upon amplitude. Finite
amplitude effects do not arise in the solution of a linear model, and are a matter of degree in the solution
of a nonlinear model. For example, in the problems studied here, the amplitude is clearly determined by
ηo , the initial thickness anomaly. If the shallow water equations were linear (they aren’t) then the
solution η (x, y, t)/ηo would be invariant to all changes in ηo . Even in a nonlinear model, we can expect
that η (x, y, t)/ηo should be invariant in the limit of very small ηo (say ηo = 1 m with H = 500 m)
simply because the nonlinear terms noted above will be very, very small compared to the dominant
linear terms. The amplitude of the base case, ηo = 50 m, has been chosen to correspond with the
observed SSH signature of mesoscale eddies (Fig. 1). At this amplitude, there are several interesting
and important finite amplitude phenomena, including very large float (material) displacements, even
though a great deal of what we can see in a solution for η (x, y, t) and V(x, y, t) can be attributed to linear
wave dynamics (more on this in Part 3).
The linear subset of the shallow water model results when D( )/Dt is replaced everywhere by
∂ ( )/∂ t, and the variable layer thickness h is replaced by the constant initial thickness, H. In Cartesian
components the linear shallow water model is
∂h ∂u ∂v
= −H( + ), (28)
∂t ∂x ∂y
∂u ∂h
= −g0 + f v, (29)
∂t ∂x
∂v ∂h
= −g0 − f u. (30)
∂t ∂y
The stratification, represented in this one layer model by g0 and H and assumed homogeneous, is a
significant factor in the setting the phase speed of waves. The f is these equation is in general
dependent upon latitude (or north distance, y). The spatial variation of f determines the kinds of waves
that are possible.
2.3 Boundary and initial conditions define the problem
The mass and momentum equations (14) and (15) could be used to model a wide range of phenomena.
The definition of a specific, solvable problem follows from the specification of an initial condition on h
and V throughout the model domain, and the definition of boundary conditions along the edges of the
domain.
The problems studied here are variants on the classical problem of geostrophic adjustment1 and
specifically, geostrophic adjustment of a mesoscale eddy-like feature at mid-latitude. The stratification
is chosen to be typical of the subtropical main thermocline,
H = 500 m and δ ρ = 2 kg m−3 , (31)
and the resulting gravity wave speed,
g(δ ρ /ρo )H ≈ 3.1 m s−1 ,

p
C=
is that of an internal or baroclinic gravity wave (Sec. 2.4.2). The initial condition is either a
one-dimensional ridge with half-width L,
η (x, t = 0) = ηo if | x |≤ L, or else η (x, t = 0) = 0, (32)
(here in Sec. 3) or a two-dimensional eddy with radius = L (in Sec. 4).14 The half-width and the
amplitude of the initial interface displacement are chosen to be comparable to observed mid-latitude
mesoscale eddies (Fig. 1),
L = 100 km and η0 = 50 m. (33)
The initial velocity can be one of several forms, the default being a state of rest,
V(x, t = 0) = 0. (34)
If the initial velocity and vorticity vanish, then the initial ridge has a potential vorticity q = f /(H + ηo)
compared with f /H in the outlying fluid. This potential vorticity anomaly is preserved during
geostrophic adjustment, the full import of which will become clear in Sec. 4.5.
The computational domain will be either one-dimensional, having a width of several thousand
kilometers, where that is appropriate (Secs. 3 and 4, where f is presumed to be a constant), or
two-dimensional and 4000 km on a side (Sec. 5). There is no attempt to compute the fluid state outside
of the computational domain, and so something other than the momentum and continuity equations will
have to be imposed on the boundaries.
The only energy in the initial state is associated with the initial eddy which is placed near the
center of the model domain. It is then reasonable that waves will be radiated outward only, i.e., that
14 The initial ridge given by Eqn.
(32) has a very sharp edge, which in (numerical) practice is smoothed over a few horizontal
grid points. Nevertheless, when released suddenly, this initial condition tends to produce energetic gravity waves that are not
highly realistic of most natural phenomenon. Wind or tidal forcing acting on the ocean are by comparison rather slowly
varying in time and space, and so less apt to produce energetic gravity waves. The initial condition may be easily changed to
something smoother, and for that matter, the models may be readily configured to allow realistic wind or tidal forcing as in
Part 3.
nothing will come into the model domain from the outside. A plausible and generally effective
representation of this one-way, outward transfer is made by imposing a radiation boundary condition
along the boundaries:
∂ψ ∂ψ
= −Urad , (35)
∂t ∂n
where ψ is thickness or a velocity component, n is the direction normal and outward from the boundary
and Urad is the appropriate velocity component normal to the boundary. This amounts to imposing a
one-dimensional advection process normal to the boundary and at the speed Urad , which is obviously
going to be very important in what follows. The appropriate value of Urad depends upon the dominant
process local to the boundary, and that may change with time. Here, during the first 20 days of an
√
experiment, Urad is taken to be the gravity wave speed, Urad = C = g0 H, the fastest wave speed in the
shallow water system (Sec. 3.1). This serves well to usher along the gravity waves that first reach a
distant boundary. But after that comes trouble ... very, very slowly ... in the form of low frequency
Rossby waves (in Part 3).
2.4 Appendix to Sec. 2: Normal modes of a two layer model ocean
The shallow water model applied to a homogeneous (unstratified) ocean stands on its own. But if the
shallow water model is applied to internal or baroclinic motions (defined below), then some discussion
of the correspondence between observed and modelled layer thicknesses and phase speeds seems
necessary. In that vein it is useful to examine briefly the wave properties of the simplest two layer
model. For this purpose rotation and nonlinearity may be ignored, and the motion presumed to be in
one horizontal dimension only (x, t). The stratification is represented by two homogeneous layers, Fig.
(5), an upper layer 1, having an undisturbed thickness H1 and density ρ1 , and a lower layer, 2, with
undisturbed thickness H2 and slightly greater density, ρ2 = ρ1 + δ ρ . If this was meant to represent the
open ocean, then the interface between the layers would correspond with the middle of the main
thermocline and typical values would be H1 = 500 m, H2 = 3500 m, ρ1 = 1030 kg m−3 and
ρ2 = ρ1 + δ ρ with δ ρ = 2 kg m−3.
The interest here is in large scale, fairly slowly varying motions for which vertical accelerations are
very, very gentle, g, and the pressure is consequently hydrostatic, i.e., due to the weight of the fluid
overhead. The bottom pressure is Pb = g(ρ1 h1 + ρ2 h2 ), and the pressure within the layers varies as
P1 (x, z, t) = gρ1 (h1 (x, t) + h2 (x, t) − z) and P2 (x, z, t) = gρ1 h1 (x, t) + gρ2 (h2 (x, t) − z),
where h(x, t) is the space and time-varying layer thickness. The pressure gradient divided by the density
Figure 5: A two layer model appropri-

ate to the stratified ocean. The dashed
lines are level surfaces. The undis-
turbed layer thicknesses are H1 and
H2, and the actual thicknesses are h1
and h2. The upper surfaces of the lay-
ers are displaced to η1 = h1 + h2 −
(H1 + H2 ) and η2 = h2 − H2 . The
pressure of the atmosphere above is
presumed horizontally uniform and so
negligible.
is then in each layer,

1 ∂ P1 ∂ h1 ∂ h2 1 ∂ P2 ρ1 ∂ h1 ∂ h2
= g( + ) and =g +g . (36)
ρ1 ∂ x ∂x ∂x ρ2 ∂ x ρ2 ∂ x ∂x
With these results in hand, the momentum and continuity equations for the velocity and layer
thicknesses are (anticipating a matrix format):
∂ u1 ∂ h1 ∂ h2
+ g +g = 0,
∂t ∂x ∂x
∂ u1 ∂ h1
H1 + = 0,
∂x ∂t
ρ1 ∂ h1 ∂ u2 ∂ h2
g + +g = 0,
ρ2 ∂ x ∂t ∂x
∂ u2 ∂ h2
H2 + = 0. (37)
∂x ∂t
Notice that the thickness of a given layer can vary only if there is divergence within that layer, and, that
the pressure gradient has a dependence upon both layer thicknesses. Thus the pressure gradient couples
the layers together.
To find the wave properties of this system, presume that a wave may exist in both layers:
u1 (x, t) = U1cos(kx − ω t), u2(x, t) = U2cos(kx − ω t), h1 (x, t) = H1 + Γ1cos(kx − ω t) and
h2 (x, t) = H2 + Γ2cos(kx − ω t). The U1, Γ1 , U2, Γ2 are constant but to this point unknown amplitudes.
Notice that the wave frequency and wavenumber are the same in the two layers, which is implicit in ’a
wave’, and are also unknown. After substitution of this wave form, the governing equations (37) may be
written in matrices as
  
−ω gk 0 gk U1
 H k −ω 0 0  Γ1 
 
 1
   = 0.

ρ1
 0 gk ρ2 −ω gk  U2

0 0 H2k −ω Γ2
The simplest and most insightful description of a multi-part, linear system such as this will often
be in terms of the normal mode frequencies and structure, also called the eigenvalues and eigenvectors.
The normal modes are independent of one another, and hence the governing equations written in terms
of the normal modes will be decoupled, which is not the case with Eqns. (37) written in the layer-wise
thickness and velocity. To find the frequencies of the normal modes it is necessary to solve the
characteristic equation of the coefficient matrix (same as saying that the determinant must vanish),
which is a fourth order polynomial in ω
ρ2 − ρ1
ω 4 − gk2 (H1 + H2)ω 2 + g2k4 H1H2 ( ) = 0.
ρ2
This characteristic equation is bi-quadratic (the cubic and linear terms vanish) and can be readily solved
as a quadratic equation for ω 2 ,
q
2
gk (H1 + H2) g2k4 (H1 + H2)2 − 4g2k4 H1H2 ( ρ2ρ−ρ1 )
2 2
ω = ± . (38)
2 2
ρ2 −ρ1
In the oceanic case in which ρ2 1 this may be simplified by factoring the term under the square
root and then using the binomial theorem that (1 + ε )1/2 ≈ 1 + ε /2 when ε 1 as applies here,
gk2 (H1 + H2)

2 2 (H1 + H2 ) H1H2 ρ2 − ρ1
ω = ± gk 1−2 ( ) . (39)
2 2 (H1 + H2)2 ρ2
2.4.1 Barotropic normal mode
The larger of the two roots of Eqn. (39) is the frequency squared of the external or barotropic mode,15
15
The terms barotropic and baroclinic describe the dependence of density upon pressure. A barotropic fluid is one in
which the density can be written as a function of the pressure only, ρ = ρ (P). This would hold always in a fluid that was
homogeneous, for example, in the shallow water model. It would also hold in a fluid that was inhomogeneous and density
stratified, provided that the density surfaces and pressure surfaces were everywhere parallel. This would likely be true in
barotropic u barotropic η
Figure 6: Normal modes of a two-
0 0 layer ocean. (upper) The barotropic
z/(H1 + H2)
normal mode; current profile (up-

−0.5 −0.5
per left, arbitrary amplitude) and the
sea surface and interface displace-
ments (upper right, blue and green
−1 −1 lines). (middle) The baroclinic nor-
−1 0 1 −5 0 5
mal mode. The currents and the sea
baroclinic u baroclinic η surface and interface displacements
have a self-consistent structure, but
0 0 the thicknesses, H1 = 500 m and H2 =
z/(H1 + H2)
1000 m, and the density difference,

−0.5 −0.5 δ ρ = 300 kg m−3, were chosen
to make this structure legible rather
than realistic of the ocean thermo-
−1 −1
−1 0 1 −5 0 5 cline. (lower) The reduced gravity
approximation of the baroclinic mode
reduced gravity u reduced gravity η uses only the density interface to com-
pute the pressure gradient in the up-
0 0 per layer on the assumption that the
z/(H1 + H2)
lower layer currents and pressure gra-

−0.5
u=0
−0.5
∇P =0
dient vanish. This would not be a
good approximation for the stratifica-
−1 −1 tion shown above because H1/H2 is
−1 0 1 −5 0 5 not 1.
current Kx
2 H1 H2 ρ2 − ρ1
ωbtr = gk2(H1 + H2) − gk2 ( ). (40)
(H1 + H2) ρ2
Since the trailing factor involving the density difference is very, very small, ≈ 0.002, an excellent
approximation of this frequency is p
ωbtr ≈ ±k g(H1 + H2 ). (41)
a resting state, but not when motion causes vertical displacements of density surfaces that will often greatly exceed the
associated displacement of pressure surfaces. In that case, the fluid would be described as baroclinic, meaning that density
and pressure surfaces intersect, and that density varied with more than the pressure, i.e., with horizontal position and time at
a given pressure. Thus, barotropic is a special case in which density surfaces are always parallel with pressure surfaces, while
baroclinic is all else. The distinction is important in that the pressure gradient can (will) generate vertical shear in a baroclinic
fluid, but not in a barotropic fluid. Thus a shallow water model will be realistic of a barotropic fluid and flow. However, a
shallow water solution will require some interpretation if, as here, it is meant to represent a baroclinic phenomenon.
If the intent is to reproduce this dispersion relationship within a shallow water model, then the
single-layer equivalent gravity and layer thickness are
ge = g, and He = H1 + H2, (42)
which might have been guessed without help from a two layer model. The same kind of result is less
obvious for the baroclinic mode coming next.
By putting the appropriate frequency (41) back into the governing equations (37) we can solve for
the ratio of the amplitude of the layer thicknesses changes,16
ρ2 ρ1 1/2

Γ1 2 2
= (H1 + (H1 + H2 − 2H1H2 + 4H1H2 ) )/H2 − 1 . (43)
Γ2 2ρ1 ρ2
Noting that ρ1 /ρ2 is very close to 1, then as a very good approximation,
Γ1 H1
≈ . (44)
Γ2 H2
In this linear model the amplitude of either Γ1 or Γ2 is arbitrary, but the ratio Γ1 /Γ2 is determined by
the dynamics. Eqn. (44) indicates that the layer thicknesses oscillate in phase and with an amplitude
that is proportional to the undisturbed layer thickness. For some purposes it is helpful to know the
displacement of the upper surfaces of the layers; the sea surface, sometimes called the free surface, is
displaced from its resting (level) height by
η1 = (h1 + h2) − (H1 + H2),
and the density interface between the layers by
η2 = h2 − H2.
In the barotropic mode, the amplitude of the vertical displacement of the free surface compared to the
displacement of the interface is then, using Eqn. (44),
Γ1 + Γ2 H1 + H2
= ≈ 1.14, (45)
Γ2 H2
for our nominal, open ocean stratification. The vertical displacement is thus a maximum at the sea
surface and decreases linearly to zero at the bottom (Fig. 6, upper). The displacements η1 (x, t) and
η2 (x, t) have a cos(kx − ω t) time- and space-dependence that is common to both layers, which is what
distinguishes a mode from an arbitrary motion. The pressure gradient is due almost entirely to the
16 The algebra required for a three layer model is a bit tedious and is a good application for symbolic mathematics; see
twolayer eig.m linked in Sec. 7.
10
barotropic Figure 7: The dispersion relation
frequency, ω, hour−1
8 −1 for gravity waves of a two layer,

C = 200 m sec
nonrotating model of a stratified
6 ocean. The solid lines are the
full dispersion relation, and the
4 dashed red line is an approxima-
tion discussed in the main text.
2 baroclinic At a common wavenumber, the
C = 3 m sec−1 barotropic and baroclinic modes
have vastly different frequencies
0
−0.1 −0.05 0 0.05 0.1 and phase speeds.
−1
wavenumber, k, km
displacement of the free surface, and the currents are the essentially the same in the two layers (uniform
with depth). The density interface η2 moves up and down exactly as does the pressure at that level, and
thus the density could be written as a function of the pressure in the initial state and during the
subsequent motion. This kind of wave motion in which ρ = ρ (P) everywhere in the fluid is, as here,
often termed barotropic.
The phase speed of a barotropic wave for a nominal, open ocean depth H1 + H2 = 4000 m is very
fast,
ωbtr p
= g(H1 + H2) ≈ 200 m sec−1 ≈ 680 km hour−1 ,
Cbtr =
k
or comparable to that of a jet transport. High frequency (5 to 20 min period) barotropic waves of this
kind (often called tsunamis, Japanese for harbor wave) are the primary oceanic response to a rapid
vertical displacement of the sea floor (a few meters in a few minutes over a large horizontal scale).
Because the fluid velocity associated with open ocean tsunami waves is comparatively gentle, a few
centimeters per second, these waves undergo very little dissipation, and, aside from two-dimensional
spreading, may arrive on a far distant shore with a significant amplitude. Lower frequency barotropic
waves (periods 1/2 to 1 day) are the principal components of the open ocean, astronomical tides. These
near-daily frequency waves gravity waves are modified substantially by Earth’s rotation and are often
called inertia-gravity waves (Sec. 3.2, also called Poincare waves).17
17 An
excellent resource for tsunami waves is http://www.tsunami.noaa.gov/ A classic paper on open ocean tides is by
http://articles.adsabs.harvard.edu//full/1944MNRAS.104..244P/0000254.000.html
2.4.2 Baroclinic normal mode
The smaller of the two roots of Eqn. (39) is labeled the internal or baroclinic normal mode, and
s
H1 H2 ρ2 − ρ1
ωbcl = ±k g ( ). (46)
(H1 + H2 ) ρ2
Baroclinic gravity waves have a much lower frequency and phase speed than does a barotropic wave of
the same wavelength, s
ωbcl Cbcl g0 H1H2 1
= ≈ ≈ ,
ωbtr Cbtr g(H1 + H2)2 70
mainly because the reduced gravity, g0 = g(ρ2 − ρ1 )/ρ2 = g(2/1000), is very much less than the full
gravity, g. If the intent is to model the baroclinic dispersion relation using a shallow water model, then
H1 H2
ge = g0 and He = (47)
(H1 + H2)
will give an appropriate phase speed.18 It is usually the case that H2 exceeds H1 by a factor of 5 - 10,
and as a fair approximation, He ≈ H1. The baroclinic, long gravity wave phase speed is then, for the
nominal stratification,
ωbcl p 0
Cbcl = = g H1 ≈ 3.1 m sec−1 ≈ 270 km day−1 ,
k
the red line of Fig. (7).
The ratio of layer thicknesses in the baroclinic mode is

ρ2 ρ1 1/2

Γ1 2 2
= (H1 − (H1 + H2 − 2H1H2 + 4H1H2 ) )/H2 − 1 . (48)
Γ2 2ρ1 ρ2
Making the same approximation noted earlier for the barotropic mode,
Γ1 ρ2
≈ − ≈ −1, (49)
Γ2 ρ1
and hence the layer thicknesses oscillate out of phase, and with nearly equal amplitude (Fig. 6, middle).
There is a small but very important displacement of the free surface,
Γ1 + Γ2 ρ2 − ρ1
≈− = −0.0017,
Γ2 ρ1
18 Most authors define the equivalent depth from phase speed only, i.e., He = Cbcl 2 /g with g the full gravity. With that
definition, the baroclinic He ≈ 1 m (!). In other words, the baroclinic gravity wave mode has the phase speed of a 1 m thick,
homogeneous (barotropic) layer. In a linear shallow water model and problem, the phase speed is all that matters and this is
a useful definition. However, if the shallow water model includes finite amplitude effects and is intended to simulate realistic
amplitudes, then the actual layer thickness is relevant, and the equivalent values of Eqn. (47) seem more apt.
for the given stratification. An interface (thermocline) displacement of 50 m will thus be accompanied
by a free surface displacement of about 50 ∗ (δ ρ /ρo ) ≈ 10 cm, which is readily detectable to
satellite-based, altimetric methods, as in Fig. (1). The currents in the upper and lower layers are exactly
out of phase, and their ratio is such that the net transport vanishes.
2.4.3 Reduced gravity approximation of the baroclinic normal mode
The transport (depth integrated velocity) of the baroclinic normal mode = u1 H1 + u2H2 = 0. Thus the
upper and lower layer currents are in the ratio u1 /u2 = −H2 /H1. In the limit that H1/H2 → 0, the lower
layer current and pressure gradient are much, much less than in the upper layer. In that case, an
approximation of Eqn. (36) is that
∂ h2 ρ1 ∂ h1
=− ,
∂x ρ2 ∂ x
which may be used to eliminate the h2 term from the upper layer pressure gradient,
∂ P1 ρ2 − ρ1 ∂ h1 ∂ h1

=g = g0 .
∂x ρ2 ∂x ∂x
The density interface can thus be used as a proxy for the sea surface insofar as the pressure gradient in
the upper layer is concerned, provided that the full gravity that multiplies the sea surface slope is
replaced by the much smaller reduced gravity, g0 , multiplying the much greater slope of the density
interface. In that way all reference to the lower layer may be eliminated from the upper layer equations.
This is sometimes referred to as a one and a half layer model, where the half refers to the deep, resting
lower layer, and also called the reduced gravity approximation, (Fig. 6, lower). This is an appropriate
interpretation of the shallow water model when applied to a simulation of baroclinic, mesoscale eddies
which are mainly upper ocean (layer one) phenomena, and plausible also for modelling some aspects of
the wind-driven circulation of subtropical gyres (Fig. 8 and more in Part 3). The reduced gravity
approximation may also be used in a multi-layered model, so long as it is appropriate to approximate
the deep flow and pressure gradient as vanishing.19
19 Some modal questions for you: 1) Any arbitrary configuration of layer thicknesses may be decomposed into the normal
modes. What would you infer is the two layer model composition of the offhand sketch of layer thicknesses in Fig. (5)?
2) In place of two active layers, suppose three layers of equal thickness. What would you guess for the modal structure
of η ? Check your answer against the (numerical) eigenvectors of twolayer eig.m (linked in Sec. 7.1). 3) The reduced
gravity approximation is counterintuitive in that the density interface displacement is used to compute the hydrostatic pressure
anomaly in the layer above the density interface. Rather than try to explain that while developing the shallow water model in
Sec. 2.1, this valuable and sensible approximation (or interpretation) was deferred to this appendix. The essential, physical
connection between pressure and mass or thickness anomaly is that a comparatively thick upper layer is a region of high
pressure anomaly, as in Fig. (8). What do you think would happen if, due to a sign error in the code of a numerical model
(we all make them) this relationship was reversed?
Figure 8: A cross section of the North Atlantic subtropical thermocline, sliced east-west along
35oN, and viewed looking toward the north. The upper panel is SSH as in Fig.1, but here
the monthly average for September over about twenty years of measurement. The lower panel
is the long-term, September average of density along 35o N from the World Ocean Atlas 2001
(http://www.nodc.noaa.gov/OC5/WOA01/pr woa01.html). The tilt of the thermocline mirrors the tilt
of the sea surface so that high SSH corresponds to a thick, warm upper layer. A result is that the pressure
gradient and geostrophic velocity are comparatively small in the abyssal ocean, roughly 1500 - 2500
m depth, suggestive of a reduced gravity approximation (Section 2.4.3). Not shown in this figure are
bottom-trapped density currents found at depths in excess of about 3000 m that make up the lower limb
of the meridional overturning circulation. This figure was kindly provided by Iam-Fei Pun of WHOI.
3 GRAVITATIONAL ADJUSTMENT, F = 0 34
3 Gravitational adjustment, f = 0
The sequence of experiments described here and in Part 3 differ mainly in the way that Earth’s shape
and rotation are represented. Of course, we know that the Earth is almost round, and that the Coriolis
parameter,
f = 2Ωsin(latitude)
varies significantly with latitude.20 Nevertheless, it can be very useful to treat this latitudinal variation
by a series of approximations. In the first experiment, here in Sec. 3, the model domain is presumed flat
and not rotating, f = 0, so that the shallow water model is at it’s simplest, gravity waves only. Then in
Sections 4 and 5, f = fo = constant, as if the Earth was rotating but flat, in which case there are gravity
waves, inertial oscillations and geostrophic motion all at once. With the results from those two
experiments in hand, we will be ready for spatially varying f (latitude) in Part 3; first a mid-latitude case
and then an equatorial case. These latter two experiments will include all of the phenomena that came
before, plus a low frequency Rossby wave motion.
3.1 Just gravity waves
The experiment starts at t = 0 when the thickness anomaly of dense water is released. Within the first
tens of minutes the ridge begins to slump under the force of gravity, releasing potential energy and
generating currents (kinetic energy). In this experiment, the motions take the form of two equal,
outgoing solitary ’wave pulses’ of amplitude ηo/2 and width 2L. These pulses move at a steady speed
√
that is very close to the (baroclinic) gravity wave speed, C = g0 H, and hence they run off of the model
domain in just a few days. Other than these two discrete wave pulses, there is nothing. If you (like the
author) had expected the response to look something like the waves excited on the surface of a pond by
an errant golf ball, or tsunami waves generated on the surface of the ocean by a moving sea floor, then
this solution will seem strange indeed. This raises a string of questions, some that are specific to this
case and others that are surprisingly general. In the first place, how can we know that this numerical
solution is trustworthy? 21 These pulses look nothing like elementary gravity waves and so is it
20 A (semi-) trick question, good for an oral exam: given a fluid velocity, V, how does the three-dimensional Coriolis force
vary with latitude?
21 Looking at the plot of a numerical solution is a little bit like looking at Saturn through a new, homemade telescope
(Fig. 11, Part 1). Barring catastrophic errors, the largest eddies and the rings will be obvious and unmistakable. However, at
some level of detail an active skepticism is healthy, even essential to avoid over-interpretation or outright error. For example,
would you notice that Saturn is somewhat flattened if you didn’t anticipate it? How could you discern genuine flatness
from an optical distortion? One way would be to observe other, better-known phenomenon, e.g., the Moon, to calibrate the
instrument’s performance and extrapolate from there. Something like that can be done for numerical solutions by verifying
appropriate to call them ’wave’ pulses?
3.2 An exact solution of the linear, one-dimensional wave equation
There is a well-known, exact analytic solution for the initial value problem of the linear
one-dimensional wave equation that can serve as a very useful reference for some aspects of this
numerical solution. The wave equation of the linear, nonrotating shallow water system (Eqns. 28 - 30
with f = 0 and v = 0) is readily found by eliminating u in favor of η (recall that η = H + h with H a
constant and so ∂ h/∂ t = ∂ η /∂ t),
∂ 2η 2
0 ∂ η
=gH 2. (50)
∂ t2 ∂x
Given that the initial data is
η (x, t = 0) = ηo(x), (51)
then the D’Alembert solution
1 1
η (x, t) = ηo (kx − ω t) + ηo (kx + ω t), (52)
2 2
solves (50) and (51) provided that p
ω =k g0 H. (53)
The relationship ω (k) given by Eqn. (53) is called the dispersion relation, and is the crucial and
distinguishing property of a wave system.22 This is a particularly simple dispersion relation, ω (k) being
√
a straight line with slope g0 H, Fig. (7). Waves of all wavenumber thus have the same phase speed,
ω p 0
Cp = = g H =C (54)
k
and the same group speed,

∂ω p 0
Cg = = g H =C (55)
∂k
Because all of the waves in this system (that is to say, waves of all wavenumbers) have the same phase
speed, the system is said to be nondispersive. As a consequence, the initial shape, ηo (x), is retained in
that they follow adequately the conservation laws for energy and potential vorticity, among others, and that they reproduce
adequately known solutions.
22
Some important things have gone by rather quickly here. You should verify that 1) any function whose argument is
√
(kx ± ω t) satisfies the elementary wave equation provided that ω /k = g0 H, and 2) the initial data is satisfied by Eqn. (52)
for any ηo (x). However, if ηo H does not hold, then the assumption that the system is linear would not be appropriate.
Figure 9: An experiment in gravita-

tional adjustment (no rotation) solved
by the numerical model geoadj 1d.m
(Sec. 7.1). (upper) The nondi-
mensional (scaled) thickness anomaly,
η (x, t)/η0 , following the release of a
rectangular ridge that had an initial
amplitude η0 = 50m. The thin red
lines have a slope given √ 0by the grav-
ity wave speed, C = g H = 3.1 m
s−1 in this experiment. This solution
was computed by the numerical model,
geoadj 1d.m (Sec. 7). (lower) A snap-
shot of the solution at t = 1.2 days. The
green line shows the initial thickness
anomaly. The current, which is shown
by the array of vectors plotted above, is
3 o
these currents are rotated 90 ccw in the x-direction only and has been ro-
2.5 tated 90o counterclockwise to be made
visible here. The current is scaled with
η/η0 and U/(C η0/H)
2
Cηo/H (= 0.3 m s−1 for the case ηo =
1.5 time, days = 1.2 50 m). The blue dots and green crosses
1 at depth = -1.1 show the present posi-
tions and the initial positions of floats.
0.5
An animation of the lower panel is at:
0 www.whoi.edu/jpweb/ga1-lat0.flv
−0.5
−1
−500 0 500
X distance, km
the propagating wave pulses. Examples of dispersive wave systems arise in later problems. Notice that
the phase and group speed in this linear model are dependent upon the stratification g0 and H only, and
are independent of the amplitude of the motion.
Gravity waves are the only possible nontrivial motion in this system. There is nothing that picks
out a favored direction, and hence it is isotropic (from Greek iso + tropos, equal in all directions). Thus
by symmetry half of the initial ridge propagates in one direction, and half goes the other.
The current that accompanies the wave pulses is in the direction of the pulse propagation, i.e.,
toward positive x within the right-traveling pulse. The velocity and the pressure gradient within the
pulses is thus colinear (parallel or anti-parallel), a relationship often termed longitudinal in this context,
and that is characteristic of pure gravity wave motion. On the other hand, geostrophic motion is
transverse insofar as the velocity and the pressure gradient are approximately normal to one another
(Part 1, Sec. 5 and other examples coming soon). The observed velocity/pressure gradient relationship
thus provides a useful qualitative clue to dynamics at the level of momentum balance.
Perhaps the two most important points are that 1) comparison with the exact D’Alembert solution
shows that the numerical solution of Fig. (9) is qualitatively correct (more details below), and 2) the
isolated pulses seen in the numerical solution are gravity waves in all respects save their appearance.23
Their structure is that of the initial condition and is retained because of the nondispersive property of
shallow water (nonrotating) gravity waves combined with the one-dimensional geometry of this specific
problem.
3.3 The choice of scales
There is usually more than one way to plot a model solution (Fig. 9). It is highly desirable to use
coordinates that will allow for use by the broadest possible audience. This generally requires the use of
nondimensional coordinates, i.e., that independent variables, distance, time, etc. be reported in units
that are intrinsic (natural) to the model system rather than in meters or seconds. A curious feature of the
√
linear, nonrotating shallow water model, Eq. (50), is that the gravity wave speed, C = g0 H, is the only
intrinsic scale; there is no intrinsic horizontal length scale or time scale, as there will be when rotation is
included in Sec. 3.1. Dimensional scales that are convenient for the specific initial condition, kilometers
23
There is no doubt that this is true for the linear D’Alembert problem, but this is not proven for the wave pulses of the nu-
merical solution. Rather, this analysis leads to the (plausible) hypothesis that the wave pulses found in the numerical solution
have the properties and parameter dependence of elementary gravity waves. How can you test this, given the opportunity to
specify the parameters of new numerical experiments (geoadj 1d.m of Sec. 7)? For now, exclude finite amplitude effects by
keeping the amplitude small, ηo /H ≤ 0.1, say, and omit friction, which will be discussed below.
1 Figure 10: The energy balance for the

gravitational adjustment problem (no ro-
tation, Sec. 3.1) evaluated over a control
0.5 volume −500 ≤ x ≤ 500 km. Energy is
nondimensionalized with the initial po-
energy/PE0
tential energy, PE0 , and time is in days.

0 The first four curves shown here have
the same meaning as the curves with the
same name (and color) in the parcel on a
KE
−0.5 slope problem of Part 1, Sec. 5. Fric-
PE
KE + PE
tional work (FW) is vanishingly small
FW since the frictional parameter r was effec-
−1 PW tively zero in this experiment.
0 1 2 3 4
time, days
for horizontal distance and days for time, were therefore chosen for Fig. (9) and will be retained in all
later plots to facilitate comparison with this first experiment.
The scales for the dependent variables, η and u, were here chosen with the goal that the scaled
(nondimensional) variables should be independent of the initial amplitude, η0 , in the limit of small
amplitude. The obvious scale for the thickness anomaly is simply η0 , the initial value. Thus the
thickness anomaly is plotted as η /ηo . (To preserve the important sign of ηo it is necessary to use the
absolute value, η / | ηo |, though the plot legends do not show this.) The appropriate scale for the fluid
√
speed is less obvious. The gravity wave speed C = g0 H seems a promising candidate at first. However,
C depends only upon the fluid properties and is independent of the amplitude, here ηo , while the fluid
speed should increase roughly linearly with ηo . Hence, C alone is evidently not an appropriate scale for
the fluid speed. As a guess let’s try instead Cηo/H, i.e., plot u/(Cηo /H) (Fig. 9). Note that the current
amplitude is about 0.5 in these nondimensional units. If this is indeed an appropriate scaling, then the
nondimensional (or normalized) fluid speed should remain about 0.5 when ηo is varied within a small
amplitude regime.
3.4 Energy and potential vorticity balances
The energy balance (Sec. 2.3.2) was evaluated over a domain −500 ≤ x ≤ 500 km (Fig. 10). Several
features of this energy balance are in common with the parcel on a slope (Part 1, Sec. 5.3). First, the
only source of energy is the potential energy, PE, stored in the initial ridge, as is true in all of the
geostrophic adjustment problems that follow. Second, after the ridge is released and begins to slump,
the decrease of PE is accompanied by a closely comparable increase of kinetic energy, KE. Thus the
total energy, KE + PE, is approximately conserved in this ideal shallow water model, aside from small
losses (a few percent per week) due to numerical viscosity and to a small but numerically necessary
horizontal diffusion (discussed below).
Other aspects of this energy balance are quite different from that of the single parcel experiment.
Notice that once the outward-going wave pulses are fully separated, t ≥ 2L/C = 0.8 days, the kinetic
energy and potential energy are thereafter equal. This kind of energy equipartition is a characteristic of
(nonrotating) linear gravity wave motion and so provides a useful check on the numerical solution in the
much more common case that an exact analytic solution is not available. A second difference is the
domain over which the energy balance was diagnosed. There was no choice about where to evaluate the
energy balance for a single parcel; of course the balance was computed for the moving parcel. In fluid
dynamics parlance that would likely be called a Lagrangian (material-following) balance. Given the
present fluid model, there are a number of other ways that an energy balance could be evaluated. The
balance could follow a wave pulse, or it could follow any specific fluid column. Here, the energy
balance is evaluated over a fixed control volume, which is oftentimes called an Eulerian balance. As
one consequence, energy could leave the control volume through the sides beginning at t ≈ 1.6 days
when the wave pulses reach the edges of the control volume. This energy flux is accounted by the
pressure work term of Eqn (23).24
Potential vorticity is interesting in quite a different way — there isn’t any! There is no planetary
vorticity since f = 0, there is no relative vorticity in the initial condition since the ridge begins at rest,
and there is no vorticity produced by the subsequent gravity wave processes because the pressure
gradient in a shallow water model has zero curl.25 Potential vorticity analysis is often invaluable for the
analysis of low frequency phenomenon (next sections), but here is a reminder that potential vorticity,
and vorticity generally, is blind to shallow water gravity waves.
24 Friction was omitted from the discussion above, but it is included in the numerical model via a linear damping of the
velocity, −rV, sometimes called Rayleigh damping. In the solution shown here r was set to a small enough value that
frictional effects were negligible. Some friction questions: 1) What value of r is required to damp an outgoing wave pulse to
half initial amplitude as it arrives at x = 500 km? 2) How does the required r vary with the nominal layer thickness, H? 3)
Can you define a nondimensional number analogous to the Ekman number, E = r/ f (Part 1, Sec. 5) that serves to predict the
amplitude of this damping effect?
25 In a shallow water model the fluid density is presumed constant in space and in time. If instead the density varies in
space, then the pressure gradient divided by the variable density may have a curl, often called a solenoidal term. This term
will produce vorticity having a component normal to the plane of the pressure and density gradients. This baroclinic process
can be important in a more general three-dimensional fluid, but can not arise here.
3.5 Finer details of the solution; finite amplitude effects
On first seeing this solution the question came up — in what ways and to what degree can we trust this
numerical solution? If we were attempting to make a realistic simulation of an observed, physical
phenomenon (we aren’t) then we would have to consider the assumptions made in formulating the
shallow water model. Here the issue is much narrower .... how can we know that the plots in front of us
represent an acceptably accurate solution to the problem posed? Comparison of the numerical solution
with the exact D’Alembert solution suggests that wave propagation aspects are simulated reasonably
well in the sense that the pulses move at nearly the phase speed C, and without significant dispersion, as
they should in a small amplitude limit. Energy is also balanced fairly closely, though there is an
unaccounted loss of a few percent in the first several days.
A close look reveals two other, small but systematic differences between the numerical solution
and the exact, linear D’Alembert solution insofar as the numerical wave pulses do not retain the exact
shape of the initial ridge. First, the numerical η (x, t) shows a small scale spatial oscillation (wavelength
several times the numerical grid interval) that is characteristic of a finite-differencing error.26 A second
and more interesting difference is that the leading edge of a numerical wave pulse is steepened slightly,
while the trailing edge is elongated or rarefied. This appears to be a genuine, finite amplitude effect due
to nonlinear processes — horizontal advection and large variations in layer thickness — that are present
in the nonlinear numerical model and solution, but not in the linear D’Alembert solution.
It is very useful to be able to make an a priori estimate of the nonlinear terms and so to form a
hypothesis regarding how finite amplitude effects may vary with the amplitude, ηo. For a longitudinal
wave in which the wave and fluid velocity are colinear, u ∝ uo cos (kx − ω t), h ∝ ho cos (kx − ω t + φ ),
the ratio of the nonlinear thickness advection term to the linear local time rate of change may be
estimated as
nonlinear u∂ h/∂ x uok uo
= = = , (56)
linear ∂ h/∂ t ω C
the ratio of the fluid speed to the wave phase speed. This ratio arises frequently and is often called the
26 This error can be made larger or smaller depending upon the smoothness of the initial condition. If allowed to grow
unchecked, this grid scale noise can readily swamp the ’physical’ content of a numerical solution, i.e., the part that is faithful
to the model equations. A simple cure is to include the smallest possible horizontal diffusion of thickness and momentum that
will serve to damp high frequency variability. That has been done in these numerical solutions, though mention of it was left
out of the previous discussion. The issue then becomes whether this ad hoc diffusion impacts significantly the aspects of the
solution that are of interest. In this experiment, a thickness diffusion D = 5 m2 sec−1 was sufficient to prevent runaway growth
√ duration of this experiment, T = 5 days, this
of the grid scale noise, though not enough to crush it completely. Over the short
diffusion causes a small but acceptable spreading of the ridge judging from DT 2 ≈ 1.5 km L.
Froude number,
uo
F= (57)
C
Given the (tentative) scale estimate of the current amplitude, uo ≈ Cηo/H (Sec. 3.1), the Froude
number can be written in external parameters and easily evaluated (recall that for the base case, η0 = 50
m and H = 500 m) as
η0
F= = 0.1.
H
The other nonlinear term in the thickness equation may be readily factored into nonlinear and linear
terms, h(∂ u/∂ x) = (η + H)(∂ u/∂ x). Under the assumption that the variable η will be no larger than
the initial value ηo , then again
nonlinear η (∂ u/∂ x) η ηo
= = ≈ .
linear H(∂ u/∂ x) H H
The important ratio, ηo/H, called simply the ’amplitude’, happens to be equal to the Froude number in
this case, but not generally. These simple estimates indicate that nonlinear terms are small but not
negligible compared to comparable linear terms in the base case in which ηo /H ≈ 0.1. It is then
plausible, but not certain from this analysis, that there should occur modest but detectable finite
amplitude effects in the corresponding numerical solution.
A complementary view of finite amplitude effects comes from considering that the gravity wave
speed likely depends upon the actual fluid thickness, h = H + η and not just the nominal thickness as in
the linear wave speed, Eqn. (54). A wave pulse, which has one sign, might then propagate at a speed
given by the average thickness taking account of the pulse amplitude,
p
C pulse = g0 (H + ηo /4) ≈ C(1 + ηo/8H), (58)
in the specific case considered here. Taking this thickness-dependent view of wave speed one step
further, the observed steepening appears to be consistent with the idea that the thickest part of a wave
pulse should propagate faster than in the undisturbed fluid, which should in consequence produce a
steepened wave pulse front.27
27 You can investigate finite amplitude effects systematically by experimenting with a much smaller and a much larger value
for the initial thickness anomaly, say ηo = 2 and then ηo = 200 m (in the setup of the numerical model geoadj 1d.m). Some
questions: 1) It is expected that finite amplitude effects should increase with ηo . But what if the nominal layer thickness,
H, is increased apace with ηo so that the ratio, ηo /H remains constant? 2) Can you verify that wave front steepening
(and rarefaction) is large or small depending upon the Froude number, and conversely is independent of, e.g., the density
difference, δ ρ , when all else is held constant? 3) Does the scale for current speed, C ηo /h, account for the change in the
actual (dimensional) current amplitude at small values of ηo /H? How about at very large values of ηo /H? 4) Test the
estimated finite amplitude wave pulse speed Eqn. (58) using numerical solutions. 5) Instead of a ridge, suppose the initial
interface displacement is a trough, η0 < 0. Use the previous result 4) to explain the steepening and rarefaction seen in this
case.
4 GEOSTROPHIC ADJUSTMENT ON A ROTATING, FLAT EARTH, F = F0 42
4 Geostrophic adjustment on a rotating, flat Earth, f = f0
The one-dimensional domain and initial condition, including the initial state of rest, are retained in the
next experiment in which the Coriolis parameter is taken as a constant. The base case latitude is taken
to be 30o N, and so
fo = 2Ωsin(latitude) = Ω = 7.292 × 10−5 sec−1,
and the inertial period is 2π / fo ≈ 1 day (more precisely, 1 day less about 4 minutes). With f constant
and non-zero, this is essentially the classical geostrophic adjustment problem treated by Rossby.1
4.1 Gravity waves plus inertial and geostrophic motions
We can directly compare the new solution with the non-rotating experiment of the previous section, and
it is clear that the effects of rotation are very, very significant (cf. Fig. 11 and Fig. 9). There are still
gravity waves that propagate away from the ridge, and the wave front expands at a speed that is close to
√
the gravity wave speed, g0 H. The current that accompanies the wave front is polarized in the
x-direction and so the velocity-pressure relationship is longitudinal, as found for pure gravity waves.
The amplitude at the wave front is, however, much less than half the initial thickness anomaly seen in
the previous, non-rotating experiment. As time runs and waves continue to arrive in the far field, say
x = 300 km, the current associated with the waves changes to a semi-circular, clockwise-rotating
motion that approximates an inertial oscillation. The frequency of rotation is at first somewhat higher
than f and then decreases with time toward f (it is easier to see this frequency shift in the animation of
Fig. 11 or in the additional plots that accompany geoadj 1d.m).
After about five days, the gravity wave radiation away from the ridge is largely completed, though
some inertial motion persists. Most of the ridge remains in place, about 75% based upon a potential
energy balance (Fig. 14). There are counter-flowing currents, or ’jets’ since they have a well-defined
maximum, along the edges of the ridge (normal to the page). , and thus the velocity-pressure
relationship is transverse. Their direction is such that the Coriolis force acting on these jets tends to
compress the ridge, while the pressure gradient associated with the sloping layer thickness, ∂ η /∂ x,
tends to spread the ridge. When these opposing forces are balanced so that 0 ≈ f v − g0 ∂ η /∂ x at every
point, the ridge is in a nearly steady, geostrophic balance. Friction aside, an exact geostrophic balance
and so an exactly steady state is possible in this one-dimensional, f -plane system.28
28 How long would you say it takes for geostrophy to hold in this experiment? Recall from Part 1 that the time scale required
for rotation to deflect a moving parcel significantly was 1/ f ≈ 1/4 days at mid-latitudes. On the other hand, energetic near-
inertial oscillations evidently persist for roughly five days. Are you thinking about an instantaneous geostrophic balance, or a
Figure 11: An experiment in

geostrophic adjustment in one-
dimension on an f -plane. (upper) The
layer thickness anomaly, η (x, t)/η0 ,
of a rotating fluid layer following the
release of a rectangular ridge having
amplitude η0 = 50 m. This solution
was computed by the numerical model
geoadj 2d.m Time and horizontal
distance are shown in days and kilo-
meters to make a direct comparison
with the nonrotating problem shown
in Fig. 9. The thin red lines have
a slope given√ by the gravity wave
speed, C = g0 H = 3.1 m s−1 in
this case. (lower) A snapshot of the
solution at t = 1.2 days. The green
line shows the initial displacement
ηo . The current, which is shown by
the array of vectors plotted above,
is scaled with Cηo /H = 0.30 m s−1
for the parameters of this experiment.
3
The blue dots and the green crosses
2.5 at depth = -1.1 show the present and
η/η0 and V/(C η0/H)
2 initial positions of floats. The floats do

1.5 Rd time, days = 1.2
not move far in this experiment, and
time, IP = 1.2 finite amplitude effects are generally
1 small and uninteresting. Note the
0.5 counter-flowing jets along the flanks
0 of the partially subsided ridge. When
the ridge reaches geostrophic balance,
−0.5
the Coriolis force acting on these jets
−1 is just sufficient to restrain the ridge
−500 0 500 from slumping further. Animation:
X distance, km www.whoi.edu/jpweb/ga1-lat30.flv
4.2 Dispersion relation for waves on the f -plane
The waves in this rotating experiment are clearly very different from the pure gravity waves of the
nonrotating case. The dispersion relationship that links the wave frequency ω and the wave number k
makes a concise and very useful characterization of their properties. To find the dispersion relation,
postulate a plane wave solution of the form u(x, y, t) = U exp(i(kx + ky − ω t)) and similarly for v and η ,
where U,V, Γ are the constant but unknown amplitudes of the velocity componentsq and η . The kx and ky
are the components of the wave number vector which has a magnitude K = kx2 + ky2 = 2π /λ .
Substitution into the linear shallow water equations (28) - (30) and rewriting the resulting algebraic
equations in a matrix format yields
−iω − f ig0 kx
  
U
 f −iω ig0 ky   V  = 0. (59)
iHkx iHky −iω Γ
This homogeneous system has nonzero solutions for [U V Γ] if and only if the coefficient matrix is
singular, and so has a vanishing determinant,
ω 3 − ω (g0 HK 2 + f 2 ) = 0, (60)
which is to say, only if ω and K are related by a dispersion relation. The roots of (60)
ω = 0 and ω 2 = g0 HK 2 + f 2 , (61)
correspond to steady, geostrophic motion, and to a gravity wave motion modified by rotation, discussed
further below.
4.3 Intrinsic scales of the f -plane
The rotating, shallow water model has just three external parameters: the Coriolis parameter, f , that
defines a natural or intrinsic time scale, 1/ f ; the reduced gravity, g0 = gδ ρ /ρo , an acceleration; and the
layer thickness, H, a degenerate length scale since there is no variation over the layer. The initial
condition adds two lengths, ηo and L. The Coriolis parameter f is an obvious (inverse) time scale for
normalizing the frequency of the dispersion relation Eqn. (61), and thus a given frequency can be said
to be high or low compared to the inertial frequency, f . By that rearrangement to nondimensional form
the dispersion relation (61) for the nonzero root then reads
ω 2 g0 H 2
= 2 K + 1, (62)
f2 f
time-mean geostrophic balance?
in which the wavenumber K is scaled by a very important derived length scale,
C
Rd =
f
called the radius of deformation (Figs 12 and 13). Just as f −1 is the intrinsic time scale of a rotating
shallow water model, Rd is the intrinsic (horizontal) length scale. It is often the case that the initial
condition or external forcing bring an external length scale to the problem, in this case the ridge width L,
and the radius of deformation is the appropriate length scale for comparison. Thus a given ridge width
is large or small as L/Rd is greater than or less than 1.29 After this very simple but significant
rearrangement to nondimensional coordinates, the dispersion relation (Figs. 12 and 13) holds for all
values of f , g0 and H, i.e., for all shallow water, f -plane models.
The value of Rd depends upon the stratification through C and the latitude through f . Over the
open, subtropical oceans, C ≈ 3 m sec−1 and not highly variable (C ≈ 2 m sec−1 at subpolar latitudes,
however). The most significant geographic variation of Rd comes from the latitudinal variation of
f −1 ∝ 1/sin(latitude). Notice that the equator is a special and an especially important case that will be
considered in Sec. 4.2. In the ocean thermocline at 30o N, Rd ≈ 40 km. The Rd of the atmosphere at this
latitude is much larger, Rd ≈ 1000 km, because the atmosphere is both thicker than the ocean
thermocline, and more strongly stratified. This very large disparity in Rd is reflected directly in the
dominant horizontal scales of the variability seen in the atmosphere and ocean (cf. Figs. 1 of Part 1 and
Part 2), i.e., much larger horizontal scales in the atmosphere.
There are three important limits in the dispersion relation that correspond with modes of the
shallow water momentum equation. The root ω / f = 0 corresponds to exactly steady, geostrophic
motion (the blue plane of Fig. 13, upper). In this f -plane model, the corresponding wavenumber can
have any orientation and any magnitude (this will not be true when the latitudinal variation of f is
acknowledged in Part 3). At first this may seem a trivial solution, but it is instead the f -plane
approximation of the slowly-varying and nearly geostrophic motion that makes up most of the
atmosphere and ocean circulation. The non-zero ω / f root (the multi-colored parabola of Fig. 13, upper
and the similar line of Fig. 13, lower) corresponds to gravity waves that are modified more or less by
rotation. There is a short and long wave limit of inertia-gravity waves: in the short wave limit, as KRd
29 Three remarks. 1) This procedure has the somewhat formal (and shallow) feeling of a dimensional analysis. A more
compelling physical interpretation of Rd as a length scale will come just below. 2) The ratio of frequencies is a nondimensional
number sometimes referred to as the temporal Rossby number, Rot = ω / f , for which there seems to be no widely accepted
symbol. The nondimensional ratio KRd is the layered model equivalent of the square root of the Burger number, B = N 2 / f 2 ,
where N 2 = (g/ρ )(∂ ρ /∂ z) is the static stability, times the aspect ratio H/λ (not used here). 3) There is no guarantee that the
value 1 will always mark the boundary between large and small values of a given nondimensional number. Factors of π or
even 2π 2 are often relevant, and in any event, size is context dependent.
Figure 12: The dispersion relation for plane

waves supported by the rotating, shallow
water model. The frequency is normal-
ized by the Coriolis parameter, f , and the
4 wavenumber (an inverse length) by the in-
verse radius of deformation, discussed in
3 the text. The (gaudy) colors are indicative
of frequency only. A wavenumber vector K
t
that has an origin at (0, 0) and that termi-
ω/f
2
nates on either the bowl shaped surface or is
1 in the plane ω = 0 is a solution, i.e., a possi-
4 ble free wave or geostrophic motion on an f-
2 plane. This dispersion relation is symmetric
0 0
−4
about the ω axis (isotropic). The frequency
−2 −2
0 2
of waves depends upon the magnitude of the
4 −4 north, k R
y d wavenumber, K, but not the direction of the
east, k R wavenumber.
x d
becomes large, inertia-gravity waves asymptote to pure gravity wave motion in which the phase and
√
group speed are C = g0 H, as in the non-rotating, shallow water model (Sec. 3.1).30 Pure gravity waves
have the fastest group speed of all, and would be expected to arrive first in the far field, just as observed
in the numerical solution. In the long wave limit, i.e., as KRd becomes very small, inertia-gravity waves
asymptote to inertial motion. In this limit the phase speed is very large (infinite at ω / f = 1). The entire
field of motion consists of clockwise-rotating velocity having uniform phase, and thus no divergence
and no associated η or pressure anomaly. Exact inertial motion, like exact geostrophic motion, amounts
to a kind of frozen state that is unable to evolve. Consistent with this, the group speed in the limit
ω / f → 0 is zero. Exact inertial motion requires an infinite horizontal scale, and hence can not be
realized in a finite ocean basin, or for any number of other reasons, e.g., the latitudinal variation of f and
the existence of other ocean currents. While exact inertial motion is not expected in the ocean, rotary
currents that turn in the direction of an inertial motion and that have frequencies within 5-10 percent of
f — near-inertial motions — are very common, for example whenever there is a rapidly changing wind
stress on the ocean surface (an example was noted in Section 5 of Part 1).31 Notice that there is a low
frequency range, 0 < ω < f , within which there are no free motions possible in this system.
30
If the wavelength is less than the layer thickness, H, then the wave motion will become a short gravity wave that is
surface-trapped. These are the high frequency, wind-generated waves that make up most of the sea state. Their frequency can
be high enough that the pressure is nonhydrostatic and so these short wavelength gravity waves, sometimes called deep-water
waves, are not in the domain of the shallow water system.
31 Go back to the linear shallow water system Eqns. (28) - (30) and identify which of the terms are relevant in each of these
three limits or modes.

Figure 13: (upper) The dispersion rela-

tion for plane waves supported by the rotat-
f−plane waves and geostrophy ing, shallow water model. This is a slice
4 through the isotropic surfaces of the pre-
vious figure. The blue line at zero fre-
inertia−gravity wave quency is the steady, geostrophic solution
gravity wave
ω/f
2 and the multi-colored parabola is the dis-

inertial motion persion relation for inertia-gravity waves.
geostrophic The straight red lines are for pure grav-
0 ity waves ( f = 0); in these nondimensional
−4 −3 −2 −1 0 1 2 3 4
coordinates their slope is ±1 for right-
and left-going waves.(lower) The phase and
4 group speeds of the inertia-gravity waves
Cp/C and Cg/C
Cp, phase speed (solid and dashed blue lines) normalized

Cg, group speed by the gravity wave speed, C. Inertia-
2
gravity waves are dispersive in that their
phase speed varies with K. The group speed
0 goes to zero for very long waves, KRd 1,
−4 −3 −2 −1 0 1 2 3 4 that have frequencies near f , i.e., for near-
K Rd inertial motions. For short waves, KRd 1,
the group speed asymptotes to that of pure
gravity waves (the red line) and in that limit
Cg = C p = C.
Intermediate scale waves fall somewhere between the limiting cases of pure gravity wave and pure
inertial motion, and are appropriately called inertia-gravity waves (or inertia-gravity or rarely,
gravity-inertia waves). Their phase speed depends upon K. Hence the inertia-gravity waves of an
f −plane are dispersive, so that the shape of η0 (x) is not, in general, preserved as these waves propagate
away from a source (cf. the nonrotating, pure gravity wave experiment of Sec. 3.1). Instead, the initial
form η0 (x) will become more or less spread out, or ’dispersed’, over time as waves of different
wavelengths propagate at different C p.32
The dispersion relation is an invaluable guide to the properties of the waves that arise in a
geostrophic adjustment experiment, but it can’t tell us everything we might want to know. For example,
the dispersion relation alone gives no hint to the amplitude of any process, i.e., it does not tell how
much of the initial ridge will disperse into waves vs. survive into a geostrophic steady state. Notice too
that a given wavenumber in Fig. (12) can be either a time-dependent wave, or, a steady geostrophic
32 The
linear initial value problem can also be solved usefully via Fourier transform. An example is carried out by the script
ftransform.m (Sec. 7), which allows the choice of a dispersion relation. In the non-rotating gravity wave case, the initial
thickness anomaly ηo all goes into propagating waves. In the rotating case, a more or less significant fraction of ηo remains
in the geostrophically balanced end-state as we will see shortly.
Figure 14: Energy balance of the f -

plane geostrophic adjustment prob-
1
lem evaluated over −500 ≤ x ≤ 500
km. The energy and work terms are
nondimensionalized with the initial
0.5
potential energy, PE0 , and time is
in days. At 30o N, 1 day ≈ 1 iner-
energy/PE0
tial period. These data may be com-

0
pared directly with that of the non-
rotating experiment (Fig. 10) and
are remarkably different. In this ro-
−0.5
KE tating experiment, very little energy
PE was radiated out the sides of the do-
KE + PE main, PW KE + PE, and was in-
FW stead retained within the geostrophi-
−1 PW cally balanced ridge as potential en-
0 1 2 3 4 ergy, mainly, plus significant kinetic
time, days energy.
motion. What distinguishes these two very different kinds of motion? These important questions lead to
an analysis of potential vorticity conservation, described next.
4.4 The steady state inferred from potential vorticity conservation
A key result of this experiment is that most of the ridge survived the adjustment process, e.g., about
75% of the initial potential energy (Fig. 14). The conservation of potential vorticity provides real
insight into why this was the case in this experiment, and why it will not be in others. In words, the
q−conservation law Eqn. (27) states that q of a parcel (column) is unchanged by the process of
geostrophic adjustment, or,
q(α ) = qo (α ), (63)
where α is a parcel tag. Assuming that qo is known, then the task is to follow the parcels (all of the
parcels) during the adjustment process. However, if the amplitude is small in the sense that parcel
displacements are small compared to the scales over which the currents vary significantly, then a very
helpful approximation of (63) is that
q(x) = q0 (x), (64)
with x the usual (Eulerian) spatial coordinate. In other words, if the horizontal displacements during the
adjustment process are small, then q will be conserved in place, approximately, and q conservation is
linear. This may be checked by comparing the horizontal displacement of floats against the width of the
adjusted portion of the ridge. In the base case (Fig. 11) there is a detectable, outward displacement of
floats during the adjustment process, but the float displacement is small compared to the width of the
adjusted portion of ridge, the radius of deformation, discussed below. As well, the geostrophic currents
that form along the edge of the adjusted ridge are in a direction that is perpendicular to the thickness
and velocity gradients and so cause no horizontal advection. The result is that q is indeed conserved in
place, approximately, (Fig. 15, upper) as presumed in Eqn. (64). Note that the changes of v and η that
occur during adjustment are significant so that linear q conservation is by no means trivial (Fig. 15).
The potential vorticity in this one-dimensional (x-dependent only) case is just
f + ∂ v/∂ x
q= .
H +η
The initial state included the ridge within which the layer thickness was H + η0, and elsewhere the
thickness was H. The ridge was presumed to be at rest, and so the initial relative vorticity vanishes. The
q distribution in the initial and in the final state is then, within the initial ridge,
f + ∂ v/∂ x f
= if | x | ≤ L, (65)
H +η H + η0
and outside of the ridge
f + ∂ v/∂ x f
= if | x | > L. (66)
H +η H
Note that both v(x) and η (x) are unknowns, but in addition, we have come to expect that the steady
state should be in geostrophic balance, and hence we have also the geostrophic relationship between v
and η , viz.
∂ ηgeo
0 = f vgeo − g0 .
∂x
When this vgeo is substituted into the q conservation law, say for the region x ≥ L, Eqn. (66), there
follows a second order, ordinary differential equation in ηgeo(x):
g0 H d 2 ηgeo
− ηgeo = 0.
f 2 dx2
The boundary conditions on η may be applied over four segments, x ≤ −L, −L < x ≤ 0, 0 < x ≤ L
and x ≥ L. For example for the right-most segment, x ≥ L, the boundary conditions are
ηgeo(L) = η0/2 and ηgeo(∞) = 0.
The solution for this segment only is then

ηo x−L
ηgeo(x) = exp (− ). (67)
2 Rd
By matching the solutions of the four segments, ηgeo(x) may be constructed across the entire ridge with
a result that compares fairly well with the the numerical solution (Fig. 15, middle). Two notable
features of this solution:
1) The horizontal length scale is Rd : Compared with the initial ridge, the adjusted ridge is ’deformed’
over the horizontal e-folding scale Rd . The width of the adjusted portion of the ridge is then about 3Rd .
This result suggests a clear physical interpretation of Rd . Recall from Part 1, Sec. 5.3 that 1/ f is the
time it takes for rotation to turn a parcel (started from rest) by one radian with respect to an
impulsively-applied force. Rd is then the distance that a gravity wave will propagate in the time 1/ f ,
and hence Rd is proportional to the distance that the ridge will spread outward (as a gravity wave)
before being arrested by the Coriolis force. Implicit in this is a definite meaning to the phrase ’large
scale’ as applied to a mass anomaly of half-width or radius L: the appropriate standard against which to
compare (or measure) a horizontal length is Rd . The radius of deformation varies quite a lot with
latitude, as noted earlier, and is very large (though not infinite) as f vanishes. The nonrotating
(gravitational) adjustment experiment of Sec. 3 can now be seen as the limit Rd /L → ∞ in which
rotational effects are expected to be negligible — literally zero in that case — compared to the
spreading of an unbalanced thickness anomaly by gravity wave processes.
2) q-conservation and the geostrophic jets: The result Eqn. (67) gives an essentially complete
account of the geostrophic jets. The jets have a relative vorticity ∂ v/∂ x (Fig. 15, lower) that follows
from the effect of column thickness change and q conservation. For example, over the right-most
segment, x > L, the layer thickness increased (the water column was stretched) from H to H + η (x, t)
during the adjustment process (Fig. 15, middle). The initial potential vorticity was just f /H since there
was no initial velocity and no relative vorticity. To maintain constant q during the adjustment process,
the relative vorticity ξ = ∂ vgeo/∂ x = f η /H > 0 of the jet must thus be positive, or cyclonic, in this
region. Stretching of the water column appears to produce positive relative vorticity by concentrating
the planetary vorticity (Fig. 4, upper). To an Earth-bound, rotating observer, who sees stars turning
overhead but not the rotation due to planetary velocity, it looks as if this relative vorticity has come out
of nowhere. To an inertial observer for whom the same stars stand still, this looks like familiar angular
momentum conservation, albeit of a deformable object.
4.5 If there is no potential vorticity anomaly
Note that if the initial q was spatially uniform (no q anomaly) the geostrophic, q-conserving solution
(67) would be η (x) = 0, i.e., vanishing geostrophic steady state. In the experiment just described the
initial ridge was a q anomaly because the initial velocity was assumed to vanish and hence the q of the
Figure 15: Elements of the poten-

1.1 q0 tial vorticity balance for the f −plane
geostrophic adjustment experiment. (up-
qfinal
q/(f/H)
1 per) Initial (green) and final (blue) q(x)

(the latter is at time = 6 days). Note that
q of the numerical solution is almost un-
0.9
changed by the process of geostrophic
−200 −100 0 100 200
adjustment, as should have been the case,
despite that there are large changes in the
1 thickness and the current. (middle) The
η
0
ηfinal initial and final thickness anomaly (green
η/η0
0 and blue lines). The red line, ηgeo, is the

ηq conserv steady state thickness anomaly computed
−1 from Eqn. (67) under the assumption of
−200 −100 0 100 200
geostrophy and q conservation. (lower)
Time-mean current (blue vectors). The
0.5 velocity maxima, or jets, at x = ± 100
km = ±L, coincide with the adjusted part
V/(C η0/H)
of the initial ridge. These jets are in a

0 near geostrophic balance as evidenced by
a very close match with the geostrophic
velocity, (g0 / f )∂ η f inal /∂ x (dashed red
−0.5
−200 −100 0 100 200 line).
distance, km
initial ridge was f /(ηo + H) compared with f /H outside the ridge. A vanishing initial velocity is the
simplest initial state to define and a reasonable place to start experimenting, but it is not necessarily the
most realistic or common initial condition. An experiment that allows a more general initial condition
that includes velocity will make clear that insofar as the possible steady (adjusted) state is concerned,
the initial ridge is defined by q its anomaly, generally, and not by η alone, as it may have seemed from
this first rotating experiment.
For this more general case it is preferable to assume an initial ηo (x) that has zero mean (though
still called a ridge), e.g.,
η (x, t = 0) = η0 sin(π x/L) if | x |< L, and η (x, t = 0) = 0 otherwise.
Now consider one of three initial velocity fields that vanish outside of the thickness anomaly:
vanishing velocity, V (x, t = 0) = 0, just as before,
geostrophic velocity, V (x, t = 0) = (g0 / f )∂ η (x, t = 0)/∂ x, or,

Figure 16: Three experiments in

geostrophic adjustment in which the
initial thickness anomaly was a sin-
gle wave η0 sin(π x/L). The initial
velocity was either (top) vanish-
ing, as in the previous experiments,
(middle) geostrophically balanced,
in which event nothing happens
when the ridge is released, or (bot-
tom) the velocity that gave zero q
anomaly of the ridge with respect to
the outlying fluid. In this latter case
there is no geostrophic, steady state;
all of the initial ridge is dispersed
by inertia-gravity waves. These
experiments were carried out with
geoadj 1d uic.m, linked in Section
7. An animation of this last case is
www.whoi.edu/jpweb/zeroq1d.mp4
5 GEOSTROPHIC ADJUSTMENT IN TWO DIMENSIONS 53
Z x
vanishing q anomaly V (x, t = 0) = f η (x, t = 0)/Hdx.
−L
The first of these experiments (Fig. 16, upper) is very much like our previous experiment in that about
two thirds of the initial ridge survives into a steady, geostrophic state. The second experiment is a
non-event (Fig. 16, middle); nothing happens when a geostrophically-balanced ridge is released onto an
f -plane (this would not be true if the latitudinal variation of f was retained, Part 3).
The third experiment is the important and perhaps surprising one: a ridge having zero q anomaly is
completely dispersed by inertia-gravity waves, with nothing surviving into a steady state (Fig. 16,
lower). This holds regardless of the width of the ridge. Evidently the width of a ridge is best defined
from it’s q anomaly, and not by width alone, insofar as the existence of a steady state is concerned. This
last experiment, vanishing initial q anomaly, may seem a special, contrived case, but in fact it is very
common: tidal waves in the open ocean have large spatial scales, wavelengths of hundreds to thousands
of kilometers, and yet tidal waves propagate as almost free waves that show no sign of adjusting to the
kind of geostrophic balance that characterizes mesoscale eddies and gyres. The reason is that tidal
waves are generated by a small gravitational/centrifugal imbalance, the tidal force, that does not exert a
torque and so do does not by itself generate a potential vorticity anomaly. Thus the difference between
two identical wavenumbers (Fig. 12), one being a steady geostrophic motion and the other being a
wave, is evidently the association (or not) of a potential vorticity anomaly.
5 Geostrophic adjustment in two dimensions
Part 3 will consider the effects of spatially-varying f , and for that purpose it will be necessary to utilize
a two-dimensional model domain (and a numerical model written in Fortran, geoadj 2d.f, Sec. 7). Here
we will briefly note the (minor) consequences of going from one to two dimensions in an f -plane
geostrophic adjustment experiment. The initial thickness anomaly is taken to be a right cylinder having
a radius L = 100 km that is typical of mesoscale eddies, and the stratification and central latitude are as
before.
The initial slumping of the thickness anomaly produces gravity wave radiation away from the
center (cover graphic) that is familiar from the one-dimensional f -plane case, e.g., the wave front
√
expands at a rate very close to the expected maximum gravity wave speed, C = g0 H, and the
amplitude at the wave front is a small fraction of η0 . A difference is that the amplitude at the wave front
decreases with time due to geometric spreading ∝ 1/r.
Another difference in detail compared to a one-dimensional experiment is that circular motion with
azimuthal speed Uθ around an eddy center imposes curvature on parcel trajectories ∝ 1/R and thus an
outward (radially) directed centrifugal force. The steady, radial component momentum equation in
polar coordinates is (Sec. 3, Part 1), moving the radial pressure gradient to the left side,
Uθ2 1 ∂P
f Uθ + + = 0,
R ρ ∂r
coriolis + centri f ugal + pgradient = 0.
The Coriolis force can be positive (anti-cyclonic flow) or negative (cyclonic flow), the centrifugal force
is always positive, and the pressure gradient can have either sign. Assuming a given pressure term, this
may be solved as a quadratic equation for Uθ (Fig. 18). The magnitude of the centrifugal force may be
assessed from the ratio of the centrifugal force to the Coriolis force,
centri f ugal Uθ2/R Uθ

= = = Ro ,
coriolis f Uθ fR
an important nondimensional number often called the Rossby number. In this application the Rossby
number has the sign given by Uθ , e.g., for an anti-cyclone, Ro ≤ 0. Using this to rewrite the radial
component equation,
1 ∂P
f Uθ (1 + Ro ) + = 0.
ρ ∂r
In the case of a high pressure anti-cyclone, Fig. (17), the centrifugal force adds constructively to
the pressure gradient force, and in steady state, this combined centrifugal plus pressure gradient force
must be balanced by a centripetal Coriolis force. The magnitude of the azimuthal velocity in a steady,
high pressure anticyclone thus differs from the geostrophic speed computed from the pressure gradient
and f alone, i.e., as if strictly geostrophic, by a factor 1/(1 + Ro ) ≈ 1.04 in this experiment. The small
effect of curvature in this experiment is not apparent in Fig. (17) but see Fig. (19). A steady flow in
which geostrophy is modified by curvature is said to be in gradient balance or sometimes just
’balanced’.33
33
What if the eddy is a low pressure, cyclone? The discussion of curvature made allusion to centrifugal force acting on a
parcel and is thus inherently Lagrangian. Our shallow water model is Eulerian, and there must be an equivalent representation
of this process in the Cartesian form of the Eulerian momentum equations, Eqns. (12) and (13). Can you discern what term(s)
it is? Hint: consider a control volume (in 2-d) set on the northern edge of the eddy and estimate the net east-west advection
of north-south momentum into this control volume.
500
time = 25 days
Figure 17: Velocity and interface
displacement (color contours) from
a two-dimensional f -plane experi-
north, km
ment. This is a snapshot at 25

0 days. This eddy is a high in thick-
ness and pressure and so the bal-
−1
0.3 m sec anced flow is anti-cyclonic. The
big vector at lower left has a mag-
nitude 0.5Cη0 /H and serves as a
speed scale.
−500
−500 0 500
east, km
Figure 18: The red line shows the lo-

cus of possible solutions for steady,
gradient balanced flow over a very
large range of the Rossby number.
The schematic force balance for a
high pressure anti-cyclone like Fig.
(17) is at lower left; P is the pressure
gradient force (outward), F is the
Coriolis force (inward) and C is the
centrifugal force, always away from
the center of a curved flow. The an-
ticyclonic eddy of Fig. (17) has a
very small Rossby number, about -
0.04, so that C F.
6 SUMMARY 56
0.25
0.2 velocity
geo vel Figure 19: The geostrophic veloc-
0.15 ity computed from the thickness
0.1 field (green line) and the actual az-
imuthal velocity (blue line) in a
−1
speed, m sec
0.05
slice through the high pressure, anti-
0 cyclonic eddy of Fig. (17). This
−0.05 snapshot came from time = 25 days,
when the eddy was nearly steady.
−0.1
The actual velocity is just slightly
−0.15 larger than the geostrophic velocity,
−0.2 about as expected from the Rossby
number of this eddy, Ro = −0.04.
−0.25
−200 −100 0 100 200
x distance, km
6 Summary
In Section 1 we asked, What circumstances lead to a near geostrophic balance? and proceeded to
look for answers by analyzing several problems in geostrophic adjustment posed in a shallow water
model. The experiments started with a thickness anomaly having a horizontal scale L = 100 km and
amplitude, η0 = 50 m, that matched observed mesoscale eddies. Once released, the anomaly was free to
evolve according to the physical processes allowed.
1) When rotation is omitted by setting f = 0, the anomaly is quickly and completely dispersed in
space by the propagation of shallow water gravity waves. These waves are nondispersive, having
√
phase speed and group speed C = g0 H ≈ 3 m sec−1 for the stratification typical of the mid-latitude
oceans. These gravity waves are the only possible nontrivial free motions when f = 0. In the
one-dimensional, gravitational adjustment problem of Sec. 3, the solution consists of discrete,
propagating pulses that retain the shape of the initial ridge and so look nothing like elementary (sine)
waves.
A given experiment is characterized also by the amplitude of the initial thickness anomaly, η0/H.
If the initial amplitude is small, ηo /H 1, the pulses propagate at very nearly the expected gravity
wave speed, C. For a larger initial amplitude, say ηo /H ≥ 0.2, finite amplitude effects include
distortions (from linear) in the shape of the pulses, and a systematic change in the propagation speed
depending upon the sign of ηo .27
2) When rotation is included by an f -plane approximation, the transient response includes gravity
6 SUMMARY 57
waves that are more or less modified by rotation, and a possible steady state. The wave frequency
depends upon wavelength compared to the radius of deformation, Rd = C/ f , the intrinsic horizontal
length scale of a rotating, shallow water system. The highest frequency and shortest waves are
nondispersive gravity waves, while longer wavelengths are near-inertial rotary motions having very
large phase speed but very small group speed. In this long wavelength limit the waves are highly
dispersive.
3) If the initial ridge is a potential vorticity anomaly with respect to the outlying fluid, and if the
initial ridge is sufficiently wide compared to Rd , then some fraction of the ridge will survive the
adjustment process and reach geostrophic balance. In the case of a one-dimensional ridge, the
geostrophic current takes the form of jets whose half-width is Rd . In an inviscid f -plane experiment,
geostrophic balance can be exact, and exactly steady, i.e., could persist forever. (As we will see in Part
3, this is not the case when the latitudinal variation of f is also considered.) In a two-dimensional
experiment, an exactly steady balance is also possible, and may include a flow curvature effect
proportional to the Rossby number, U / f L, with U the current speed.
If a calculation of the final, geostrophically balanced state of an essentially linear experiment had
been the only goal, then the best route would surely be the analytic solution of Section 4.4 built upon q
conservation and geostrophy. On the other hand, if the need was to know the inertia-gravity waves that
accompany geostrophic adjustment or if the experiment was in a large amplitude regime, then it would
be necessary to solve the full shallow water momentum and continuity equations that contain
q-conservation plus gravity wave dynamics and nonlinear processes.
This concludes our discussion of the classical ( f -plane) geostrophic adjustment problem, and it is
fair to ask if we can now claim to understand the circumstances that lead to abnear geostrophic balance
of ocean gyres and mesoscale eddies. The rotating, shallow water geostrophic adjustment experiments
of Secs. 4 and 5 are a promising result in the sense that a mass and potential vorticity anomaly having
the width, stratification and rotation of a mid-latitude oceanic mesoscale eddy was found to adjust
spontaneously to a steady, geostrophic balance. This is a very robust result of the numerical model and
experiments insofar as it is not sensitive to numerical or physical details. Hence, we are lead to expect
that mid-latitude, mesoscale (or larger) mass anomalies of the kind seen in Fig. 2 should indeed be in
nearly geostrophic balance given their size and nominal stratification.34 But to understand geostrophic
balance in the fullest sense requires knowing the limits in parameter space (at what L/Rd , ηo/H etc.)
this result will not hold, and to know at least a little bit of what happens then. This is going to require
34 This leaves open why there are mass anomalies in the first place. This very important question goes well beyond the
present discussion, but see, for example, Tulloch, R., J. Marshall, C. Hill and K. S. Smith, ’Scales, growth rates and spectral
fluxes of baroclinic instability in the oceans’, J. Phys. Oceanogr., 41, pp 1057-1076, 2011. DOI: 10.1175/2011JPO4404.1.
7 SUPPLEMENTARY MATERIAL 58
some further experimentation and analysis.35 Running a new experiment is as easy as changing one of
the parameters of the physical problem, e.g., the latitude or the initial ridge width. To gain some
confidence in the numerical solution one might experiment with some details of the numerical model as
well, e.g., the grid interval, time step, etc. With some experience it will become evident that the
subsequent cataloging and interpretation of the result is best couched in terms of one of the derived
parameters — the radius of deformation and the potential vorticity — which offer a very significant
economy. To this point then, it appears that we have a solvable and relevant model of geostrophic
adjustment, and a few key concepts that will help interpret new model solutions and help us develop an
intuition for geostrophic balance.
What’s next? An understanding of geostrophy is an essential first or second step in a study of the
atmosphere and ocean, but be assured (or warned?) that there is much, much more to learn. In Part 3 we
will take up a fundamental question that comes from a second look at the large scale mass (pressure)
and circulation fields of the ocean and atmosphere, How do small departures from geostrophic
balance lead to time-dependence and to the marked east-west asymmetry of low frequency
phenomenon? The plan is to solve several further experiments in geostrophic adjustment that take
account of the often very important variation of f with latitude.
7 Supplementary material
The most up-to-date version of this essay plus the Matlab and Fortran source codes noted in the text
may be downloaded from the author’s public access web site:
www.whoi.edu/jpweb/aCt.update.zip
35 Using the Matlab scripts geoadj 1d.m and geoadj 1d uic.m to generate new solutions, 1) find and explain how the width
of a geostrophic jet varies with latitude and with H. What do you expect if the latitude is from the southern hemisphere?
Before trying these experiments in the numerical model, do a thought experiment and test your developing intuition against
the numerical solution. 2) Show that a plausible velocity scale for a geostrophic jet is C ηo /H. How/why does the jet speed
vary with H? Show that a reasonable a priori estimate of parcel displacement is C ηo /H f and that the condition for linearity
with respect to horizontal advection is then the familiar one, ηo /H 1. 3) The power of the potential vorticity analysis
becomes evident when one considers different kinds of initial conditions, say ηo = 0, but with an unbalanced current, u = uo
if | x |≤ L, and otherwise u = 0 (with u the velocity component in the x-direction). Then assume the same profile but with the
initial velocity directed in the y direction. Which of these initial states has vorticity, and thus a potential vorticity anomaly that
should remain in the adjusted state? Solutions for these initial velocity distributions can be generated by geoadj 1d uic.m. 4)
The solution shown here in Sec. 4 was made with a vanishingly small value of the linear damping coefficient. What happens
to the wave response and to the geostrophic currents when damping is made significant? Connect the results you obtain here
with the Ekman number developed in Part 1, Sec. 5.
7 SUPPLEMENTARY MATERIAL 59
Matlab and Fortran source codes include the following, all of which are public domain for all
educational purposes.
twowaves.m shows the result of superposing two sine waves whose dispersion relation may be
specified arbitrarily.
twolayer eig.m solves for eigenmodes of two and three layer models using symbolic math.
ftransform.m solves an initial value problem in one-dimension and without rotation. Shows the wave
forms that result from dispersion that is normal, anomalous or none.
geoadj 1d.m is a shallow water model used for the one-dimensional geostrophic adjustment
experiments of Secs. 3.1 and 3.2. The latitude, thickness anomaly, friction, etc. may be easily varied
from experiment to experiment. These all assume that the initial velocity vanishes.
geoadj 1d uic.m as above, but in this script the initial condition includes several choices for the initial
current, including vanishing potential vorticity.
geoadj 2d.for is a two-dimensional shallow water model used to do the geostrophic adjustment
experiment of Sec. 5. Written in a very simple Fortran. Solutions like those shown in Sec. 5 can be
generated on a reasonably capable PC in a few minutes. Model parameters including latitude, the kind
of f model, eddy amplitude, etc. are entered from the keyboard. Output is saved to disk.
galook.m is a Matlab script used to plot the data generated and saved by geoadj 2d.for.
The additional animations are:
ga2d eta latxx.mpg; animations of thickness anomaly, η , for latitude xx = 10, 20, 40, and 60oN.
ga2d u latxx.mpg; the currents that go with the η animations above.

A Coriolis Tutorial: Part 2: A Rotating Shallow Water Model, Potential Vorticity and Geostrophic Adjustment

Uploaded by

Copyright:

Available Formats

A Coriolis Tutorial: Part 2: A Rotating Shallow Water Model, Potential Vorticity and Geostrophic Adjustment

Uploaded by

Document Information

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

A Coriolis Tutorial: Part 2: A Rotating Shallow Water Model, Potential Vorticity and Geostrophic Adjustment

Uploaded by

Copyright:

Available Formats

a Coriolis tutorial

Part 2: a rotating shallow water model, potential vorticity and

time = 1.5 days

Version 7 May 30, 2013

2 An ideal shallow water model 7

2.2.4 Linear and nonlinear; finite amplitude effects . . . . . . . . . . . . . . . . . . . 22

4 Geostrophic adjustment on a rotating, flat Earth, f = f0 42

5 Geostrophic adjustment in two dimensions 53

1 Large-scale flows of the atmosphere and ocean

Aviso, SSH, 2007, week 40

100oW 80oW 60oW 40oW 20oW 0o

Figure 2: Schematic of a single fluid

2 An ideal shallow water model

2.1 Mass and momentum balance of an ideal fluid layer

The start is a prescription for

Figure 3: A control volume defined

and the hydrostatic pressure within the layer (z < h) is

P(x, y, z) = Pa (z) − gρ0 z + gδ ρ (h(x, y, t) − z). (6)

∂ (hu) ∂ (hu2 + 12 g0 h2 ) (huv)

∂ (hv) ∂ (huv) ∂ (hv2 + 21 g0 h2)

2.2 Solving and diagnosing the shallow water system

2.2.1 Wave velocities and fluid velocities

2.2.2 Energy balance

2.2.3 Potential vorticity balance

Figure 4: Some illustrations of q-conservation,

The terms in potential vorticity are:

f = ∇×VΩ is the planetary vorticity of the planetary velocity, VΩ = Ω × X, due to Earth’s

f + ξ = ∇×(VΩ + V) is appropriately termed the absolute vorticity, since it is indeed the

h is the layer (or column) thickness.

2.2.4 Linear and nonlinear; finite amplitude effects

2.3 Boundary and initial conditions define the problem

H = 500 m and δ ρ = 2 kg m−3 , (31)

and the resulting gravity wave speed,

g(δ ρ /ρo )H ≈ 3.1 m s−1 ,

η (x, t = 0) = ηo if | x |≤ L, or else η (x, t = 0) = 0, (32)

2.4 Appendix to Sec. 2: Normal modes of a two layer model ocean

Figure 5: A two layer model appropri-

is then in each layer,

gk2 (H1 + H2)

2.4.1 Barotropic normal mode

normal mode; current profile (up-

1000 m, and the density difference,

lower layer currents and pressure gra-

ge = g, and He = H1 + H2, (42)

η1 = (h1 + h2) − (H1 + H2),

and the density interface between the layers by

8 −1 for gravity waves of a two layer,

2.4.2 Baroclinic normal mode

The ratio of layer thicknesses in the baroclinic mode is

2.4.3 Reduced gravity approximation of the baroclinic normal mode

3.1 Just gravity waves

appropriate to call them ’wave’ pulses?

3.2 An exact solution of the linear, one-dimensional wave equation

and the same group speed,

Figure 9: An experiment in gravita-

3.3 The choice of scales

1 Figure 10: The energy balance for the

tential energy, PE0 , and time is in days.

3.4 Energy and potential vorticity balances