Coastal Ocean Processes
Hans Burchard, Ulf Gräwe and Lars Umlauf
October 14, 2015
Contents
1 Basic equations
1.1 Dynamic equations . . . . . . .
1.1.1 Reynolds decomposition
1.1.2 Shallow-water equations
1.1.3 Boundary condition . . .
1.2 Tracer equations . . . . . . . .
.
.
.
.
.
1
1
2
3
6
7
2 Frictional effects and Ekman theory
2.1 Pressure-driven, non-rotating flows . . . . . . . . . . . . . . . . .
2.2 Ekman theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
11
11
15
3 Vertically integrated shallow-water equations
21
4 Bottom gravity currents
4.1 Introduction . . . . . . . . . . . . . . . . .
4.2 Governing equations . . . . . . . . . . . .
4.3 Non-rotating flows . . . . . . . . . . . . .
4.3.1 Similarity solutions . . . . . . . . .
4.3.2 Formulation as a diffusion problem
4.4 The effect of rotation . . . . . . . . . . . .
4.5 Gravity current in a rotating channel . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
23
23
24
26
26
28
29
33
5 Entrainment
5.1 Introduction . . . . . . . . . . . . . . . . . . . .
5.2 Basic framework . . . . . . . . . . . . . . . . .
5.3 Evolution of the surface mixed layer . . . . . . .
5.3.1 Energetics . . . . . . . . . . . . . . . . .
5.4 Entrainment models . . . . . . . . . . . . . . .
5.4.1 Models based on energy partitioning . .
5.4.2 Models based on stability considerations
5.5 Further reading . . . . . . . . . . . . . . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
36
36
37
38
39
41
42
43
45
i
.
.
.
.
.
.
.
.
.
.
.
.
.
.
6 Surface waves in the coastal ocean
6.1 Observational evidence . . . . . . . . . . . .
6.1.1 Wave height and period . . . . . . .
6.1.2 The random-phase/amplitude model
6.1.3 The variance density spectrum . . . .
6.2 Small-amplitude wave theory . . . . . . . . .
6.3 Wave kinematics and pressure . . . . . . . .
6.3.1 The dispersion relationship . . . . . .
6.3.2 Phase velocity and group velocity . .
6.3.3 Wave induced pressure . . . . . . . .
6.4 Wave energy . . . . . . . . . . . . . . . . . .
6.5 Wave propagation . . . . . . . . . . . . . . .
6.5.1 Shoaling . . . . . . . . . . . . . . . .
6.5.2 Refraction . . . . . . . . . . . . . . .
6.5.3 Currents . . . . . . . . . . . . . . . .
6.6 Radiation stress and wave setup . . . . . . .
6.6.1 Radiation stress . . . . . . . . . . . .
6.6.2 Wave setup . . . . . . . . . . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
47
47
48
49
49
51
53
54
54
55
56
57
57
58
58
59
59
60
7 Shallow-water waves and tides
7.1 Tidal current profiles . . . . . . . . . . . . . . . . . . . . . . . . .
63
63
8 Estuarine circulation
8.1 Observational evidence . . . . . . . . . . . . . . . . . . . . . . . .
8.2 The Knudsen relation . . . . . . . . . . . . . . . . . . . . . . . . .
8.3 Idealised dynamics . . . . . . . . . . . . . . . . . . . . . . . . . .
69
69
74
75
ii
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
Chapter 1
Basic equations
1.1
Dynamic equations
The hydrodynamical processes discussed in this lecture will be based on a wellknown mathematical framework usually referred to as the Navier-Stokes equations. Not only in the ocean, these equations have been shown to successfully
describe the flow of so-called viscous Newtonian fluids under virtually all circumstances. Their derviation from material properties of viscous fluids, and the
basic laws of classical mechanics and thermodynamics is discussed in detail in the
lecture on “Hydrodyamics”, and in the following only some of the main results
are briefly summarized.
Here, the Navier-Stokes equations are formulated in Cartesian coordinates,
with the x1 -axis (or x-axis) pointing eastwards, the x2 -axis (or y-axis) pointing
northwards and the x3 -axis (or z-axis) pointing upwards. The origin of the z
axis is at the mean (or reference) sea level such that the bulk of the ocean is
associated with negative z-values. The components of the velocity vector in this
local coordinate system are denoted as u1 = u, u2 = v and u3 = w. Under the socalled Boussinesq approximation (density variations are taken into account only
in the buoyancy term), the Navier-Stokes equations on the rotating Earth can be
written as
1 ∂p
g
∂ 2 ui
∂ui ∂(ui uj )
+ 2ǫijk Ωj uk = −
− ρδi3 + ν 2 ,
+
∂t
∂xj
ρ0 ∂xi ρ0
∂xj
(1.1)
with the pressure p, the gravitational acceleration g, the density ρ, the constant
reference density ρ0 (typically set to 1000 kg m−3 ), the kinematic viscosity ν
(typically of the order of 10−6 m s−2 ), and the Earth rotation Ωj with components
Ω1 = 0, Ω2 = Ω cos φ and Ω3 = Ω sin φ, where φ is the local latitude, and
Ω = 2π/86164 s−1 the rotation rate of the Earth. The third term in (1.1) is thus
recognized as the Coriolis force.
1
The balance of mass is identical to that of an imcompressible fluid:
∂ui
= 0,
∂xi
(1.2)
where it should be clear that, with the Boussinesq assumption, this does not
imply that density is constant, as evident already from the appearence of the
buyoancy term in (1.1).
To close this system of 4 equations for the 5 state variables u1 , u2 , u3 , p and ρ,
the density is usually calculated from a so-called equation of state as a function of
pressure, temperature, and salinity. This requires additional transport equations
for the latter two quantities as discussed in more detail below.
1.1.1
Reynolds decomposition
In virtually all geophysical problems, the solution of the Navier-Stokes equations
is further complicated by the fact that the flows are generally turbulent, i.e.
involving stochastic fluctuations that are impossible to predict on longer time
scales. The standard procedure to attack this problem is the decomposition of
any flow variable X into a mean component, hXi, obtained from the ensemble
average of many (infinitely many, to be precise) identical flow experiments, and
a fluctuating component X ′ such that X = hXi + X ′ . After Obsborn Reynolds,
who first introduced this important concept, this decomposition into mean and
fluctuating parts is often referred to as the Reynolds decomposition.
By definition, the ensemble average denoted by the angular brackets has the
properties hX ′ i = 0 and hhXii = hXi, i.e. the average of fluctuations is zero, and
ensemble averaging does not change a mean quantity. Further, it can be shown
that the averaging operator obeys the following important properties:
1. Linearity:
hX + λY i = hXi + λhY i,
(1.3)
where λ is a constant.
2. Derivatives and averages commute:
∂X
∂hXi
=
,
∂x
∂x
∂X
∂t
=
∂hXi
,
∂t
(1.4)
3. Product average:
hXhY ii = hXihY i.
(1.5)
Inserting the decomposition into mean and fluctuating components for all
variables in (1.1) and (1.2), averaging the resulting expressions, and using the
2
averaging rules summarized in (1.3) - (1.5), the so-called Reynolds Averaged
Navier-Stokes (RANS) equations are obtained:
∂hui i ∂(hui ihuj i)
+
+ 2ǫijk Ωj huk i =
∂t
∂xj
(1.6)
1 ∂hpi
g
∂
−
− hρiδi3 +
ρ0 ∂xi
ρ0
∂xj
∂hui i
′ ′
ν
− hui uj i ,
∂xj
∂hui i
=0 .
(1.7)
∂xi
Since all terms except the advection terms are linear, only the appearence of
the latter is changed due to the following transformation: hui uj i = hui ihuj i +
hu′i u′j i. Thus, an additional divergence term appears, which has been merged
together with the existing diffusion term in (1.6). The quantity hu′i u′j i is called the
Reynolds stress tensor. It is symmetric, and hence introduces 6 new unknowns,
such that the mathematically closed Navier-Stokes equations have been converted
into an underdetermined set of equations for the mean quantities. The “closure”
of this problem is a complex topic in itself, and not discussed here in any detail
(interested students are referred to our lecture on “Marine Turbulence”). Some
very basic ideas for closing the system will be discussed below.
1.1.2
Shallow-water equations
The Reynolds averaged Boussinesq equations in (1.6) describe most geophysical
flow problems - but in many cases nearly identical results can be obtained also
with much simpler equations. The most important simplification in geophysical
applications is related to the fact that geophysical flows (oceans, lakes, atmospheres) are often very shallow, i.e. their vertical scale H is much smaller than
their horizontal scale L. Assuming H/L ≪ 1 (small aspect ratio), it has been
3
shown in the lecture on Hydrodynamics that (1.6) become
∂hui
∂hui
∂hui
∂hui
+ hui
+ hvi
+ hwi
− f hvi =
∂t
∂x
∂y
∂z
∂hui
∂
1 ∂hpi
′ ′
ν
+
− hu w i ,
−
ρ0 ∂x
∂z
∂z
∂hvi
∂hvi
∂hvi
∂hvi
+ hui
+ hvi
+ hwi
+ f hui =
∂t
∂x
∂y
∂z
∂hvi
∂
1 ∂hpi
′ ′
ν
+
− hv w i ,
−
ρ0 ∂y
∂z
∂z
(1.8)
∂hpi
= −ghρi .
∂z
where we have introduced the Coriolis parameter f = 2Ω sin φ for convenience.
These equation are usually called the shallow water equation. The most important change is identified in the vertical component, which has degenerated into
the so-called hydrostatic balance for the pressure. Note that in the derivation of
(1.8), we have also assumed that the horizontal Reynolds stresses are neglible.
This additional assumption does not directly follow from the shallow-water assumption because without knowledge about the horizontal and vertical turbulent
scales nothing can be said about the relative importance of vertical and horizontal Reynolds stresses. For the idealized flows considered in this lecture this
assumption does not impose any constraints but in more complex shallow-water
flows it is not guaranteed to be justified.
In analogy with the gradient relations for the molecular fluxes, the vertical
turbulent fluxes of momentum (or the vertical Reynolds stresses) appearing in
(1.8) are often expressed as
hu′ w ′i = −νt
∂hui
,
∂z
hv ′ w ′i = −νt
∂hvi
,
∂z
(1.9)
where it should be noted that the vertical turbulent diffusivity, νt , usually is
several orders of magnitude larger than its molecular equivalent ν. If νt would
be known, the equations would be closed. At the most basic level, one could
for example assume that νt corresponds to a prescribed constant, from which
some interesting solutions of (1.8) can be obtained (see below). It is important
to realize, however, that in essentially all relevant flows the turbulent diffusivity
exhibits are strong spatial and temporal variability that is inconsistent with the
assumption of constant νt . More complex and more generally valid models for
the computation of νt are available from turbulence theory but none of them has
so far proved to be universally applicable to all flows.
4
Inserting (1.9) into (1.8) and omitting the brackets denoting the ensemble averages (this convention will be used in all of the following, except noted otherwise)
leads to the following form of the three-dimensional shallow water equations:
∂u
∂u
∂u
∂u
+u
+v
+w
− fv =
∂t
∂x
∂y
∂z
∂u
∂
1 ∂p
(νt + ν)
+
−
ρ0 ∂x ∂z
∂z
∂v
∂v
∂v
∂v
+u
+v
+w
+ fu =
∂t
∂x
∂y
∂z
1 ∂p
∂v
∂
−
(νt + ν)
+
ρ0 ∂y ∂z
∂z
(1.10)
∂p
= −gρ .
∂z
The pressure p is typically eliminated from the shallow water equations by
vertically integrating the hydrostatic equilibrium from a position z in the water
column up to the surface η:
∂p
= −gρ
∂z
⇒ p(η) − p(z) = −g
Z
η
ρ dz
z
Z η
∂p(η) ∂p(z)
∂
⇒
ρ dz
−
= −g
∂x
∂x
∂x z
Z η
∂η
∂ρ
∂p(η) ∂p(z)
−
= −g
dz − gρ(η)
⇒
∂x
∂x
∂x
z ∂x
(1.11)
where in the last step Leibniz’ rule of differential has been used. With the atmospheric pressure at the sea surface denoted as pa = p(η) and the approximation
ρ(η) ≈ ρ0 (consistent with the Boussinesq assumption), the pressure gradient can
be formulated as
Z
1 ∂p
∂η
g η ∂ρ
1 ∂pa
−
= −g
−
dz −
,
(1.12)
ρ0 ∂x
∂x ρ0 z ∂x
ρ0 ∂x
where we have divided by ρ0 . Thus, the horizontal pressure gradient includes
contributions from (a) the sea surface slope (external or barotropic pressure gradient), (b) the integrated density gradient (internal or baroclinic pressure gradient), and (c) the atmospheric pressure gradient. The treatment of the pressure
5
gradient in the y-direction is completely analogous. After substitution of these
results into (1.10), the shallow water equations can be expressed in the following
form:
∂u
∂u
∂u
∂u
+u
+v
+w
− fv =
∂t
∂x
∂y
∂z
Z
∂η
∂u
g η ∂ρ
1 ∂pa
∂
−g
(νt + ν)
−
dz −
+
∂x ρ0 z ∂x
ρ0 ∂x
∂z
∂z
∂v
∂v
∂v
∂v
+u
+v
+w
+ fu =
∂t
∂x
∂y
∂z
Z
∂η
∂v
g η ∂ρ
1 ∂pa
∂
−g
(νt + ν)
−
dz −
+
∂y ρ0 z ∂y
ρ0 ∂y
∂z
∂z
1.1.3
(1.13)
.
Boundary condition
As discussed in more detail in the lecture on “Hydrodynamics”, at the boundaries
different types of boundary conditions have to be satisfied, depending on both
the flow kinematics and the material properties of the fluid. Kinematic boundary conditions are derived from purely kinematical arguments, and have to be
obeyed by all fluids, irrespective of their material properties. They follow from
the assumption that the boundary of the fluid consists of the same material fluid
particles for all times. For the geometries considered here, this is mathematically
expressed as
∂η
∂η
∂η
+ u(η)
+ v(η)
,
w(η) =
∂t
∂x
∂y
(1.14)
∂H
∂H
− v(−H)
.
w(−H) = −u(−H)
∂x
∂y
where z = −H(x, y) denotes the bottom, and z = η(x, y, t) the free surface.
These important relations have to be satisfied for any flow, and for any type of
fluid (e.g., viscous or inviscid).
The expressions in (1.14) are useful for obtaining an evolution equation for
the free surface η, starting from the continuity equation (1.2) in component form:
∂u ∂v ∂w
+
+
=0 .
∂x ∂y
∂z
(1.15)
To this end, (1.15) is integrated term-by-term from the bottom z = −H(x, y)
to the surface z = η(x, y, t). Then, the horizontal derivatives are moved out of
the integrals using Leibniz’ rule, and the terms involving the vertical velocity are
6
replaced with the help of (1.14). The yields
Z η
Z η
∂
∂
∂η
u(z) dz −
v(z) dz
= −
∂t
∂x −H
∂y −H
(1.16)
∂V
∂U
−
,
= −
∂x
∂y
where the integrated velocities U and V are usually referred to as the horizontal
transports. This relation, which will prove to be useful in later chapters, therefore
states that any horizontal divergence in the transports results in a local rise or
fall of the water level.
For the special case of a viscous fluid, so-called dynamic boundary conditions
are derived from the assumption that the velocity of a fluid particle directly
adjacent to a rigid wall moves with the speed of the wall (e.g., it is zero if the
wall doesn’t move). For the bottom boundary, this yields the simple relation
u = 0 at z = −H(x, y) .
(1.17)
Equivalent to these Dirichlet-type boundary conditions, von Neumann boundary
conditions involving the momentum flux, instead of momentum itself, can be
formulated for the lower boundary:
ν
∂u
τb
=
∂n
ρ0
at z = −H(x, y) ,
(1.18)
where τ b is the momentum flux into the fluid at the bottom, and n denotes the
local coordinate normal to the bottom (directed away from the fluid into the
bottom). Note that the bottom momentum flux is equal to the stress exerted by
the bottom on the fluid. Similarly, the boundary conditions at the surface is
∂u
τs
ν
=
∂n
ρ0
at z = η(x, y, t) ,
(1.19)
where n is again the outward unit normal vector, and τ s denotes the momentum
flux into the fluid at the surface, or the stress exerted by the atmosphere on the
fluid.
1.2
Tracer equations
From the appearance of ρ in the Boussinesq equations (1.1), it is evident that
both θ and S may be essential for the dynamics of the system because these
quantities have a strong effect on the density of the fluid. It is known from
the thermodynamics of viscous fluids that the dependency of density on these
7
quantities is given by an empirical equation of state that is of the following general
form:
ρ = ρ(θ, S, p) .
(1.20)
For water at surface pressure, this relation is displayed in Figure 1.1, illustrating
that density generally increases with salinity but decreases with temperature, except for temperatures below the temperature of maximum density (approximately
4 °C for fresh water).
1025
1000.2
1020
25
−3
−3
ρ (kg m )
20
1015
ρ (kg m )
1000
15
1010
10
999.8
999.6
1005
5
1000
995
999.4
0
0
5
10
T (°C)
999.2
15
0
5
10
T (°C)
15
Figure 1.1: Left: Density as function of temperature for salinities ranging from 0
to 25 g/kg. Right: same as left, but enlarged for zero salinity (freshwater). The
square dots denote the freezing point, while the bullets denote the maximum
density temperature.
In many idealised process-oriented applications, the full non-linearity of the
equation of state is not needed. In this cases, a linearised form is frequently used:
ρ
= 1 − α(θ − θ0 ) + β(S − S0 )
(1.21)
ρ0
with constant reference values ρ0 , θ0 and S0 , thermal expansion coefficient α, and
haline contraction coefficient β. Typical values for temperate ocean water are α =
2 · 10−4 K−1 and β = 8 · 10−4 kg g−1 . Note that in (1.21), the pressure dependency
is entirely ignored. This is justified for the idealized problems discussed in the
following, as well in situations with small water depth and thus small pressure
variations. In many other cases, however, the compressibility of sea water does
matter (e.g. by adiabatically changing the temperature), and more sophisticated
models for the equation of state are required.
8
It is often convenient to use the buoyancy
b = −g
ρ − ρ0
ρ0
(1.22)
instead of the density, which simplifies some expressions due to the scaling by
g/ρ0 .
The budget equation for temperature is derived from the first law of thermodynamics, assuming, consistent with the Boussinesq assummption, that compressibility effects can be ignored (see lecture on “Hydrodynamics”). Furthermore, the
same Reynolds averaging procedure as for the Navier-Stokes equations is carried
out. Additionally, we allow for short-wave solar radiation as the only external
energy source. With this, the temperature equation is of the following form:
∂θ
∂θ
∂
∂θ
∂θ
∂θ
∂θ
∂
Kh
−
Kh
+u
+v
+w
−
∂t
∂x
∂y
∂z ∂x
∂x
∂y
∂y
(1.23)
∂
∂θ
1
∂I
−
(νt′ + ν θ )
=
,
∂z
∂z
ρ0 Cp ∂z
with the eddy diffusivity νt′ , the molecular diffusivity ν θ , the horizontal eddy
diffusivity Kh , the heat capacity at constant pressure Cp (≈ 3850 J kg−1 ) and the
solar radiation in the water I (in W m−2 ).
The salinity budget equation can be derived from simple mass conservation
and subsequent Reynolds decomposition:
∂S
∂
∂S
∂S
∂S
∂S
∂
∂S
Kh
−
Kh
+u
+v
+w
−
∂t
∂x
∂y
∂z
∂x
∂x
∂y
∂y
(1.24)
∂
∂S
−
(νt′ + ν S )
=0,
∂z
∂z
where ν S denotes the molecular diffusivity of salt.
Bottom boundary conditions for temperature and salinity are typically zeroflux conditions, whereas at the surface freshwater and heat fluxes are needed. The
freshwater flux is given by the net effect of evaporation and precipitation (this
also needs to be considered for the surface elevation equation (1.16)), and the heat
fluxes typically consist of sensitive (temperature difference between sea surface
and atmosphere), latent (heat loss due to evaporation) and long-wave (emission
of long wave radiation) contributions. As we will see below, dynamically relevant
is only the surface buoyancy flux,
B0 = −νt′
∂b
, for z = η ,
∂z
resulting from the effective heat and salinity fluxes at the free surface.
9
(1.25)
In many coastal problems, so-called suspended particulate matter (SPM),
i.e., particles with a certain settling velocity ws > 0 is of great importance,
mainly because it (i) influences the light climate of coastal waters, (ii) consists of
organic, degradable material and (iii) may on long time scales change the bottom
bathymetry of coastal waters. SPM is typically quantified as a concentration C
and may therefore be calculated by means of a conservative tracer equation with
an additional settling velocity:
∂C
∂C
∂
∂C
∂C
∂C
∂
∂C
Kh
−
Kh
+u
+v
+ (w − ws )
−
∂t
∂x
∂y
∂z
∂x
∂x
∂y
∂y
(1.26)
∂
∂C
−
νt′
=0,
∂z
∂z
In simple cases, it is assumed that there are no net fluxes of SPM through the
surface and the bottom.
10
Chapter 2
Frictional effects and Ekman
theory
2.1
Pressure-driven, non-rotating flows
Let’s consider an unstratified (ρ = ρ0 ), stationary, horizontally homogeneous flow
over flat bottom. If we assume in addition that the flow is non-rotational (f = 0)
with a constant (prescribed) surface elevation gradient, the dynamic equations
(1.13) reduce to a simple balance between the pressure gradient and the vertical
divergence of the turbulent momentum flux (the molecular flux can usually be
ignored):
∂
∂u
∂η
νt
=g
,
(2.1)
∂z
∂z
∂x
where we have assumed, without loss of generality, that the pressure gradient is
aligned with the x-axis. The momentum balance in (2.1) can be integrated over
the whole water column to yield
∂η
τxs τxb
+
= gH
,
ρ0 ρ0
∂x
(2.2)
where τxs and τxb denote the surface and bottom stresses, respectively, that have
already been defined in (1.18) and (1.19). These stresses can also be interpreted
as the surface and bottom momemtum fluxes using the convention that fluxes
entering the fluid are counted positive. We ignore the surface stress in this simple
example, and note that the bottom stress is often expressed
in terms of the
b
b
b b
friction velocity u∗ implicitly defined as τx /ρ0 = −u∗ u∗ . Assuming here and in
the following, for simplicity, that τxb ≤ 0 for ∂η/∂x ≤ 0 (bottom stress opposes
the mean flow), (2.2) becomes
− (ub∗ )2 = gH
11
∂η
,
∂x
(2.3)
indicating that the imposed pressure gradient due to the surface tilt is compensated by the bottom stress.
Integration of (2.1) from the bottom to an arbitrary position in the water
column, again assuming zero surface stress and inserting (2.3), results in
νt
z
∂u
= −(ub∗ )2 ,
∂z
H
(2.4)
which, introducing the local momentum flux,
τx = −ρ0 νt
∂u
,
∂z
(2.5)
can also be written as
z b
τ .
(2.6)
H x
This expression indicates that for a purely pressure-driven flow the momentum
flux is a linear function of depth. Division of (2.4) by νt and vertical integration
leads, after including the no-slip bottom boundary condition u(−H) = 0, to an
expression of the form
Z
(ub∗ )2 z z ′ ′
dz .
(2.7)
u(z) = −
H
−H νt
τx = −
The solution of (2.7) for constant νt is the following parabolic velocity profile (a
special form of the Hagen-Poiseuille flow):
z 2 (ub∗)2 H
1−
.
(2.8)
u(z) =
2νt
H
Alternatively, introducing the vertically averaged flow velocity,
Z
1 0
u(z) dz ,
ū =
H −H
the parabolic profile in (2.8) may also be written as
z 2 3
u(z) = ū 1 −
.
2
H
(2.9)
(2.10)
This solution holds with perfect accuracy for laminar flow (i.e., replacing νt by
the molecular viscosity ν). The structure of the velocity profile for a typical
value of the turbulent diffusivity is shown in figure 2.1. However, for realistic
turbulent flows, the assumption of a constant eddy viscosity is only a rather poor
approximation as discussed in the following.
As most turbulent wall-bounded shear flows, marine bottom boundary layers
are generally characterized by a thin near-bottom region with particularly strong
12
vertical gradients. In this region, the assumption of a constant eddy diffusivity
fails, and there are two frequently used approaches to overcome this problem.
In the first approach, it is attempted to parameterize the thin, unresolved nearbottom layer with the help of an expression relating the flow “just above the
bottom “ to the bottom stress. In the most simple case, the following linear
expression is used:
(ub∗ )2 = F u(−H)
(2.11)
where F denotes the linear drag coefficient. Note that νt is still assumed to
be constant but the no-slip condition (u = 0) at the bottom is longer imposed.
Under these conditions, the integration of (2.4) yields
z 2 FH
u(z) = u(−H)
1−
+1 ,
(2.12)
2νt
H
which, for given vertical mean velocity ū, can be re-expressed as
z 2 FH
1−
+1
2νt
H
u(z) = ū
.
2 FH
+1
3 2νt
(2.13)
The second approach relies on the more realistic assumption of a parabolic
profile for the turbulent diffusivity:
z
b
,
(2.14)
νt = κ|u∗ |(z + z0 + H) −
H
where κ = 0.4 is a universal model constant (the so-called von Kármán constant),
and z0 the bottom roughness, both found from flow experiments. The physical
motivation of (2.14) can be understood from considering the linear limiting case
very close to the bottom:
νt = κ|ub∗ |(z + z0 + H) .
(2.15)
This expression states that the viscosity very close to the lower boundary is the
product of a turbulent velocity scale, set proportional to the only available turbulent velocity scale ub∗ in this problem, and a turbulent length scale, set proportional to the distance from the wall, which is the only relevant length scale here.
The diffusivity of the form (2.15) then follows from dimensional arguments. It is
worth noting that this is completely analogous to the kinematic theory of gases,
predicting that the molecular viscosity is the product of a velocity scale (speed
of sound) and a length scale (mean free molecular path). The parabolic profile in
(2.14) provides nothing more than a simple extrapolation of two boundary layers
with linear diffusivities as in (2.15) at the top and bottom, respectively, for the
whole water column. More details on this topic are discussed in the lecture on
Marine Turbulence.
13
Inserting (2.14) into (2.7), and integrating leads to a logarithmic variation of
the velocity:
z + H + z0
ub∗
.
(2.16)
ln
u(z) =
κ
z0
Alternatively, if the vertically averaged current velocity ū is given instead of the
bottom friction velocity ub∗ , the logarithmic law in (2.16) can be expressed as
ū
z + H + z0
u(z) =
ln
.
(2.17)
z0
1 + zH0 ln zH0 + 1 − 1
Figure 2.1 shows velocity profiles for the constant eddy viscosity, following (2.10),
and two profiles with the more realistic parabolic eddy viscosity, following (2.17),
all profiles being normalised with the depth-mean velocity. Clearly, the resulting
shape for the constant eddy viscosity differs strongly from the logarithmic profiles
resulting from the parabolic eddy viscosity, particularly in view of the much
stronger near-bottom shear of the former.
0
z/H
constant νt
parabolic νt , z0 = 10−4 m
parabolic νt , z0 = 10−2 m
-0.2
νt = 0.01 m2 s−1 , F = 0.001 m s−1
-0.4
-0.6
-0.8
-1
0
0.2
0.4
0.6
0.8
u / ū
1
1.2
1.4
1.6
Figure 2.1: Stationary velocity profiles with constant and parabolic eddy viscosity
profiles.
It is worth noting that very close to the bottom, the linear stress profile can
be approximated by a constant-stress layer, implying that the local momentum
flux is equal to the bottom stress:
νt
∂u
= (ub∗ )2 = const.
∂z
14
(2.18)
Inserting (2.15), valid only in the near-bottom region, into (2.18), we find again
the logarithmic profile in (2.16). This profile, often referred to as the logarithmic
wall layer, is one of the most well-known results from turbulence theory for wallbounded shear flows. It can alternatively be derived from purely dimensional
arguments (not discussed here), and therefore corresponds to an exact solution
of the Navier-Stokes equations. If extended to the whole water column, however,
additional modeling assumptions are involved that may not always hold.
Finally, it should be noted that the logarithmic law (2.16) can also be rewritten
as
(ub∗ )2 = CD u2 ,
(2.19)
with the quadratic drag coefficient defined as


CD = 

ln
2

κ

z + H + z0 
z0
.
(2.20)
For the range of applicability of this relation (i.e. very close to the wall), typical
values are of the order of CD ≈ 10−3 . Such a quadratic relation between the
stress and the speed is characteristic for many other turbulent flows (e.g., for
flows in pipes), whereas laminar flows typpically exhibit a linear relationship as
in (2.8).
2.2
Ekman theory
A more general view of the water column structure is obtained by including the
effects of Coriolis accelerations but otherwise retaining the assumption from the
previous section. According to (1.13), the horizontal momentum budget is then
of the following form:
∂
∂u
∂η
∂v
∂η
∂
νt
+ fv = g
νt
− fu = g
,
,
(2.21)
∂z
∂z
∂x
∂z
∂z
∂y
which are often referred to as the Ekman equations. If for a moment frictional
effects are neglected in (2.21), we obtain a balance between the pressure gradient and the Coriolis force, the so-called geostrophic balance with the resulting
geostrophic velocities defined as
ug = −
g ∂η
,
f ∂y
vg =
g ∂η
f ∂x
.
(2.22)
Note that the geostrophic flow occurs at a right angle with respect to the pressure
gradient, i.e. currents follow lines of constant pressure (isobars). Inserting the
15
definitions (2.22) into (2.21), we obtain the relation
∂u
∂v
∂
∂
νt
+ f (v − vg ) = 0 ,
νt
− f (u − ug ) = 0 ,
∂z
∂z
∂z
∂z
(2.23)
illustrating that frictional effects lead to an ageostrophic component of the flow.
In one of the classical examples, originally devised by the Swedish oceanographer Vagn Walfrid Ekman, the flow near the surface of an infinitely deep ocean
with negligible surface pressure gradient (ug = vg = 0) is considered. Assuming,
without loss of generality, that the wind stress is aligned with the x-axis, the
horizontal transports,
Z 0
Z 0
U=
u dz , V =
v dz ,
(2.24)
−∞
−∞
for this flow are obtained by vertically integrating (2.23) from z = −∞ to z = 0.
This yields:
τs
U =0,
V =− x ,
(2.25)
f ρ0
with the striking interpretation that the so-called Ekman transport is exactly to
the right of the imposed wind stress.
In order to obtain managable analytical expressions for the vertical variability
of the Ekman transport, we assume in the following a constant viscosity, recalling
from the above discussion that this represents the laminar case but should be
treated with caution if applied to turbulent flows. For this case, the Ekman
equations in (2.23) become
νt
∂2u
+ f (v − vg ) = 0 ,
∂z 2
νt
∂2v
− f (u − ug ) = 0 ,
∂z 2
(2.26)
with boundary conditions (1.17) and (1.19). Introducing the complex velocity
W = (u − ug ) + i(v − vg ) ,
(2.27)
we are able to transform (2.26) into the more compact relation
νt
∂2W
− if W = 0 .
∂z 2
(2.28)
The no-slip condition (u = v = 0) at the bottom leads to
W = −ug − ivg
for z = −H ,
(2.29)
and, following (1.19), the surface stress is included in the form
νt
τys
τs
∂W
= x +i
∂z
ρ0
ρ0
16
for z = 0 .
(2.30)
Using an exponential ansatz of the form
W (z) ∝ eiλz ,
(2.31)
the linear eigenvalue problem in (2.28) leads to
(νt λ2 + if )eiλz = 0 .
(2.32)
For arbitrary z, this relation can only be satisfied if λ2 = −if /νt , or, taking the
root, for
1−i
(2.33)
λ=±
∆
where we have introduced the Ekman depth
r
2νt
∆=
.
(2.34)
f
The general solution is therefore of the form
W (z) = Ae−(i+1)ζ + Be(i+1)ζ ,
(2.35)
where A and B are two complex constants determined by the boundary conditions, and ζ = z/∆ denotes the non-dimensional coordinate.
The expressions for A and B are obtained by inserting (2.35) into the bottom
and surface boundary conditions (2.29) and (2.30), which yields
Ae(i+1)δ + Be−(i+1)δ = −ug − ivg ,
(2.36)
where δ = H/∆ is the ratio of the water depth to the Ekman depth, and
B−A=
τxs + iτys
ρ0 νt (i + 1)/∆
.
(2.37)
The relations in (2.37) and (2.36) form two equations for the determination of
the two unknowns A and B.
The re-transformation of the complex solution to real velocity components
u and v is straightforward but quite tedious, and here we only discuss the final
result as presented by Simons (1980):
α
β
γ+ε s γ−ε s
1
u = 1−
ug − vg +
τ +
τ ,
c
c
∆ρ0 f
c x
c y
(2.38)
α
β
1
γ+ε s γ−ε s
v = 1−
vg + ug +
τ −
τ ,
c
c
∆ρ0 f
c y
c x
17
with c = cosh(2δ) + cos(2δ), and a number of parameters depending on the
non-dimensional vertical coordinate:
α(ζ) = cosh(δ + ζ) cos(δ − ζ) + cosh(δ − ζ) cos(δ + ζ) ,
β(ζ) = sinh(δ + ζ) sin(δ − ζ) + sinh(δ − ζ) sin(δ + ζ) ,
(2.39)
γ(ζ) = sinh(2δ + ζ) cos(ζ) + sinh(ζ) cos(2δ + ζ) ,
ε(ζ) = cosh(2δ + ζ) sin(ζ) + cosh(ζ) sin(2δ + ζ) .
For large depths compared to the Ekman layer thickness (δ → ∞), the surface
and bottom boundary layers are dynamically decoupled, and greatly simplified
solutions can be obtained. For the example of an infinitely deep ocean with
vanishing surface pressure gradient (ug = vg = 0), and the surface stress directed
into the x-direction (τys = 0), the solution in (2.38) becomes:
τxs
u = eζ cos(ζ) + sin(ζ)
,
∆f ρ0
v = eζ [sin(ζ) − cos(ζ)]
τxs
∆f ρ0
(2.40)
.
The following observation can be made from (2.40). The velocity vector the
decreases exponentially away from the surface with a superimposed clockwise
rotation. This type of motion has become famous in physical oceanography under
the name Ekman spiral. From the first term on the right hand side, the Ekman
depth ∆ is identified as the exponential decay rate with depth. Quite peculiar is
also the fact that, directly at the surface, the velocity exhibits an angle π/4 with
respect to the direction of the surface stress, i.e. it is not aligned with the wind.
Finally, it is worth noting that from the vertical integral of (2.40), the Ekman
transport in (2.25) can be recovered. An example of the surface Ekman spiral is
provided in Figure 2.2.
To derive a solution for the bottom Ekman layer for deep water, we assume
that the surface momentum flux is zero (τxs = τys = 0) and that the geostrophic
velocity has a vanishing y-component (so the surface tilt is purely in the ydirection). Using similar simplifications as for the near-surface solution, (2.38)
for this case yields
u = 1 − e−(δ+ζ) cos(δ + ζ) ug ,
v = e−(δ+ζ) sin(δ + ζ)ug
18
.
(2.41)
To derive (2.41) we have used that near the bottom of deep water δ + ζ ≪ δ − ζ.
For the interpretation of (2.41), it is helpful to note that the shifted vertical
coordinate ξ = δ + ζ appearing in the function arguments starts at the bottom,
pointing upwards.
Figure 2.2 reveals that this solution corresponds to an Ekman spiral as well,
now however with the velocity vector veering in the anti-clockwise direction towards the bottom with exponentially decreasing speed. If the distance from the
bottom is much larger than the Ekman layer thickness (ξ = δ + ζ ≫ 1) the solution converges towards the interior inviscid, geostrophically balanced velocity
ug . The Ekman transport found from the integral of (2.41) exhibits an angle π/2
to the left of the interior flow ug . Directly at the bottom, the velocity is zero
(no-slip) and the bottom stress according to (1.18) is
∂u
νt
τxb
= −νt
= − ug ,
ρ0
∂z
∆
τyb
∂u
νt
= −νt
= − ug ,
ρ0
∂z
∆
(2.42)
so the bottom stress is not opposite to the interior flow (as in the non-rotating
case discussed in the previous section) but rotated at an angle of π/4.
19
Surface Ekman spiral
Surface Ekman spiral
0
1
-1
v / us
z/∆
0.5
-2
-3
0
-0.5
-4
u
v
-5
-1
-0.5
0
-1
0.5
1
(u,v)
-1
(u, v) / us
Bottom Ekman spiral
-0.5
0
0.5
1
u / us
Bottom Ekman spiral
5
0.5
v / ug
(z + H) / ∆
1
4
3
2
0
-0.5
1
0
-0.2
u
v
0
0.2 0.4 0.6 0.8
-1
1
1.2
(u,v)
-1
(u, v) / ug
-0.5
0
0.5
1
u / ug
Figure 2.2: Surface (upper panels) and bottom (lower panels) Ekman solutions
shown as profiles (left panels) and hodographs (right panels). The green dots in
the hodographs are solutions of the Ekman dynamics for each step of ∆/2.
20
Chapter 3
Vertically integrated
shallow-water equations
In many applications the vertically integrated currents and the resulting surface
elevation are of major interest, for example for storm surge and tsunami predictions, and for many other wave phenomena with wave length much larger than
the water depth. Therefore, it is often attractive to first vertically integrate the
three-dimensional shallow water equations from the bottom at z = −H(x, y) to
the surface at η(x, y, t). Carrying out this vertical integration and inserting the
dynamic boundary conditions (1.18) and (1.19) and the kinematic boundary conditions (1.14), the u-component of the shallow water equations (1.13) is of the
following form:
Z
Z η
Z η
Z η
∂
∂
∂ η
u dz +
uu dz +
uv dz − f
v dz =
∂t −H
∂x −H
∂y −H
−H
(3.1)
Z η Z η
∂η
∂ρ
g
D ∂pa τxs τxb
−gD
−
dz̄dz −
+
+ .
∂x ρ0 −H z ∂x
ρ0 ∂x
ρ0 ρ0
Inserting the vertically integrated transports U and V defined in (1.16) and using
the approximations
Z η
Z η
UU
UV
,
,
(3.2)
uu dz ≈
uv dz ≈
D
D
−H
−H
(3.1) is of the following form:
∂ UV
∂ UU
∂U
+
− fV =
+
∂t
∂x D
∂y
D
∂η
g
−gD
−
∂x ρ0
Z
η
−H
Z
z
η
D ∂pa τxs τxb
∂ρ
+ .
dz̄dz −
+
∂x
ρ0 ∂x
ρ0 ρ0
21
(3.3)
The v-component of (1.13) may be transformed and approximated accordingly:
∂ VV
∂ UV
∂V
+
+ fU =
+
∂t
∂x D
∂y
D
(3.4)
Z Z η
∂ρ
g η
D ∂pa τys τyb
∂η
−
dz̄dz −
+
+ .
−gD
∂y ρ0 −H z ∂y
ρ0 ∂y
ρ0 ρ0
22
Chapter 4
Bottom gravity currents
4.1
Introduction
An interesting class of geophysical flows that we have ignored so far is driven by
density differences inside the fluid that, in a gravitational field, lead to additional
body force terms in the momentum budget (1.8). If the density is smaller than
that of the ambient fluid, the gravity current is buoyant, and thus moves near the
surface. So-called river plumes, in which fresh (and thus light) river water moves
on top of saline ocean water form a well-known example that can be observed at
many places in the coastal ocean.
Contrary, if the density is larger than that of the ambient, the fluid sinks
and propagates as a bottom gravity current or dense bottom plume along the
topography. Such flows may be driven by higher salinity, lower temperature, or
higher sediment concentration, which all result in a positive density contrast with
respect to the ambient fluid. Bottom gravity currents in the ocean occur on a
wide range of scales, from small-scale flows of a few meters width that propagate
in small topographic corrugations, up to large-scale “overflows” that may have a
width of up to 100 km, and a thickness of the order of 100 m. The latter form
an important link in the global ocean circulation, and no useful ocean climate
predictions can be made without their proper description. The same is true
for the Baltic Sea, in which dense bottom gravity currents are formed by saline
(dense) waters from the North Sea that, driven by their negative buoyancy, slowly
travel along the bottom towards the deeper basins of the central Baltic Sea.
Here, we aim at a simple dynamical description of bottom gravity currents
in different configurations of increasing complexity. Our major goal will be the
understanding of the basic physical principles that govern the flow, which requires
the introduction of a number of approximations that simplify the mathematical
description but may not always be justified. We focus on frictionally controlled
bottom gravity currents, which contrasts the frictionless (hydraulic) treatment
that can be found in many classical textbooks on buoyancy-driven flows (e.g.
23
Baines, 1995).
4.2
Governing equations
In our simple model, we assume that a bottom gravity current consists of a
homogeneous water mass with constant density ρ that is larger than the constant
density ρ0 of the ambient fluid. We will see below that the dynamically relevant
quantity related to this density difference is the relative buoyancy
ρ − ρ0
B = −g
,
(4.1)
ρ0
which is constant and negative, consistent with the idea the dense bottom gravity
currents are negatively buoyant. The gravity current is bounded from below at
z = −H by bottom topography, and from above at z = η by a moving interface
separating the gravity current from the ambient fluid (Figure 4.1). This interface
is assumed to be material, i.e. no ambient fluid is allowed to enter the gravity
current. With this assumption, the possibility of mixing with ambient fluid is
ignored, which is consistent with our assumption of constant B, and greatly
simplifies the analysis. Clearly, however, the assumption of negligible mixing is
rather serious, and may not always be justified.
Figure 4.1: Geometry and coordinate system for the description of dense bottom
currents on a slope with variable topography.
The velocity components in the x and y directions are denoted by u and
v, respectively, which may vary horizontally but not in the vertical. The latter
implies that the transports U and V introduced in (1.16) can be expressed as U =
uD and V = vD, where D denotes the local thickness of the gravity current such
that η = −H + D. The key for the following investigations are the assumptions
(a) that the flow is shallow (i.e. lateral scales are much larger than changes in
D and η), and (b) that the ambient fluid above the gravity current is infinitely
deep and stagnant (u = v = 0) at all times. This geometric simplification is often
referred to as the 1.5-layer approximation.
24
Then, the gravity current will be described by exactly the same vertically
integrated equations that also describe an unstratified flow with a free-surface
(see Chapter 3), except for the fact that the acceleration of gravity, g, appearing
in (3.3) and (3.4) has to be replaced by the buoyancy, −B. This last step can
be understood from the following argument. Let’s assume that the pressure at
a certain reference level zr above the gravity current corresponds to an arbitrary
but constant reference pressure pr . Then, the hydrostatic pressure at a point
inside the gravity current is given by
p = pr + (zr − η)ρ0 g + (η − z)ρg
.
(4.2)
The pressure gradient appearing in (1.10) is therefore
−
1 ∂p
∂η
∂η ρ
ρ − ρ0 ∂η
∂η
=
g−
g = −g
=B
,
ρ0 ∂x
∂x
∂x ρ0
ρ0 ∂x
∂x
(4.3)
where in the last step the definition of the buoyancy in (4.1) has been used. A
similar equation follows for the y-direction. Vertical integration of the pressure
gradient in (4.3) than results in BD∂η/∂x and BD∂η/∂y for the two components
of the interfacial pressure gradient, which proves the statement made above.
The geometry of the problem can be further simplified by assuming that the
flow is homogeneous in the x-direction. This implies that all x-derivatives can be
ignored, and that the flow can be described in a two-dimensional plane spanned
by the y- and z-axes (see Figure 4.1). Under these conditions, the integrated
continuity equation in (1.16) becomes
∂vD
∂η
=−
,
∂t
∂y
(4.4)
expressing the simple fact that any lateral flow convergence results in a local
change in interface height if volume is to be conserved.
The integrated momentum budget in (3.3) and (3.4) can be re-written as
∂uD ∂uvD
τb
+
− f vD = x ,
∂t
∂y
ρ0
∂η τyb
∂vD ∂v 2 D
+
+ f uD = BD
+
,
∂t
∂y
∂y ρ0
(4.5)
where we have replaced the factor g by −B as discussed above, and used the facts
that (i) the density is constant inside the gravity current (see above), (ii) the
atmospheric pressure pa is irrelevant, and (iii) there is no stress at the interface.
The latter is consistent with the assumption of negligible mixing (why?). The
momentum budget in (4.5) expresses a balance between acceleration, Coriolis
force, pressure gradient due to the interface tilt, and bottom friction.
25
Equations (4.4) and (4.5) are only closed if the bottom stresses can be expressed in terms of other known flow properties. Such models generally relate
the stresses to the velocities, resulting in the most simple case in a relation of the
form
τxb
= −Ku ,
ρ0
(4.6)
τyb
= −Kv ,
ρ0
where K denotes the linear drag coefficient. The advantage of this model is its
linearity that makes it often possible to find closed analytical solutions of the
problem. The disadvantage is that a linear relationship between the velocity and
the stress is generally not observed in turbulent flows. As discussed in the context
of (2.19), a more realistic but also more complex alternative is a quadratic drag
law of the form
τxb
= −CD |u| u ,
ρ0
(4.7)
τyb
= −CD |u| v ,
ρ0
where CD denotes the (dimensionless) quadratic drag coefficient, and |u| = (u2 +
v 2 )1/2 the total speed.
4.3
Non-rotating flows
We first consider the non-rotating case with f = 0. Under these conditions the
flow is fully two-dimensional in the y-z plane (u = 0), and the momentum budget
in (4.5) reduces to
∂vD ∂v 2 D
∂η
+
= BD
− CD |v| v ,
∂t
∂y
∂y
(4.8)
where we have used (4.7) for the bottom stress.
4.3.1
Similarity solutions
Due to the strong non-linearities of the coupled system (4.4) and (4.8), a generally
valid analytical solution of the problem cannot be found. For some special cases,
however, similarity solutions exist, revealing a number of interesting features of
the system. One solution of this type can be constructed by assuming that both
the bottom slope, Sy = ∂H/∂y, and the down-slope speed, v = v0 , are constant.
Assuming v0 > 0 (downslope flow), it is then possible to introduce a similarity
variable,
ζ = y − v0 t ,
(4.9)
26
hoping that self-similar solutions of the form η(ζ), D(ζ), and v(ζ) exist. Using (4.9), it is easy to see that the partial derivatives transform into ordinary
derivatives according to
d
∂
d
∂
= −v0
,
=
,
(4.10)
∂t
dζ
∂y
dζ
which strongly simplifies the problem in two ways. First, using (4.10) and the fact
that v = v0 is constant, it is easy to show that (4.4) is automatically satisfied.
Second, with exactly the same argument, it can be shown that the rate and
advection terms on the left hand side of (4.8) cancel (although they do not vanish
individually). Then, (4.8) can be written in self-similar form as
CD v02
,
(4.11)
B
where we have introduced the constant parameter R, and used the convention
that derivatives with respect to ζ are denoted by a prime. Note that if we had
used the linear drag parameterization in (4.6), we would have arrived at exactly
the same expression, however, with R = Kv0 /B. Recall that B and thus R are
negative.
In spite of its simple appearance, (4.11) is a non-linear differential equation for
which solutions are not obvious. Combining tabulated solutions for differential
equations, or using a symbolic algebra package, it can, however, be shown that
the full solution is given by the expression
"
#!
2ζ
Sy
− 1+ R
R
D(ζ) = −
,
(4.12)
1 + W −e
Sy
D(D ′ − Sy ) = R ,
R=
which for v0 > 0 is valid for −∞ < ζ ≤ 0. Note that the symbol W in (4.12)
denotes the so-called Lambert-W function; here, we don’t bother too much about
the properties of this function except for noting that W (−e−1 ) = −1 and W (0) =
0. The former implies that at ζ = 0 the plume thickness is zero; the latter
indicates that for ζ → −∞ (e.g. for large times) we have D → −R/Sy, indicating
that the thickness of the gravity current approaches a maximum value.
For some parameters that are typical for the coastal ocean, the evolution of
a gravity current on a plane slope according to (4.12) is shown in Figure 4.2.
The gravity current is seen to follow a steep front (it can be shown that the
interface slope in fact becomes infinite at the front), moving downward with constant velocity v0 as required by (4.9). The thickness increases continuously, and
slowly approaches the limiting value −R/Sy = 5 m that is reached for t → ∞.
It is worth noting that the shallow-water assumption inherent in the vertically
integrated equations (4.5) cannot hold in the vicinity of the front, where vertical and horizontal variations become comparable. The real shape of the nose of
the gravity current is therefore expected to be somewhat different from the solution derived here, which is confirmed by available laboratory experiments and
numerical simulations of the full (non-hydrostatic) equations.
27
z (m)
5
t=1h
0
−5
z (m)
5
t=2h
0
−5
z (m)
5
t=3h
0
−5
z (m)
5
t=4h
0
−5
0
2
4
6
8
10
y (km)
Figure 4.2: Solution (4.12) for different times as indicated. The parameters used
are Sy = 5 × 10−4 , v0 = 0.5 m s−1 , CD = 10−3 , and B = −0.1 m s−2 .
4.3.2
Formulation as a diffusion problem
While the solution in (4.12) reveals a number of important properties, it is valid
only under quite restrictive conditions (e.g., only for constant velocity). To derive
a more generally valid framework, we first consider (4.8) under the additional
assumption that the rate and advection terms on the left hand side cancel, and
that friction is linear as in (4.6). From this, we obtain an expression of the form
BD
∂η
= Kv ,
∂y
(4.13)
which is a simple balance between the down-slope pressure gradient and bottom
friction. Note that (4.13) includes the similarity solution discussed above, for
which the rate and advection terms were shown to cancel as well. Re-arranging
(4.13) as an expression for the down-slope transport,
vD =
B 2 ∂η
D
,
K
∂y
(4.14)
and inserting into (4.4) yields
∂η
∂
=
∂t
∂y
η ∂η
ν
,
∂y
where
νη = −
B 2
D
K
28
.
(4.15)
(4.16)
This suggests that the motion of the gravity current is governed by a diffusion
equation with a non-constant diffusivity ν η . According to (4.16), we find ν η → 0
as D → 0 near the boundaries of the gravity current, such that ν η remains
strictly positive and bounded everywhere, as required for any well-posed diffusion
problem. It is interesting to note that the non-linearity of this diffusion problem
may lead to the development of steep fronts, which is somewhat counter-intuitive
in view of the usual expectation that diffusion has a tendency to smooth gradients.
It is easy to show (see assignments) that for the more complex case of quadratic
friction in (4.7), the expression for the interface diffusivity becomes
η
ν =
−B
CD
21 − 21
∂η D 23
∂y .
(4.17)
This quadratic solution retains the crucial properties of positivity and boundedness even near the boundaries, where ν η → 0.
4.4
The effect of rotation
A measure for the relative importance of rotational effects (that we have ignored
so far) and frictional effects is the so-called Ekman number. In the context of
(4.5), this quantity is defined as
Ek =
|τ b |
,
ρ0 f D |u|
(4.18)
which, using the quadratic drag law (4.7), can be re-written as
Ek =
CD |u|
fD
.
(4.19)
Using typical mid-latitude values (CD = 10−3 , f = 10−4 s−1 , u = 0.1 − 1 m
s−1 ), it is clear that rotational effects will generally be comparable to, or even
dominating over, frictional effects, unless the gravity current is very thin in which
case one finds Ek ≫ 1 (rotational effects can be ignored).
Searching for a simple framework for the analysis of the effects of rotation, we
start from the general momentum budget in (4.5). As above, we assume that the
rate and advection terms cancel, which removes so-called internal interface waves
from the solution. Since, here, we are not concerned about wave phenomena, this
assumption is convenient and greatly simplifies the analysis. However, it should
be kept in mind that it may not always be justified. Note that assuming that
the rate and advection terms cancel still allows for the evolution of fronts at the
plume boundaries, as we have already seen above.
29
With these assumptions, (4.5) can be written as
f vD = Ku ,
f uD = BD
∂η
− Kv ,
∂y
(4.20)
where we have again used the linear drag law in (4.6) for simplicity.
For K = 0 (no friction), (4.20) is recognized as a geostrophic balance between
the cross-slope pressure gradient due to the interface tilt, and the along-slope
flow, similar to the flow discussed in the context of (2.22). In this case, the
interface of the gravity currents remains stationary for all times, which contrasts
the non-rotating case discussed above. This is one manifestation of the fact
that rotation generally inhibits the conversion of potential to kinetic energy in
geophysical flows. For K 6= 0, however, the relative influence of rotation becomes
smaller, and it is clear from (4.20) that a cross-slope transport will be induced.
This leads to a deformation of the interface according to (4.4), and eventually to
a downslope motion of the gravity current. Friction is therefore essential for the
conversion of potential to kinetic energy associated with the downward migration
of dense bottom currents in rotating reference systems.
Eliminating the along-slope velocity u from (4.20), we obtain
B
∂η
f2 2
,
(4.21)
vD 1 + 2 D = D 2
K
K
∂y
which can be written more compactly as an expression for the down-slope transport
BδL D 2 ∂η
vD =
.
(4.22)
f δL2 + D 2 ∂y
where we have introduced the bulk Ekman layer thickness δL = K/f for linear
bottom friction. Qualitatively, the interpretation of δL as the near-bottom layer
affected by frictional effects is similar to that of the Ekman layer thickness defined
in (2.34) above. Our assumption of vertically constant velocities, however, is
clearly not consistent with the frictionally-induced veering of the velocity vector
(Ekman spiral) discussed in Chapter 2. This is the price that has to be paid for
the advantage of working with the vertically integrated equations. The physical
interpretation δL becomes clear from constructing, analogous to the case with
quadratic friction in (4.19), the Ekman number as the ratio of friction and Coriolis
terms in (4.20). In the linear case, we find Ek = K/(f D), which is equivalent
to Ek = δL /D, indicating that the ratio of the Ekman layer thickness to the
local plume depth describes the relative importance of friction in the momentum
budget.
As for the non-rotating case, it is evident that the transport in (4.22) can
be expressed as the product of the interface gradient and a turbulent diffusivity,
30
defined here as
νη = −
BδL D 2
,
f δL2 + D 2
(4.23)
which, similar to (4.16), is strictly positive, bounded, and converges to zero near
the boundaries of the gravity current. It is easy to show (see assignment) that
for f → 0, the diffusivity in (4.23) becomes identical to the expression for nonrotating flows in (4.16). Finally, inserting (4.23) into the continuity equation
equation (4.4) reveals that also in the rotating case the evolution of the interface
is governed by a diffusion equation of the form (4.15). It can be shown that also
for the case with quadratic friction, a diffusion equation can be derived. The
corresponding expression for the diffusivity is, however, much more complex and
outside the scope of this introductory text.
0
−200
t=0h
z (m)
−400
−600
−800
v (m s−1)
u (m s−1)
−1000
1
0.5
0
0.4
0.2
0
D (m)
−0.2
100
50
0
0
50
100
y (km)
150
200
Figure 4.3: Solution of (4.20) for a gravity current on a plane with constant slope
H ′ = 5 × 10−3 . Parameters used here are B = −0.01 m s−2 , K = 10−3 m s−1 ,
and f = 10−4 s−1 . The initial distribution of the thickness is parabolic (see last
panel), and large compared to the Ekman layer thickness δL (gray-shaded area
in last panel).
For any given interface elevation η and bottom topography H, the plume
thickness D = η + H is known, and (4.22) can be used to compute the downslope
speed v. The along-slope velocity u then follows from the second relation in (4.20).
As an example, we discuss some results for a gravity current with (initially)
parabolic thickness on a sloping plane. The results shown in Figure 4.3 correspond
31
to typical parameters for large scale gravity currents on the continental slope. The
along-slope velocity is seen to reach values above 1 m s−1 near the lower edge
of the gravity currents, where the interface gradient is largest, and, according to
the second relation in (4.20), largest geostrophically balanced speeds have to be
expected. Note the reversal of the along-slope velocity at the upper edge, where
the interface elevation changes sign.
0
−200
t = 72 h
z (m)
−400
−600
−800
v (m s−1)
u (m s−1)
−1000
1
0.5
0
0.4
0.2
0
D (m)
−0.2
100
50
0
0
50
100
y (km)
150
200
Figure 4.4: Same as Figure 4.3 but now for t = 72 h.
As discussed above, the down-slope velocity evident in Figure 4.3 causes a
deformation of the interface that is predicted by the diffusion equation (4.15),
using the diffusivity in (4.23). These equations have been solved numerically,
and the situation after 3 days of integration is displayed in Figure 4.4. The most
striking feature visible in this figure is a thin layer of thickness comparable to
the Ekman layer thickness δL following a sharp front that quickly travels in the
downward direction. This process, sometimes referred to as Ekman drainage,
leads to gradual loss of fluid in the main body of the plume that is much thicker
than the Ekman layer thickness, and only weakly affected directly by friction.
In contrast to the swift descent of the lower edge of the gravity current, the
upper edge remains almost stationary, which leads to an overall widening of the
flow. Near the thickest part of the gravity current (the “vein”), we find δL ≪ D,
indicating that the Ekman number is small, and the momentum budget is close
to a geostrophic balance. The results for the quadratic drag law (not shown) are
similar if the drag coefficient CD is adjusted accordingly.
32
4.5
Gravity current in a rotating channel
An interesting extension of the two-dimensional geometry in (4.5) can be obtained
by allowing for the possibility of an interface tilt in the x-direction. Clearly, to
be consistent with our assumption of homogeneity in this direction, this interface
tilt has to be constant. As an example, imagine an infinitely long channel with
arbitrary cross-section that is tilted into the x-direction such that both bottom
slope and interface tilt are identical: S = −∂η/∂x = ∂H/∂x. With the additional
pressure-gradient term BD∂η/∂x = −BDS introduced by this, (4.5) can be
written as
∂uD ∂uvD
+
− f vD = −BDS − CD |u| u ,
∂t
∂y
(4.24)
∂η
∂vD ∂v 2 D
+
+ f uD = BD
− CD |u| v ,
∂t
∂y
∂y
where we have parameterized the bottom stress with the quadratic drag law
(4.7). The momentum budget in (4.24) is solved simultaneously with our usual
continuity equation (4.4) for two-dimensional flows.
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
z (m)
0
−5
−10
−15
0.0 h
z (m)
0
−5
−10
−15
2.0 h
z (m)
0
−5
−10
−15
4.0 h
z (m)
0
−5
−10
−15
6.0 h
z (m)
0
−5
−10
−15
−10
8.0 h
−5
0
y (km)
5
10
Figure 4.5: Numerical solution of (4.4) and (4.24) for different times as indicated.
Model parameters are: S = 5×10−4 , B = −0.10 m s−2 , CD = 10−3, and f = 10−4
s−1 . Color shading indicates the down-channel speed in m s−1 .
A numerical solution of (4.4) and (4.24) for some parameters that are typical
for bottom gravity currents passing through a system of submarine channels in
the Western Baltic Sea is shown in Figure 4.5. The flow is seen to accelerate
to nearly stationary conditions within a few hours, with higher speeds observed
33
in the deep part of the channel, where the gravity current is thicker and the
relative importance of frictional effects is smaller. The increasing down-channel
speeds are accompanied by a growing cross-channel interface tilt, which provides
a clear indication for the effect of rotation. No interface waves are evident (a
closer examination reveals the existence of small-amplitude waves that are quickly
damped due to friction).
Given the quick adjustment time of this flow, it may be interesting to look
for stationary solutions of (4.24). In this case, the acceleration terms vanish by
definition, and (4.4) requires v = 0 everywhere for the cross-channel velocity.
Without any further assumptions, (4.24) can thus be written as
BDS = −CD u2 ,
f u = B(D ′ − H ′ ) ,
(4.25)
where the primes denote derivatives with respect to y (in the cross-channel direction). Equation (4.25) is recognized as (a) a balance between bottom friction and
the down-channel pressure gradient, and (b) a purely geostrophic balance due to
the cross-channel tilt of the interface. Eliminating u from (4.25), we obtain
DSf 2
= (D ′ − H ′ )2 ,
(4.26)
CD B
which can be rewritten as a first-order differential equation for the plume thickness
12
1
S
′
D =−
D 2 + H′ ,
(4.27)
δG
where, for convenience, we have introduced the frictional length scale
CD B
δG = − 2
.
(4.28)
f
−
For the linear drag law (4.6), as shown in the assignments, we would have obtained
an equation that is rather similar to (4.27) but somewhat easier to handle because
the non-linearity in D disappears.
It is evident from (4.27) that for any given cross-channel bathymetry H the
solution only depends on the parameter S/δG that involves a combination of
all other parameters of the problem. Numerical solutions of (4.27) for some
values of this parameter are shown in Figure 4.6. In the first panel, the value
of S/δG corresponds to the parameter set given the caption of Figure 4.5. The
same cosine-shaped topography has been used, and, as expected, the stationary
solution in Figure 4.6 and the time-dependent solution for large times shown
in Figure 4.5 closely correspond. Larger values of S/δG (e.g. due to steeper
slopes, smaller friction, smaller buoyancy contrast) lead to thinner and faster
plumes with a larger cross-channel tilt. Stationary solutions of this type have
been compared with good success to a number of real oceanic gravity currents of
largely different scales and parameter ranges.
34
z (m)
0
−5
−10
8
−1
S/δ × 10 = 5 m
−15
z (m)
0
−5
−10
8
−1
8
−1
S/δ × 10 = 20 m
−15
z (m)
0
−5
−10
−15
−10
S/δ × 10 = 80 m
−5
0
y (km)
5
10
Figure 4.6: Upper panel: solution of (4.27) for the parameters given in the caption
of Figure 4.5. Lower panels correspond to different values of the parameter S/δG
as indicated. Note that in contrast to Figure 4.5 the gray shaded areas only
indicate the fluid inside the gravity current, and not the down-channel speed.
35
Chapter 5
Entrainment
5.1
Introduction
It is frequently observed that an unconfined turbulent region has a tendency to
expand into quiescent ambient fluid as long as there is a sufficient energy supply
for turbulence. On the small scale, turbulent eddies acting near the boundary
between the turbulent and ambient fluids intermittently ”entrain” small amounts
of ambient fluid, thereby increasing the net volume of the turbulent region. The
growth of clouds provides a good example where these entrainment processes
become directly observable due to the presence of microscopic water droplets that
serve as an optical tracer (Figure 5.1). In this case, entrainment is energized by
turbulent motions associated with plumes of light fluid rising in the interior of the
cloud. Entrainment is also important in many engineering flows, as for example
in the propulsion jet from a rocket engine shown in Figure 5.1, where shear
instabilities near the boundary of the jet form the energy source for turbulence
and entrainment. The expansion of the jet due to entrainment is clearly evident
in this example.
While the description of the complex geometry of the sharp, instantaneous
interface between the turbulent and non-turbulent regions shown in the above
examples is mathematically challenging, it has turned out that the motion of the
averaged interface often follows relatively simple scaling laws. The key component
of such models is a relation between the entrainment velocity, i.e. the average
velocity of the interface with respect to the ambient fluid, and a characteristic
velocity of the turbulent flow (e.g. the mean speed of the jet shown in Figure
5.1). An excellent and exhausting discussion of the entrainment concept including
numerous applications to geophysical flows can be found in paper by Turner
(1986).
In the following, we concentrate on a class of entraining flows that is of particular interest in oceanographic applications. In the upper ocean, turbulence
caused by the wind stress at the water surface is known to lead to a strongly
36
Figure 5.1: Examples for entraining interfaces separating a highly turbulent region from weakly turbulent ambient fluid. Left panel: cumulus cloud. Right
panel: turbulent jet created by a rocket engine.
turbulent, essentially well-mixed surface layer that deepens as a result of entrainment into the stratified, weakly turbulent interior. Due to its importance for the
oceanic heat and momentum budget, this process has been investigated in detail
in the oceanographic literature. In the following, we provide a short introduction
into the most relevant concepts.
5.2
Basic framework
The essential mechanisms of mixed-layer entrainment in the ocean can be understood from simplified versions of the transport equations for momentum (1.8),
temperature (1.23), and salinity (1.24). Since for most conditions mixed-layer
entrainment can be described as a purely vertical mixing process, the major simplifications arise from assuming horizontal homogeneity in all equations. Further,
although Earth’s rotation may be quite important for entrainment in some cases,
it is ignored here because the basic mechanisms and concepts can be illustrated
more clearly without rotational effects. Thus, ignoring all horizontal gradients
and rotation, and assuming that the molecular vertical flux of momentum is much
smaller than the turbulent flux hu′ w ′ i, the momentum budget in (1.8) reduces to
∂u
∂hu′ w ′ i
=−
∂t
∂z
.
(5.1)
Here, and in the following, we omit the angular brackets denoting mean flow
quantities (as in hui), and assume that the flow direction coincides with the
x-direction.
With analogous assumptions, simplified one-dimensional equations can be derived also for temperature and salinity from (1.23) and (1.24), respectively. It is
easy to show (see assignments) that the resulting two equations can be combined
37
into a single equation for the mean density ρ, provided the equation of state is
assumed to be linear as in (1.21). From this, using the definition of buoyancy
given in (1.22), the transport equation for buoyancy is given by
∂hb′ w ′ i
∂b
=−
,
∂t
∂z
(5.2)
where we have ignored the effect of radiative heating such that the right hand side
corresponds to the vertical divergence of the turbulent buoyancy flux. Working
with buoyancy instead of temperature and salinity is often preferable because
buoyancy is the dynamically relevant quantity.
5.3
Evolution of the surface mixed layer
Experience has shown that turbulence associated with the surface wind stress
creates a nearly well-mixed layer of thickness h near the surface. Thus, inside
this layer, buoyancy can only be a function of time: b = bm (t). At the base of
the surface mixed layer, density rapidly increases until the ”interior” density just
below the mixed layer is reached. Such regions with strong vertical density gradients, often denoted as ”pycnoclines”, are a typical phenomenon at the interface
between turbulent and non-turbulent regions in stratified flows. In the following,
the pycnocline at the base of the mixed layer is described in the most simple way:
as a discontinuity in buoyancy denoted by ∆b.
To investigate the effect of stratification below the mixed layer, in the following
we consider two idealized models that will serve as a simple representation of
typical oceanic conditions. The first (Figure 5.2a) is a simple 2-layer fluid with
zero stratification below the mixed layer:
(
bm for z > −h
(5.3)
b=
0 otherwise ,
where z points upwards with z = 0 at the surface. Here, we have arbitrarily
assumed that b = 0 below the mixed layer, which is dynamically irrelevant but
convenient for the following investigations.
In the second case (Figure 5.2b), we assume a constant vertical buoyancy
gradient N 2 = ∂b/∂z below the mixed layer:
(
bm
for z > −h
(5.4)
b=
N 2 z otherwise ,
such that without a mixed layer (h = 0) the whole water column is linearly
stratified, and b = 0 at the surface.
38
Figure 5.2: Model stratification for mixed-layer deepening (a) without and (b)
with stratification below the mixed layer.
For simplicity, we assume that the fluxes of heat and salinity, and thus the
buoyancy flux, vanish at the surface and below the mixed layer. Vertical integration of (5.2) then yields
Z
d 0
b dz = 0 ,
(5.5)
dt −∞
indicating that the total buoyancy in the domain remains constant. This result
can be used to obtain information about the temporal evolution of bm and ∆b
during mixed-layer deepening (ḣ > 0) due to entrainment (see assignments). For
the 2-layer stratification in (5.3), it is easy to show that (5.5) is satisfied only if
bm h = ∆b h = const.
(5.6)
for all times. For the case of linear stratification, (5.4), assuming that h = 0
initially, we find that
1
(5.7)
bm = − N 2 h = −∆b
2
is required to satisfy (5.5). Thus, while in the 2-layer case the buoyancy jump at
the base of the mixed layer decreases in time, it increases for linear stratification
below the mixed layer. This difference has considerable implications for the
energetics of entrainment as discussed in the following.
5.3.1
Energetics
During the entrainment process, dense fluid from the mixed-layer base is lifted
up and mixed with fluid in the surface layer by wind-generated turbulence. This
lifting of dense water implies that work against gravity is performed, or, in other
words, that potential energy is increased at the expense of (turbulent) kinetic
energy in the mixed layer. The ultimate energy source for this process is the
work done by the surface wind stress.
39
To express this energy conversion in mathematical terms, we first note that
the potential energy per unit area is defined as
Z 0
Z 0
g(ρ − ρref )
Ep =
z dz = −
(b − bref ) z dz ,
(5.8)
ρ0
−∞
−∞
where ρref (z) and bref (z) refer to a given reference state with respect to which
Ep is defined. In (5.8), we have also assumed that Ep is referenced with respect
to the surface, z = 0, which does not affect any of the following results but is
convenient for the analysis.
For the 2-layer case in (5.3), we identify the reference state with the initial
stratification at t = 0, corresponding to bref = b0 inside a mixed-layer of thickness
h0 , and bref = 0 below. Then, from (5.6) and (5.8), it can be shown (see assignments) that the increase of potential energy with respect to the reference state
is Ep = b0 h0 (h − h0 )/2. From this, using (5.6), it follows that
∆bh
ḣ ,
(5.9)
2
illustrating that the rate at which Ep increases is proportional to the entrainment
velocity we = ḣ. Similarly, for the linearly stratified case, taking the undisturbed
initial stratification bref = N 2 z as the reference, it can be shown that the potential
energy follows from Ep = N 2 h3 /12. Using (5.7), a short computation illustrates
that the rate of change of Ep is identical to that of the 2-layer case in (5.9). In
spite of this formal similarity, however, it is important to realize that the temporal
evolution of ∆b is very different for the two cases. In the linearly stratified case
the buoyancy jump ∆b increases in time, whereas for a 2-layer fluid, it decreases.
For any given rate of mixed layer deepening ḣ, the rate Ėp at which turbulent
motions have to work against gravity thus increases more rapidly in the linearly
stratified case compared to the 2-layer stratification.
A general expression for the change of potential energy for arbitrary stratification can be obtained from multiplying (5.2) with z, and integrating by parts:
Z 0
Ėp = −
hb′ w ′i dz ,
(5.10)
Ėp =
−∞
where we have exploited the fact that the buoyancy flux hb′ w ′i vanishes at both
the upper and lower integration limits. Expression (5.10) illustrates that any
increase in potential energy (Ėp > 0) is associated with a net downward buoyancy
flux (hb′ w ′i < 0). This mirrors the work done by turbulent eddies in moving
buoyant fluid down (or, vice-versa, dense fluid up). Of key interest for describing
this process with the help of a theoretical model is therefore the question about
the energy source for this mixing process.
In order to find a relation for the kinetic energy as the second important form
of mechanical energy, we vertically integrate (5.1) resulting in
d(uh)
= u2∗ ,
dt
40
(5.11)
where the surface wind stress τs has been expressed in terms of the (squared)
friction velocity u2∗ = τs /ρ0 , which is a frequently used convention in studies of
boundary-layer turbulence. The stress below the mixed layer has been assumed
to be negligible. In this simple model, the surface layer moves in a slab-like
fashion with vertically constant velocity u. Multiplying (5.11) by u, it can be
shown that
dEkin
u2
= u2∗ u − ḣ ,
(5.12)
dt
2
where the vertically integrated kinetic energy is denoted by Ekin = u2 h/2. The
right hand side of (5.12) reveals that part of the surface energy flux u2∗ u is used
to accelerate newly entrained fluid at the bottom of the mixed layer (last term).
In this context, it is very important to note that our assumption of a vertically
constant velocity u eliminates any transfer from mean to turbulent kinetic energy
by turbulence shear production, and ultimately to dissipation. In this type of
mixed-layer models, dissipation is thought to take place in the strongly sheared
boundary layers at the top and bottom of the mixed layer that are not explicitly
considered (although they are crucial for the total energy budget). The term u2∗u
in (5.12) thus has to be interpreted as the energy flux into the mean flow rather
than the total energy input into the whole mixed layer.
5.4
Entrainment models
If the wind stress is assumed to be constant in time and the flow in the mixed-layer
has reached a quasi-steady state, the entrainment velocity we = ḣ is thought to
primarily depend on only 3 parameters: the friction velocity u∗ , the mixed layer
depth h, and the buoyancy jump ∆b at the base of the mixed layer. In some
cases, it may be necessary to include additional parameters describing the stratification below the mixed layer because internal wave motions in that region may
drain energy from the mixed layer and thus affect the entrainment process. This
effect, however, is ignored in most studies. Also ignored is the direct impact of
the molecular diffusivities for momentum, heat, and salinity because their effect
is generally believed to be negligible in the strongly turbulent mixed layer. Note
that the assumption of negligible molecular effects is not justified in most laboratory experiments of entrainment, and caution is always advised when translating
such experiments to oceanic conditions.
Under these assumption, 2 non-dimensional products can be constructed for
the 4 dimensional key parameters. The most frequently adopted non-dimensional
parameters are the entrainment parameter Eτ and the frictional Richardson number Rτ :
∆bh
ḣ
(5.13)
, Rτ = 2 ,
Eτ =
u∗
u∗
41
for which we expect to find a non-dimensional relationship of the form
Eτ = Êτ (Rτ ) ,
(5.14)
If the wind stress is variable, or if the system has not yet reached a steady state,
the function Ê has to include information about the history of the processes (e.g.,
about the history of the wind stress). As discussed in the following, modeling
assumptions and observations have led to different suggestions for the functional
relationship between Eτ and Rτ .
5.4.1
Models based on energy partitioning
Early models for wind-driven entrainment were based on the simple idea that
a (fixed) portion of the energy imparted by the wind is converted to potential
energy by mixing up fluid from the base of the mixed layer. The rest of the energy
is either dissipated, or used to accelerate the mixed layer as described in (5.12).
While it is straightforward to compute the increase in potential energy due to
entrainment as shown above, a model for the energy input due to the wind is less
evident because the partitioning between Ekin and the amount of energy that is
dissipated is unclear. Most authors have assumed that the wind energy input
is proportional to u3∗ , which implicitly assumes that most of the wind energy is
dissipated in a strongly sheared near-surface layer, the rest being available for
mixing.
Accepting this scaling, the assumed proportionality between the wind energy
input and the gain in potential energy can be expressed as
u3∗ = c∆bhḣ ,
(5.15)
where c is a model constant, and the right hand side follows from (5.9). Rewriting
(5.15) in non-dimensional form yields
1
Eτ = Rτ−1 ,
c
(5.16)
which is recognized as a special incarnation of (5.14) with an inverse dependency
on the frictional Richardson number.
In the 2-layer case with constant ∆bh, (5.15) can be integrated to yield
h=
u3∗
t + h0 ,
cb0 h0
(5.17)
where h0 and b0 are the mixed-layer depth and buoyancy, respectively, for t = 0.
The linear increase of h reflects the constant entrainment velocity that follows
directly from (5.16) and the fact that both u∗ and Rτ are constant in this case
(for constant wind stress).
42
In the linearly stratified case, we have ∆bh = N 2 h2 /2 as shown above, and
(5.15) integrates to
31
2 1
6
(5.18)
u∗ N − 3 t 3 ,
h=
c
for which Kato and Phillips (1969) have suggest the value c = 2/5. Compared to
the constant entrainment velocity in (5.17), the linearly stratified case exhibits
decreasing entrainment in time, which mirrors the rapidly increasing rate of work
against gravity required for entrainment.
5.4.2
Models based on stability considerations
A different class of entrainment models is based on the idea that the interplay
between the destabilizing effect of the shear across the mixed-layer base, and the
stabilizing effect of the pycnocline represented by ∆b, governs the entrainment
process. Instead of scaling the non-dimensional parameters with the friction
velocity u∗ as in (5.13), in this approach non-dimensional parameters are constructed with the help of the velocity jump ∆u = u (the lower layer is assumed
to be stagnant), while the rest of the parameters remains identical. Analogous
to (5.13), the entrainment parameter and the bulk Richardson number can be
defined according to
ḣ
∆bh
E = , Rv = 2 ,
(5.19)
u
u
suggesting that a non-dimensional relationship of the form
E = Ê(Rv ) ,
(5.20)
quantifies the entrainment process. Often, the dependency on Rv is replaced
−1/2
by an (equivalent) dependency on the bulk Froude number F r = Rv . The
relationship in (5.20) should not be understood as a replacement for the generally
valid relation (5.14). Rather, as shown in the following, (5.20) corresponds to a
class of models that can ultimately be used for deriving expressions of the form
given in (5.14).
One advantage of (5.20) is its applicability to a broad range of entraining
flows that do not necessarily involve the surface friction velocity as an important
parameter. Examples include surface and bottom mixed layers, dense bottom
gravity currents, and buoyancy-driven intrusions. However, in spite of numerous
laboratory and field investigations, a generally valid functional form for (5.20) is
still lacking. Most authors seem to agree on a power-law dependency of the form
E ∝ Rvp ,
(5.21)
with negative exponent p < 0. In any case, if the functional form of (5.20) is
prescribed, the problem is closed, and the temporal evolution of momentum in
43
(5.11) and buoyancy in (5.6) or (5.7), respectively, can be determined with the
help of (5.20).
Some interesting insights into the general behavior of stability-based entrainment laws can be gained by investigating a particularly simple form (5.20). This
model is based on the experimental observation that the mixed layer velocity u
adjusts, as a result of entrainment, exactly such that the bulk Richardson number
remains close to a limiting value Rv ≈ 0.6. This closure assumption is mathematically expressed as Ṙv = 0, or, using the definition of the bulk Richardson
number in (5.19) as
Ṙv = 0
⇒
(∆bh)˙ 2∆bhu̇
=
u2
u3
.
(5.22)
Replacing the time derivative of ∆bh by the expressions corresponding to the
2-layer and linearly stratified cases, and using
1 2
u∗ − uḣ ,
(5.23)
u̇ =
h
derived from the momentum balance in (5.11), (5.22) can be re-written as
ḣ = n
u2∗
,
u
(5.24)
where n = 1 corresponds to the 2-layer case, and n = 1/2 to the linearly stratified
case. One important observation is that (5.24) can be re-written as
1
− 12
E∗ = nRv2 Rτ
,
(5.25)
which is recognized as a special form of (5.14), derived from the simple model
assumption that the interface is marginally stable for a fixed value of Rv .
Assuming Ṙv = 0 and using the different expressions for ∆bh in (5.6) and
(5.7), respectively, (5.25) can be integrated to yield the temporal evolution of the
mixed-layer depth. For the 2-layer case, this results in
1
Rv 2 2
h=
u∗ t + h0 ,
(5.26)
b0 h0
which, similar to the energy-based model in (5.17), predicts a linear increase in
time. Note, however, that the dependency on the wind stress is different. For
the linearly stratified case, we obtain
1
1
h = (2Rv ) 4 u∗ (t/N) 2
.
(5.27)
with a different temporal behavior compared to (5.18). It has been shown that
(5.26) and (5.27) predict the correct evolution of h, at least during the initial
phase of mixed-layer deepening.
44
5.5
Further reading
The key paper describing the entrainment concept in a broad geophysical context
is Turner (1986). The classical references for energy-based entrainment models
are Turner and Kraus (1967) and Kraus and Turner (1967), after which this class
of models is often referred to as the “Krauss-Turner” model. Pollard et al. (1973)
have pointed out that shear instability at the mixed layer base is likely to be the
primary reason for entrainment, similar to entrainment in dense bottom gravity
currents discussed by Turner (1986). Pollard et al. (1973) have also suggested to
use a constant value for Rv as the most simple stability-based entrainment model,
which was later confirmed by Price (1979) from a reanalysis of laboratory experiments. The importance of rotational effects has been emphasized by Pollard et al.
(1973) and Niiler (1975). In the review paper by Umlauf and Burchard (2005),
stability-based models for surface layer entrainment are compared with the predictions from state-of-the-art turbulence models.
45
46
Chapter 6
Surface waves in the coastal
ocean
6.1
Observational evidence
When the surface of a body of water is disturbed in the vertical direction, the
force of gravity will act to return the surface to its equilibrium position. The
returning surface water has inertia that causes it to pass its equilibrium position
and establish a surface oscillation. This oscillation disturbs the adjacent water
surface, causing the forward propagation of a wave. A wave on the water surface is thus generated by some disturbing force which may typically be caused
by the wind, a moving vessel, a seismic disturbance of the shallow sea floor, or
the gravitational attraction of the sun and moon. These forces impart energy to
the wave which, in turn, transmits the energy across the water surface until it
reaches some obstacle such as a structure or the shoreline which causes the energy
to be reflected and dissipated. The wave also transmits a signal in the form of
the oscillating surface time history at a point. As a wave propagates, the oscillatory water motion in the wave continues because of the interaction of gravity
and inertia. Since water particles in the wave are continuously accelerating and
decelerating as the wave propagates, dynamic pressure gradients develop in the
water column. These dynamic pressure gradients are superimposed on the vertical hydrostatic pressure gradient. As the wave propagates energy is dissipated,
primarily at the airwater boundary and, in shallower water, at the boundary between the water and the sea floor. The different wave generating forces produce
waves with different periods. Wind-generated waves have a range of periods from
about 1 to 30 s with the dominant periods for ocean storm waves being between
5 and 15 s. Vessel-generated waves have shorter periods, typically between 1 and
3 s. Seismically generated waves (tsunamis) have longer periods from about 5
minutes to an hour and the dominant periods of the tide are around 12 and 24
hours.
47
The simplest and often most useful theory (considering the effort required in its
use) is the two-dimensional small-amplitude or linear wave theory first presented
by Airy (1845). This theory provides equations that define most of the kinematic
and dynamic properties of surface gravity waves and predicts these properties
within useful limits for most practical circumstances. The assumptions required
deriving the small-amplitude theory, an outline of its derivation, the pertinent
equations that result, and the important characteristics of waves described by
these equations are presented in this chapter.
6.1.1
Wave height and period
In the description of wind waves it is common to define the wave height H as the
vertical distance between the highest and the lowest surface elevation (crest to
trough) in a wave. A wave will thus have only one height, with the definition of
the mean wave height H as
N
1 X
Hi
(6.1)
H=
N i=1
where i is the sequence of the wave in a record.
This definition if rather straight forward, but not often used, because H deviates
from visual estimated wave heights. Instead, another wave height, called the
significant wave height HS (just as in visual observations is used. It is defined as
the mean of the highest one-third of waves in a wave record:
N/3
HS = H1/3
1 X
=
Hj
N/3 j=1
(6.2)
where j is the rank number of the sorted wave heights (j=1 the highest wave,
j=2 the second-highest, etc.). One can also define the amplitude a as H1/3 /2.
It is equally natural to define the period T of a wave as the time interval between
the start and the end on the wave (time between two crests/troughs, two downward/upward zero crossings). If the wave period is defined as the zero-crossings,
it is called T0 and in analogy tho the mean wave height H, one can define the
mean wave period T , as
N
1 X
T =
T0,i
(6.3)
N i=1
Mostly, only the significant wave period TS is used:
N/3
TS = T1/3
1 X
T0,j
=
N/3 j=1
48
(6.4)
6.1.2
The random-phase/amplitude model
The basic model for describing the moving surface elevation η(t) is the randomphase/amplitude model, in which the surface elevation is considered as the sum
of a large number of harmonic (monochromatic) waves, each with a constant
amplitude ai and phase ψi :
η(t) =
N
X
ai cos(2πfi t + ψi )
(6.5)
i=1
with fi the frequency. It is important to note, that ai and ψi are random variables, which are characterised by there respective probability distributions p. In
the random-phase/amplitude model, the phase at each frequency fi is uniformly
distributed between 0 and 2π:
p(ψi ) =
1
2π
for 0 < ψi ≤ 2π
and the amplitude ai follows a Rayleigh distribution:
πai 2
π ai
exp − 2
for ai ≥ 0
p(ai ) =
2π µ2
4µ
(6.6)
(6.7)
with µ the expected value of the amplitude µ = hai i.
For typical oceanic applications fi vary between 0.05-1.0 Hz.
Although the random-phase/amplitude model is widely used, the following remarks should be made:
ˆ The random-phase/amplitude model generates a stationary Gaussian process. But this is only true for durations of less than 30 minutes. In additions,
the model rely on independed wave components (a sum over independed
monochromatic waves), which is not fully true, because they interact to
some degree. However, if the waves are not to steep and the water is not to
shallow, the random-phase/amplitude model can be used as a good basic
model to describe ocean waves.
ˆ In Eq. 6.5 the summation is done of discrete frequencies, but in reality this
would be a continuum of frequencies.
6.1.3
The variance density spectrum
The amplitude spectrum provides enough information to describe the sea-surface
elevation realistically (the phase is not important). However, it is more relevant
to present the spectral information in a different way: consider the variance
h 12 ai 2 i, rather than the mean hai i. Thus, we consider the variance spectrum and
not the amplitude spectrum. Netherveless, we still consider discrete frequencies,
49
whereas in Nature this concept does not exist. Thus the discrete model needs
to be modified. This is done by distributing the variance h 21 ai 2 i over frequency
intervals ∆fi at frequency fi . The resulting variance density spectrum E ∗ (fi ) is
then,
1 1 2
h ai i
(6.8)
E ∗ (fi ) =
∆fi 2
This spectrum is defined for all frequencies, but still varies discontinuously from
one frequency band to the next one. A continuous version is obtained by having
the frequency band ∆fi approaching zero:
1 1 2
h ai
∆f →0 ∆f 2
E(fi ) = lim
(6.9)
The variance density spectrum E(f ) gives a complete description of the surface
elevation in a statistical sense. Thus all statistics of the wave field can be expressed in terms of the spectrum.
To show how useful this concept is, one needs to recall that the variance of the
surface elevation η(t) by definition the average of the squared surface elevation
η 2 is (assuming zero mean). For a harmonic wave with amplitude a, the variance
is η 2 = 21 a2 . Thus for a wavefield, with a large number of individual waves, the
variance η 2 is:
N
X
1
2
η =
h ai 2 i for hηi = 0
2
i=1
Therefore it follows that the total variance η 2 in the continuos spectrum is
Z ∞
2
η =
E(f )df
(6.10)
0
Thus by knowing the spectrum, simple integration will give us the variance. A
quantity which is directly related to the variance is the energy Etotal .
The energy can be expressed as
1
Etotal = ρg a2
2
(6.11)
with ρ the water density and g the gravitational acceleration. Because we defined
above η 2 = 21 a2 , the energy of a wave field simply be computed as:
Z ∞
2
Etotal = ρgη = ρg
E(f )df
(6.12)
0
For simplicity, in the following sections we will use only monochromatic waves,
thus the variance density spectrum can be expressed with a delta peak. Although,
this is a simplification, by considering the superposition of a large number of
monochromatic waves, we can recover the random-phase/amplitude model.
50
6.2
Small-amplitude wave theory
The small-amplitude theory for two-dimensional, freely propagating, periodic
gravity waves is developed by linearising the equations that define the free surface
boundary conditions. With these and the bottom boundary condition, a periodic
velocity potential is sought that satisfies the requirements for irrotational flow.
This velocity potential, which is essentially valid throughout the water column
except at the thin boundary layers at the airwater interface and at the bottom,
is then used to derive the equations that define the various wave characteristics (e.g., surface profile, wave speed, pressure field, and particle kinematics).
Specifically, the required assumptions are:
1. The water is homogeneous and incompressible, and surface tension forces
are negligible. Thus, there are no internal pressure or gravity waves affecting
the flow, and the surface waves are longer than the length where surface
tension effects are important (i.e., wave lengths are greater than about 3
cm).
2. Flow is irrotational. Thus there is no shear stress at the airsea interface
or at the bottom. Waves under the effects of wind (being generated or
diminished) are not considered and the fluid slips freely at the bottom and
other solid fixed surfaces. Thus the velocity potential φ must satisfy the
Laplace equation for two-dimensional flow:
∂2φ ∂2φ
+ 2 =0
∂x2
∂z
and
∂φ
=u ,
∂x
∂φ
=w
∂z
(6.13)
(6.14)
3. The bottom is stationary, impermeable, and horizontal. Thus, the bottom is not adding or removing energy from the flow or reflecting wave
energy. Waves propagating over a sloping bottom, as for example when
waves propagate toward the shore in the nearshore region, can generally be
accommodated by the assumption of a horizontal bottom if the slope is not
too steep.
4. The pressure along the airsea interface is constant. Thus, no pressure is
exerted by the wind and the aerostatic pressure difference between the
wave crest and trough is negligible.
5. The wave height H is small compared to the wave length L and water depth
D.
H
L
H
≪1 ,
≪1 ,
≪1
(6.15)
L
D
D
51
is also called the wave steepness. If this ratio is bigger than
The ratio H
L
1/8, waves usually start to break.
This assumption allows one to linearise the higher order free surface boundary conditions and to apply these boundary conditions at the still water line
rather than at the water surface, to obtain an easier solution. This assumption means that the small-amplitude wave theory is most limited for waves
in deep water and in shallow water far away from the breaking point.
6.
The small-amplitude wave theory is developed by solving Eq. 6.13 for the a 2D domain (we neglect all variations in the y-direction), with the appropriate
boundary conditions for the free surface and the bottom.
At the bottom there is no flow perpendicular to the bottom which yields the
bottom boundary condition:
w=
∂φ
= 0 at z = −D
∂z
(6.16)
At the free surface there is a kinematic boundary condition that relates the vertical component of the water particle velocity at the surface to the surface position:
w=
∂η
∂η
+u
∂t
∂x
at z = η
(6.17)
The Bernoulli equation for unsteady irrotational flow may be written
p
1 2
∂φ
u + w 2 + + gz +
=0
2
ρ
∂t
(6.18)
where g is the acceleration of gravity, p is the pressure, and ρ is the fluid density.
At the surface where the pressure is zero the dynamic boundary condition becomes:
∂φ
1 2
u + w 2 + gz +
= 0 at z = η
(6.19)
2
∂t
Eq. 6.20 and the Eq. 6.19 have to be linearised and applied at the still water line
rather than at the a priori unknown water surface. This yields for the kinematic
boundary condition (Eq. 6.20)
w=
∂η
∂t
at z = 0
(6.20)
and for the dynamic boundary condition becomes (Eq. 6.19)
1 2
∂φ
u + w 2 + gz +
= 0 at z = η
2
∂t
52
(6.21)
and for the dynamic boundary condition
gη +
∂φ
= 0 at z = 0
∂t
(6.22)
Employing the Laplace equation (Eq. 6.13), the bottom boundary condition (Eq.
6.16), and the linearised dynamic boundary condition (Eq. 6.22), one can derive
the velocity potential for the small-amplitude wave theory (see Ippen, 1966 or
Dean and Dalrymple, 1984). The most useful form of this velocity potential is
φ=
ga cosh (k(D + z))
sin (kx − ωt)
ω
cosh(kD)
(6.23)
with a the wave amplitude a = H/2.
Inserting the velocity potential into the linearised dynamic boundary condition
(Eq. 6.22) with z=0 to directly determine the equation for the wave surface
profile:
η = a cos (kx − ωt)
(6.24)
Thus, the small amplitude wave theory yields a cosine surface profile. This is
reasonable for low amplitude waves, but with increasing wave amplitude the
surface profile becomes vertically asymmetric with a more peaked wave crest and
a flatter wave trough.
Combining Eq. 6.20 and Eq. 6.22 by eliminating the water surface elevation yields
∂2φ
∂φ
+g
= 0 at z = 0
2
∂t
∂z
(6.25)
Then, inserting the velocity potential, differentiating, and rearranging we have
ω 2 = gk tanh (kD)
(6.26)
This equation is commonly known as the dispersion equation. For a spectrum of
waves having different periods (or lengths), the longer waves will propagate at a
higher speed and move ahead while the shorter waves will lag behind.
6.3
Wave kinematics and pressure
The particle velocities can readily be obtained form the velocity potential φ (Eq.
6.23) and the definition of the velocities (Eq. 6.14):
u=
∂φ
cosh (k(D + z))
= ωa
sin(ωt − kx)
∂x
sinh(kD)
(6.27)
w=
sinh (k(D + z))
∂φ
= ωa
cos(ωt − kx)
∂z
sinh(kD)
(6.28)
and
53
In deep water, i.e. when kD → ∞, the expression for the wave velocity component
u and w reduces to
u = ωaekz sin(ωt − kx) ,
w = ωaekz cos(ωt − kx)
(6.29)
These expressions show that, for deep water, the wave induced velocities decrease
exponentially with the distance from the surface.
In shallow water, i.e. when kH → 0, the expression for the wave velocity component u and w reduces to
z
ωa
cos(ωt − kx)
(6.30)
sin(ωt − kx) , w = ωa 1 +
u=
kD
D
The first of these two expressions show that, in very shallow water, the amplitude
of the horizontal velocity is constant over the vertical, whereas the second expression shows that the amplitude of the vertical varies linearly along the vertical.
6.3.1
The dispersion relationship
As given above the dispersion equation (Eq. 6.26) reads as:
ω 2 = gk tanh (kD)
or
gT 2
tanh
L=
2π
2πD
L
(6.31)
(6.32)
For deep water (tanh (kD) → 1 for kD → ∞) the dispersion relationship approaches:
p
gT0 2
ω = gk0 or L0 =
(6.33)
2π
with k0 and L0 the deep water wave number and wave length respectively.
For shallow water (tanh (kD) → kD for kD → 0) the dispersion relationship
approaches:
p
p
ω = k gD or L0 = T gD
(6.34)
6.3.2
Phase velocity and group velocity
The propagation speed of the surface wave profile, i.e. the phase speed, is readily
obtained from the dispersion relationship
r
g
g
tanh (kD)
(6.35)
c = tanh (kD) =
ω
k
54
This shows that, in general the phase speed depends on the wave number and
therefore on frequency. Such waves are called dispersive waves.
In deep water (tanh (kD) → 1 for kD → ∞) the expression reduces to:
r
g
g
or c0 =
T0
(6.36)
c0 =
k0
2π
In shallow water (tanh (kD) → kD for kD → 0) the expression reduces to:
p
c = gD
(6.37)
This equation shows that, in shallow water, the phase speed does not depend on
wave length of frequency. Thus the waves are said to be non-dispersive.
Note: In engineering applications, the phase velocity is often called celerity.
A second important quantity to describe surface waves is the group-velocity
cgroup = cg =
∂ω
= nc
∂k
(6.38)
where c is the phase speed (Eq. 6.35) and n is (from the dispersion relationship
Eq. 6.26)
2kD
1
1+
(6.39)
n=
2
sinh(2kD)
Since 0 ≤ kH ≤ ∞ and therefore 0 ≤ 2kD/ sinh(2kD) ≤ 1, the expression for
n shows that n varies between n = 0.5 (deep water) and n = 1 (shallow water).
This implies that the speed of an individual wave (phase speed) is always larger
than the speed of the group:c ≥ cg .
6.3.3
Wave induced pressure
An analytical expression for the pressure is readily derived by substituting the
solution for the velocity potential (Eq. 6.23) into the linearised Bernoulli equation
p
∂φ
+ gz +
=0
ρ
∂t
(6.40)
which yields
p = −ρgz + ρga
cosh (k(D + z))
sin (ωt − kx)
cosh(kD)
(6.41)
The first term on the right hand side is the hydrostatic pressure. It is obviously
independent of the presence of a wave (at least in first order theory). The second
55
term is due to the wave and therefore represents the wave induced pressure, pwave
In deep water the wave induced pressure is
pwave = pw sin (ωt − kx)
with pw = ρga
cosh (k(D + z))
cosh(kD)
(6.42)
This expression is a propagating pressure wave in phase with the surface elevation.
In deep the amplitude of the wave induced pressure is
pw = ρgaekz
(6.43)
These expressions show that, for deep water, the wave induced pressure decrease
exponentially with the distance from the surface.
In shallow the amplitude of the wave induced pressure amplitude is constant
along the vertical:
pw = ρga
(6.44)
6.4
Wave energy
The total mechanical energy in a surface gravity wave is the sum of the kinetic
and potential energies. The instantaneous potential energy (i.e. mass times
elevation, at a given moment of time) of a slice of water, relative to z = 0, is
therefore ρgz∆x∆z∆y. The corresponding wave induced potential energy EP in
the entire water column, is equal to the potential energy presence of the wave
minus the potential energy in the absence of the wave over a wave period. Thus
Z TZ η
Z TZ 0
Z TZ η
EP =
ρgzdzdt −
ρgzdzdt =
ρgzdzdt
(6.45)
0
−D
0
−D
0
0
For a harmonic wave with amplitude a, this integral reads as
Z T
1 2 1
ρη = ρga2
EP =
2
4
0
(6.46)
where η is replaced by Eq. 6.24.
Doing the same exercise for the kinetic energy EK , EK can be written as
Z TZ η
1 2
EK =
ρu dzdt
(6.47)
0
−D 2
Replacing u by Eq. 6.14, the integral can be expressed as
1
(6.48)
EK = ρga2
4
Therefore, in the small amplitude wave theory, EK = EP . The total wave energy
E over on wave period is then
1
E = EK + EP = ρga2
2
56
(6.49)
6.5
6.5.1
Wave propagation
Shoaling
A harmonic wave, propagating over slowly changing topography, retains its frequency, but, since the dispersion relationship is still valid
ω 2 = gk tanh(kD)
(6.50)
it wavelength will decrease, if the depth decreases and the phase speed also decrease
r
g
tanh(kD)
(6.51)
c=
k
The same holds for the group velocity cg
1
cg = nc with n =
2
1+
2kD
sinh(2kD)
(6.52)
As the wave propagates into shallower water, the phase speed approaches the
group velocity and the wave will be less and less dispersive.
Assume a beach with parallel depth contours and normal incidence (perpendicular
to the coastline). Further assume that the whole situation is stationary. This
enables a simple energy balance. At first we have to define the energy flux
P = Ecg . As the wave propagates forward the energy flux at one point must
equal the energy flux at a subsequent point minus the energy added, and plus
the energy dissipated and reflected per unit time between the two points. In
first-order theory and over reasonably short distances it is common to neglect the
energy added, dissipated, or reflected, giving
1
1
P1 = P2 → [Ecg ]1 = [Ecg ]2 → ρga1 2 cg,1 = ρga2 2 cg,2
2
2
so that
a2 =
r
cg,1
a1
cg,2
(6.53)
(6.54)
If one takes the up-wave boundary in deep water and replace the index 1 with ∞
and drop the index 2, one can define the shoaling coefficient Ksh as
r
cg,∞
Ksh =
(6.55)
cg
Therefore, as the wave approaches the shore, the shoaling coefficient will start
to grow and therefore the wave height. However, one can also see as the wave
reaches the shore and cg is 0, the shoaling coefficient shows a singularity and tends
to infinity. Obviously, the linear wave theory breaks down long before that.
57
6.5.2
Refraction
As a wave approaches the same straight coast as above, but now with an angle
(oblique incidence), the wave will slowly change its direction as its approaches
the coast. This is due to the depth variations along the wave crest, with a
corresponding variation of phase speed along the crest
r
g
c=
tanh(kD)
(6.56)
k
Because refraction is a 2 dimensional problem, the wave number turns into a
wave number vector k = (kx , ky ).
Using the analogy to geometrical optics, one can apply Snel’s Law
sinθ
d
=0
(6.57)
dm
c
where m is the normal of the wave crest and theta the angle between m and the
normal of the shore line. This can also be written as
sinθ
= constant
c
(6.58)
Thus the wave direction at any point can be written as
sin θ =
c
sin θ∞
c∞
(6.59)
where c∞ is the deep water phase speed and θ∞ is the deep water wave direction.
This expression shows that, when the wave approaches the shore line, where
c = 0, the wave direction will be θ = 0. Thus all waves always reach the shore at
a right angle, independently of the deep water direction.
6.5.3
Currents
If a wave propagates in an area with constant depth and a constant current, the
linear theory is still valid in a frame of reference moving with the current.
The frequency of the wave moving in this frame of reference is called relative of
intrinsic frequency, denoted as σ and is expressed as
σ 2 = gk tanh(kD)
(6.60)
In a fixed frame of reference, the frequency of the wave is called absolute frequency,
denoted as ω and is expressed as
ω = σ + kUn
58
(6.61)
where Un is the component of the current in the wave direction (normal to the
crest). The propagation velocity of the wave energy cg,absolute in this fixed frame
→
−
of reference, is obtained by adding as vectors the current velocity U to the group
velocity relative to the current, cg,relative
→
−
cg,absolute = U + cg,relative
(6.62)
The direction of wave energy transport is therefore generally not normal to the
wave crest in the presence of an ambient current.
6.6
Radiation stress and wave setup
In fluid flow problems, some analyses are best carried out by energy considerations
(e.g., head loss along a length of pipe) and some by momentum considerations
(e.g., force exerted by a water jet hitting a wall). Similarly, for waves it is
better to consider the flux of momentum for some problem analyses. For wave
analyses, the flux of momentum is commonly referred to as the wave radiation
stress which may be defined as the excess flow of momentum due to the presence
of waves Longuet-Higgins and Stewart (1964). Problems commonly addressed
by the application of radiation stress include the lowering (setdown) and raising
(setup) of the mean water level that is induced by waves as they propagate into
the nearshore zone, the interaction of waves and currents, and the alongshore
current in the surf zone induced by waves obliquely approaching the shore.
6.6.1
Radiation stress
The horizontal flux of momentum at a given location consists of the pressure
force acting on a vertical plane normal to the flow plus the transfer of momentum
through that vertical plane. The latter is the product of the momentum in the
flow and the flow rate across the plane. From classical fluid mechanics, the
momentum flux from one location to another will remain constant unless there
is force acting on the fluid in the flow direction to change the flux of momentum.
If we divide the momentum flux by the area of the vertical plane through which
flow passes, we have for the x direction
p + ρu2
(6.63)
For a wave, we want the excess momentum flux owing to the wave, so the radiation
stress Sxx for a wave propagating in the x direction becomes
Z 0
Z η
2
(p + ρu ) dz −
ρgdz
(6.64)
Sxx =
−D
−D
59
where the subscript xx denotes the x-directed momentum flux across a plane
defined by x =constant. In Eq. 6.64 p is the total pressure given by Eq. 6.41 so
the static pressure must be subtracted to obtain the radiation stress for only the
wave. The overbar denotes that the first term on the right must be averaged over
the wave period. Inserting the pressure (Eq. 6.41 and the particle velocity from
Eq. 6.14 leads to (Longuet-Higgins and Stewart (1964))
1
2kD
ρga2 1
= E 2n −
(6.65)
+
Sxx =
2
2 sinh(kD)
2
For a wave traveling in the x-direction there also is a y-directed momentum flux
across a plane defined by y =constant. This is
ρga2
kD
1
=E n−
(6.66)
Sxx =
2
sinh(kD)
2
The cross components Sx y = Sy x are both zero.
In deep water Eqs. 6.65 and 6.66 become
Sxx =
E
2
,
Syy = 0
(6.67)
In shallow water Eqs. 6.65 and 6.66 become
Sxx =
3E
2
,
Syy =
E
2
(6.68)
so, like wave energy E, the radiation stress changes as a wave propagates through
water of changing depth (as well as when a force is applied).
If a wave is propagating in a direction that is situated at an angle θ to the specified
x-direction, the radiation stress components become
Sxx = E n (cos2 θ + 1) − 12 Syy = E n sin2 θ + 1 − 21
(6.69)
2
E
Sxy = 2 n sin θ = Syx
6.6.2
Wave setup
When a train of waves propagates toward the shore, at some point, depending
on the wave characteristics and nearshore bottom slope, the waves will break.
Landward of the point of wave breaking a surf zone will form where the waves
dissipate their energy as they decay across the surf zone.
As the waves approach the breaking point there will be a small progressive set
down of the mean water level below the still water level. This setdown is caused
by an increase in the radiation stress owing to the decreasing water depth as the
waves propagate toward the shore. The setdown is maximum just seaward of the
60
breaking point. In the surf zone, there is a decrease in radiation stress as wave
energy is dissipated. This effect is stronger than the radiation stress increase owing to continued decrease in the water depth. The result is a progressive increase
or setup of the mean water level above the still water level in the direction of the
shore. This surf zone setup typically is significantly larger than the setdown that
occurs seaward of the breaking point.
The equations that predict the wave-induced nearshore setdown and setup can be
developed by considering the horizontal momentum balance for two-dimensional
waves approaching the shore (Longuet-Higgins and Stewart (1964)). The net
force caused by the cyclic bottom shear stress is reasonably neglected.
∂Sxx
∂η
∂Qx
=−
− ρgD
(6.70)
∂t
∂x
∂x
with D the water depth D = d + η.
For stationary situations, thus all time derivatives are zero; this can be reduced
to an ordinary differential equation
dη
dSxx
+ ρg (d + η)
=0
dx
dx
(6.71)
If one further assumes that η ≪ d,
dη
1 dSxx
=−
dx
ρgD dx
(6.72)
This implies that, if the radiation stress is positive, dSxx /dx > 0, the slope of the
mean surface is negative, dη/dx < 0 (set-down), and if the if the radiation stress
is negative, dSxx /dx < 0, the slope of the mean surface is positive, dη/dx > 0
(set-up).
If no energy dissipation occurs (one might discuss about this assumption), than
the variation of the wave amplitude is due solely to shoaling. Integrating Eq. 6.71,
taking η = 0 in the deep water, using the shoaling equation for the amplitude
(Eq. 6.54), one can find (Longuet-Higgins and Stewart (1964))
η=−
a2 k
1
2 sinh(2kD)
(6.73)
The minus sign in this expression shows that shoaling lowers the mean sea level,
thus a set-down occurs. In the shallow water limit, this can be approximated
with
1 H2
(6.74)
η≈−
16 D
This equation indicates that the set-down is proportional to the square of the wave
height (which is generally increasing towards the shore, if no breaking occurs),
so that the setdown increases as the wave propagates towards the shore. At
61
the location where the wave starts to beak, the wave height-to-depth ration is
typically Hbr /dbr = γ ≈ 0.8, so that the set-down at this point is 4-5% of the
local water depth or wave height. At the point of breaking, the amplitude starts
to decrease. This implies that dSxx /dx < 0 and the slope of the mean water
surface changes sign and a set-down starts to turn into a set-up.
3
= 16
ρgH 2 and also
If one again assumes shallow water, and this Sxx = 3E
2
that the wave height remains equal to a fixed fraction of the local water depth
H = γ(d + η), it follows that
Sxx =
3
ρgγ 2 (d + η)2
16
(6.75)
Substituting this result into Eq. 6.71, gives
dη
dd
= −K
dx
dx
(6.76)
with K = 83 γ 2 /(1 + 38 γ 2 ).
Since the water depth decreases towards the shore, the mean surface elevation in
the surf zone tilts up towards the shore (dη/dx > 0)
62
Chapter 7
Shallow-water waves and tides
7.1
Tidal current profiles
Before oscillating velocity profiles are considered, the simple oscillation of forced
and damped vertically averaged tidal flow ū is briefly investigated:
∂ ū
∂η
F
= −g
− ū
∂t
∂x H
(7.1)
with the linear friction velocity F , which can be considered as the product of a
constant drag coefficient, cd and a constant reference velocity, uref . When the
forcing is harmonic, i.e.,
∂η
= η̄x cos(ωt),
(7.2)
∂x
then the solution is easily found as
−g η̄x
F
ū = F 2
cos(ωt) + ω sin(ωt) .
(7.3)
2
H
+
ω
2
H
For a frictionally dominated situation with F/H ≫ ω, surface slope and velocity
are in phase (both proportional to cos). For a frictionless situation with ω ≫
F/H, surface slope and velocity are out of phase.
In the following, velocity profiles in a tidal flow with linearised friction are
considered. Without the effect of Earth rotation and density differences and with
constant eddy viscosity the linearised dynamic equation for the eastward velocity
profile reads as
∂u
∂η
∂2u
= −g
+ νt 2 .
(7.4)
∂t
∂x
∂z
At the surface, the dynamic boundary condition (1.19) with zero surface momentum flux will be used, for the bottom the linearised bottom boundary condition
(2.11) is applied. The ansatz for solving (7.4) is to separate the vertical coordinate
63
and the time:
u(z, t) = R [û(z) exp(iωt)] ,
(7.5)
η(t) = R [η̂ exp(iωt)] ,
i.e., the real part of the solution is considered as the physically relevant solution.
Inserting (7.5) into (7.4) gives the following relation:
iω û = −g
∂ η̂
∂ 2 û
+ νt 2 .
∂x
∂z
(7.6)
A solution ansatz for (7.6) is
û = A exp(dz) + B exp(−dz) + C,
(7.7)
where we directly see:
d=
iω
νt
1/2
;
C=−
g ∂ η̂
.
iω ∂x
(7.8)
Inserting the zero-flux surface boundary condition into (7.7) gives
A = B,
(7.9)
and insertion of the linearised bottom boundary condition (2.11) gives
dAνt (exp(−dH) − exp(dH)) = F [A (exp(−dH) + exp(dH)) + C] ,
such that the complete solution reads as
(exp(dz) + exp(−dz))
û(z) = C
+1
E
(7.10)
(7.11)
with
E=
dνt
(exp(−dH) − exp(dH)) − (exp(−dH) + exp(dH)) .
F
C can be eliminated by defining the vertical mean velocity amplitude:
Z
1 0
û(z) dz
ûmax =
H −H
(7.12)
(7.13)
and inserting this into (7.7):
C=
ûmax dHE
.
exp(dH) − exp(−dH) + dHE
64
(7.14)
With this, the following solution is obtained:
exp(dz) + exp(−dz) + E
û(z)
= dH
ûmax
exp(dH) − exp(−dH) + dHE
(7.15)
Since û(z) is a complex velocity profile, the real part of its multiplication with
exp(iωt) will result in a linear combination of a sin and a cos function, i.e., a
profile of harmonic oscillations where the phase shift is a function of the vertical
position.
Figure 7.1 shows solutions for the velocity field u(z, t) from (7.5) with the
z-dependend velocity profile û(z) from (7.15) for two different water depths.
Clearly, for shallow water frictional effects dominate over the whole water column
whereas for deeper water frictional effects are constrained to the lower part of
the water column.
It can be clearly seen that the solution depends on two non-dimensional parameters,
1/2
ω
−1/2
J =i
dH = H
,
(7.16)
νt
which can be interpreted as the non-dimensional depth and
Y =
(ωνt )1/2
,
F
(7.17)
which reflects the effect of bottom friction. For details, see Prandle (2009). When
scaling the eddy viscosity with the water depth H and the maximum bottom
friction velocity u∗max , compare to (2.14),
νt ∝ u∗max H
(7.18)
and scaling the linear bottom friction coefficient F with the bottom friction velocity,
F ∝ u∗max
(7.19)
we see that
J2 ∝ Y 2 ∝
ωH
= St ,
u∗max
(7.20)
with the inverse Strouhal number, St . This can be interpreted in a way that
the fundamental external parameter of the tidal flow is here the inverse Strouhal
number. Another external parameter however is given by the way how exactly
the bottom friction is applied. In realistic scenarios (not using a constant eddy
viscosity) this parameter will be the non-dimensional bottom roughness, z0 /H.
It has been shown by Baumert and Radach (1992) that also for a spatially
and temporally variable eddy viscosity the simple tidal oscillation depends on one
65
Figure 7.1: Velocity profiles (upper panels) and velocity time series at z/H =
−0.9 (red) and z/H = 0 (green) for shallow water (H = 1 m, left) and relatively
deep water (H = 10 m, right). For these simulations the following other parameter have been used: ω = 2π/44714 s−1 , νt = 10−4 m2 s−1 and F = 10−4 m s−1 .
The solutions have been calculated as the velocity field u(z, t) from (7.5) with
the z-dependend velocity profile û(z) from (7.15).
non-dimensional parameter only. Assuming a spatially variable eddy viscosity
and a harmonic pressure gradient (surface slope) forcing, then (7.4)
∂u
∂u
∂
νt
,
(7.21)
= −gSmax cos(ωt) +
∂t
∂z
∂z
with the maximum surface slope Smax . When introducing non-dimensional parameters
z
u
νt
, t′ = tω, z ′ = , νt′ = √
,
(7.22)
u′ = √
H
gHSmax
H gHSmax
equation (7.21) can be cast into the following non-dimensional form:
′
∂u′
∂
′
′ ∂u
S̃t ′ = − cos(t ) + ′ νt ′ ,
∂t
∂z
∂z
66
(7.23)
with
S̃t = √
ωH
,
gHSmax
(7.24)
which can also be interpreted as an inverse Strouhal number. Since the only
parameter included in (7.23) is S̃t , it is clear that the dynamics does only depend
on this parameter. However, to calculate the eddy viscosity, typically a turbulence
closure model is used, which introduces the non-dimensional bed roughness as
a further parameter. In a generalisation of the work by Baumert and Radach
(1992), Burchard (2009) demonstrated how a number of other non-dimensional
parameters come in when effects of Earth rotation, wind and horizontal density
gradients are considered.
It should be noted that there are many definitions for the Strouhal number,
which originally has been defined by Strouhal (1878) as the ratio of the velocity
at which a wire is moved through the air and the product of its thickness and
the frequency of the generated sound. In oceanography, it has been defined as
clockwise and anti-clockwise Strouhal numbers for investigating rotational tidal
flow, as the ratio of the clockwise or anti-clockwise tidal velocity components
to the product of the water depth with the sum or difference between tidal and
inertial frequencies (multiplied by π, see Prandle (1982)).
Examples of flows for two different inverse Strouhal numbers are shown in
figure 7.2. The results have been obtained by applying the General Ocean Turbulence Model (GOTM, see www.gotm.net and Umlauf et al. (2005); Umlauf and Burchard
(2005)). It should be noted that the surface gradient forcing was proportional to
− sin(t′ ). It can be clearly seen that the velocity profile has more vertical shear
for the low inverse Strouhal number (relatively large friction velocity) and is also
significant less than a quarter tidal period out of phase with the surface gradient
forcing, which is in accordance with the simple result shown in (7.3) for frictional
tidal currents. In contrast, for the case of a large Strouhal number, the tidal
current is out of phase by a quarter tidal period compared to the surface gradient
forcing almost over the entire water column, which is in accordance with (7.3)
for frictionless tidal currents.
67
S̃t = 4.5 · 10−2
20
S̃t = 4.5 · 10−1
15
10
u′
u′
5
0
-5
-10
-15
z ′ = −0.9
z′ = 0
-20
0
0.2
0.4
0.6
0.8
1
t′ /(2π)
2.5
2
1.5
1
0.5
0
-0.5
-1
-1.5
′
-2 z =′ −0.9
z =0
-2.5
0
0.2
0.4
0.6
0.8
1
t′ /(2π)
Figure 7.2:
Non-dimensionalised tidal currents calculated with GOTM
(www.gotm.net) using a two-equation turbulence closure model and a harmonic
forcing. Shown are time series of currents of a periodic state at two depths, near
the bottom at z ′ = −0.9 (red) and near the surface at z ′ = 0 (green). Left: small
inverse Strouhal number. The flow is clearly asymmetric due to the effect of
friction. Right: large inverse Strouhal number. Here the flow is fairly symmetric
with the currents being fully out of phase with the surface gradient forcing, which
was proportional to sin(t′ ). For both simulations, the water depth was H = 10 m,
and the maximum surface slope was Smax = 10−5 . For the small inverse Strouhal
number, the tidal period was 44714 s, and for the other case it was 4471.4 s.
68
Chapter 8
Estuarine circulation
8.1
Observational evidence
Tidal estuaries are the lower reaches of rivers where the discharged freshwater
is mixed with oceanic water by tidal motions. A horizontal density gradient
due to the horizontal salinity gradient is an essential dynamic feature of tidal
estuaries. As one of the most remarkable phenomena, distinct estuarine turbidity
maxima (ETMs) are observed at the upper limit of the salt intrusion. These are
established by high concentrations of suspended particulate matter (SPM) which
is accumulated in this part of the estuary. This is shown for the Elbe Estuary
in figure 8.1. For two different river run-off situations, high SPM concentrations
are observed at salinities below 4 g kg−1 .
Observations of tidal mean (residual) profiles of velocity and SPM concentration explain the accumulation of SPM in the upper limit of the salt intrusion, see
figure 8.2: The lower part of the residual is directed upstream, i.e., against the
residual river run-off. The SPM, due to its tendency to sink down to the bottom
against the vertical mixing, is therefore subject to a net upstream transport in the
entire region of horizontal density gradients. At the upper limit of the salt intrusion where the horizontal density gradient is vanishing, this transport mechanism
is absent, such that SPM accumulated due to the divergence of horizontal net
transport. Mechanisms which drive the upstream near bottom residual velocity
will be investigated in section 8.3.
Another key process in tidal estuaries and coastal seas is tidal straining: During flood dense off-shore waters are sheared over less dense onshore waters, leading
to static instability and subsequent full vertical mixing. Thus, during flood there
is the tendency for the mixing the water column from the bottom (where the
shear is strongest) upwards to the surface. During ebb the opposite happens, less
dense onshore waters are sheared over denser offshore waters, leading to stable
stratification, suppression of vertical mixing and stratification. This process of
tidal straining has been first described by Simpson et al. (1990).
69
54°10'
WGE
Be 10.96
Trischen
Großer Vogelsand
Friedrichskoog
bsan
ee
-K
d
730
Neufelder
Sand
Neuwerk
Medemsand
720
No
rd
-O
sts
Scharhörn
54°00'
an
al
Gel
740
Brunsbüttel
Stör
700
710
690
Cuxhaven
53°50'
680
Medem
Medemgrund
Gr. Knechtsand
Glückstadt
isc
W
a
ck
e
lb
ere
üd
.S
hh
Oste
u
670
Krü
Pagensand
Pinnau
660
53°40'
Stade
650
640
inge
Watt
630 Hamburg
Lü
he
Schw
Es
te
Strom-km 700
0
10km
8°10'
20'
30'
40'
50'
9°00'
10'
20'
30'
40'
Mühlenberger Loch
50'
53°30'
10°00'
Figure 8.1: Observations of estuarine turbidity maxima in the Elbe estuary.
Upper panel: Map of the Elbe Estuary with river kilometres. Lower left: Synoptic
observations of surface salinity and suspended particulate matter concentrations
along the Elbe Estuary at a river runoff of 707 m3 s−1 at Neu Darchau (upstream
of Hamburg). Lower right: Same as lower left panel, but for a weak river runoff of
464 m3 s−1 . The two lower panels have been kindly provided by Jens Kappenberg
(Helmholtz-Zentrum Geesthacht).
70
Figure 8.2: Observations of tidal residual velocity, SPM concentration, horizontal
SPM flux and salinity in the Elbe Estuary (Kappenberg et al. (1995)).
Depending on the relative force of the horizontal density gradient (see section
8.3 for details), Simpson et al. (1990) could observe three different regimes: Permanently mixed waters (weak density gradient forcing), permanently stratified
waters (strong density gradient forcing) and an intermediate state which they
called Strain-Induced Periodic Stratification (SIPS).
In coastal waters such as Liverpool Bay, all three situation may occur subsequently at a certain position, depending on the external forcing conditions, see
Verspecht et al. (2009).
Detailed direct observations of SIPS have been made by Rippeth et al. (2001),
see figure 8.4.
71
Fig 2.a Temperature (Degrees C)
0
16
Depth (m)
1 5 .5
-20
15
1 4 .5
14
-40
1 3 .5
LB2
13
1 2 .5
-60
12
1 1 .5
11
100
90
80
70
60
50
40
30
20
10
0
Fig 2.b Salinity (PSU)
0
3 4 .4
Depth (m)
3 4 .2
34
-20
3 3 .8
3 3 .6
3 3 .4
-40
3 3 .2
33
3 2 .8
-60
3 2 .6
3 2 .4
3 2 .2
100
90
80
70
60
50
40
30
20
10
0
Fig 2.c Sigma T (kg/m3)
0
2 6 .6
2 6 .4
2 6 .2
Depth (m)
26
-20
2 5 .8
2 5 .6
2 5 .4
2 5 .2
-40
25
2 4 .8
2 4 .6
2 4 .4
2 4 .2
-60
24
2 3 .8
2 3 .6
100
90
80
70
60
50
40
30
20
10
0
Distance Along Transect (km)
Figure 8.3: Hydrographic transect across Liverpool Bay on July 7, 1999, showing
the horizontal gradients of temperature (upper panel), salinity (middle panel)
and density anomaly (lower panel), see Rippeth et al. (2001) and Simpson et al.
(2002) for details. The arrow shows where the water column observations by
Rippeth et al. (2001) have been carried out.
72
Liverpool Bay, observed temperature
Liverpool Bay, observed salinity
0
0
T /◦ C
16.0
-5
-10
15.5
-10
-15
15.0
14.5
-20
S/psu
33.6
33.4
z/m
z/m
-5
33.2
-15
33.0
32.8
-20
32.6
14.0
-25
-25
-30
-30
186.8
187.0
187.2
187.4
187.6
186.8
187.0
Julian Day 1999
187.2
187.4
187.6
Julian Day 1999
Liverpool Bay, observed u-component
Liverpool Bay, observed v-component
0
0
u/m s−1
0.8
0.6
0.4
0.2
0.0
-0.2
-0.4
-0.6
-0.8
z/m
-10
-15
-20
v/m s−1
-5
0.4
-10
z/m
-5
-25
0.2
0.0
-15
-0.2
-20
-0.4
-25
-30
-30
186.8
187.0
187.2
187.4
187.6
186.8
187.0
Julian Day 1999
187.2
187.4
187.6
Julian Day 1999
Liverpool Bay, observed dissipation rate
0
log10 (ε/m2 s−3 )
-5
-3
-4
z/m
-10
-5
-15
-6
-7
-20
-8
-25
-30
186.8
187.0
187.2
187.4
187.6
Julian Day 1999
Figure 8.4: Strain-Induced Periodic Stratification (SIPS) in Liverpool Bay: Temperature T̄ , salinity S̄ and velocity components ū and v̄ and turbulent dissipation
rate from observations by Rippeth et al. (2001) and simulations without nudging
to observed T̄ and S̄ for two tidal cycles. See figure 8.3 for the location of the
observed water column within the horizontal density gradient.
73
Figure 8.5: Sketch of an estuary to explain the Knudsen relation. Graphics taken
from MacCready and Geyer (2010).
8.2
The Knudsen relation
The intensity of an estuarine exchange flow has already been described by Knudsen
(1900) for the example of the Darss Sill which is the main sill into the Baltic Sea.
With QR being the river volume run-off (plus net precipitation) into the estuary,
Q1 being the outflow at the mouth of the estuary and Q2 being the inflow at the
mouth of the estuary (see figure 8.5), then the steady state volume balance is
expressed as
Q1 = Q2 + QR .
(8.1)
With s1 denoting the outflow salinity and s2 denoting the inflow salinity, the salt
budget of the estuary can be calculated as (note that the salinity of the river is
sR = 0)
V ṡa = Q2 s2 − Q1 s1 ,
(8.2)
where st is the volume averaged salinity. Assuming a stationary salinity in the
estuary (which a good approximation on long time scales), i.e., ṡa = 0, the inflow
and outflows can be estmated on the basis of their salinities (which are much
easier to measure):
Q1 =
s2
QR ;
s2 − s1
Q2 =
s1
QR ,
s2 − s1
(8.3)
relations which can be used to get a simple relation for the ratio of the inflow to
the outflow:
s1
Q2
=
(8.4)
Q1
s2
In his historical paper, Knudsen (1900) gives for the Darss Sill values of s1 = 8.7
and s2 = 17.4, such that Q2 /Q1 =0.5, i.e., the inflow is half of the outflow (as well
as half of the total river run-off + net precipitation).
74
8.3
Idealised dynamics
To reproduce the basic dynamics leading to estuarine circulation, Earth rotation,
horizontal velocity gradients and atmospheric pressure gradients are neglected in
the dynamic equations (1.13), leading to the following one-dimensional momentum equation:
Z
∂u
∂η
∂
g 0 ∂ρ
∂u
(νt + ν)
= −g
−
−
dz.
(8.5)
∂t
∂z
∂z
∂x ρ0 z ∂x
If, additionally, a temporally and spatially constant horizontal density gradient
is assumed, then the equations simplify to the following form:
∂u
∂η
∂
∂u
(νt + ν)
= −g
−
− zbx
(8.6)
∂t
∂z
∂z
∂x
with
bx =
∂b
g ∂ρ
=−
,
∂x
ρ0 ∂x
(8.7)
see (1.22). If furthermore the eddy viscosity is temporally constant and the molecular viscosity is neglected, and additionally the temporal average of a periodically
oscillating solution ū is considered, then the dynamic equation is of the following
form:
∂ ū
∂
νt
= −g η̄x − zbx ,
(8.8)
−
∂z
∂z
with the temporally averaged surface elevation gradient, η̄x . By assuming periods of Tp = 44714 s, the semi-diurnal principal lunar tide is considered in the
following.
Integration from an arbitrary position z in the water column to the surface
(where we assume zero stress, i.e., no wind stress), results in
∂ ū
z 2 bx
+ zE,
=
∂z
2 νt
(8.9)
where for simplicity E = g η̄x /νt has been defined. Subsequent integration from
the bottom at z = −H (using ū(−H) = 0) to an arbitrary position in the water
column leads to
3
2
z
z
H 3 bx
H2
ū(z) =
+
E.
(8.10)
+
−
6
6
νt
2
2
To eliminate E, we define the residual run-off ur as the vertically and tidally
averaged velocity:
Z
1 0
u(z) dz.
(8.11)
ur =
H −H
75
When taking the vertical average of (8.10), and using the result for eliminating
E in (8.10), the following classical solution for estuarine circulation is obtained
(Hansen and Rattray (1965); MacCready and Geyer (2010)):
z 2
3 z 2
bx H 3 z 3
8
+9
−1 −
− 1 ur
ū(z) =
48νt
H
H
2 H
z 2
z 3
3 z 2
= 8
+9
− 1 ue −
− 1 ur ,
H
H
2 H
(8.12)
with the exchange flow intensity, ue = (bx H 3 )/(48νt ). The deviation of the
resulting vertical salinity anomaly
Z
1 0
′
s =s−
dz
(8.13)
H −H
(deviation from the vertical mean salinity) can be obtained from the following
balance between salinity advection and vertical mixing:
ūsx = νt′
∂ 2 s′
,
∂z 2
(8.14)
see assignment 1. The salt deviation profile will come out as a fifth-order polynomial (due to the double integration of the third-order polynomial for ū) with
zero-gradients at the surface and the bottom. In (8.14), the constant in time and
space horizontal salinity gradient is defined as
sx =
bx
gβ
(8.15)
with the haline contraction coefficient β from (1.21).
The suspended particulate matter equation (1.26) simplifies for stationary
and horizontally homogeneous conditions to
∂C
∂
′ ∂C
− ws
νt
= 0.
(8.16)
−
∂z
∂z
∂z
Vertical integration from the bottom to an arbitrary position z in the water
column leads to
′ ∂C ′ ∂C − ws C(z) − νt
= −ws C(−H) − νt
,
(8.17)
∂z z
∂z z=−H
where the right hand side is zero due to boundary condition that the net SPM
flux through the bottom vanishes. With this, the vertical net SPM flux is zero for
76
all z, such for a spatially and temporally constant eddy diffusivity the following
condition is obtained:
ws
∂C
= − ′ C,
(8.18)
∂z
νt
the solution of which is
ws
C = C0 exp − ′ z .
νt
With Cb = C(−H), this can be reformulated to
C
ws H z
= exp − ′
+1 .
Cb
νt
H
(8.19)
(8.20)
Figure 8.6 shows the residual velocity, salinity anomaly and particulate suspended matter profiles for constant in time and space eddy viscosity and diffusivity, resulting from equations (8.12) and (8.19). The resulting SPM flux profile for
q = ūC/Cb as well as the analytical solution of (8.14). The similarity of these four
profiles with the Elbe Estuary measurements shown in figure 8.2 are remarkable.
With this, the classical estuarine circulation profile can be explained by the
presence of the horizontal density gradient, such that the circulation is density
driven. Therefore this type of circulation is also called gravitational circulation.
This gravitational circulation has been identified by Postma and Kalle (1955) as
the major driver of estuarine turbidity maxima, as shown in figure 8.1.
The dynamics of the up-estuary SPM flux as theoretically derived here is
governed by two dimensionless parameters. As shown in (8.12) the exchange flow
velocity can be expressed as
bx H 3
,
(8.21)
ue =
48νt
such that with (7.18) the dimensionless exchange flow can be formulated as
bx H 2
ue
∝
=: Si,
u∗max
(u∗max )2
(8.22)
with the Simpson number Si, which quantifies the non-dimensional horizontal
density gradient. It can also be interpreted as the ratio of the stratifying forcing
due to the density gradient forcing and the de-stratifying forces due to vertical
mixing. The importance of Si for estuarine and coastal dynamics has first been
acknowledged by Simpson et al. (1990). Later this number has been phrased as
horizontal Richardson number Rx (Monismith et al. (1996)) and just recently it
was renamed to Si to honour John Simpson for his substantial scientific contributions to the Physical Oceanography of coastal, estuarine and shelf sea waters.
The other relevant non-dimensional parameter governs the SPM profile:
ws H
ws
≡ ∗ := R,
′
νt
umax
77
(8.23)
with the Rouse number R, which characterises the balance between setting and
vertical mixing of SPM.
However, gravitational circulation is not the only process leading to nearbottom upstream residual velocities.
The role of tidal straining for the generation of estuarine circulation has been
first proposed by Jay and Musiak (1994), see figure 8.7: The flood profile is denoted as Uf and shows in comparison to a symmetric flood profile +A a negative
velocity defect in the upper and a positive velocity defect in the lower part of the
water column, caused by a relatively large flood eddy viscosity Kzf which mixes
large amounts of momentum down to the bottom. This is because during flood
denser, saltier offshore water is sheared over less dense onshore waters, leading
to strong enhancement of vertical mixing due to statically unstable stratification, see also the enhanced turbulent dissipation rate during flood as shown in
the observations presented in figure 8.4. During ebb, the opposite mechanism is
active, leading to suppressed downward transport of down-estuary momentum.
Averaged over a complete tidal cycle, an excess of up-estuary transport is mixed
downwards, leading to up-estuary near-bottom residual circulation, in addition
to the classical gravitational circulation shown in figure 8.6. This tidal asymmetry had been described before as tidal straining by Simpson et al. (1990), who
concentrated on the strain-induced periodic stratification (SIPS, see figure 8.4)
caused by this mixing asymmetry. Using a one-dimensional numerical model
and observations, Stacey et al. (2008) could show the importance of tidal mixing
asymmetry for the generation of estuarine circulation. It was recently shown by
means of a one-dimensional numerical model study only considering longitudinal flow processes that for a SIPS situation tidal straining accounts for about
2/3 of the residual circulation (Burchard and Hetland (2010)). This confirmed
an earlier result by Burchard and Baumert (1998) who used a two-dimensional
vertical-longitudinal model for a tidal channel to prove the dominant role of tidal
straining for the generation of estuarine circulation.
78
Residual SPM profile
z/H
Residual velocity profile
0
0
-0.2
-0.2
-0.4
-0.4
-0.6
-0.6
-0.8
-0.8
-1
-0.1 -0.08 -0.06 -0.04 -0.02
-1
0
0.02 0.04
0
0.2
0.4
0
-0.2
-0.2
-0.4
-0.4
-0.6
-0.6
-0.8
-0.8
0
0.005
0.8
1
Residual salinity anomaly profile
Residual SPM flux profile
0
-1
-0.01 -0.005
0.6
C/Cb
ū / m s−1
0.01
0.015
-1
-0.3
0.02
q / m s−1
-0.2
-0.1
0
0.1
0.2
0.3
s′
Figure 8.6: Residual velocity, ū (upper left), SPM normalised by bottom SPM,
C/Cb (upper right), horizontal SPM flux normalised by bottom SPM, q = ūC/Cb
(lower left), and salinity anomaly, s′ (lower right) for constant in time and space
eddy viscosity and diffusivity. The following parameters have been used here:
H = 20 m, νt = 3·10−3 m2 s−1 , νt′ = 1.5·10−3 m2 s−1 , bx = 10−6 s−2 , sx = −1.49·10−4
m−1 , ur = −0.02 m s−1 and ws = 2 · 10−4 m s−1 .
79
Figure 8.7: Sketch explaining the process of tidal straining. This graphic
has been taken from MacCready and Geyer (2010) who have redrawn it from
Jay and Musiak (1994).
80
Bibliography
Baines, P. G., 1995: Topographic Effects in Stratified Flows. Cambridge Monographs on Mechanics, Cambridge University Press, Cambridge, UK, ISBN:
0–551–43501–3.
Baumert, H. and G. Radach, 1992: Hysteresis of turbulent kinetic energy in
nonrotational tidal flows. J. Geophys. Res., 97 (C?), 3669–3677.
Burchard, H., 2009: Combined effects of wind, tide and horizontal density gradients on stratification in estuaries and coastal seas. J. Phys. Oceanogr., 39,
2117–2136.
Burchard, H. and H. Baumert, 1998: The formation of estuarine turbidity maxima due to density effects in the salt wedge. A hydrodynamic process study. J.
Phys. Oceanogr., 28, 309–321.
Burchard, H. and R. D. Hetland, 2010: Quantifying the contributions of tidal
straining and gravitational circulation to residual circulation in periodically
stratified tidal estuaries. J. Phys. Oceanogr., 40, 1243–1262.
Hansen, D. V. and M. Rattray, 1965: Gravitational circulation in straits and
estuaries. J. Mar. Res., 23, 104–122.
Jay, D. A. and J. D. Musiak, 1994: Particle trapping in estuarine tidal flows. J.
Geophys. Res., 99 (C10), 20 445–20 461.
Kappenberg, J., G. Schymura, and H.-U. Fanger, 1995: Sediment dynamics and
estuarine circulation in the turbidity maximum of the Elbe River. Netherlands
Journal of Equatic Ecology, 29, 229–237.
Kato, H. and O. M. Phillips, 1969: On the penetration of a turbulent layer into
stratified fluid. J. Fluid Mech., 37 (4), 643–655.
Knudsen, M., 1900: Ein hydrographischer Lehrsatz. Ann. Hydrogr. Mar. Meterol., 28, 316–320.
Kraus, E. B. and J. S. Turner, 1967: A one-dimensional model of the seasonal
thermocline - II. The general theory and its consequences. Tellus, 19, 98–105.
81
Longuet-Higgins, M. S. and R. W. Stewart, 1964: Radiation stresses in water
waves; a physical discussion, with applications. Deep-Sea Res., 11 (4), 529–
562.
MacCready, P. and W. R. Geyer, 2010: Advances in estuarine physics. Annu.
Rev. Marine. Sci., 2, 35–58.
Monismith, S. G., J. R. Burau, and M. Stacey, 1996: Stratification dynamics and
gravitational circulation in Northern San Francisco Bay. San Francisco Bay:
the ecosystem, J. T. Hollibaugh, Ed., American Association for the Advancement of Science, Pacific Division, San Francisco, 123–153.
Niiler, P. P., 1975: Deepening of the wind mixed layer. J. Mar. Res., 33, 405–421.
Pollard, R. T., P. B. Rhines, and R. O. R. Y. Thompson, 1973: The deepening
of the wind mixed layer. Geophys. Fluid Dyn., 3, 381–404.
Postma, H. and K. Kalle, 1955: Die Entstehung von Trübungszonen im Unterlauf der Flüsse, speziell im Hinblick auf die Verhältnisse in der Unterelbe. Dt.
Hydrogr. Z., 8, 138–144.
Prandle, D., 1982: The vertical structure of tidal currents. Geophysical and Astrophysical Fluid Dynamics, 22, 29–49.
Prandle, D., 2009: Estuaries: Dynamics, mixing, sedimentation amd morphology.
Cambridge University Press, 236 pp. pp.
Price, J. F., 1979: On the scaling of stress driven entrainment experiments. J.
Fluid Mech., 90 (3), 509–529.
Rippeth, T. P., N. R. Fisher, and J. H. Simpson, 2001: The cycle of turbulent
dissipation in the presence of tidal straining. J. Phys. Oceanogr., 31 (8), 2458–
2471.
Simons, T. J., 1980: Circulation models of lakes and inland seas. Canadian Bulletin of Fisheries and Aquatic Sciences, 203.
Simpson, J. H., J. Brown, J. Matthews, and G. Allen, 1990: Tidal straining,
density currents, and stirring in the control of estuarine stratification. Estuaries
and Coasts, 13, 125–132.
Simpson, J. H., H. Burchard, N. R. Fisher, and T. P. Rippeth, 2002: The semidiurnal cycle of dissipation in a ROFI: model-measurement comparisons. Cont.
Shelf Res., 22, 1615–1628.
Stacey, M. T., J. P. Fram, and F. K. Chow, 2008: Role of tidally periodic density
stratification in the creation of estuarine subtidal circulation. J. Geophys. Res.,
113, doi: 10.1029/2007JC004 581.
82
Strouhal, V., 1878: Über eine besondere Art der Tonerregung. Annalen der
Physik, 241, 216–251.
Turner, J. S., 1986: Turbulent entrainment: the development of the entrainment
assumption, and its applications to geophysical flows. J. Fluid Mech., 173,
431–471.
Turner, J. S. and E. B. Kraus, 1967: A one-dimensional model of the seasonal
thermocline – I. A laboratory experiment and its interpretation. Tellus, 19,
88–97.
Umlauf, L. and H. Burchard, 2005: Second-order turbulence closure models for
geophysical boundary layers. A review of recent work. Cont. Shelf Res., 25,
795–827.
Umlauf, L., H. Burchard, and K. Bolding, 2005: GOTM – Scientific Documentation. Version 3.2. Marine Science Reports 63, Leibniz-Institute for Baltic Sea
Research, Warnemünde, Germany.
Verspecht, F. I., H. Burchard, T. P. Rippeth, M. J. Howarth, and J. H. Simpson, 2009: Processes impacting on stratification in a region of freshwater influence: Application to Liverpool Bay. J. Geophys. Res., 114, 11 022, doi:
10.1029/2009JC005475.
83

The First Integral Method for (3+1) dimensional modified KdV

TOURNAMENT INFORMATION Lil Big South Date: January 17, 18

Document

Electromagnetism: Example Sheet 3

Optimization of an Oscillating Wave Surge Converter

Florian Kogelbauer, MSc BSc Some Topics in Mathematical Water

Plenary 4 - Laura Carse

流体力学に現れる偏微分方程式系の数学解析 Mathematical Analysis

our Exam1 - Bama.ua.edu