CRWR Online Report 0409
Hydrostatic and Nonhydrostatic Internal Wave Models
by
Bridget M. Wadzuk, B.S., M.S
Graduate Research Assistant
and
Ben R. Hodges, Ph.D.
Principal Investigator
The University of Texas at Austin
December 4, 2004
Final Report to the Office of Naval Research
under Contract No. N000140110574
This document is available online via World Wide Web at
http://www.crwr.utexas.edu/online.html
Copyright
by
Bridget Marie Wadzuk
2004
Abstract
Numerical models have become an indispensable tool for ocean and inland flow
modeling. Such models typically use the hydrostatic approximation based on the
argument that their horizontal length scales are longer than the vertical length
scales. There are a wide variety of physical processes in oceans and inland water
systems, and many of these processes are adequately modeled with the hydrostatic
approximation. However, internal waves contribute to the physics that influence
mixing in a density stratified system and have been previously shown to be non
hydrostatic. The neglect of nonhydrostatic pressure in a hydrostatic model is
problematic since nonhydrostatic pressure plays a significant role in internal wave
evolution balancing nonlinear wave steepening. Where nonhydrostatic pressure is
neglected in a model, the governing equations are missing a piece of the physics
that control the internal wave evolution, so it should not be surprising that the
evolution may be poorly predicted.
Despite the knowledge that the nonhydrostatic pressure is necessary for
correctly modeling the physics of a steepening internal wave, the high
computational cost of solving the nonhydrostatic pressure has limited its use in
largescale systems. Furthermore, the errors associated with hydrostatic modeling
of internal waves have not been quantified. This research quantifies the differences
between hydrostatic and nonhydrostatic simulations of internal wave evolution
and develops a method to a priori determine regions with nonhydrostatic behavior.
In quantifying the errors and differences between the two models this research
provides the characteristics of model error with grid refinement. Additionally, it is
shown that hydrostatic models may develop high wavenumber ?solitonlike?
features that are purely a construct of model error, but may seem to mimic physical
behaviors of the nonhydrostatic system. Finally, it is shown that regions of
significant nonhydrostatic pressure gradients can be identified from a hydrostatic
model. This latter finding is a building block towards coupling local non
hydrostatic solutions with global hydrostatic solutions for more efficient
computational methods. The work presented here provides the foundations for
future nonhydrostatic model development and application.
i
Chapter 1. Introduction
Abstract ............................................................ i
Chapter 1. Introduction............................1
1.1 Motivation................................................1
1.2 Objectives ................................................3
1.3 Background..............................................3
1.3.1 Internal Waves.................................. 3
1.3.2 Modeling........................................... 5
1.4 Approach..................................................6
Chapter 2. Methods...................................7
2.1 Numerical Model .....................................7
2.1.1 Governing Equations........................ 7
2.1.2 Solution Methods .............................. 8
2.1.3 Iterative Methods ? Pressure Poisson
Equation 10
2.1.4 Boundary Conditions...................... 10
2.2 Analysis Methods ..................................10
2.2.1 Background..................................... 11
2.2.2 Numerical Diffusion ....................... 13
2.2.3 Numerical Dissipation.................... 16
2.2.4 Spectral Analysis ............................ 17
Chapter 3. Verification and validation of
the nonhydrostatic solver ............................19
3.1 Verification by the Method of
Manufactured Solutions...............................19
3.1.1 Setup ............................................... 19
3.1.2 Error Analysis................................. 19
3.2 Validation and Convergence for an
Unsteady Internal Wave...............................22
3.2.1 Setup ............................................... 22
3.2.2 Results............................................. 23
3.3 Validation Against Regime Theory and
Laboratory Experiment................................25
3.3.1 Regime Theory................................ 25
3.3.2 Horn?s Laboratory Experiment and
Model Simulations.......................................... 26
3.3.3 Results............................................. 28
3.4 Summary................................................32
Chapter 4. Comparison of Hydrostatic
and Nonhydrostatic models ..........................33
4.1 Laboratory Scale Comparison ...............33
4.1.1 Setup ............................................... 33
4.1.2 Scenario Simulations ...................... 34
4.1.3 Changing the Grid Resolution ........ 39
4.2 LakeScale Comparison.........................47
4.2.1 Setup ............................................... 47
4.2.2 Discussion....................................... 52
4.3 Summary................................................53
Chapter 5. Isolation of nonhydrostatic
effects .............. ...............................................59
5.1 Numerical Theory..................................59
5.2 Application of the Nonhydrostatic
Pressure Isolation Method...........................62
5.2.1 Setup................................................62
5.2.2 Results.............................................62
5.3 Summary ...............................................64
Chapter 6. Conclusions and
Recommendations .........................................71
6.1 Summary Discussion.............................71
6.2 Conclusions ...........................................72
6.3 Recommendations .................................73
Acknowledgements .......................................73
Bibliography ..................................................75
iii
iv
Chapter 1. Introduction
Chapter 1. Introduction
1.1 Motivation
Biogeochemical processes in oceans, lakes and
estuaries evolve against the background of vertical
density stratification. Internal waves are one
mechanism which carry energy and momentum
through a basin and contribute to mixing events in
littoral regions that prompt ecosystem changes.
Internal waves are motions that occur beneath the
freesurface of a densitystratified waterbody. As an
external force (i.e. wind or river inflow) moves the
density layers from their equilibrium position, an
internal wave is initiated to restore the system to
equilibrium (Figure 1.1). Internal waves propagate
through the basin and interact with the basin
boundaries; this internal wave ? slope interaction is
an important source of energy that transports
nutrients, biota and contaminants through the water
column (Imberger and Ivey, 1993; Javam, et al.,
1999). Hydrodynamic processes (e.g. internal
waves, tides and eddies) are transport mechanisms in
coastal oceans and lakes and their scope and
magnitude depends on the processes? speed and
length scales. Table 1.1 shows a large range of
scales for different hydrodynamic processes and
internal waves are considered a middlescale process.
Quantitative assessment of nutrient and
constituent transport in a stratified basin requires
modeling of internal wave evolution (Imberger,
1994). To have an unambiguous water quality
model, the complementary hydrodynamic model,
which provides the flow field, should simulate all the
significant physics (Gross, et al., 1999).
Hydrodynamic models often neglect physics that are
irrelevant to the focus of an investigation or can be a
priori scaled as small compared to dominant
(Marshall, et al., 1997). Largescale ocean and
littoral modeling have traditionally made the
hydrostatic approximation, which neglects non
hydrostatic pressure and subsequently vertical
momentum (e.g. POM: Blumberg and Mellor, 1987;
EFDC: Hamrick, 1992; ROMS: Haidvogel, et al.,
2000; Ezer, et al., 2002). The hydrostatic
approximation fails at open boundaries (Mahadevan,
et al., 1996a) and at steep slopes with strong vertical
velocities (Weilbeer and Jankowski, 2000; Horn, et
al., 2001; Horn, et al., 2002). Linear waves that are
damped by viscous effects may be considered to
b)a)
c) d)
Figure 1.1: Internal wave evolution in a twolayer system. a) System at rest, where the equilibrium
position for the freesurface (thick dashed blue line) and pycnocline (thick dashed red line) is flat. b)
Wind is applied over the free surface (large red arrow) which moves the freesurface (thin blue line) and
pycnocline (thin red line) from their respective equilibrium positions. c) Continued application of wind
over the freesurface, so the equilibrium positions for the freesurface and pycnocline are tilted. The
freesurface and pycnocline oscillate around the equilibrium position. d) The wind force stops and the
freesurface and pycnocline equilibrium position is flat. The freesurface and pycnocline oscillate around
the equilibrium position, setting up a basinscale seiche.
1
Wadzuk & Hodges: Hydrostatic and Nonhydrostatic Internal Wave Models
2
Process/Activity Speed (m/s) Length
Oceanic Turbulence 0.1  100 1 cm  100 m
Internal Waves 1  3 100 m  1 km
Diurnal Tides 0.1  10 1 km  1000 km
Mesoscale Eddy 1 10 km  1000 km
Rossby Waves 0.1  9 1000 km  10000 km
Surface Gravity Waves 1  20 1 m  100 m
Sound waves 1400 200 km
behave hydrostatically, while internal waves that
steepen nonlinearly are inherently nonhydrostatic
(Long, 1972). It is suggested that nonhydrostatic
pressure is essential for proper modeling of internal
wave development and propagation (Laval, et al.,
2003). When numerical models neglect non
hydrostatic pressure, they are neglecting a physical
attribute of internal wave evolution. This neglect has
repercussions on vertical transport at the topographic
boundaries of a basin. Internal waves that
nonlinearly steepen and degenerate into solitons can
propagate to a basin?s boundary where the wave may
shoal and mix, energizing a turbulent benthic
boundary layer (Boegman, et al., 2003); the
interaction of internal waves and boundary layer
mixing is a significant mechanism for nutrient
transport and other biogeochemical processes (De
Silva, et al., 1997; Nishri, et al., 2000). Modeling the
boundary layer is not of interest to the present
research, but properly modeling the internal wave
evolution which contributes to boundary layer is.
During the past decade there has been interest
simulating internal waves in numerical models,
marked by the development of several non
hydrostatic models (e.g. Mahadevan, et al., 1996a,b;
Marshall, et al., 1997; Stansby and Zhou, 1998;
Casulli, 1999; Weilbeer and Jankowski, 2000;
Fringer and Street, 2003; Daily and Imberger, 2003;
Yuan and Wu, 2004). However, there has been little
work done on quantifying the differences between
hydrostatic and nonhydrostatic models or the scales
for which nonhydrostatic pressure impacts internal
wave evolution. If the proper scales (i.e. sufficient
grid resolutions) are not used in a nonhydrostatic
model, then the results of the nonhydrostatic model
are no different than a hydrostatic model?s solution,
as described in ?4.2
Presently, computational power is insufficient
for application of nonhydrostatic models to large
domains (i.e. coastal oceans) with sufficient grid
resolution. By example, one computational ?work
unit? may be considered the processing power
required for one operation across the entire domain.
Thus, a work unit scales on the number of cells
within a domain (e.g. a grid of N ? M size: 1 work
unit = N ? M). A hydrostatic model?s basic
operations are solving for the three velocity
components, salinity transport, temperature transport
and the freesurface; the computational effort for
each of these operations scales on one work unit,
thus totaling six work units for each timestep. A
nonhydrostatic model solves for these operations
twice (once before the nonhydrostatic solver and
once after) and a pressure solution. The work
associated with the pressure solution is N ? M ? the
number of smoothing iterations per timestep; the
number of iterations per timestep for the non
hydrostatic pressure solver may be O(10
2
) (?3.2 ff.).
If the nonhydrostatic model is applied with the same
spatial grid and timestep as the hydrostatic model,
the nonhydrostatic model requires an additional 10
6
work units, that is, about two orders of magnitude
more computational effort. However, as discussed in
?4, the nonhydrostatic model also requires a finer
spatial grid and timestep than a hydrostatic model. A
typical hydrostatic coastal ocean model applies a
spatial grid with an aspect ratio (?z/?x) of O(10
3
)
(Marshall, et al., 1997; Tartinville, et al., 1998),
while ?4 shows that a nonhydrostatic model requires
a grid with an aspect ratio of at least O(10
2
). This
refined spatial grid increases the required work units
by another order of magnitude, to ~10
3
.
Furthermore, ?4 discusses the need for a non
hydrostatic model to be applied to a smaller timestep,
which increases the number of work units by another
order of magnitude. Thus, O(10
4
) work units are
needed for the nonhydrostatic model, while the
hydrostatic computational effort scales on O(10)
work units. According to Moore?s law (Moore,
1965), computational ability doubles about every two
years. Assuming that we are reaching the maximum
capacity of present computers with the hydrostatic
model, the computational power needed to run a
coastal ocean nonhydrostatic model with a
sufficiently fine spatial grid and timestep will not be
available for another 26 years! It is computational
power that limits our ability to use nonhydrostatic
models. This motivates the present research to
quantify the spatial grid and timestep necessary to
apply a nonhydrostatic model and develop the
theoretical foundations for future methods that may
circumvent the constraint of computational power.
Table 1.1: Scales of oceanic processes (Gill, 1982;
CushmanRoisin, 1994; Kantha & Clayson, 2000a;
Miropol?sky, 2001)
Chapter 1. Introduction
1.2 Objectives
The present research has two main objectives to
address the motivational issues within internal wave
modeling:
? quantify the differences between hydrostatic
and nonhydrostatic simulations of internal
wave evolution, and
? develop a method to a priori determine
regions with nonhydrostatic behavior.
The first objective is achieved by: 1)
quantifying accumulation of error within both
hydrostatic and nonhydrostatic models of internal
waves; and 2) comparing the spacetime evolution of
internal wave characteristics in both models.
Comparing the models? internal wave evolution and
error accumulation provides a means of defining the
conditions for which a hydrostatic model is
acceptable and those conditions requiring a non
hydrostatic model to capture an internal wave?s
propagation. It is shown that a model?s error
accumulation and representation of internal wave
evolution are dependent on spatial and temporal
resolution. Delineation of the spatial and temporal
resolution necessary to adequately resolve non
hydrostatic processes is included within the first
objective.
The second objective is achieved by using a
hydrostatic model to estimate nonhydrostatic
pressure effects. Using this estimation, local regions
where the nonhydrostatic pressure significantly
contributes to internal wave evolution are identified.
Through reaching this objective, the present work
provides the theoretical basis for future model
development for a nonhydrostatic solution in
isolated, ?nonhydrostatic? areas, concurrent with a
hydrostatic solution elsewhere.
3
1.3 Background
1.3.1 Internal Waves
Internal waves occur in stably stratified fluids when
a water parcel is displaced by some external force
(wind, inflow, etc.) and is restored by buoyancy
forces; the restoration motion may overshoot the
equilibrium position and set up an oscillation thereby
forming an internal wave (Turner, 1973; Kantha and
Clayson, 2000a). The buoyancy frequency (N) is the
upper limit of wave frequencies (?) that can
propagate through a system; thus, the only internal
waves of interest for the present research are those
where N?< , as waves where remain local
and do not significantly contribute to the system?s
overall transport and energy (Turner, 1973; Lighthill,
1978; Javam, et al., 1999). The lower limit of
internal wave frequencies is the Coriolis or inertial
wave frequency, f. The inertial frequency is defined
as f = 2?sin?, where ? is the angular velocity of the
earth?s rotation and ? is the latitude. At f, wave
motions become inertial oscillations, where fluid
parcels have horizontally circular trajectories
(CushmanRoisin, 1994). While Coriolis forces are
known to affect internal wave evolution, their impact
is principally threedimensional. To limit the
number of processes under consideration, this project
is restricted to vertical stratification in two
dimensional box models where Coriolis effects
cannot be represented. Thus, this project is focused
primarily on internal waves dominated by gravity,
which is typically the case for waves breaking on
sloping boundaries (De Silva, et al., 1997).
N?>
Field studies by Garrett and Munk (1979)
indicated that most internal waves follow a
consistent energywave spectrum (Figure 1.2).
There is a full range of wave frequencies and
associated energy, where waves of lowfrequency are
associated with high energy, while those of high
frequency have low energy. Internal waves also
evolve over a wide spectrum of spatial scales. An
initial lowfrequency, long wavelength internal wave
may degenerate into a train of highfrequency, short
Figure 1.2: EnergyFrequency Spectrum. f is the
Coriolis frequency, N is the BruntV?is?l?
frequency and ? is the frequency, cph (cycles per
hour) is a frequency scale and m
2
/cph is an energy
scale. (Garrett and Munk, 1979)
Wadzuk & Hodges: Hydrostatic and Nonhydrostatic Internal Wave Models
wavelength waves (Kao, et al., 1985; Horn, et al.,
2001; Boegman, et al., 2003), as seen in Figure 1.3.
This degeneration is connected to the nonlinearity of
internal waves. Several field studies (Farmer, 1978;
Osborne and Burch, 1980; Apel, 1981; Wiegand and
Carmack, 1986) and model results (Segur and
Hammack, 1982; Hutter, et al., 1998; Horn, et al.,
2002) have verified that an internal wave will behave
nonlinearly if the initial wave is sufficiently steep.
Steepness is defined as the ratio of a wave?s
amplitude to its wavelength. A ?sufficiently? steep
wave is one whose timescale for steepening is
shorter than its timescale for damping (Farmer, 1978;
Horn, et al., 2001). The relationship between the
total water depth and thickness of the upper layer can
also contribute to a wave?s nonlinearity.
The following description of wave evolution is
adapted from Kinsman (1965), Whitham (1974) and
Lighthill (1978). A wave may be viewed as a
recognizable feature of a disturbance that moves with
a finite velocity. An infinite number of sinusoidal
waves, each with a discrete frequency and discrete
wavelength, may be superposed to create the
characteristic material surface recognizable as a
wave (Figure 1.4 A). Wave speed (c):
c
k
?
= (1.1)
where ? is the wave frequency and k is the total
wavenumber (i.e.
22 2
kjlm=++; j, l and m are the
wavenumbers in the x, y and zdirections,
respectively) is dependent on the dispersion
relationship for an internal wave:
( )
()
22
22
22 2
jl
N
jlm
+
?=
++
(1.2)
where the buoyancy frequency, N
2
is:
2
o
gd
N
dz
?
=
?
(1.3)
It follows that wave speed is dependent upon
the density gradient of fluid. Waves on different
density structures propagate at different speeds (i.e.
the crest of the wave, with one density gradient, is
moving faster than the trough of the wave, with
another density gradient, as seen in Figure 1.4 B&C).
This density gradient dependency produces a
nonlinearly steepening wave. Wave speed is also
dependent on wavelength (the inverse of the
wavenumber). As the different sinusoidal waves
comprising the material wave move at different
speeds, the material wave spreads out over time, or
disperses, into a series of waves (Figure 1.4 D).
The wave speed dependence on the dispersion
relationship shows there is a relationship between
Figure 1.3: Laboratory experiment, conducted by Horn, et al. (2001), of a steepening, initial
basinscale wave that degenerates into a train of solitons. The first panel is the initial wave setup.
The second panel shows the wave front steepening. The third panel shows the beginning of
smaller waves developing behind the lead wave. The last two panels are the train of solitons
traveling through the domain. The basin is 6m long and 0.29m high.
4
Chapter 1. Introduction
nonlinear acceleration and a dispersive force. The
dispersive force that opposes nonlinear wave
deformation is the nonhydrostatic pressure gradient.
The dispersion of a wave leads to an energy cascade
from an initial, lowfrequency wave to a resulting,
high frequency wave series that is more susceptible
to wave breaking (van Haren, 2004). When non
hydrostatic pressure is not modeled while
nonlinearities are modeled (i.e. typical of a
hydrostatic model), there are two possible results.
The first possibility is a wave steepens unabated and
creates artificial mixing events through density
overturns. The second possibility is a build up of
numerical error that balances nonlinear steepening.
A
The two dominate forms of numerical error are
numerical diffusion of mass (from herein called
numerical diffusion) and numerical diffusion of
momentum. The latter results in dissipation of
energy and is henceforth called numerical
dissipation. Numerical diffusion artificially weakens
sharp gradients across the wave front. Thus, as a
wave steepens, numerical diffusion increases and
smoothes the wave front. This smoothing impedes
wave steepening and reduces the available potential
energy in wave, but increases the background
potential energy of the system (Figure 1.5 A). The
increase in background potential energy is seen in
the thickening of the layer that contains the density
gradient between the constant density upper and
lower layer, which is called the pycnocline (see
?2.2.1 for a description of available and background
potential energy). Numerical dissipation has a
similar effect in damping wave steepening. As a
wave steepens, the velocity shear increases across
the wave front, thereby increasing the artificial
spreading of momentum through model error. The
artificial spreading of momentum reduces the kinetic
energy of the wave and inhibits steepening (Figure
1.5 B).
B
C
D
Figure 1.4: Cartoon of wave steepening
and dispersion. A) An initial wave with a
unique wave speed and amplitude. B&C)
Wave steepening where the crest is
moving faster than trough. D) Wave
degeneration into a train of waves.
1.3.2 Modeling
Several hydrodynamic models for largescale
applications use the hydrostatic approximation (e.g.
POM: Blumberg and Mellor, 1987; EFDC: Hamrick,
1992; ROMS: Haidvogel, et al., 2000). The
hydrostatic approximation neglects nonhydrostatic
pressure and vertical acceleration, but models
B A
5
Figure 1.5: Cartoon of wave steepening. A) Nonlinear steepening wave with numerical
diffusion. As the wave steepens, diffusion acts across the wave front and reduces the
density gradient, available potential energy and steepness. B) Nonlinear wave
steepening with numerical dissipation. The initial wave is primarily moving in a
horizontal direction. As the wave steepens, the vertical velocities in front and behind the
wave increase, increasing the shear. The increased diffusion of momentum (dissipation)
reduces the wave?s kinetic energy and limits the wave steepness.
Wadzuk & Hodges: Hydrostatic and Nonhydrostatic Internal Wave Models
6
nonlinear horizontal momentum. The hydrostatic
approximation is generally applicable to processes
with small aspect ratios (depth: length ? O[10
3
]); for
example, most largescale ocean and atmospheric
problems. However, the hydrostatic approximation
breaks down for mesoscale systems [aspect ratio ~
O(10
1
? 10
2
)] and smallscale events, such as the
steepening and breaking of nonlinear internal waves
(Kantha and Clayson, 2000a). Hydrostatic models
often predict excessive steepening (Wadzuk and
Hodges, 2003), which induces artificial breaking and
diffusion, thereby developing a weakened pycnocline
(Laval, et al., 2003). Hydrostatic models also show
increased dissipation of internal wave energy which
artificially reduces internal wave amplitude. Internal
wave dispersion is not physically modeled in
hydrostatic models (as it requires nonhydrostatic
pressure), but may still appear as a numerical
phenomenon in a simulated evolution (Hodges and
Delavan, 2004). The aforementioned errors, which
are all numerical, result in an incorrect evolution of
an internal wave and its energy. Thus, the need for
nonhydrostatic models is evident (Horn, et al., 2001;
Laval, et al., 2003).
Several quasihydrostatic and nonhydrostatic
models for largescale systems have been developed
(e.g. Mahadevan, et al., 1996a,b; Casulli and
Stelling, 1998; Casulli, 1999; Fringer and Street,
2003). Despite this work, there has been little
research on quantifying the timestep and spatial grid
scales necessary to model nonhydrostatic behavior.
The present research shows that a fine spatial and
temporal resolution is needed to capture non
hydrostatic behavior (?4 ff.), thus significantly
increasing computational time over hydrostatic
models. The aforementioned models applied the
nonhydrostatic model over the entire domain, but
Stansby and Zhou (1998) recognized that non
hydrostatic effects are local and suggested
application of the nonhydrostatic solution only in
these limited areas. However, their work provides
neither a method to locally solve the pressure, nor a
method to locate regions where the nonhydrostatic
pressure is important. The present research has
developed a method to numerically identify regions
with significant nonhydrostatic effect.
1.4 Approach
Achieving the objectives of this study required both
hydrostatic and nonhydrostatic models. The
University of Western Australia Centre for Water
Research?s Estuary and Lake Computer Model
(CWRELCOM: Hodges, 2000) was used as the
hydrostatic model. As part of this project, a non
hydrostatic pressure solver was developed as a
modification to CWRELCOM. CWRELCOM is a
threedimensional hydrodynamic model that applies
the incompressible, Boussinesq and hydrostatic
approximations to solve the Euler equations on an
Arakawa C staggered grid with a semiimplicit
method and a moving freesurface. The model
temporal accuracy is firstorder for the flow field,
thirdorder ULTIMATEQUICKEST (Leonard,
1991) for scalar transport and firstorder for the free
surface. Surface thermodynamics are solved by bulk
transfer models that are not used or investigated in
this work. The EulerLagrange method for
momentum is thirdorder spatially accurate, while
spatial gradients in other terms are secondorder.
Details of the nonhydrostatic solver developed for
the present research are described in ?2 ff.
The existing and modified models were used to
perform the analysis fulfilling the two research
objectives. The first chapter of this document
provides the reader with the motivation, objectives
and background of this project. The second chapter
provides the numerical methods applied in the
hydrostatic model and the design of the non
hydrostatic solver, as well as the analytical methods
used to quantify the performance of hydrostatic and
nonhydrostatic models. The third chapter
demonstrates the validation and verification of the
nonhydrostatic solver. Chapter four uses the tools
developed in chapter two to examine the first
objective (comparing the differences between
hydrostatic and nonhydrostatic models), including
an evaluation of model skill on predicting internal
wave evolution and quantification of numerical error.
Chapter five develops and demonstrates a method to
isolate nonhydrostatic effects. Chapter six
concludes this work and makes suggestions for
future work in this area of hydrodynamic modeling.
Chapter 2. Methods
Chapter 2. Methods
Several different numerical and analytical methods
are used in the present work to develop the models
used in analysis and to provide a means to quantify
and compare model results. Governing equations
(?2.1.1) are used to model the physics of internal
wave evolution and are applied to the numerical
model with a pressure Poisson solution implemented
by the fractionalstep method described in ?2.1.2 and
?2.1.3. Along with the governing equations and
numerical methods, a set of approximations and
conditions are made to simplify the problem and
provide boundaries for the solution space. In
assessment of model skill, both qualitative and
quantitative techniques are used. Section 2.2
provides qualitative methods for assessing model
skill by examining the pycnocline displacement and
quantitative methods for evaluating model skill
through computing numerical dissipation and
diffusion.
2.1 Numerical Model
This sections discusses a description of the
governing equations of CWRELCOM and the non
hydrostatic solver and a justification for the inviscid
approximation (?2.1.1). The methods used to
develop the nonhydrostatics solver (?2.1.2 and
?2.1.3) and the applicable boundary conditions
(?2.1.4) are also described.
2.1.1 Governing Equations
The governing equations for the present research are
the Euler equations, the freesurface equation with
the kinematic boundary condition and the linear
equation of state for density. In Cartesian tensor
form, these are:
Incompressible momentum:
()
i
i
o
z'
nh
oi
DU 1
g
Dt x x
P1
x
?
?
??
?
?? ?
?
=?? + ?
???
?
?
?
??
?
dz
?
?
?
(2.1)
Incompressible continuity:
i
i
U
0
x
?
=
?
(2.2)
Freesurface ? Kinematic boundary condition
integrated over depth with continuity applied:
b
Udz
tx
?
?
?
?? ?
=?
??
?
(2.3)
Density:
( )
o
1S?=? +? (2.4)
In Equations (2.1)  (2.4), i = 1, 2, 3; ? = 1, 2; U
i
is
the velocity; ? is the freesurface elevation; z? is the
vertical position in the water column; P
nh
is the non
hydrostatic pressure; ? is the density; ?
o
is the
reference density; g is gravity. The haline expansion
coefficient, ?, is constant for a linear equation of
state and S is salinity. The density can be
decomposed as:
o
?=? +? (2.5)
where
? is a small departure from the reference
density. When the hydrostatic approximation is
used, as in CWRELCOM, the nonhydrostatic
pressure gradient is neglected in Equation (2.1).
Within the hydrostatic limit, vertical motions are
assumed small (horizontal velocity ~ 10
1
ms
1
,
vertical velocity ~ 10
4
to 10
5
ms
1
); therefore the
vertical momentum equation is neglected (Cushman
Roisin, 1994). However, Long (1972) showed that
the inclusion of vertical acceleration and non
hydrostatic pressure is necessary to properly model
internal wave evolution. Thus, CWRELCOM has
been modified as described in this chapter to include
a nonhydrostatic pressure component based on
Equation (2.1).
The present research?s focus is on quantifying
hydrostatic and nonhydrostatic models? ability to
represent internal wave evolution and error
accumulation. All analysis for the present research
neglects viscous effects and diffusion terms in the
momentum and transport equations for two reasons:
1) the effects of viscous damping and viscous effects
in wave breaking are at the two extremes of wave
phenomena, which are of lesser interest than the
evolution of solitons; and 2) by neglecting viscosity
and diffusion the dissipation of energy and diffusion
of mass can be used to quantify numerical error.
Viscous damping acts primarily through the wall
boundary layers, which are not modeled in the
present research, and the shear between density
7
Wadzuk & Hodges: Hydrostatic and Nonhydrostatic Internal Wave Models
layers. The viscous effect between
density layers is immaterial in the
present work as the timescale for
damping is longer than the timescale
for steepening, which is a phenomenon
of interest for an evolving nonlinear
wave. The modulus of decay provides
the timescale for significance of
viscous effects. This viscous timescale
is the time it takes for a wave?s
amplitude to decrease to e
1
of its
initial amplitude (Lamb, 1932):
? (m) ? (s) ? (d)
10000 1.267E+12 1.466E+07
100 1.267E+08 1.466E+03
1 1.267E+04 1.466E01
0.01 1.267E+00 1.466E05
1.0E05
1.0E03
1.0E01
1.0E+01
1.0E+03
1.0E+05
1.0E+07
1.0E+09
0.01 0.1 1 10 100 1000 10000
Wavelength (m)
M
odul
us of
D
ecay (
d
)
2
2
8
?
?=
??
(2.6)
where ? is the modulus of decay, ? is
wavelength and ? is viscosity. Figure
2.1 shows the modulus of decay for
several different wavelengths. This
viscous timescale is small for very short waves (e.g.
capillary waves, Kinsman, 1965), so viscosity is
important in their evolution and attenuation.
Conversely, the viscous timescale for gravity waves
(e.g. ? > O[1m]) is long compared to the timescales
of forcing and wave evolution (see discussion of
wave steepening timescale, ?3.3.3 ff.), so viscosity is
insignificant for the evolution of longer waves.
Viscous and diffusive effects in wave breaking and
mixing are not of interest in the present work as
largescale coastal oceans and lake models do not
have sufficient grid resolution to capture the details
of these processes. The present research is focused
on wave behavior prior to breaking, when the wave
can be considered a smooth material surface that
effectively suppresses turbulent diapycnal diffusion
of mass and momentum.
Figure 2.1: Modulus of decay for as a function of wavelength,
determined from Equation (2.6). Viscosity is 10
6
m
2
/s.
2.1.2 Solution Methods
The fractional step method (Kim and Moin,
1985) was used to incorporate the nonhydrostatic
solver into the existing hydrostatic CWRELCOM.
The fractional step method is mass conservative and
has secondorder temporal accuracy for 3D
incompressible flows with fine spatial discretization
and coarse temporal discretization (Casulli and
Stelling, 1998; Casulli, 1999; Armfield and Street,
1999). The fractional step method can be outlined as
follows:
? Solution of hydrostatic velocity field and
free surface
? Solution of nonhydrostatic pressure
? Update of velocity field and free surface to
reflect nonhydrostatic pressure.
The velocity field is initially approximated via
a hydrostatic solution, providing U* and W*, as well
as a hydrostatic free surface (?*).
( )
()
11
22
1
2
1
2
11
22
11
22
**
*n i1 i
i,k i,k
i,k
nn
i1,k i,k
o i,k
b
*n
i,k i,k
*n * *
ii i,ki,k i,k
i,k
b
UGUgt
x
1
gt z
x
WGW
t
UU z
x
+
++
+
?
+
+
++
?
+?
??
???
=????
?
??
??
?? ? ?? ???
??
??
=
??
?
?=?? ? ???
?
??
?
?
(2.7)
where subscripts indicate the cell center (integers)
and cell faces (fractions), as seen in Figure 2.2. The
xdirection cell length is ?x
i
, while ?x
i+1/2
is the cell
center to center distance. The Lagrangian
discretization of velocity from the EulerLagrange
method is represented by GU
n
, which uses a particle
pathline in the time ?n? velocity field to estimate the
Lagrangian momentum term at time ?n+1.? This
method reduces the artificial damping associated
with most loworder methods and can be applied to
flows that have a CourantFriedrichLewy (CFL)
condition up to two (Hodges, 2000).
8
Chapter 2. Methods
9
The hydrostatic estimated velocity [Equation
(2.7)] is subtracted from the discretized non
hydrostatic momentum equation, written as:
( )
11
22
1
2
1
2
1
2
11
22
1
2
n1 n1
n1 n i1 i
i,k i,k
i,k
nn
i1,k i,k
o i,k
b
n1 n1
nh,i 1,k nh,i,k
o i,k
n1 n1
nh,i,k 1 nh,i,kn1 n
i,k i,k
o i,k
UGUgt
x
1
gt z
x
PP
t
x
PP
t
WGW
z
++
+ +
++
+
?
+
+
++
+
+
++
++
++
+
??
???
=????
?
??
?? ? ?? ?
??
??
??
?
?
?? ?
??
??
??
??
?
?
=??
?
??
??
?
?
?
(2.8)
to obtain velocity correction equations:
11
22
1
2
n1 n1
nh,i 1,k nh,i,kn1 *
i,k i,k
o i,k
PP
t
UU
x
++
++
++
+
??
?
?
=??
??
??
? (2.9)
11
22
1
2
n1 n1
nh,i,k 1 nh,i,kn1 *
i,k i,k
o i,k
PP
t
WW
z
++
++
++
+
??
?
?
=??
??
??
? (2.10)
The discrete form of continuity [Equation (2.2)] is:
11
22
n1 n1 n1 n1
i 1,k i,k i,k 1 i,k
i,k i,k
UUWW
0
xz
++++
+
++
??
+
??
= (2.11)
which can be applied to the divergence of Equations
(2.9) and (2.10) to yield a pressure Poisson equation
(in Einstein summation):
2*
onh i
ii i
PU
xx t x
???
=
?? ??
(2.12)
Equation (2.12) can be written discretely as:
i+
1
/
2
i
1
/
2
k+
1
/
2
k
1
/
2
i,k
k+1
k1
i+1i1
11
22
11
22
11 1
22 2
n1 n1 n1 n1
nh,i 1,k nh,i,k nh,i,k nh,i 1,k
i,k i,k i,k
n1 n1 n1 n1
nh,i,k 1 nh,i,k nh,i,k nh,i,k 1
i,k i,k i,k
** * *
i,ki,k i,k i,k
o
i,k i
PPPP
1
xx x
PPPP
1
zz z
UU WW
tx z
++++
+?
++++
+?
+? +
??
??
??
??
?? ?
??
??
?
+ ??
?? ?
??
?
=+
?? ?
?
,k
??
??
??
(2.13)
Solution of Equation (2.13) provides a non
hydrostatic pressure, P , which is used to update
the horizontal velocity field using Equation (2.9).
The updated horizontal velocity field is used to
update the freesurface:
n1
nh,i,k
Figure 2.2: Schematic of grid field
+
()11
22
n1 n
n1 n1
,k,k ,k
kb
t
UUz
x
+
??
?
++
??+ ??
?
=
?=?
??
?
????
?
??
?
?
(2.14)
If Equation (2.10) is used to update the vertical
velocity, the resulting velocity field may not be
solenoidal due to the residual in the nonhydrostatic
pressure solver. To ensure a solenoidal velocity field
for mass transport, the vertical velocity is updated
diagnostically from the continuity equation:
11
22
11
22
n1 n1
i,k i,k
n1 n1
i,ki,k i,k
i,k
UU
WW z
x
++
+?
++
+?
???
?=? ?
?
??
? (2.15)
The freesurface update [Equation (2.14)] has
been treated differently in various models. The
quasihydrostatic models of Mahadevan, et al.
(1996a) and Casulli and Stelling (1998) solve only
for the hydrostatic freesurface [Equation (2.7)] and
do not account for nonhydrostatic effects on surface
evolution. Casulli (1999) addressed this problem by
correcting the surface elevation for the intermediate
step [Equation (2.7)] and after the pressure Poisson
equation was applied [Equation (2.14)]. The
approach of Casulli (1999) was to advance the
velocities, pressure and freesurface at time level
?n+1? while still applying the vertical discretization
(?z) of the previous time step (n). His approach
removes a nonlinearity that arises when using
inconsistent time discretization of ??z? and the free
surface, which Hodges (2004) showed exists, but is
secondorder in time. Thus, the firstorder temporal
accuracy of CWRELCOM makes the inconsistency
irrelevant and Casulli?s (1999) treatment appropriate
Wadzuk & Hodges: Hydrostatic and Nonhydrostatic Internal Wave Models
for the updated freesurface. A more complex
approach by Chen (2003) applies a double predictor
corrector method, updating the velocity first after the
nonhydrostatic pressure is resolved and again after
the freesurface is updated. The double predictor
corrector allows the nonhydrostatic change in the
freesurface to alter the velocity field. While this
method is secondorder temporally accurate, Chen
(2003) states that a comparison between the results
with his method and those of Casulli (1999) show no
significant improvement in model performance.
10
2.1.3 Iterative Methods ? Pressure Poisson
Equation
The pressure Poisson equation [Equation (2.13)] is
an elliptic equation, which is inefficient for direct
solution by matrix inversion (Mahadevan, et al.,
1996a) and is therefore solved with iterative
methods. Jacobi, GaussSeidel (GS), and successive
overrelaxation (SOR) iterative methods were used
in the development of the nonhydrostatic solver as a
learning exercise for the author. The Jacobi and GS
methods were impractical because of long
convergence times (i.e. redblack point GS: ~3500
iterations for the L
?
norm to converge to 10
8
in the
manufactured solution case described in ?3.1 ff).
SOR, with a redblack point iteration scheme, has
proved to have adequate convergence times (i.e. red
black point SOR: ~1000 iterations for the L
?
norm to
converge to 10
8
; this reduced the run time by about
70% for the manufactured solution case in ?3.1 ff.).
The SOR scheme used is presented in Equation
(2.16) with the optimal overrelaxation factor (?)
determined by Equation (2.17).
()
i,kn1 n
nh,i,k nh,i,k
i,k
n1 n1
i 1,k nh,i 1,k i 1,k nh,i 1,k
i,k
n1 n1
i,k 1 nh,i,k 1 i,k 1 nh,i,k 1
i,k
Q
P1P
a
aP aP
a
aP aP
a
+
++
++ ??
++
++ ??
??
=? + ??
??
??
?+
??
?
?
+
??
?
?
?
?
?
(2.16)
g
2
1sin
N
?=
??
?
??+
??
(2.17)
where Q is the source term [i.e. the right hand side of
Equation (2.13)], a is the coefficient that represents
that spatial discretization and N
g
is the number of
grid cells in the computational domain.
Convergence was measured using by the residual
(R):
(2.18)
n1 n1
i,k i,k i 1,k nh,i 1,k i 1,k nh,i 1,k
n1 n1 n1
i,k 1 nh,i,k 1 i,k 1 nh,i,k 1 i,k nh,i,k
RQaP aP
aP aP aP
++
++ ??
++ ??
=? ?
??
+
When the residual is reduced to a specified
convergence criteria (?3.1 ff), the nonhydrostatic
pressure is considered resolved and the velocity field
and free surface are updated with the new non
hydrostatic pressure. Other methods may be used to
solve the pressure Poisson equation more efficiently
(e.g. Mahadevan, et al., 1996b; Marshall, et al.,
1997; Casulli and Stelling, 1999; Casulli, 1999;
Fringer and Street, 2003). Some possibilities include
the SOR applied in linerelaxation form, conjugate
gradient method or multigrid methods. The
multigrid method is well suited to handle stiff
problems, such as the pressure Poisson equation, and
it may use a GaussSeidel iterative scheme, so the
nonhydrostatic solver could be adapted into multi
grid form. However, the multigrid method may be
difficult to implement over complex boundaries (He,
et al., 1996).
2.1.4 Boundary Conditions
The model boundary conditions were assigned so
that the vertical solid boundaries perfectly reflect
wave propagation and neither vertical nor horizontal
boundaries have viscous boundary layers. To
achieve these requirements, the velocity is assigned
Dirichlet conditions normal to solid boundaries (U =
0 or W = 0) and Neumann conditions tangential to
solid boundaries (dU/dx = 0 or dW/dz = 0). Scalar
quantities, including the nonhydrostatic pressure,
have Neumann conditions (zero gradient) at all solid
boundaries. The velocity and scalar conditions are
implemented in CWRELCOM under the
designation of ?freeslip? boundaries. The non
hydrostatic pressure at the free surface is a Dirichlet
condition requiring P
nh
= 0. The boundary
conditions for the nonhydrostatic pressure are built
directly into the nonhydrostatic solver in the SOR
coefficients, ?a,? in Equation (2.16).
2.2 Analysis Methods
Both qualitative and quantitative analysis methods
are used to compare and contrast hydrostatic and
nonhydrostatic model results. The simplest
qualitative analysis is visual examination of internal
wave evolution over several wave periods in both
models and theory. The human eye/brain
Chapter 2. Methods
combination is capable of distinguishing the overall
scope of phase/amplitude errors in an intuitive way
that provides the first level of screening for model
results. Quantitative analysis methods developed
herein are focused on the energy changes within the
system. A principle reason for conducting model
experiments using the inviscid/diffusionless
equations is that any change in the modeled total
energy (E
T
) is directly attributable to numerical error.
Therefore, energy changes are an integrator of
numerical dissipation and diffusion errors in the
model results. Furthermore, there is an energy shift
to smaller wavelengths associated with the
degeneration of an internal wave, which is used to
identify numerical error.
2.2.1 Background
Total energy (E
T
) in a system is comprised of kinetic
(E
K
) and potential energy (E
P
), written below for a
2D system with a constant breadth (B):
L
2
Ki
0b
1
EUBd
2
?
=?
??
zx
P
B
(2.19)
11
(2.20)
L
P
0b
E g zBdzdx
?
=?
??
so that E
T
may be written as:
(2.21)
TK
EEE=+
The potential energy can further be separated
into background potential energy (E
B
) and available
potential energy (E
A
) (Lorenz, 1955). The E
B
is the
potential energy when the system is at its lowest
possible energy state. That is, if a system is
adiabatically brought to rest and the density field
settles without mixing, then the system is considered
to have its lowest possible potential energy, which is
E
B
. The E
A
is the potential energy of a system
displaced from its lowest possible energy field, or
also described as the E
P
available to be transferred to
E
K
. Thus, the EA is defined as the difference
between the E
P
and E
B
:
AP
EEE? ? (2.22)
A simple monochromatic standing wave is an
oscillatory exchange of E
A
(the potential energy in
the inclined wave) for E
K
(kinetic energy when the
wave is flat), and back to potential energy again. In
an inviscid/diffusionless system, the exchange of E
A
for E
K
must be conservative, so it is convenient to
define their sum as dynamic energy (E
D
):
(2.23)
DA
EEE?+
K
The relationship between energies is seen in Figure
2.3
Mass diffusion is the spreading, or weakening,
of a stratified density field and is characterized by
the increase in the background potential energy (E
B
)
as heavier particles are diffused upwards and lighter
particles diffused downwards. Since E
B
is the
En
e
r
g
y
Wave Period
01/4 13/41/2
E
D
E
K
E
A
E
P
E
B
En
e
r
g
y
En
e
r
g
y
En
e
r
g
y
Figure 2.3: Schematic of energy behavior over a wave period, for a inviscid/diffusionless
monochromatic standing wave. The black lines above the graph show the evolution of the wave
shape. The purple ( ) line is the potential energy (E
P
), the red ( ) line is the
background potential energy (E
B
), the blue ( ) line is the available potential energy (E
A
), the
yellow ( ) line is the kinetic energy (E
K
) and the green ( ) line is dynamic energy (E
D
).
The E
P
is the sum of E
B
and E
A
and the E
D
is the sum of E
A
and E
K
.
Wadzuk & Hodges: Hydrostatic and Nonhydrostatic Internal Wave Models
density field when restored to its lowest energy state,
it can only change through physical diffusion of the
density gradient (Winters, et al., 1995). That is,
advection and wave propagation can not directly
increase EB. However, in the physical world they
may cause shear instabilities and wave breaking
which will more readily allow physical diffusion to
occur.
Resort all water parcels (grid cells) by
density, keeping track of the volume of each
grid cell.
Fill the basin, starting at the bottom with
the heaviest parcel, where the volume
associated with this parcel is spread evenly
over the entire horizontal area. Continue
filling the basin from most to least dense, as
in Figure 2.5.
12
Numerical diffusion of mass is
indistinguishable from physical diffusion of mass in
the resulting smoothing and spreading of gradients,
which makes it difficult to quantify the numerical
diffusion error when the governing equations include
physical diffusion. In a model that sets physical
diffusion to zero, any change in E
B
is model error
(Laval, et al., 2003), which is commonly referred to
as numerical diffusion. Numerical diffusion is
dependent on the coarseness of the grid; that is, fine
grids reduce numerical diffusion (Hodges, et al.,
2000). Thus, as a grid is refined and numerical
diffusion is inhibited, modelers typically presume
that model representation of E
B
and internal wave
evolution are improved; however, as shall be seen,
this is not always true (?4 ff.).
Calculate the E
P
for this resorted system
using Equation (2.20), where ??? is the resorted
density, ?z? is the height of the layer in the
resorted basin, ?dz? is the thickness of the
resorted layer and ?dx? is the length of the
basin. This is the E
B
.
In the physical world, dynamic energy (E
D
) is
decreased by: 1) transfer of available potential
energy (E
A
) and kinetic energy (E
K
) to background
potential energy (E
B
) through turbulent mixing or 2)
loss of E
K
to heat through the action of viscosity.
The latter is called ?dissipation,? while the former
(mixing) is where waves and advection enhance the
rate of physical diffusion. Shear instabilities and
wave breaking move a fluid around, (i.e. ?stirring?)
which raises some of the heavy fluid, increasing E
A
and decreasing E
K
. If physical diffusion is present, a
portion of the stirred fluid is irreversibly mixed.
Thus, a portion of the E
K
that went into E
A
to stir the
fluid is, in a sense, converted into E
B
. Therefore,
physical diffusion increases the E
B
of a system while
decreasing the E
D
of the system. The portion of E
K
lost due to diffusion is not considered dissipation of
energy because the E
K
was converted to E
B
rather
than heat. Without physical diffusion, stirred fluids
cannot irreversibly ?mix? so the heavier fluid would
eventually fall back to its lowest energy state,
thereby converting the E
A
back to E
K
with no impact
on E
B
or E
D
. Hence, in a diffusionless, inviscid
model any decrease in E
D
represents energy lost in
the form of dissipation. As dissipation decreases E
K
,
there is less energy available for conversion to E
A
and consequently a reduction in the E
D
of a system.
Figure 2.4: Flow chart of sorting algorithm for
background potential energy (E
B
).
The evolution of background potential energy
(E
B
) and dynamic energy (E
D
) provide a measure of
numerical diffusion and numerical dissipation.
While the total potential energy for any particular
density distribution is a relatively simple
computation, the calculation of background potential
energy is more difficult. Following the example of
Winters, et al. (1995), an algorithm was developed to
determine E
B
, as described in Figure 2.4. The model
approximates reality by assigning each cell in a
system a unique value for density. By example, a
5?5 random domain would result in 25 layers in the
resorted system, as in Figure 2.5.
Figure 2.5 provides a simple example, but a
larger domain results in significantly more layers
(i.e. 4380 for a 60?73 grid). The resorted density
field yields a vertical resolution that has a much finer
scale than the model?s representation of the vertical
processes. When the full resorted density field is
used in calculations, the results are ?noisy? and are a
misleading representation of the model?s resolvable
vertical structure. To return the data to the
appropriate analysis scale, the resorted density
profile is ?binned? into groups of densities to yield a
resorted density domain that has the same number of
layers as the original vertical resolution (Figure 2.6).
Chapter 2. Methods
13
The binned, resorted density field yields a domain
that produces smoother calculations (Figure 2.7).
2.2.2 Numerical Diffusion
Physical diffusion changes stratification by
smoothing density gradients transferring more dense
fluid into lower density regions and viceversa.
Typically, changes in the background potential
energy provide a measure of the physical diffusion in
a system. In a diffusionless system, any changes in
the stratification can be attributed to numerical
diffusion. Therefore, changes in the background
potential energy for a diffusionless system provide
an estimate of numerical diffusion. Physical
diffusion is modeled by diffusivity for each of the
scalars, such as salinity or temperature. Thus, the
concept of diffusivity is used to develop a numerical
diffusion coefficient based on changes in the
background potential energy. The numerical
diffusion coefficient acts as a global parameter to
compare model performance and error accumulation.
The resorted density field allows calculation of
the E
B
per unit area (A):
()
()
B
B
b
Et
t,z gz dz
A
?
=?? ?
? ?
?
(2.24)
where ?
B
is the binned, resorted density. The
evolution of E
B
is:
()
()
B
B
b
Et
dd
t,z gzdz
dt A dt
?
??
=?
??
??
?
(2.25)
The free surface for the background density
field is constant in a system without inflows or
outflows, so the derivative passes through the
integral in Equation (2.25). Discretizing Equation
(2.25) for ?N? layers yields:
( )
()
N
B
B,k k k k
k1
Et
dd
t,z z gz
dt A dt
=
??
? ?
?=??
??
?
? ?
?
??
??
? ?
(2.26)
Figure 2.5: Schematic of a random 5?5 density field resorted into 25 layers from most to least dense.
The blue indicates least dense and red indicates most dense.
The density term for layer ?k? can be represented as:
Figure 2.6: Schematic of a random 5?5 density field resorted into 25 layers from most to least dense
and then binned into 5 layers. The blue indicates least dense and red indicates most dense.
Wadzuk & Hodges: Hydrostatic and Nonhydrostatic Internal Wave Models
0 2 4 6 8 10 12 14 16 18
0
0.05
0.1
0.15
0.2
0.25
t/T
E
B

E
B
(0
)
not binned
binned
E
B

E
B
(0
)
Figure 2.7: Typical change in background potential energy (units: Joules) versus the time
normalized by a wave period (T). The green line is before the density field is binned, the blue
line is after the density field is binned, with the errorbars representing the standard deviation.
()
()
()
()
Bk
B,k k k k
k
Bk
Bk
mt,z
dd
t,z z z
dt dt A z
1d
mt,z
Adt
mt,z
A
??
??
????= ?
??
??
?
??
=? ?
??
=
(2.27)
where m?= ? and with ?zAz?= ?
k
constant in
time. The vertical mass flux, , represents the rate
of change in the resorted density field; any change in
E
B
m
B
is due only to the diffusion of mass in the density
stratified system. The relationship in Equation (2.27)
can be substituted into Equation (2.26) to obtain:
() ()
N
BBk
k
k1
Et mt,z
d
gz
dt A A
=
??
=
??
??
?
(2.28)
Fick?s law states:
C
qD
z
?
=?
?
(2.29)
where D is the diffusion coefficient, q is the mass
flux per unit area and C is the mass concentration.
The vertical mass flux per unit area in Equation
(2.28) may be expressed in the form of Fick?s law:
()
()
()
B
mz z
Dz
Az
??
=?
?
(2.30)
Making the substitution for the vertical mass flux per
unit area into Equation (2.28), yields
( ) ()
N
B
kk
k1
Et t,z
d
Dgz
dt A z
=
????
=?
??
?
??
?
Bk
D
(2.31)
which can be redefined as
N
kk
k1
F
=
=?
?
(2.32)
where
( )
()
B
Bk
kk
Et
d
F
dt A
t,z
gz
z
??
=
??
??
??
?=?
?
(2.33)
F is computed from model results using the rate
of change of E
B
and ? is computed from the resorted
density profiles. Local diffusion coefficients (D
k
)
cannot be directly computed from the global change
in E
B
in Equation (2.32). However, we can estimate
an approximate global diffusivity based on a
simplified treatment of the density field. In all the
present investigations, hyperbolic tangents are used
to construct the density field:
14
Chapter 2. Methods
1
py
base py
py
x
ez cos z
L1
tanh
h
2
2
????
??
+? ? ?
??????
????
?=? + ??
??
??
(2.34)
where ?
base
is the density at the center of the
pycnocline, ??
py
is the density change across the
pycnocline, ? is the amplitude, x is the horizontal
grid location, L is the length of the basin, z is the
vertical grid location and h
py
is the thickness of the
pycnocline. The height at the center of the
pycnocline is z
py
and the initial thickness (h
py
) is split
equally on each side of the pycnocline center line to
designate the initial pycnocline position. Equation
(2.34) yields a representative density profile (Figure
2.8). The hyperbolic tangent profile models a
continuous stratification in a real world system with
three distinct layers; two approximately uniform
density layers separated by a region of rapidly
changing density. Therefore, dividing a stratified
system into three layers: 1 = upper, 2 = pycnocline
and 3 = lower (Figure 2.8); Equation (2.32) can be
written
(2.35)
11 2 2 3 3
DDD?+?+?=
15
F
The upper and lower layers have approximately
constant density, so ?
1
= ?
2
= 0, which leaves
()
B
2
2
2
k
Et
d
dt A
F
D
gz
z
??
??
?
==?
??
?
?
?
(2.36)
Equation (2.36) can be approximated as
()
B
2
13
py
py
Et
d
dt A
D
gz
h
??
??
?
??
??
?
? ??
??
??
(2.37)
where ?
1
and ?
3
are the densities of the upper and
lower layers, respectively. Equation (2.37) yields a
global approximation for D across the pycnocline
derived from the change in the background potential
energy system and changes in the pycnocline
thickness and density field. This approximation
provides a tool with which numerical diffusion can
be estimated for the simple threelayer system. Since
this approximation is taken at each time step there is
the possibility of oscillations due to smallscale anti
diffusive fluxes. The smallscale antidiffusive
fluxes arise from nonmonotonic advection
(Leonard, 1991). Oscillations form around sharp
wave fronts (i.e. a discontinuity) and are referred to
as Gibb?s phenomenon (Kreyszig, 1999). The
transport algorithm used in the present research is
ULTIMATE QUICKEST; ULTIMATE (Leonard,
1991) smoothes the oscillations in QUICKEST
(Leonard, 1979) so that the scalar concentration
decreases monotonically in a onedimensional
system. Lin and Falconer (1997) found that
ULTIMATE QUICKEST is not monotonic for a
twodimensional system, such as used in the present
research.
998 999 1000 1001 1002
0
1
2
3
4
5
6
7
Density (kg/m
3
)
H
e
i
ght
(
m
)
H
e
i
ght
(
m
)
Layer 1
Model performance over multiple wave periods
is the focus of this work; thus, the change in E
B
and
the pycnocline density gradient are binned over a
wave period to smooth antidiffusive fluxes and
provide a clear indication of numerical diffusion.
The average form of Equation (2.37) is
Layer 2
()
B
T
2
13
py
py
T
Et
1d
dt
Tdt A
D
1
gz dt
Th
??
??
?
??
=?
??
???
??
?
??
?
?
(2.38)
Layer 3
Figure 2.8: Typical density profile constructed by a
hyperbolictangent function. The dashed lines
separate layers 1, 2 and 3.
where T is the wave period. The change in h
py
is
slow if the upper and lower layer densities remain
Wadzuk & Hodges: Hydrostatic and Nonhydrostatic Internal Wave Models
constant and the pycnocline has not diffused over the
entire water column (Figure 2.9). With this
assumption, Equation (2.38) can be discretized as
()
() ()
{}
1
2
1
2
Bn Bn1
n
2 n
13
py
py
1
ET ET
AT
DT
gz
h
?
?
?
?
?
=?
??
???
??
??
(2.39)
where the over bar represents the binned data over a
wave period.
2.2.3 Numerical Dissipation
Energy dissipation in a physical fluid is controlled by
molecular viscosity acting at the smallest scales of
velocity shear. However, the energy at the smallest
scales of motion is fed by the larger scales of motion
(Kantha and Clayson, 2000b; Kundu and Cohen,
2002), so the larger scales are considered to be the
energycontaining scales and control the rate at
which viscosity can work. In modeling, turbulence
is commonly discussed in terms of the eddy viscosity
associated with a largescale velocity shear, so that
the eddy viscosity characterizes the enhancement of
turbulent dissipation for the largescale features. The
concept of eddy viscosity is wellentrenched in our
understanding of dissipative phenomenon, so it is
useful to develop a measure of a numerical viscosity
associated with the dissipation by numerical error.
The mechanical energy dissipation rate,
typically presented on a unit mass basis (?),
(Batchelor, 1967) in a 2D, incompressible fluid is
22
uw u w
zx x z
??
?? ? ??????
?=? + + +
??
?????
?? ? ?
?????
??
2
?
?
?
(2.40)
In a diffusionless system, the dissipation rate,
integrated over a volume, is equivalent to the rate of
change of the E
D
.
()
D
2
o
22
dE
d
dt
uw
zx
uw
22
xz
?
?
=???
?
????
=?+?? +
?
??
??
??
?
?
?
???? ??
d+ +?
?
?? ??
??
?? ??
?
?
?
?
(2.41)
Since
? ? , the density perturbation from the
reference density (
? ) can be neglected, yielding
2
D
o
22
dE uw
dt z x
uw
22
xz
?
?
????
??? +
?
??
??
??
?
?
?
???? ??
d+ +?
?
?? ??
?? ??
?
?
?
(2.42)
The principle velocity shear (and hence
numerical dissipation) occurs across the pycnocline.
A characteristic numerical viscosity associated with
numerical dissipation is developed that is consistent
with the shear characteristics in the pycnocline.
Integrating Equation (2.42) over the volume of the
pycnocline and substituting the characteristic
numerical viscosity of the pycnocline (?
py
) for the
molecular viscosity (?) yields:
2
D
opy
22
py
dE uw
dt z x
uw
22
xz
?
????
??? +?
??
??
???
?
?
???? ??
+ +??
?? ??
?? ??
?
?
(2.43)
where represents a spatial mean. If the
pycnocline volume (
py
? ) is characterized by
horizontal area (A) and pycnocline thickness (h
py
)
0 1 2 3 4 5 6 7 8 9 10
0
0.1
0.2
t/T
h
py
/
H
h
py
/
H
Figure 2.9: Typical change in pycnocline thickness (normalized by total depth, H) over time
(normalized by the wave period, T). This case was scenario 2 for the 30?73 grid (?3.3).
16
Chapter 2. Methods
then numerical viscosity, binned over a wave period,
may be approximated as
() ()
1
2
py D n 1 D n
n
2
opy
1
22
1
ET ET
T
uw
Ah
zx
uw
22
xz
+
?
?
??
?= ?
??
?
? ?
??? ??
?
?? +
?
??
?? ?
??
?
??
??
????? ??
?
++
?
?? ??
?? ??
?
??
(2.44)
where an overbar represents the binned data over the
wave period. For internal waves, the dominant
velocity gradient is ?u/?z, so Equation (2.44) can be
reasonably reduced to
() ()
1
2
py D n 1 D n
n
1
2
opy
1
ET ET
T
u
Ah
z
+
?
?
?
?= ?
?
?
????
?????
??
??
??
??
?? ?
??
????
?
?
(2.45)
Equation (2.45) provides an approximation for the
numerical viscosity in the model, which is used to
compare the relative importance of numerical
viscosity with molecular viscosity.
2.2.4 Spectral Analysis
As modeled internal waves evolve, dynamic energy
is transferred to smaller scale waves or lost to
numerical error. The computation of numerical
viscosity (?2.2.3) allows quantification of the energy
dissipated. In this section we discuss how spectral
analysis is used to quantify the transfer of energy
into shorter wavelength features.
Spectral analysis describes the distribution of
signal power over wave frequencies (or
wavenumbers). For this work, the spectral analysis
decomposes the spatial structure of pycnocline
displacement into the power spectral density
associated with different wavenumbers (k = 2?/?). A
Matlab
?
toolbox function (i.e. ?pwelch?) was used to
obtain the power spectral density. Welch?s method
is used, which is outlined in the following steps: 1)
the data is divided into overlapping segments, 2) a
Hamming window is used to calculate the modified
periodogram for each segment and 3) the modified
periodograms are averaged to obtain the estimated
power spectral density for the entire data set (Signal
Processing Toolbox User?s Guide, 1998). A small
data set may reduce the resolution of the power
spectral density estimation, thus the data set is
replicated to create a longer data set (i.e. in this
work, the replicated data set is used to include basin
scale waves). For all analysis in this work, the signal
data is divided into eight segments with 50% overlap
between adjacent segments and a Hamming window
applied on each segment, 32 replicates of the initial
half wavelength wave are used (to total 16
wavelengths). The power spectral density is
computed at each binned period for all wavelengths,
thus the data set is sampled at 2?/?x, which
corresponds to wavenumber.
As an internal wave steepens and forms
solitons, the initial longwavelength wave evolves
into a train of shortwavelength waves. The power
spectral analysis provides the power distribution of
the internal wave over different wavenumbers. As a
wave evolves, an energy shift may occur,
transferring energy from the longwavelength wave
to smallerwavelength waves. For instance, the peak
energy at time = 0 coincides with the initial, long
wavelength wave. However, as the wave evolves
and dynamic energy (E
D
) decreases, there is an
energy transfer through different wavelength waves,
which appears as a shift in the peak power
wavenumber through time. The shift in peak power
should coincide with the development of solitons. If
a wave does not degenerate into solitons, the peak
power should remain with the initial, long
wavelength wave. If a bore develops, then there may
be a shift in peak power as the bore provides higher
wavenumber components since the wave is no longer
sinusoidal. Comparing expected wave behavior with
the spectral analysis provides an assessment of
model skill, in terms of the model?s ability to
represent soliton formation and the evolution of
model error.
17
Wadzuk & Hodges: Hydrostatic and Nonhydrostatic Internal Wave Models
18
Chapter 3. Verification and validation
Chapter 3. Verification and validation of
the nonhydrostatic solver
pi/4 pi/2 3pi/4 pi 5pi/4 3pi/2 7pi/4 2pi
pi/4
pi/2
3pi/4
pi
5pi/4
3pi/2
7pi/4
2pi
9pi/4
1
0.8
0.6
0.4
0.2
0
0.2
0.4
0.6
0.8
1
Verification and validation of a computer code are
two necessary and basic steps to demonstrate that the
code does what it is intended to do. According to
Roache (2002), code verification is an evaluation of
error from a known solution, a purely mathematical
exercise, while code validation demonstrates the
accuracy with which the mathematical model
captures the physical phenomena based on theory or
measured in the field or laboratory. This chapter
verifies the nonhydrostatic pressure solver with an
analytical solution, validates the solver?s
convergence with a simple test case and validates the
model?s solution against theory and a laboratory
experiment.
Figure 3.1: Solution space for the manufactured
solution. The colorbar represents pressure in
kg/ms
2
.
3.1 Verification by the Method of
Manufactured Solutions
The purpose of the manufactured solution is to verify
the model?s accuracy and establish the convergence
criterion for the numerical model used in this work.
3.1.1 Setup
Roache (2002) suggests using the ?Method of
Manufactured Solutions? to verify a numerical code.
This method uses a continuous mathematical
solution independent of the code. The manufactured
solution must be nontrivial and should exercise all
terms in the numerical code and the corresponding
boundary conditions. The manufactured solution is
applied to the numerical model?s homogeneous
governing equations, producing a nonhomogeneous
source term that is discretized and added to the
numerical model?s governing equations; the resulting
discrete nonhomogeneous governing equations are
used to approximate the exact manufactured solution.
This study uses a manufactured twodimensional
analytical solution:
(3.1) Pcos(x)cos(z)=
Figure 3.1 shows the solution for Equation (3.1).
The pressure Poisson equation has the form
22
22
PP
Q
xz
??
+=
??
(3.2)
The second derivative of Equation (3.1) is:
2
2
2
2
P
cos(x) cos(z)
x
P
cos(x) cos(z)
z
?
=?
?
?
=?
?
(3.3)
This produces a manufactured source term:
(3.4) Q 2 cos(x) cos(z)=?
which is used as the discrete source term for the non
hydrostatic solver [e.g. the right hand side of
Equation (2.13)] to compute the discrete pressure
field of the manufactured solution. Thus, the
manufactured solution provides a known solution
space which can be compared to the model?s solution
driven by the source term, Equation (3.4). This
approach allows verification of model performance
for different spatial grids and time steps.
The verification test domain is a square box, 2?
m in length and 7.5 m in height, represented by four
different grids, as described in Table 3.1. The
boundary conditions of the manufactured solution
are the same pressure boundary conditions used in
other simulations within the present research
(?2.1.4); the domain top (freesurface) has a
Dirichlet boundary condition, while all other sides
(solid boundaries) have Neumann boundary
conditions.
3.1.2 Error Analysis
It is necessary to examine the grid convergence error
to assess the model?s accuracy at different grid
19
Wadzuk & Hodges: Hydrostatic and Nonhydrostatic Internal Wave Models
20
)
resolutions. Following Roy (2003), grid
convergence error is analyzed using the L
1
and L
2
spatial error norms. The spatial error norms consider
the entire domain to estimate the order of accuracy of
the model.
As defined by Roy (2003), the discretization
error at grid point (x,z):
(3.5) () () (
D grid exact
x,z P x,z P x,z?= ?
is used in the L
1
and L
2
spatial error norms:
()
()
grid
grid
N
D
n1
1
grid
1/2
N
2
D
n1
2
grid
x,z
Lnorm
N
x,z
Lnorm
N
=
?
=
?
?
=
?
?
??
=
??
??
?
?
?
(3.6)
Table 3.1: Grids used in the
manufactured solution (?3.1).
where N
grid
is the number of grid points in the
computational domain. Figure 3.2 shows the L
1
and
L
2
spatial error norms. The three different grids used
(Table 3.1 cases A  C) were refined in both the
horizontal and vertical direction. As the grid is
refined, the spatial error norms decrease with
secondorder behavior, verifying that the non
hydrostatic solver is spatially secondorder accurate.
In real world problems, analytical solutions
typically do not exist. Another way to look at error
is the residual [from Equation (3.2)]:
22
grid grid
2
P (x,z) P (x,z)
R(x,z) Q(x,z)
xz
??
?? ?
2
Case Grid (nx ? nz)
A 10 ? 10
B 20 ? 20
C 40 ? 40
D 41 ? 30
(3.7)
Figure 3.2: The L
1
and L
2
spatial error norms of the manufactured solution for grid
refinement (Table 3.1 cases A ? C).
10
1
10
2
10
3
10
2
10
1
# grid points in one direction
No
r
m
L
1
norm
L
2
norm
1
st
order slope
2
nd
order slope
No
r
m
No
r
m
Chapter 3. Verification and validation
The L
1
, L
2
and L
?
residual norms are then defined:
()
()
()
grid
grid
N
grid
n1
1R
grid
1/2
N
2
grid
n1
2R
grid
Rgrid
Rx,z
Lnorm
N
Rx,z
Lnorm
N
L norm max R x, z
=
=
?
=
??
??
?
=
??
??
=
?
?
?
(3.8)
Figure 3.3 shows the L
1
, L
2
and L
?
residual
norms for increasing number of iterations that the
nonhydrostatic solver performs for the
manufactured solution. This iteration count is
equivalent to the number of iterations per timestep of
the pressure solution in an unsteady problem (?3.2).
The change in residual norms with number of
iterations measures the model?s evolution towards
convergence for a specific amount of computational
effort. The norms continue to decrease with
increasing number of iterations until the residual
calculation [Equation (3.7)] reaches machine
accuracy. However, the solution at machine
accuracy will differ from the analytical solution by
the truncation error in the discrete equations (Hirsch,
1988). Convergence is defined when the norms
decrease by a specified order of magnitude. Hirsch
(1988) gives the example of, ?a secondorder
accurate space discretization, with ?x = 10
2
, will
produce an error of the order of 10
4
on the solution,
which cannot be reduced further even if the residual
equals 10
14
.? Thus, for secondorder discretization
of the nonhydrostatic solver, the solution is
considered converged when the L
?
residual norm
decreases by four orders of magnitude. The
1E12
1E11
1E10
1E09
1E08
1E07
1E06
1E05
0.0001
0.001
0.01
0.1
1
10
100
1 10 100 1000 10000
# Iterations
Nor
m
norm inf
Norm 2
Norm 1
Nor
m
Nor
m
Nor
m
Figure 3.3: Residual norms for manufactured solution (case D, Table 3.1). Solid blue line is L
?
norm,
dotted blue line is L
2
norm and dashed red line is L
1
norm. The number of iterations are the number of
iterations the nonhydrostatic solver performs per timestep.
21
Wadzuk & Hodges: Hydrostatic and Nonhydrostatic Internal Wave Models
manufactured solution drops by four orders of
magnitude (Figure 3.3) in about 50 iterations. Figure
3.4 shows the comparison of L
?
spatial error and
error norms. The L
?
spatial error norm flattens at the
same number of iterations (~50) where the L
?
residual norm has decreased by four orders of
magnitude. Thus, the L
?
residual norm can be used
to reasonably represent the error when it is not
possible to calculate the L
?
spatial error norm (i.e.
where there is no analytical solution).
In summary, applying the nonhydrostatic
solver to the manufactured solution shows that as the
grid is refined, the error decreases and the computed
solution converges to the manufactured solution.
The spatial error norms verify that the method is
secondorder accurate. The residual norms are a
good proxy for the spatial error norms when an
analytical solution is not available. Thus, the L
?
residual norm is used to establish the convergence
criterion in the above verification of the non
hydrostatic solver; which is the reduction of the L
?
residual norm by four orders of magnitude. This
convergence criterion is used in all further model
simulations.
3.2 Validation and Convergence for an
Unsteady Internal Wave
The purpose of this test case is to determine the
number of iterations to reach a converged solution
(i.e. based on ?3.1) in an unsteady internal wave.
This test case is sufficiently simple, yet exemplifies
many of the nuances of internal wave modeling.
3.2.1 Setup
The validation test basin used was 10 m long and 7.5
m deep. The initial wave was a cosine wave with
amplitude of 1.125 m, where the upper layer depth to
total depth ratio was 0.3. A hyperbolic tangent
function [Equation (2.34)] was used to construct the
0.0000001
0.000001
0.00001
0.0001
0.001
0.01
0.1
1
10
100
1 10 100 1000
# iterations
No
rm
residual
spatial error
No
rm
Figure 3.4: Spatial error and residual L
?
norms for manufactured solution for case D. The red
dotted line is where the L
?
residual norm decreases by four orders of magnitude (i.e. the
convergence criterion). The number of iterations are the number of iterations the nonhydrostatic
solver performs per timestep.
22
Chapter 3. Verification and validation
density profile (Figure 2.8). The test basin?s density
profile has an initial pycnocline that is 5m thick with
a density change of 4 kg/m
3
. A square grid of 0.25
m * 0.25 m was applied to this domain.
3.2.2 Results
Akin to the manufactured solution in ?3.1, the L
1
, L
2
and L
?
residual norms for the internal wave test case
decrease with increasing number of iterations that the
nonhydrostatic solver performs within each timestep
of the simulation (Figure 3.5). However, the
decrease in convergence is slower than the
manufactured solution (Figure 3.6). The L
?
residual
norm decreases by four orders of magnitude in about
50 iterations for the manufactured solution case,
while the L
?
residual norm decreases by the same
amount in about 180 iterations for the internal wave
test case. The longer convergence time for the test
case compared to that of the manufactured solution is
discussed in ?6.
The estimated root mean square error is used to
evaluate the density field for the different number of
iterations that the nonhydrostatic solver performs
per timestep. The estimated RMS error compares the
difference in density at each grid cell within the
pycnocline between a simulation with 1000
iterations, which is considered the converged case,
against simulations with a lesser number of iterations
(i.e. 10, 50, 75, 100, 200 and 500). The estimated
root mean square error (?
RMS
) of the density field is:
()
grid
1/ 2
N
2
k,n 1000,n
RMS n1
py grid
k
N
=
??
???
??
?
?
=
??
??
??
?
?
rapid change in the density profile (e.g. Figure 2.8,
(3.9)
where k is different numbers of iterations the non
hydrostatic solver performs per timestep and ??
py
is
the density change across the pycnocline. The
estimated RMS error is computed only within the
pycnocline (e.g. layer 2 in Figure 2.8) because layers
1 and 3 are essentially uniform density. The
pycnocline is defined by the region where there is
1.E08
1.E07
1.E06
1.E05
1.E04
1.E03
1.E02
1.E01
1.E+00
1.E+01
1.E+02
1 10 100 1000
# Iterations
No
rm
Norm 1
Norm 2
Norm inf
Figure 3.5: Residual norms for the internal wave test case. The solid purple line is the L
?
norm, the
dotted blue line is the L
2
norm and the dashed red line is the L
1
norm. The number of iterations are the
number of iterations the nonhydrostatic solver performs per timestep.
23
Wadzuk & Hodges: Hydrostatic and Nonhydrostatic Internal Wave Models
24
n
00
or
s
r
by
in
ror
drostatic pressure solution requires
more
n
1.E08
1.E07
1.E06
1.E05
1.E04
1.E03
1.E02
1.E01
1.E+00
1.E+01
1.E+02
1 10 100 1000
# Iterations
No
rm
Internal wave test case
Manufactured solution
No
rm
No
rm
No
rm
No
rm
No
rm
Figure 3.6: L
?
norm for the manufactured solution (case D, dashed blue line) and the internal
wave test case (solid red line). The number of iterations are the number of iterations the non
hydrostatic solver performs per timestep.
the pycnocline is demarcated between 998.1 kg/m
3
and 1001.9 kg/m
3
). The RMS error is normalized by
the density change across the pycnocline (e.g. in this
case, 3.8 kg/m
3
). Figure 3.7 shows the RMS error
values for the six different nonhydrostatic solver
iterations per timestep (i.e. 10, 50, 75, 100, 200 and
500). Using 200 iterations per timestep, the RMS
error is one or more order of magnitude smaller tha
the 10, 50, 75 and 100 iteration cases. After the
initial growth in the first 200 timesteps, the error
growth rate for the 200 iteration case [O(10
5
) in 8
timesteps] is slower than the growth for cases with
less iterations [i.e. the 100 iteration case grows to
O(10
3
) in 800 timesteps]. Of course, the RMS err
at 200 iterations is larger than the 500 iteration case,
but using 500 iterations increases the computational
time (200 iterations: 0.61 s/timestep; 500 iterations:
0.75 s/timestep) without necessarily improving the
model results. That is, as discussed in ?3.1.2, the
solution cannot reduce the error below truncation a
iterations are increased past the convergence
criterion (decrease in L
?
residual norm by fou
orders of magnitude) (Hirsch, 1988). At 200
iterations, the L
?
residual norm has decreased
four orders of magnitude, so any further decrease
the L
?
residual norm provided by increasing the
number of nonhydrostatic solver iterations per
timestep is unnecessary. This estimated RMS er
analysis supports the findings from the norm analysis
in Figure 3.5.
The nonhy
computational time. For example, on a Dell
Workstation PWS530 Xeon Processor with a 2.4
GHz CPU and 3.5 GB of RAM, the hydrostatic
model takes 0.31 seconds/timestep, while the no
hydrostatic model takes 0.66 seconds/timestep for
Chapter 3. Verification and validation
25
,
e
ame
l,
on is for
on equation estimates the
press
ls
d
s
3.3 Validation Against Regime Theory and
This section compares the numerical model?s
ves
3.3.1 Regime Theory
in referred to as Horn)
ry
, 2)
that
oretical
imescales
prov ning
awn
the same timestep and grid resolution where the
convergence criteria of the nonhydrostatic solve
reducing the L
r is
?
norm by four orders of magnitude;
for a grid with 60 horizontal grid cells ? 73 vertical
grid cells (0.01 m horizontal resolution ? 0.004 m
vertical resolution) and a timestep of 0.012 seconds
simulating 10 wave periods requires 90000 timestep
iterations, the hydrostatic model takes 7.9 hours,
while the nonhydrostatic model takes 16.4 hours.
However, the above mechanics understates the scop
of the problem: as discussed in ?4 ff., the non
hydrostatic model requires a finer timestep and
spatial resolution than the hydrostatic model to
capture the resolvable processes, therefore a more
meaningful comparison is the computational time
required for a coarseresolution hydrostatic model to
a fine resolution nonhydrostatic model. While the
choice of the coarseresolution grid is somewhat
arbitrary, for illustration we might consider a
hydrostatic solution on a 30 ? 73 grid on the s
problem, with a 0.2 s time step and 1800 timestep
iterations. The resulting solution time is 0.09 hours
(5.5 minutes). With this example, the difference in
computational expense between the hydrostatic and
nonhydrostatic models is evident; the spatial and
temporal requirements of the nonhydrostatic mode
combined with the pressure Poisson solution, render
it a huge computational undertaking. The solution of
the pressure Poisson equation accounts for 89% of
the increase in simulation time. Updating the
velocity and freesurface field increases the
simulation time by 10%. The above discussi
a twodimensional solution; a threedimensional
solution may further increase the difference in
computational effort.
The pressure Poiss
0 200 400 600 800 1000
10
8
10
7
10
6
10
5
10
4
10
3
10
2
10
1
10
0
# Timesteps
10
50
75
100
200
500
?
RMS
ure at a point using the pressure from the
surrounding cells. Other nonhydrostatic mode
(e.g. Marshall, et al., 1997; Casulli, 1999) use the
pressure from time ?n? as the initial estimate of the
pressure field for time ?n+1,? however the present
research found that the time ?n? pressure may be a
poor approximation for the time ?n+1? pressure fiel
due to the internal wave propagation (?6.1 ff.).
Setting the pressure field to zero at each timestep
provides a pressure Poisson solution that converge
faster than solution with the previous timestep?s
pressure field (Figure 3.8).
Figure 3.7: Root mean square error (?
RMS
) for
density field in the test case. Each line
represents the ?
RMS
for a different number of
iterations per timestep. The density range is
3.8 kg/m
3
.
Laboratory Experiment
representation of internal wave evolution to the
evolution of several laboratoryscale internal wa
and theoretical wave evolution.
Horn, et al. (2001; here
compared internal wave evolution in a laborato
experiment with four theoretical timescales of
different internal wave phenomena: 1) damping
steepening, 3) shear instability (KelvinHelmholtz
billow formation) and 4) internal bore formations.
The theoretical timescales are a function of the basin
and internal wave dimensions, such as length, height
and amplitude; Figure 3.9 shows a schematic of a
basin. From the theoretical timescales, Horn
developed a regime diagram (e.g. Figure 3.10)
uses the theoretical timescales to determine the
prevailing wave phenomena under different
conditions. Horn tested and validated the the
regimes in a laboratory experiment.
Regime diagrams based on the t
ide the boundaries between damping, steepe
into solitons, and formation of KelvinHelmholtz
billows or undular bores. A single ?universal?
regime diagram (i.e. for all waves) cannot be dr
as the regime boundaries are functions of wave
amplitude, upper layer thickness, total depth, basin
Wadzuk & Hodges: Hydrostatic and Nonhydrostatic Internal Wave Models
Figure 3.8: L
?
norm for the nonhydrostatic solver when the pressure at
the present timestep is approximated from zero (solid blue line) and the
previous timesteps solution (dashed red line).
26
ata
in
rn, Regime 1 in Figure 3.10
is the
ar
n
rnal
d,
while the wave evolution in regimes 1 and 2
ct the
d for
n?s Laboratory Experiment and Model
Simulations
series of model simulations: length
l
length, pycnocline thickness, and the density change
across the pycnocline. While there are many
possible ways to graph this multidimensional d
set, Horn showed that a graph of wave amplitude to
upper layer depth ratio (a/h) and upper layer depth to
total depth ratio (h/H) provides a clear delineation of
the timescales and regimes. Figure 3.10 is an
example of a regime diagram for the basin used
Horn?s experiments.
As discussed in Ho
region where viscous damping dominates
nonlinear steepening, typically for internal waves
with small amplitudes. Regime 2 is where nonline
steepening is dominant and solitons are formed.
Regime 2 flows have been observed many times i
different field studies (e.g. Boegman, et al., 2003;
Farmer, 1978; Wiegand and Carmack, 1986).
Regime 3 is supercritical flow, with large wave
amplitudes and a shallow pycnocline, which are
associated with internal hydraulic jumps and inte
bores. Regime 4 physics are dominated by Kelvin
Helmholtz billows, which predominantly occur in
shallow systems with deep pycnoclines; they are a
result of shear instabilities at the interface between
layers. Regime 5 physics are characterized by large
amplitude internal waves that develop into bores and
may include KelvinHelmholtz billows. Waves
developing in the highmixing regimes 3, 4 and 5
(supercritical flow, KelvinHelmholtz billows and
bores) have the onset of principle phenomena
occurring within onequarter of the wave perio
(damping and soliton formation) typically requires
one or more wave periods to significantly affe
wave characteristics. Horn?s regime diagram
provides a framework to validate the nonhydrostatic
solver. Simulation experiments were performe
several scenarios corresponding to the different
regimes.
3.3.2 Hor
The dimensions of Horn?s laboratory experiment
were used for a
(L) = 6 m, height (H) = 0.29 m and width = 0.3 m
(Figure 3.9). The regime diagram in Figure 3.10 is
defined for the above scales. Horn?s experimental
basin is a closed system, with solid boundaries on al
sides and the top. The present model is designed
0.0001
0.001
0.01
0.1
1
10
100
1 10 100 1000
# iterations
Nor
m
P = Previous Pressure
P = 0
Hh
h
H
a
?
L
?
2
?
1
Figure 3.9: Schematic of basin.
Chapter 3. Verification and validation
27
1 m)
inic
f
rn?s
tank
diffe
3
hown in Table 3.2;
these
a /
h
h / H
Figure 3.10: The regime boundaries from laboratory experiments; ordinate: amplitudeupper
layer depth ratio, abscissa: upper layer depthtotal depth ratio. A typical interface thickness, 1
cm, was used to determine the timescales, T
KH
and Td. The star represents KelvinHelmholtz
billows and bore, the diamond is the broken undular bore, the triangle is the solitons, the
square is steepening, and the circle is damped linear waves. (from Horn, et al, 2001)
only for freesurface simulations, so this is one area
where the model and experimental conditions
diverge. However, the freesurface remains
essentially flat (maximum displacement of 0.00
as the initial conditions are an entirely barocl
flow, which has no significant coupling to the free
surface (Kantha and Clayson, 2000b). The model
boundary conditions are as described in ?2.1.4. The
model viscosity is set to zero, which is a second area
where the model and experimental conditions
diverge. However, as discussed in ?2.1.1, this only
effects the damping rate of modeled waves at fine
grid scales when the numerical dissipation is less
than the physical dissipation. Thus, it can be
expected that the model may show development o
solitons (regime 2) under conditions when Ho
experiments show damped linear waves. Modeling
the viscosity in Horn?s experiment is impractical as
the principle viscous damping in the experiment is
not the shear across the wave interface, but is the
boundary layers on the sides and top of the tank.
Modeling the drag of these boundary layers is
impractical with the present model and of little
relevance to the present subject of interest.
Horn?s experiments were set up in a tilted
with two fluids of different density. The density
rence between the two layers was 20 kg/m
across a pycnocline thickness of 0.01 to 0.02 m. The
tilted tank was quickly moved to a horizontal
position to initialize the wave. This impulsive start
creates an initial condition that is approximately a
linear tilt. Model simulations of Horn?s experiments
are initialized with a cosine wave having an
amplitude that corresponds to the tank?s initial angle
of tilt in Horn?s experiment. The difference in initial
wave setup is the third area of divergence between
the model simulation and the laboratory experiment.
However, this difference is not expected to affect the
wave evolution. A three layer hyperbolic tangent
salinity profile [Equation (2.34)] was used to
establish the initial density profile that is close to
Horn?s condition (Figure 3.11).
Nine scenarios were modeled to reproduce
some of Horn?s experiments, as s
scenarios were chosen based on Horn?s
reported results (see Horn, et al., 2001 Figures 5 and
6). These scenarios exemplify the influence of
nonlinearity, and the subsequent development non
Wadzuk & Hodges: Hydrostatic and Nonhydrostatic Internal Wave Models
28
o
?s not available for numerical
qualitative comparisons are made to
f
n
aboratory
essary to review the principal
diffe
osity in
e
ary
el
the
e
hydrostatic pressure gradients, on internal wave
evolution. Group A maintains the same depth ratio
(h/H) and varies the amplitude ratio (a/h) to illust
the relationship between nonlinearity and initial
wave amplitude. Group B changes the depth ratio
and amplitude ratio. Changing the depth ratio als
influences the nonlinearity of the system; that is, as
the depth ratio decreases, the nonlinearity increases.
Scenarios 6 and 9 have a depth ratio of 0.5. At this
depth ratio, the wave is considered weakly nonlinear.
Therefore nonlinear steepening and soliton formation
are not significant wave evolutions. Scenario 6 has a
small amplitude ratio (a/h = 0.18), so the expected
phenomenon is a damped wave. Scenario 9 has a
large amplitude ratio (a/h = 1.0) and Kelvin
Helmholtz billows are expected to form.
3.3.3 Results
rate
Horn raw data was
comparisons, so
figures from Horn, et al. (2001). Specifically, the
model simulations are compared qualitatively against
Horn?s laboratory experiment to demonstrate the
model?s ability to capture soliton development in
Figures 3.12  3.14. The model is compared to the
experimental results by comparing the evolution o
pycnocline displacement at a point in the tank; the
amplitude of the wave and solitons, time of
emergence of solitons and number of solitons withi
a wave series are the three major points of
comparison.
Before a discussion of the model and l
results, it is nec
rence between the model simulations and
Horn?s experiments; that is the neglect of visc
the model simulations. The experimental tank is
narrow and has a lid, so there is significant area (i.e.
the surface area of the interface is 1.8 m
2
and the
surface area of the tank is 3.8 m
2
) that creates a
viscous boundary layer which will dominate the
wave damping. The computational domain in th
model simulations is the same size, but the bound
layers are not modeled. Without viscosity, the mod
simulations have no physical damping to restrain
wave steepening. The model does have numerical
dissipation which may cause some damping,
however, the numerical viscosity due to numerical
dissipation is less than molecular viscosity for
wave periods up to the emergence of solitons (?4.1
ff.). Therefore, as steepening is occurring, the wav
experiences negligible damping, so the modeled
wave tends to develop solitons sooner and with
greater amplitudes than in the experiments with
viscous effects.
Wave Evolution
Interfacial displace
wave evolution. T
ment can be used as a measure of
he interfacial vertical
Table 3.2: Scenarios used in the laboratoryscale
case (?3.3)
Scenario h/H a/h
Group A 10.300.15
20 3
30.300.45
40 6
50.300.9
Group B 60.500.18
70423
20.300.
80245
90.501.0
10 15 20 25 30
0
5
10
15
20
25
He
i
g
h
t
(
c
m)
Model
Horn
He
i
g
h
t
(
c
m)
He
i
g
h
t
(
c
m)
He
i
g
h
t
(
c
m)
He
i
g
h
t
(
c
m)
He
i
g
h
t
(
c
m)
He
i
g
h
t
(
c
m)
He
i
g
h
t
(
c
m)
Figure 3.11: Initial density profile for Horn?s
laboratory experiment (red, dashed line) and the
modeled hyperbolic tangent density profile.
Density anomaly
3
Chapter 3. Verification and validation
29
tory data to
l
ped waves.
Phys e
ent;
ing
the le each
.
teepen and develop solitons, however
Horn ing
ed
orted
o
rved by Horn:
the n
ut the
e
billo
is 13
elvin
the
displacement recorded by the central wavegauge in
Horn?s tank is used to compare the labora
the isopycnal displacement in the center of the mode
domain. Figure 3.12 shows the model and Horn?s
results for group A (?3.3.2., Table 3.2). In general,
the model matches the qualitative evolution of wave
phenomena and the period of the waves, but over
predicts the wave amplitude. This is consistent with
comparing an inviscid model to a laboratory
experiment with viscous fluid. There tends to be one
less soliton in the wave series in the model
simulations than in the laboratory experiments; this
will be discussed later in this section.
Scenario 1 is within regime 1 where the
dominating wave behavior is viscous dam
ical viscous damping is not present and th
damping due to numerical dissipation is less than
physical damping. Therefore, it is expected that
modeled waves within regime 1 will have less
damping (i.e. these wave will have slightly larger
?bumps? develop) than in the laboratory experim
this is seen in a comparison of scenario 1 between
the model simulation and Horn?s experiment.
Scenarios 24 are in regime 2 with soliton
development, which the model represents, show
ading wave with the largest amplitude and
successive wave in the train with smaller amplitudes
The model results qualitatively follow Horn?s trend
of more rapid soliton formation with increasing
nonlinearity.
Scenario 5 is in regime 2 and should
theoretically s
observed a broken undular bore in the lead
wave (Figure 3.12). The wave evolution, as
described by the model simulation, shows bore
development (Figure 3.15), but the breaking
described by Horn is not observed. While detail
data on the breaking phenomena were not rep
by Horn, it is likely their spatial scales are too fine t
be captured in the present model. Behind the bore,
the model shows soliton formation, which is also
indicated in Horn?s results. The model simulation
for scenario 5 shows slightly fewer solitons in the
wave series than Horn?s experiment.
The model results for group B (?3.3.2., Table
3.2) show the same characteristics obse
onlinearity of internal wave evolution increases
as the amplitude ratio (a/h) increases and as the depth
ratio (h/H) decreases (Figure 3.13). With increased
nonlinearity, nonhydrostatic pressure effects
increase, which allows the nonhydrostatic model
simulations to depict soliton formation. Witho
nonhydrostatic pressure, the modeled evolution of
the wave is quite different (?4 ff.). Similar to the
results for group A, the model simulations of group
B generally match the qualitative evolution of wav
phenomena and the period of the waves seen in
Horn?s laboratory experiment. The wave amplitude
is overpredicted in the model simulations.
Scenario 9 has a large amplitude ratio (1.0) and
is on the boundary between KelvinHelmholtz
ws (regime 4) and bore formation (regime 5).
The timescale for KelvinHelmholtz formation
s and the timescale for bore formation is exactly one
quarter of the wave period (25 s). Bores and Kelvin
Helmholtz billows form when their respective
timescales are less than one quarter of the wave
period (Horn, et al. 2001). For this scenario, K
Helmholtz billows were observed by Horn, while
model simulation showed the wave degenerate into a
system of higher mode waves (Figure 3.16). Kelvin
Helmholtz billows are a finescale phenomenon,
which require a small grid resolution to model them.
The grid scale (60?73) used may be too coarse
horizontally to capture this event.
Emergence of Solitons
A critical point of c
of solitons. Horn define
omparison is the emergence
s the emergence of solitons
as, ??
s
oliton
the time when the waves are sufficiently well
separated that the depth (measured from the crest of
the leading wave) of the trough between the leading
solitons (measured from the crest of the leading
wave) is 25% of the amplitude of the leading wave?
(Horn, et al., 2001). Horn used three wavegauge
spaced approximately 1.5 m apart. As solitons may
have emerged between wavegauges, Horn reports the
time it took for the wave to move between the
wavegauge that showed solitons and the upstream
wavegauge. Horn?s definition and method of s
emergence was used to identify the emergence of
solitons for the model simulations. The period of the
emergence of solitons for Horn?s experiments and
model simulations is seen in Figures 3.12 and 3.13.
Horn defines the timescale of steepening as:
s
L
T =
??
where L is the basin length, ? is the amplitude and
(3.10)
12
hh3
c
12
2hh
?
?= (3.11)
Wadzuk & Hodges: Hydrostatic and Nonhydrostatic Internal Wave Models
30
Scenario 1: h/H = 0.30, a/h = 0.15
Scenario 2: h/H = 0.30, a/h = 0.30
Scenario 3: h/H = 0.30, a/h = 0.45
Scenario 4: h/H = 0.30, a/h = 0.60
Scenario 5: h/H = 0.30, a/h = 0.90
0 50 100 150 200 250 300 350 400
5
5
?
(c
m
)
?
(c
m
)
5
5
?
(c
m
)
?
(c
m
)
5
5
?
(c
m
)
5
5
?
(c
m
)
5
5
?
(c
m
)
Time (s)
Figure 3.12: Interfacial displacements for Group A at the center of the tank. The figures on the left are
Horn?s experiment (Horn, et al., 2001). The figures on the right are the model simulations. In both sets
of figures, the dotted lines indicate the range in which solitons emerged.
Scenario 8: h/H = 0.20, a/h = 0.45
Scenario 2: h/H = 0.30, a/h = 0.30
Scenario 7: h/H = 0.40, a/h = 0.23
Scenario 6: h/H = 0.50, a/h = 0.18
2
2
?
(c
m)
2
2
?
(c
m
)
2
2
?
(c
m
)
0 50 100 150 200 250 300 350 400
2
2
?
(c
m)
Time (s)
Figure 3.13: Interfacial displacements for Group B at the center of the tank. The figures on the left are
Horn?s experiment (Horn, et al., 2001). The figures on the right are the model simulations. In both sets of
figures, the dotted lines indicate the range in which solitons emerged.
Chapter 3. Verification and validation
31
where h
1
is the upper layer thickness, h
2
is the lower
layer thickness, H is the total depth, c is the wave
speed:
12
hh
cg'
H
= (3.12)
g? is reduced gravity. T ng
was calcu rted in
Table 3.3. The model simulations have soliton
formation slightly earlier than the timescale of
steepening and Horn?s observations (Figures 3.14).
This discr scid
approxim
Numb
he timescale of steepeni
lated for each scenario and are repo
epancy can be attributed to the invi
ation in the model.
er of Solitons
Scenarios 35 had less solitons emerge in the model
simulation than in Horn?s laboratory experiment.
The reason for this is unknown. The number of
solito the
difference in wave am time of soliton
be
.1 ff.).
?73)
ngth (12 m) by 120 grid
cells. The author?s hypothesis that this resolution
would be sufficient is incorrect. There needs to be
more investigation into the lower limit of resolution
needed to capture all of the solitons in the wave
series. However, the grid resolution must be
balanced with viscosity since horizontal grid
refinement decreases the numerical viscosity to the
same order of magnitude or smaller than molecular
viscosity (?4.1 ff). With decreased numerical
Figure 3.14: Comparison of timescale of
steepening with Horn?s observation of the
emergence of solitons (blue lines with circle)
and model simulation observation of the
emergence of solitons (red lines with square).
ns is not a function of viscosity, like
plitude and
emergence. The number of solitons appears to
dependent on the horizontal grid resolution (?4
The grid resolution examined in this section (60
represents the initial wavele
0 50 100 150 200 250 300 350 400 450
0
50
100
150
200
250
300
350
400
450
T
S
(s)
So
l
i
ton em
ergence (s
)
0 50 100 150 200 250 300 350 400
5
5
?
(c
m
)
?
Figure 3.16: Isopycnal displacement for scenario 9 at the center of th
(c
m
)
e tank. h/H = 0.5 and a/h = 1.0.
Table 3.3: Timescale of steepening for scenarios
in internal wave test case (?3.3).
Scenario Ts (s)
1 427
2 213
3 142
7 446
8 124
9
4 107
571
6
1 2 3 4 5 6
0.29
Length (m)
He
i
g
h
t
(
m
)
He
i
g
h
t
(
m
)
He
i
g
h
t
(
m
)
enario 5, bore development at time 60s. Figure 3.15: Sc
Wadzuk & Hodges: Hydrostatic and Nonhydrostatic Internal Wave Models
32
scosity, the modeled wave amplitude and time of
Summary
The manufactured solution (?3.1) verified that the
in
e
solver iterations per timestep needed for
co s
ex on
h hysics
of internal wave evolution.
Nonhydrostatic pressure has a dispersive effect
on internal wave evolution. The nonhydrostatic
model clearly demonstrates the emergence of
solitons for a wave that lies within regime 2. The
major difference between Horn?s laboratory
experiment and the model simulations is the viscous
orn?s experiment are not included in the
model simulations. Horn?s experiment has a large
boundary layer area compared to volume of the
basin, thus viscous effects are prevalent in the
laboratory results. The model simulations did not
model viscosity, so solitons emerged earlier and
soliton amplitudes were larger. However, the basic
characteristics of the wave train were captured by the
model. The model was orrectly
represent the number of solitons in the wave series at
t on
h tested
a e
n of
i
vi
emergence of solitons may become significantly
different than what is predicted by theory and shown
in laboratory experiments.
3.4
effects in H
nonhydrostatic solver is accurate and convergent.
The convergence criterion was established to be
when the L
?
residual norm is decreased by four
orders of magnitude. The internal wave test bas
(?3.2) validated the convergence criterion
established by the manufactured solution and
determined the number of nonhydrostatic pressur
not able to c
he grid resolution used. In summary, the n
ydrostatic model is shown to be valid when
gainst theory and a laboratory experiment. Th
onhydrostatic model captures the physics
nternal wave evolution.
nvergence. Comparison between Horn?
periments and model results validates the n
ydrostatic solver by showing it captures the p
Chapter 4. Comparison of models
33
odel the physics of
(a nonhydrostatic process), yet
van (2004) observed significant
id
c model
nd
ares the laboratory
expe rostatic and nonhydrostatic
model simulations. One laboratoryscale internal
wave ase is examined on different grid meshes to
examine the effect of grid resolution on model skill
for the hydrostatic and nonhydrostatic models. The
lakescale simulation is also analyzed on different
grid meshes. The analytical techniques in ?2.2 are
used to quantitatively evaluate changes in energy,
which is used to compare differences between
hydrostatic and nonhydrostatic models on the
different grid meshes.
4.1 Laboratory Scale Comparison
Results of hydrostatic and nonhydrostatic models
have been compared with the same laboratory
experiments (Horn, et al., 2001) used in ?3.3. Non
hydrostatic pressure is a smallscale effect; thus, it is
advantageous to model a laboratoryscale internal
wave because the smallscale grid (e.g. 0.1 m ?
0.004 m) allows the nonhydrostatic pressure to
affect wave evolution and not be damped out due to
the size of the grid. The hydrostatic and non
hydrostatic models used eight different grid
resolutions to compare model skill for one scenario
(scenario 2, Table 3.2). The grid aspect ratios
(?z/?x) for all resolutions were of O(10
2
) or greater.
4.1.1 Setup
Basin dimensions and the scenarios are identical to
those listed in ?3.3.1 and Table 3.2. Each scenario
was separately simulated using the both hydrostatic
The timestep is dependent on the physics of the
wave and the grid resolution. Due to the
characteristics of the hydrostatic and nonhydrostatic
models, the timestep required to yield a stable
solution is different. For both models, the timestep
selection is based on the CourantFriedrichLewy
(CFL) condition for baroclinic motions. The
hydrostatic model uses the form of the CFL
condition:
Chapter 4. Comparison of Hydrostatic
and Nonhydrostatic models
Hydrostatic models do not m
wave dispersion
Hodges and Dela
wave dispersion in a hydrostatic model which
produces solitonlike formations. Such solitonlike
formations may mimic the expected internal wave
evolution, however this coincidence is a false
positive. Without nonhydrostatic pressure, the
solitonlike formations can only be an artifact of the
numerical method, and is thus dependent on gr
resolution and numerical truncation error. Previous
literature has not examined the effect of grid
resolution on hydrostatic and nonhydrostati
performance. This chapter examines laboratory a
lake scale simulations. The laboratoryscale
simulations examine nine different internal waves
(?3.3) and qualitatively comp
riments to the hyd
c
and nonhydrostatic models (?4.1.2). For all
scenarios, a model grid of 60?73 (?x = 0.1 m and ?z
= 0.004 m) was used. The pycnocline is represented
by five cells (h
py
/?z = 5). Scenario 2 is modeled on
several different grid meshes, as described in Table
4.1, to compare model error (?4.1.3). Each
simulation for the changing grid resolution was run
for ten wave periods as this time allows the internal
wave to evolve and develop into solitons. All
simulations used the inviscid/diffusionless
approximation.
x
t
CFL c
x
?
?
?
# Cells in Zdir
15 29 73
?
o
/?x
30 x x x 60
# Cells in Xdir 60 x x x 120
600 x x 1200
h
py
/?z 125
Table 4.1: Grid meshes for scenario 2 simulations
Those meshes with an ?x? were performed.
.
(1.13)
where ?t is the timestep for a specific, initial internal
wave speed (c) and grid resolution (?x). The CFL
condition is a nondimensional number which limits
the timestep and grid discretization needed to
stabilize the conditionally stable explicit Euler
scheme (Ferziger and Peri?, 2002) used in CWR
ELCOM. The baroclinic CFL condition is the most
restrictive condition, as opposed to the barotropic or
advective CFL condition, for densitystratified flows
(Hodges, 2000). In CWRELCOM, the maximum
allowed baroclinic CFL
x
= ?2 (Hodges, 2000), but
this work uses a conservative CFL
x
condition of 1/3.
The hydrostatic model is dissipative (?4.1.3 ff.), so
strong vertical motions are inhibited. The CFL
x
condition in Equation (1.13) provides a sufficient
timestep for the hydrostatic model. However, the
nonhydrostatic model is generally less dissipative,
especially for finer horizontal grid resolutions
(?4.1.3 ff.), so strong vertical motions may occur.
Wadzuk & Hodges: Hydrostatic and Nonhydrostatic Internal Wave Models
34
hus, the timestep that was appropriate for the
hydrostatic model is not suitable for the non
e
odel
l
FL
T
hydrostatic model (i.e. strong vertical motions caus
instabilities, Figure 4.1) and the timestep must be
smaller. The timestep for the nonhydrostatic m
is determined by the vertical velocity and vertica
grid resolution, yielding a new definition for the C
condition:
z
ac t
CFL
xz
?
?
??
(1.14)
where a is the wave amplitude and ac/?x is a
measure of vertical velocity of a bore. The
maximum CFL
the
theoretical number of solitons (N) that will evolve in
z
condition allowed is 1/3. The
smaller timestep simulation does show a marked
improvement in the internal wave evolution in
nonhydrostatic model (Figure 4.1), and is thus used
for further analysis. The hydrostatic model is not
significantly altered by the different timesteps.
4.1.2 Scenario Simulations
This section uses qualitative comparisons between
the laboratory experiment and the hydrostatic and
nonhydrostatic model simulations. The number of
solitons that emerge in a wave train for the
laboratory experiment and the model simulations is
compared to the theoretical number of solitons that
will evolve. Kao, et al. (1985) defined the
a train as:
M
N1? +
?
(1.15)
where M is defined for a twolayer system in a
model with the Boussinesq approximation:
1
2
12
o22
12
2 hh
hh3
ML
? ??
=?
? ?
? ?
The length of the basin is L, h
(1.16)
and the
1
is the upper layer
thickness, h
2
is the lower layer thickness and ?
o
is the
wave amplitude. The results are separated into the
simulations for group A, group B, scenario 9
number of solitons that emerge in a wave train.
Group A
Group A has the same interface depth ratio
0.3) and varies the a
(h/H =
mplitude of the wave for each
scena
olitons
ic model
ormation. Scenario 5 has the greatest
initial amplitude and forms an undular bore. The
boratory experiment shows
signal, indicating that the bore was initially turbulent
(Horn, et al., 2001). The nonhydrostatic model
iment
he
rio (Figure 4.2). At small amplitudes (scenario
1), there is some small steepening but no s
emerge; both models show similar evolution. For
scenarios with larger amplitudes (regime 2), the
nonlinearity of the wave increases causing the wave
to steepen and evolve into a train of solitons with
decreasing amplitude. The nonhydrostatic model
represents this evolution, while the hydrostat
showed bore f
la a highfrequency
shows the bore formation, but does not show a high
frequency signal or any other indication of
turbulence in the wave evolution. The non
hydrostatic model was able to characterize the
soliton train that formed in the laboratory exper
after the bore was damped by turbulent mixing. T
hydrostatic model developed the initial bore with a
slightly smaller amplitude than the nonhydrostatic
model. As time progresses, the hydrostatic model
retains the bore shape in the leading wave and
develops solitonlike features behind it.
Group B
Group B varies the interface depth ratio for each
scenario (Figure 4.3). For scenarios with larger
depth ratios, that is a thick upper layer, the
nonlinearity of the wave is small. When the depth
ratio is at middepth (h/H = 0.5, scenario 6), the
Figure 4.1: Scenario 6 pycnocline displacement
simulated with a coarse timestep, determined by
CFL , and a fine timestep, determined by CFL .
x z
0 50 100 150 200 250 300 350 400
4
2
0
2
4
Time (s)
?
(c
m
)
coarse ?t
fine ?t
?
(c
m
)
?
(c
m
)
?
(c
m
)
Chapter 4. Comparison of models
Figur l disp he lef erim 2001).
The f ddle hydro ulations. Th d lines in
Horn nd the d.
35
lacements for Group A measured at the center of the tank. The figures on t
are the nonhydrostatic model simulations. The figures on the right are the
nonhydrostatic model indicate the range of time in which solitons emerge
e 4.2: Interfacia
igures in the mi
?s experiment a
t are Horn?s exp
static model sim
ent (Horn, et al.,
e dotte
Scenario 1: h/H = 0.30, a/h = 0.15
2: h/H = 0.30Scenario , a/h = 0.30
3: h/H = 0.30Scenario , a/h = 0.45
4: h/H = 0.30Scenario , a/h = 0.60
0.Scenario 5: h/H = 30, a/h = 0.90
0 50 100 150 200 250 300 350 400
5
5
?
(c
m
)
?
(c
m
)
5
5
?
(c
m
)
?
(c
m
)
5
5
?
(c
m
)
5
5
?
(c
m
)
5
5
?
(c
m
)
Time
5
5
?
(c
m
)
?
(c
m
)
5
5
?
(c
m
)
?
(c
m
)
5
5
?
(c
m
)
?
(c
m
)
?
(c
m
)
5
5
?
(c
m
)
?
(c
m
)
?
(c
m
)
?
(c
m
)
0 50 100 150 200 250 300 350 400
5
5
?
(c
m
)
?
(c
m
)
?
(c
m
)
?
(c
m
)
?
(c
m
)
?
(c
Tim
Horn?s experiment
Nonhy
ydrostatic
drostatic model
H model
m
)
e (s)
Wadzuk & Hodges: Hydrostatic and Nonhydrostatic Internal Wave Models
system has no appreciable nonlinearity and the
hydrostatic and nonhydrostatic model produce
similar results. A thin upper layer (i.e scenario 8)
causes the wave to steepen quickly and results in a
train of solitary waves with larger amplitudes (a ~
2.5 cm) than a thicker upper layer (i.e. scenario 7, a ~
1.7 cm),
36
h is seen in the nonhydrostatic model.
The hydrostatic model is able to predict the onset of
steepenin ut the wave evolution is entirely
unphysic
Scenario
whic
g, b
al.
9
Horn, et al. (2001) observed shear instabilities
developing into KelvinHelmholtz billows in only
one experiment (scenario 9). There is no figure of
laboratory observations of this experiment, and it is
not specified how the wave evolves after the
billowing event, so no direct comparison can be
made between the laboratory experiments and model
simulations. Neither the hydrostatic, nor the non
hydrostatic models were able to capture Kelvin
Helmholtz billows. Horn, et al. (2001) observed two
bores form at both ends of the tank and propagate
towards the center, which the nonhydrostatic model
did simu The nonhydrostatic model shows the
bore pro tion, with a series of solitons behind the
bore (Fig 4.4). The hydrostatic model has neither
bore form n, nor KelvinHelmholtz billows.
Number olitons in Wave Train
late.
paga
ure
atio
of S
In gener e nonhydrostatic model underestimates
the number of observed solitons in the wave series
from the laboratory experiment (?3.3.3). The
theoretical number of solitons, as defined by Kao, et
al. (1985), for each scenario is in Table 4.2. Figure
4.5 shows how the number of solitons in a train for
the model simulations and laboratory experiments
compare with the theoretical number of solitons.
The theoretical number of solitons is always larger
than what was observed in the laboratory experiment
and the model simulations.
Summa
al, th
ry
hyd
rator
s (Fi
hows
ry ex
The non tic mode oduces results of
the labo eriment o differen
scenario s 4.2 and . The n ydrostati
model s tons emerge sooner th n the
laborato ments an
slightly grea plitude
isopycnal dis ement. T
can be attrib to the
application o viscid
approximation as discusse
in ?3.3.3. The hydrostatic model captures the
evolution for weakly nonlinear waves (scenarios 1
and 6), but presents a poor ntation of
nonlinear wave evolution ( 25, 7 and 8;
Figures 4.2 and 4.3). In sc 25, 7 and 8, the
hydrostatic model exactly m the non
hydrostatic model ste causes the
hydrostatic model vel a bore while the
nonhydrostatic model evol o a train of
solitons. Bore dev ment in the hydrostatic model
is typically followed by so onlike
formations, but these are d from observations
in the laboratory experimen nonhydrostatic
model simulations. The hy ic solitonlike
formations are near the cre bore front and
their amplitudes are small c d to the
amplitudes of the soliton tr the laboratory
experiments. A compariso ernal wave
evolution for scenario 2 is ure 4.6. This
confirms observations of h and non
hydrostatic internal wave e om the
isopycnal displacement.
Discussion
represe
scenarios
enarios
atches
epening
op into
ves int
me solit
ifferent
ts and
drostat
st of the
ompare
ains in
n of int
seen in Fig
ydrostatic
volution fr
until
to de
elop
rosta
y exp
gure
soli
peri
ter am
plac
uted
f the in
l repr
quite well f
4.3)
d
in
his
d
the
r the
onh
an i
t
c
Observations of the isopycnal displacement
demonstrate that the nonh model performs
better than hydrostatic mod spect
reproducing the int wa on of the
laboratory experim by . (2001). The
nonhydrostatic model sho fferences in
internal wave evolution am fferent
scenarios, while the hydros l either shows a
damped wave (scenario 1 a bore with small
solitonlike formations (sce 7 and 8). The
solitonlike formations beh re do not
disperse or change in ampl e progresses;
this observation is the same narios where
bore formation occurs. Th ostatic model
shows the train of solitons nd the
leading wave and decrease de with time, as
in the laboratory experimen drostatic
pressu a dispersive pro g, 1972); the
nonh static model?s ab ulate
disper verifies that the static model
captures the overarching ph onlinear wave
evolu
Tab 2: Theoretical number olitons (N) that w p for a 2layer,
Bou sq mode
Simus12345 89
N 6.2 8.3 9.9 11.3 13.6 6 12.7 1.0
ydrostatic
el, with re
ve evoluti
Horn, et al
ws clear di
ong the di
tatic mode
nd 6) or a
nario 26,
ind the bo
itude as tim
for all sce
e nonhydr
disperses behi
in amplitu
t. Nonhy
perty (Lon
ility to sim
nonhydro
ysics of n
ill develo
67
1.0 5.
to
ernal
ents
re is
ydro
sion
tion.
of sle 4.
ssine
lation
l.
Chapter 4. Comparison of models
37
Figure 4.3: Int ial displ figu
The figures in t middle ar ght a
Horn?s experi and the tons
erfac
he
ment
acements for Group B measured at the center of the tank. The
e the nonhydrostatic model simulations. The figures on the ri
nonhydrostatic model indicate the range of time in which soli
res on the left are Horn?s experiment (Horn, et al., 2001).
re the hydrostatic model simulations. The dotted lines in
emerged.
Scenario 8 = 0.2
Scenari = 0. , a/h = 0.30
Scenari 0.23
Scenari 0.18
: h/H
o 2: h/H
o 7: h/H
0, a/h = 0.45
30
= 0.40, a/h =
o 6: h/H = 0.50, a/h =
2
2
?
(c
m)
2
2
?
(c
m
)
2
2
?
(c
m
)
0 50 100 150 200 250 300 350 400
2
2
?
(c
m)
Time (s)
Horn?s experiment ro drostatiNonhyd static model Hy c model
2
2
?
(c
m
)
?
(c
m
)
?
(c
m
)
?
(c
m
)
2
2
?
(c
m
)
?
(c
m
)
?
(c
m
)
2
2
?
(c
m
)
Time
0 50 100 150 200 250 300 350 400
2
2
?
(c
m
)
?
(c
m
)
Wadzuk & Hodges: Hydrostatic and Nonhydrostatic Internal Wave Models
38
0 50 100 150 200 250 300 35
0
2
4
6
8
10
12
14
16
12345678
Simulation
#
o
f
S
o
li
t
o
n
s
Theory
Laboratory
Model
#
o
f
S
o
li
t
o
n
s
Figure 4.5: Number of solitons that develop according to theory (Ka
al., 198 d in the laboratory experiment and model simulations.
1 2 3 4 5 6
0.05
0.15
0.25
Length (m)
He
i
g
h
t
(
m
)
Figure 4.6: a) Hydrostat d b) non rostatic el of sce 2 internal wa on.
1 2 3 4 5 6
0.05
0.15
0.25
Le (m)
He
i
g
h
t
(
m
)
a)
b)
ngth
ve evolutinario modhydic an
5) an
et o,
0 400
5
5
?
(c
m
)
?
(c
m
)
Figure 4.4: Isopycna placement for scenario 9 at the center of the tank. h/H = 0.5 and a/h = 1.0. l dis
Chapter 4. Comparison of models
Wave evolution in scenario 9 includes the
appearance of shear instabilities in a Kelvin
Helmholtz billow (Horn, et al, 2001). While neither
the hydrostatic nor the nonhydrostatic models were
able to capture the billows. However, the non
hydrostatic model was able to simulate the bore
forma
39
Horn, et al. (2
evolves after the KelvinHelmholtz billow collapses,
but it is speculated that the wave will degenerate into
a train of solitons. The nonhydrostatic model shows
a series of solitons develop behind the leading bore.
If the speculation of the wave evolution is correct,
than the nonhydrostatic model may be a reasonable
estimation of the wave evolution after the initial
KelvinHelmholtz billow formation and collapse.
4.1.3 Changing the Grid Resolution
dynamic energy ev
the power spectral
Isopycnal Displacement
tion propagating from the end of the tank.
001) does not specify how the wave
Scenario 2 was examined on several different grids
(Table 4.1) to asses the effect of grid resolution on
the models? performance. Scenario 2 lies within
regime 2 and is expected to have steepening and
soliton formation. The analysis techniques discussed
in ?2.2 are used to quantify differences between the
models and grid resolutions. Specifically, the
different analysis methods used to quantify model
results are the isopycnal displacement, background
potential energy evolution, numerical diffusivity,
olution, numerical viscosity and
density.
The isopycnal displacement of scenario 2 for the
different grid resolutions is seen in Figures 4.7 and
4.8. The wave evolution for a grid refined vertically
(Figure 4.7) shows that the leading wave moves
slightly faster for finer vertical grids. The
amplitudes of all the solitons within the wave train
are smaller for finer vertical grids (~ 1 cm difference
7).
hydrostatic models. Horizontal refinement
substantially changes the pycnocline displacement
(Figure 4.8). Horizontal grid refinement also
develops differences in the wave period. As the grid
is horizontally refined, the hydrostatic and non
hydrostatic models predict a smaller wave period;
that is in Figure 4.8, the finer horizontal grid has an
average wave period of 120 s, while the coarser grid
has an average wave period of 230 s. Furthermore,
the effect of grid scale on the wave influences the
formation of soliton trains. The nonhydrostatic
model has more solitons in the wave train for the
finer grid than the coarser grid. Finally, the leading
soliton has approximately the same wavelength (~1
m) through time and for the different grids in the
nonhydrostatic model. In the hydrostatic model, the
wave evolves into a bore followed by solitonlike
formations; the wavelengths of the solitonlike
formations change in size as the grid is refined.
between the coarsest and finest grids in Figure 4.
This is seen in both the hydrostatic and non
Background Potential Energy
Changes in the background potential energy (E
B
)
represent numerical diffusion in a diffusionless
model (?2.2). The nonhydrostatic and hydrostatic
models show an increase in the background potential
0 100 200 300 400 500 600 700 800
10
10
?
(c
m
)
Time (s)
30x15
30x29
30x73
0 100 200 300 400 500 600 700
10
800
10
Time (s)
?
(c
m
)
Nonhydrostatic
a)
Figure 4.7: Isopycnal displacement of scenario 2 for vertical grid refinement, measured at the center of
basin.
Hydrostatic b)
a) Nonhydrostatic model, b) hydrostatic model.
Wadzuk & Hodges: Hydrostatic and Nonhydrostatic Internal Wave Models
40
at
(Table 4.3). In the
odels, refinement
ows less change in E
73
and 60?73) is about an order of magnitude less than
rse 15 and 60?15).
The system can be considered unphysical when
the normalized background potential energy (Figure
4.9) goes above unity. In the nonhydrostatic model,
the E
B
goes above unity around three wave periods
he
ve
ins physical for only
nt
vertical grid has more effective numerical diffusion
than a fine vertical grid, so grid refinement reduces
the numerical diffusion.
Numerical Diffusivity
energy (E
B
), with the hydrostatic model increasing
a slightly lower rate (Figure 4.9). In both models,
grid refinement in the horizontal direction (Figure
4.9 c and d) shows more change in E
B
over time. In
the nonhydrostatic model, the finest horizontal grid
(600?29) has a change in E
B
of about 92% over nine
wave periods, while the coarsest grid (30?29) has a
change in E
B
by 80% (Table 4.3). The hydrostatic
model has a change in E
B
by 89% for the finest grid
and 77% for the coarsest grid
hydrostatic and nonhydrostatic m
in the vertical direction sh
B
(Figure 4.9 a and b). The finest vertical grid (30?
for the coarsest grid (30?15). Refining the grid
horizontally for this vertical resolution decreases t
time when the E
the coa st vertical grid (30?
B
goes unphysical to 2.2 wave
periods. For the 600?29 grid, the E
B
increases abo
one around five wave periods, while the 30?29 grid
is physical until about nine wave periods. The
hydrostatic model has slightly less diffusion than the
nonhydrostatic model and rema
a short time longer. Conversely, vertical refineme
significantly improves the numerical diffusion in
both models (Figure 4.9); at ten wave periods the E
B
is well below unity for the finest vertical resolution,
irrespective of the horizontal resolution. A coarse
Nonhydrostatic
Hydrostatic
0 50 100 150
The numerical diffusion coefficient (Figure 4.10)
derived in ?2.2 confirms what was observed in the
background potential energy (E
B
); vertical
refinement reduces the system?s diffusivity, while
refinement in the horizontal increases the system?s
diffusivity (Figure 4.10). The hydrostatic model has
lower diffusivities than the nonhydrostatic model
for all grids, except 60?73. As time progresses,
there is an initial growth in the diffusivity, followed
by a continual decline. This trend matches the
behavior of E
B
, which increases rapidly at first and
flattens with time (Figure 4.9). The numerical
diffusivity is orders of magnitude larger than the
usion
is significant in the models and physical diffusion (in
diffusivity of salt, indicating that numerical diff
200 250 300 350 400
10
10
Time (s)
?
(c
m
)
0 50 100 150 200 250 300 350 400
10
10
Time (s)
?
(c
m
)
30x73
60x73
a)
b)
Figure 4.8: Isopycnal displacement of scenario 2 for h
ic mod
Grid
Hydrostatic Nonhydrostatic
30*15 90 91
60*15 79 90
600*15 77 80
30*29 89 92
% Increase
orizontal grid refinement, measured at the center
el. of basin. a) Nonhydrostatic model, b) hydrostat
60*29 78 85
600*29 77 80
Table 4.3: Percent increase in background potential
energy for scenario 2.
Chapter 4. Comparison of models
e 4.9: Background potential energy (E
41
Figur cenario al grid , B) Horizontal grid = 60, C) Vertical grid
= 15, esent hydrosta rostati odel results (h). The difference between E
B
and t
B
i mic ener
B
) of s
D) Vertical grid = 29. Lines repr non
he initial E is normalized by the init al dyna
2 for different grid resolutions. A) Horizont
tic model results (nh), markers represent hyd
gy (E
= 30
c m
D
).
0 2 4 6 8
10
2
10
1
10
0
10
1
t/T
3
3
3
3
3
3
0x15 nh
0x15 hyd
0x29 nh
0x29 hyd
0x73 nh
0x73 hyd
0 2 4 6 8
10
1
10
0
10
1
t/T
60x15 nh
60x15 hyd
60x29 nh
60x29 hyd
60x73 nh
60x73 hyd
a)
b)
[E
B
?
E
B
(0
)
]
/
E
D
(0)
0 2 4 6 8
10
1
10
0
10
1
t/T
30x29 nh
30x29 hyd
60x29 nh
60x29 hyd
600x29 nh
600x29 hyd
d)
[E
B
?
E
B
(0
)
]
/
E
D
(0)
[E
B
?
E
B
(0
)
]
/
E
D
(0)
0 2 4 6 8 10
10
1
10
0
10
1
30x15 nh
30x15 hyd
60x15 nh
60x15 nh
600x15 hyd
600x15 nh
t/T
c)
[E
B
?
E
B
(0)
]
/
E
D
(0
)
Wadzuk & Hodges: Hydrostatic and Nonhydrostatic Internal Wave Models
42
Fig 10: Numerical diffusivity (?) of scenario 2 ifferent grid resolutions. a) Horizontal grid = 3
b) horizontal grid = 60, c) vertical grid = 15, d) vertical grid = 29. Lines represent nonhydrostatic mod
results (nh), markers represent hydrostatic model results (h). ?
s
= 10
9
m
2
/s.
ure 4. for d 0,
el
0 2 4 6 8
10
0
10
1
10
2
10
3
10
4
t/T
?
/
?
s
?
/
?
s
a)
30x15 nh
30x15 hyd
30x29 nh
30x29 hyd
30x73 nh
30x73 hyd
0 2 4 6 8
10
0
10
1
10
2
10
3
10
4
t/T
?
/
?
s
?
/
?
s
30x15 nh
30x15 hyd
60x15 nh
60x15 hyd
600x15 hyd
600x15 nh
c)
0 2 4 6 8
10
0
10
1
10
2
10
3
10
4
t/T
?
/
?
s
?
/
?
s
b)
60x15 nh
60x15 hyd
60x29 nh
60x29 hyd
60x73 nh
60x73 hyd
?
/
?
s
?
/
?
s
0 2 4 6 8
10
0
10
1
10
2
10
3
10
4
30x29 nh
30x29 hyd
60x29 nh
60x29 hyd
600x29 nh
600x29 hyd
t/T
d)
Chapter 4. Comparison of models
43
e absence of turbulence and wave breaking) would
ve a negligible effect on internal wave evolution.
Dynamic Energy
th
ha
Decreases in dynamic energy (E
D
) represent the
dissipation by numerical error (?2.2). The non
hydrostatic model has lower dissipation rates than
the hydrostatic model for all grid resolutions (Figure
4.11). Horizontal grid refinement (Figure 4.11 c and
d) decreases the rate at which E
D
is dissipated. For
the finest horizontal grid resolution (600?29), the
nonhydrostatic model dissipates E
D
by 8% over nine
wave periods, whereas the hydrostatic model
dissipates E
D
by 85%. The coarsest grid (30?29)
dissipates E
D
by 80% in the nonhydrostatic model
and the hydrostatic model dissipates E
D
by 90%
(Table 4.4). Vertical grid refinement (Figure 4.11 a
and b) generally increases the dissipation rate.
However, the 600?29 grid has less dissipation than
the coarser 600?15 grid, which is counter to other
observations of vertical grid refinement (Figure
4.12).
Numerical Viscosity
Numerical viscosity (?2.2) tends to be smaller for
r resolution grids (Figure 4.13); the non
rostatic model produces numerical viscosities
that are smaller than those produced by the
hydrostatic model. Generally, the numerical
viscosity increases and then decreases as time
progresses. The decrease in numerical viscosity with
time corresponds to an overall decrease in dynamic
energy (Figure 4.11). Once the dynamic energy
begins to decrease, the velocity shears that drive the
numerical viscosity are being reduced. That is, any
system with a lower dynamic energy should
generally lose dynamic energy at a lower rate than a
system with a higher dynamic energy, and thus a
lower numerical viscosity is expected, which is seen
in Figures 4.11 a and b and 4.13 a and b. The
numerical viscosity is slightly larger than molecular
viscosity; as the grid is refined, the numerical
viscosity decreases and nears the molecular
viscosity. Numerical viscosities computed for the
hydrostatic model are slightly greater than those for
the nonhydrostatic model, which is expected as the
hydrostatic model has greater numerical dissipation
than the nonhydrostatic model.
Power Spectral Density
fine
hyd
evolution is associated with the wavelength that has
the maximum pycnocline displacement (i.e.
amplitude). For all grid resolutions, the hydrostatic
model had no shift in peak power from the initial
wavelength within ten wave periods (Figure 4.15).
However, there is still some energy transferred to
different wavelength waves. Timeslices of the
power spectral density show as the wavelength
decreases there is energy associated with certain
wavelengths that follow a harmonic pattern (Figure
4.16). The nonhydrostatic model does show a shift
in peak power between four and seven wave periods
depending on the grid resolution (Figure 4.15). Finer
horizontal grids shift the peak power sooner than
coarser horizontal grids. The peak power shifts to
different wavelengths, although the range is small (~
.12 ? 0.17 wavelengths). Once the peak power shift
has taken place, all grid resolutions for the non
hydrostatic model show similar peak
power/wavelength evolution. For example, as seen
in Figure 4.15d, the 600?29 grid shifts the peak
power at t/T = 4 to ?/?
o
= 0.12. At t/T = 4.5, the
600?29 mesh shifts to ?/?
o
= 0.15, which coincides
with the peak power shift of the 60?29 grid. Prior to
the peak power shift, the nonhydrostatic model
shows a harmonic energy shift to ?/?
o
= 0.5 and
some smaller wavelengths (Figure 4.17). However,
after the peak power shift, the peak power is located
at a wavelength ratio (?/?
o
) of 0.2 for all grid
resolutions. The spectral analysis shows that the
peak power shift is to 20% of the initial basinscale
wavelength, which is consistent with the
degeneration of the wave into a train of five solitons
(Figure 4.2). The increase in power at small
wavelengths shows that energy is being directly
transferred from the basinscale wave to the train of
solitons.
Grid
Hydrostatic Nonhydrostatic
30*15 85 77
60*15 85 47
600*15 80 20
30*29 91 80
60*29 90 52
600*29 89 8
% Increase
Table 4.4: Percent decrease in dynamic energy for
scenario 2.
0
The power spectral density (?2.2) quantifies the
energy transfer of different wavelength waves over
time (Figure 4.14). The peak power in a wave
Wadzuk & Hodges: Hydrostatic and Nonhydrostatic Internal Wave Models
44
Figure 4.11: Dy = 60, c) Vertical
grid = 15, d) V l results (h). E
D
is normalized b
namic energy (E
D
) of scenario 2 for different grid resolutions, a) Horizontal grid = 30, b) Horizontal grid
ertical grid = 29. Lines represent nonhydrostatic model results (nh), markers represent hydrostatic mode
y the initial E
D
.
0 1 2 3 4 5 6 7 8 9
10
1
10
0
30x15 n
30x15 h
60x15 n
60x15 h
600x15 h
600x15 n
h
yd
h
yd
yd
h
t/T
0 1 2 3 4 5 6 7 8 9 10
10
1
10
0
t/T
60x15
60x15
60x29
60x29
60x73
60x73
nh
hyd
nh
hyd
nh
hyd
0 1 2 3 4 5 6 7 8 9
10
1
10
0
30x15 nh
30x15 hyd
29 nh
29 hyd
73 nh
73 hyd
30x
30x
30x
30x
t/T
c)
E
D
/ E
D
(0)
a)
b)
E
D
/ E
D
(0)
E
D
/ E
D
(0)
d)
E
D
/ E
D
(0)
0 1 2 3 4 5 6 7 8 9
10
10
1
0
t/T
30x29 nh
30x29 hyd
60x29 nh
60x29 hyd
600x29 nh
600x29 hyd
Chapter 4. Comparison of models
Discussion
45
The key findings of the laboratoryscale model
simulation are: 1) in both the hydrostatic an n
hydrostatic models horizontal grid refinement always
leads to more growth in the background potential
energy and larger numerical diffusivities, indicating
the horizontal grid refinement increases the
numerical diffusion of the system, 2) the hydrostatic
model is more dissipative and somewhat less
diffusive than the nonhydrostatic model and 3) the
nonhydrostatic model is nearly free of numerical
dissipation for a sufficiently refined horizontal grid
and the controlling mechanism for model skill is
numerical diffusion. The increase in background
potential energy with horizontal grid refinement is
opposite of the conventional wisdom that reducing
grid size must improve model skill and reduce error.
This finding corroborates the finding of Hodges and
Delavan (2004) and may be explained by examining
the concept of energy exchange from ?2.2.
Dissipation decreases the dynamic energy, which is
the kinetic and available potential energy. The
available potential energy is the energy that activates
diffusion (when mixing occurs). Thus, if there is less
available potential energy to be mixed, then there is
less possible diffusion. As the grid is horizontally
refined, the dissipation rate decreases. Re ing
merical dissipation allows the backgrou
ntial energy to continue to grow, thereby
increasing numerical diffusion. The nonhydrostatic
model has less dissipation than the hydrostatic model
(Figure 4.11), and thus has more diffusion (Figure
4.9). Vertical grid refinement generally has less
numerical diffusion and more numerical dissipation
in both models, which is tied to bove idea that
the more dynamic energy is dam the less energy
available to be diffused.
The results indicate that nonhydrostatic
model numerical diffusion is s dependent on
the vertical grid resolution (F 9 a and b and
Figure 4.10 a and b) and num on is
strongly dependent on the hor grid (Figure
4.11 c and d and Figure 4.13 c . That is, as a
grid is vertically refined, num iffusion
decreases and horizontal refin decreases
numerical dissipation. Indeed merical
dissipation in fine horizontal such that the
numerical viscosity is on the o molecular
viscosity. It is an expensive t a nearly
numerically dissipative free m .g. a grid with
1200 cells per wavelength; fo sin a horizontal
grid of 600). However, consi urbulent eddy
viscosities where it is suppose ddy viscosity is
100 times the molecular viscosity, all of the grids
examined would have smaller ical viscosities
than the eddy viscosity (Figur . Ev the
eddy viscosity is 10 times the lar viscosity,
several grids (30?73, 60?73, 00?15 and
600?29) would have numerical viscosities less than
the eddy viscosity. This shows that coarser grid
resolutions may still produce l that performs
well with respect to numerical ity when
compared with eddy viscosity peculated that
with numerical viscosities low eddy viscosity,
then the numerical dissipation e 4.11) could be
considered small. Following, l with low
numerical dissipation would h ? skill controlled
by numerical diffusion.
The nonhydrostatic mo ntation of
internal wave evolution closel the
evolution seen in the laborato ent. The
wavelength and amplitude is tely the same
for different resolutions. The odel
also shows a shift in peak pow haracteristic
wavelength of the solitons. T ccurs at all
grid resolutions and is an indi he non
hydrostatic model?s ability to e physics of
internal wave evolution. The model
was not able to always captur of
observed solitons occurring in train. This
difference is most likely due t with
numerical dispersion. Numerical dispersion is
Figure 4.12: Dynamic energy of scenario 2 for the
nonhydrostatic model. Six grids are displayed
where the horizontal refinement (horizontal grid =
30, 60 and 600) is the same for two different vertical
grids (vertical grid = 15 and 29). E
D
is normalized
by the initial E
D
.
0 2 4 6 8
10
1
10
0
D
D
()
30x15
60x15
600x15
600x29
60x29
30x29
t/T
D
D
()
E
D
/ E
D
(0)
the a
ped,
for the
trongly
igure 4.
erical dissipati
izontal
and d)
erical d
ement
, the nu
grids is
rder of
ask to have
odel (e
r this ba
dering t
d the e
numer
e 4.13)
molecu
60?15, 6
a mode
viscos
. It is s
er than
(Figur
a mode
ave its
del?s represe
y matched
ry experim
approxima
nonhydrostatic m
er to the c
his shift o
cation of t
capture th
nonhydrostatic
e the numbers
the wave
o an issue
en if
d no
duc
nd nu
pote
Wadzuk & Hodges: Hydrostatic and Nonhydrostatic Internal Wave Models
46
0 2 4 6 8
10
1
10
0
10
1
10
2
t/T
?
/
?
m
?
/
?
m
30x15 nh
30x15 hyd
30x29 nh
30x29 hyd
30x73 nh
30x73 hyd
0 2 4 6 8
10
1
10
0
10
1
10
2
t/T
?
/
?
m
?
/
?
m
60x15 nh
60x15 hyd
60x29 nh
60x29 hyd
60x73 nh
60x73 hyd
0 2 4 6 8
10
2
10
1
10
0
10
1
10
2
t/T
m
?
/
?
m
?
/
?
30x15 nh
30x15 hyd
60x15 nh
60x15 hyd
600x15 hyd
600x15 nh
0 5 10
2
0
2
4
6
8
10
12
14
16
18
20
t/T
?
/
?
m
?
/
?
m
30x29 nh
30x29 hyd
60x29 nh
60x29 hyd
600x29 nh
600x29 hyd
Figure 4.13: Numerical viscosity (?) of scenario 2
grid = 60, c) vertical grid = 15, d) vertical grid =
markers represent hydrostatic model results (h).
viscosity includes negative values. ?
for
29. L ,
Note,
2/s.
a)
different grids. a) Horizontal grid = 30, b) horizontal
ines represent nonhydrostatic model results (nh)
d) the ordinate is not log scale because the numerical
m
= 10
6
m
b)
c)
d)
Chapter 4. Comparison of models
affected by discretization methods in time, which
were not analyzed in the present research; it remains
a po
10
1
10
0
10
0
normalized wavelength (? / 2L)
n
o
r
m
al
i
z
ed
t
i
m
e
(
t
/
T
)
7 6 5 4 3 2
Figure 4.14: Representative power spectral density for nonhydrostatic model of scenario 2.
The colorbar represents the log
10
relative PSD.
47
interest for future research. The
hydr odel did not show the soliton
development observed in the laboratory experiment.
Analysis of the power spectral density shows the
hydrostatic model retains most of the wave?s energy
within the bore (represented as the initial wavelength
wave in Figure 4.15). Small solitonlike formations
were observed behind the bore, but the peak power is
not shifted to them (Figures 4.15 and 4.16). The
solitonlike formations vary for different horizontal
grid resolutions in both amplitude and wavelength.
The solitonlike formations are considered an artifact
of the model and grid resolution because the soliton
like formations differ significantly from the
laboratory experiment, vary with the grid resolution
and the peak power is not transferred to them. It is
thought that the hydrostatic model will continue to
change the soliton train characteristics with grid
refinement.
4.2 LakeScale Comparison
The hydrostatic and nonhydrostatic models were
app
mod
labo
dev
Zhou
scale problems. Those that do (Marshall, et al.,
1997; Casulli, 1999), do not model internal waves.
Th testing the behavior of a nonhydrostatic
m on a largerscale internal wave has not
previously been done. Similar to ?4.1.3, the grid is
varied and the analytical methods of ?2.2 are used to
quantify the effect of grid resolution on the model?s
skill.
4.2.1 Setup
A 2D closed basin with length (L) of 12500 m and
height (H) of 50 m was chosen as representative of a
mediumsized lake (Delavan, 2003). The initial
wave shape was a cosine and the density field was
constructed with a threelayer hyperbolic tangent
function [Equation (2.34)]. The density difference
across the pycnocline was varied to assess the effect
of density gradient on model performance. The
amplitude ratio (a/h) was 0.3 and the depth ratio
(h/H) was 0.3. The wave used in this simulation has
amplitude and depth ratios that indicate steepening
and soliton formation should occur (i.e. regime 2).
Several different grid meshes were examined, as
listed in Table 4.5. The timestep used was based on
Equation (1.13); this will be discussed further in
int of
ostatic m
us,
odel
lied to a lakescale comparison to demonstrate
el performance on a larger scale than the
ratory experiment. Most of the previously
eloped nonhydrostatic models (e.g. Stansby and
, 1998; Chen, 2003) do not attempt realworld
?4.2.2 ff. Each simulation was conducted over 100
wave periods. The inviscid/diffusionless
approximation was used.
Wadzuk & Hodges: Hydrostatic and Nonhydrostatic Internal Wave Models
48
0 0.5 1
0
1
2
3
4
5
6
7
8
?/?
o
t/
T
t/
T
30x15 nh
30x15 hyd
30x29 nh
30x29 hyd
30x73 nh
30x73 hyd
0 0.2 0.4 0.6 0.8 1
0
1
2
3
4
5
6
7
8
9
?/?
o
t/
T
t/
T
60x15 nh
60x15 hyd
60x29 nh
yd60x29 h
0 0.5 1
0
1
2
3
4
5
6
7
8
9
?/?
o
t/
T
t/
T
30x15 nh
30x15 hyd
60x15 nh
60x15 hyd
600x15 hyd
600x15 nh
0 0.5 1
0
1
2
3
4
5
6
7
8
9
10
?/?
o
t
/
T
t
/
T
30x29 nh
30x29 hyd
60x29 nh
60x29 hyd
600x29 nh
600x29 hyd
Figure 4.15: Wavelength with peak
Horizontal grid = 30, b) horizon
grid = 29. Lines represent nonhydro
represent hydrostatic model results (h
PSD o
tal grid = 6
static s
yd).
f scenario 2 for different grids. a)
0, c) vertical grid = 15, d) vertical
model results (nh), marker
Chapter 4. Comparison of models
49
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
0.45
0.5
PS
D
?/?
o
30x29, t/T=2.5
60x29, t/T=2.5
30x29, t/T=7.5
60x29, t/T=7.5
600x29, t/T=2.4
600x29, t/T=7.3
PS
D
Figure 4.
and t/T =
Figure ulation of scenario 2 at t/T = 2.5
and t/T = 7.5.
16: Typical timeslices of PSD for hydrostatic model simulation of scenario 2 at t/T = 2.5
7.5.
4.
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
17: Typical timeslices of PSD for nonhydrostatic sim
?/?
o
PS
D
30x29, t/T=2.5
60x29, t/T=2.5
30x29, t/T=7.5
60x29, t/T=7.5
600x29, t/T=2.4
600x29, t/T=7.3
PS
D
Wadzuk & Hodges: Hydrostatic and Nonhydrostatic Internal Wave Models
Results
Results are analyzed for nine simulations using the
loss of dynamic energy (Figure 4.19), computed
numerical viscosity (Figure 4.20), change in
background potential energy (Figure 4.21), computed
numerical diffusivity (Figure 4.22) and the peak
power spectral density (Figure 4.23). As a note to
the results, the finest horizontal grid went completely
unstable for all vertical resolutions (81?9, 81?27,
81?81), as seen in the E
B
, E
D
, numerical diffusivity
and
adve
50
viscosity in Figures 4.19  4.22. The vertical
ctive CFL:
advective z
t
CFL w
z
?
?
?
?
(1.17)
was greater than one, indicating the presence of
strong vertical motions. The EulerLagrange Method
is stable for hydrostatic models (Hodges, 2000), but
the stability has not been investigated in this research
or elsewhere in literature for nonhydrostatic models.
Reducing the timestep to adhere to Equation (1.14)
did not mitigate the instability. The addition of
viscosity had no effect on the instability either. This
issue is discussed further in ?4.2.3.
The internal wave initially steepens with small,
solitonlike features developing (Figure 4.18) and is
damped out after ten wave periods. The system is
dissipative, losing an order of magnitude of dynamic
energy in less than ten wave periods for all grid
resolutions (Figure 4.19). The large numerical
dissipation accounts for the damping effect observed
on the wave?s evolution. For vertical and horizontal
grid refinement, there is nearly the same loss in
dynamic energy for the first three wave periods.
After this time, grid refinement causes more dynamic
energy to be dissipated. The most significant
differences in dynamic energy amongst grid
resolutions were due to refinement in the vertical
tory
scale hydrostatic simulations (i.e. the dynamic
energy decreases about an order of magnitude in ten
wave periods). However, in the laboratoryscale
# Cells in Zdir
92781
?
o
/?x
9xxx18
# Cells in Xdir 27 x x x 54
81 nh nh x 160
h
py
/?z 139
grid. The decrease in dynamic energy in the lake
scale simulations is comparable to the labora
Table 4.5: Grid meshes for lake scale comparison.
Those meshes with an ?x? were performed; nh
indicates that only the nonhydrostatic model was
used.
2000 4000 6000 8000 10000 12000
10
20
30
40
Length (m)
He
i
g
h
t
(
m
)
t/T = 0.9
He
i
g
h
t
(
m
)
2000 4000 6000 8000 10000 12000
10
20
30
40
Length (m)
He
i
g
h
t
(
m
)
t/T = 2.9
Figure 4.18: Internal wave evolution for lakescale case at t/T = 0.9 and t/T = 2.9.
Chapter 4. Comparison of models
51
Figure 4.19: Dynamic energy of lakescale
for different grid resolutions. a)Vertical grid
id
rid
= 9, b) vertical gr
81, d)horizontal g
27. Lines represent n
results (nh), markers
model results (h).
= 27, c) vertical grid =
= 9, e) horizontal grid =
onhydrostatic model
represent hydrostatic
10
1
10
0
10
1
10
2
10
3
10
2
10
1
10
0
10
1
E
D
/ E
D
(0
)
t/T
9 x 81 nh
9 x 81 h
27 x 81 nh
27 x 81 h
81 x 81 nh
81 x 81 h
E
D
/ E
D
(0
)
E
D
/ E
D
(0
)
10
1
10
0
10
1
10
2
10
3
10
2
10
1
10
0
10
1
E
D
/ E
D
(0
)
t/T
9 x 9
9 x 9
9 x 2
9 x 2
9 x 8
9 x 8
nh
h
7 nh
7 h
1 nh
1 h
E
D
/ E
D
(0
)
E
D
/ E
D
(0
)
10
1
10
0
10
1
10
2
10
3
10
2
10
1
10
0
10
1
E
D
/ E
D
(0
)
t/T
27 x 9 nh
27 x 9 h
27 x 27 nh
27 x 27 h
27 x 81 nh
27 x 81 h
E
D
/ E
D
(0
)
E
D
/ E
D
(0
)
10
1
10
0
10
1
10
2
10
3
10
2
10
1
10
0
10
1
E
D
/ E
D
(0
)
t/T
9 x
9 x
27 x
27 x
81 x
9 nh
9 h
9 nh
9 h
9 nh
E
D
/ E
D
(0
)
E
D
/ E
D
(0
)
10
1
10
0
10
1
10
2
10
3
10
2
10
1
10
0
10
1
E
D
/
E
D
(0
)
Tt/
9 x 27 nh
9 x 27 h
27 x 27 nh
27 x 27 h
81 x 27 nh
E
D
/
E
D
(0
)
E
D
/
E
D
(0
)
d)
e)
c)
b)
a)
Wadzuk & Hodges: Hydrostatic and Nonhydrostatic Internal Wave Models
nonhydrostatic model simulations had less
numerical dissipation than the hydrostatic model and
horizontal grid refinement dominated changes in
numerical dissipation, which was not seen in the
lakescale simulations. After about ten wave
periods, the numerical viscosity (Figure 4.20)
becomes scattered. At ten wave periods the wave is
damped, so there are small velocities and
subsequently low numerical dissipation. Therefore,
after ten wave periods, small fluctuations in the
dynamic energy are amplified resulting in scatter.
This scatter is irrelevant to the present research as we
are interested in wave evolution before damping. All
grids show the numerical viscosity initially increase
and then decrease after two or three wave periods
(Figure 4.20). This is similar to the laboratoryscale
results (?4.1.3). Vertically refined grids (Figure 4.20
d
52
e) showed lower values for the numerical
viscosity, similar to the laboratoryscale results; there
was not a significant difference in the numerical
viscosity for horizontally refined grids, which was
not seen in the laboratoryscale results. One other
difference between the lakescale and laboratory
scale results for numerical viscosity is that the
numerical viscosity for the lakescale simulations
were over 100 times greater than molecular
viscosity, while the laboratoryscale simulations?
numerical viscosity never exceeded 30 times
molecular viscosity.
The growth in background potential energy is
greater for horizontal grid refinement and less for
vertical grid refinement (Figure 4.21) similar to the
laboratoryscale. Vertical grid refinement dominates
changes in numerical diffusion over different grid
resolutions. The normalized change in E
B
increases
above unity and is considered unphysical around
three wave periods for the coarsest vertical grid. The
finest vertical grid has a change in E
B
from ~ 0.01 to
0.4 over 100 wave periods (Figure 4.21 d and e).
Changes in the horizontal grid have a relatively small
effect on the growth of E
B
. Likewise, the numerical
diffusivities (Figure 4.22) decrease by an order of
magnitude for vertical grid refinement, while
horizontal refinement produces relatively similar
numerical diffusivities (Figure 4.22). This confirms
the vertical grid as the controlling mechanism of
numerical diffusion, which is the same result in the
laboratoryscale simulations (?4.1.3).
The power spectral density (Figure 4.23) shows
a shift in peak power around ten wave periods; the
wavelength with peak power varies with the grid
resolution. A coarse horizontal grid (e.g. Figure
4.23a) shows all power shifting to the same
wavele (?/?
o
~ 0.5) for different vertical
finer horizontal resolution, as the
efined (Figure 4.23 b), the peak
aller wavelengths (?/?
o
= 0.2) and
ack towards larger wavelengths
e progresses. However, unlike the
imulations where the peak power
e wavelengths irrespective of the
e peak power shifts to different
the grid is horizontally refined
the nonhydrostatic model
tical results to the hydrostatic
o grid resolutions that showed
en the hydrostatic and non
ere the 27?9 and 9?27 grids.
Figures 4.21 and 4.22 show that the 27?9 and 9?27
grids had less diffusion in the nonhydrostatic model
than the hydrostatic model. After ten wave periods,
the background potential energy in the non
hydrostatic model converges to the background
potential energy in the hydrostatic model. The
numerical diffusivities in the nonhydrostatic model
are less than the numerical diffusivities for the
hydrostatic model for the entire 100 wave period
simulation.
The density gradient across the pycnocline was
varied (?? = 0.1, 1.0 and 10.0 kg/m
3
) to asses its
effect on the modeling of internal waves. There is no
significant difference between hydrostatic and non
hydrostatic simulations for the identical grid
resolution (Figure 2.24).
4.2.2 Discussion
Long (1972) states that there is a relationship
between nonlinearity and nonhydrostatic pressure; a
sufficiently steep wave has significant non
hydrostatic effects affecting the internal wave
evolution. The wave simulated in this section is
nonlinear and lies within regime 2 (a/h = 0.3, h/H =
0.3). Thus, nonhydrostatic pressure effects should
disperse the steepened wave and have different
numerical error behavior than the hydrostatic model.
However, the results show the nonhydrostatic model
is nearly identical in all grid resolutions to the
hydrostatic model. The similarity between the
hydrostatic and nonhydrostatic models can be
explained by grid resolution. Grid resolution plays a
key role in how the nonhydrostatic pressure affects
internal wave evolution. The model grid can only
represent a single, average nonhydrostatic pressure
for a single cell. Thus, for cells with small grid
aspect ratios (i.e. 9?9: aspect ratio = 4[10
3
]; 9?81:
ngth
resolutions. For a
grid is vertically r
power shifts to sm
gradually moves b
(?/?
o
= 0.5) as tim
laboratoryscale s
shifted to the sam
grid resolution, th
wavelengths when
(Figure 4.23 ce).
Application of
generates nearly iden
model. The only tw
any difference betwe
hydrostatic models w
and
Chapter 4. Comparison of models
53
essure
n
ion (in
his author,
y
ith
ted
stability may be due to
the d
r
l.,
If
only
grid mesh.
Stabi
vant to realworld scale basins.
nt
and
sure
tic
rse linearly steepening initial
wave
odel
model is numerically
e lake
tio
the
is
)]
aspect ratio = 4[10
4
]), the nonhydrostatic pr
effect is muted. The 27?9 grid has the largest grid
aspect ratio (0.012) of all the grids successfully
applied to the lakescale wave. This grid resolutio
did show small differences in numerical diffus
terms of the background potential energy and
numerical diffusivity) between the hydrostatic and
nonhydrostatic models. It is believed by t
if finer grid resolutions with larger aspect ratios had
remained stable (i.e. the 81?9 and 81?27 cases),
these cases would have shown effects from the
inclusion of nonhydrostatic pressure.
The mechanism causing the model instability
for fine horizontal grid resolutions is presently
unknown. The timestep is not controlling vertical
motions, unlike the laboratory case (?4.1.1), so this
is not a simple CFL issue. There are several possible
sources of error for this instability. There may be an
error in the nonhydrostatic model that was only
detected in this situation. The most likely possibilit
is the use of the explicit Euler discretization for the
baroclinic term in CWRELCOM. While explicit
Euler discretizations are unstable, they have been
successfully applied as part of models, such as
CWRELCOM where truncation error serves as
stabilization. Stability is clearly only a problem w
the nonhydrostatic pressure solution as the
hydrostatic model remains stable under all tes
conditions. The issue with
ecrease in numerical dissipation seen in the
laboratoryscale simulation. Another possibility fo
this instability could be a more global problem with
the application of the fractionalstep method for
internal wave evolution. While several other non
hydrostatic models have used the fractionalstep
method, there has been no reported work on grid
refinement and very little work on internal wave
modeling. Most of the known nonhydrostatic
models (Mahadevan, et al. 1996a, b; Marshall, et a
1997; Stansby and Zhou, 1998; Casulli, 1999; Chen,
2003) do not explicitly model internal waves and
present results only for laboratory scale problems.
lake or ocean scale problems are investigated,
the general circulation patterns and freesurface
elevations are reported. Daily and Imberger (2003)
model internal waves, however only on the
laboratory scale. None of the aforementioned
models present results for more than one
lity at fine resolutions in a large domain
remains an open issue for future research.
4.3 Summary
The nonhydrostatic model laboratoryscale
simulations in the present work compared well with
the theory and laboratory experiments of Horn.
Horn, et al. (2001) examined data from several
published field studies and confirmed the
applicability of the theoretical timescales used to
develop a regime diagram. Therefore, it is inferred
that the results of the laboratoryscale model
simulations are rele
The nonhydrostatic model provides a
substantially better representation of internal wave
evolution than the hydrostatic model for laboratory
scales. The features of internal wave developme
(such as soliton formation) are physical in the non
hydrostatic model, while a function of the grid
resolution in the hydrostatic model. There is no
appreciable difference between the hydrostatic and
nonhydrostatic model for weakly nonlinear,
subsequently weakly nonhydrostatic waves. In
cases with a strongly nonlinear, nonhydrostatic
wave on a sufficiently fine grid, such as in the
laboratory simulations, the nonhydrostatic pres
effect is significant. The addition of nonhydrosta
pressure dispe s the non
into a train of solitons, as observed in the
laboratory, and captures the energy cascade through
the wave train. Adding nonhydrostatic pressure
does not eliminate numerical error, but it does
change the characteristics of model error; numerical
diffusion is the mechanism that determines m
skill, whereas the hydrostatic
dissipative.
The largest grid aspect ratio applied to th
scale case for a stable nonhydrostatic simulation is
0.012 (27?9 grid). The largest grid aspect ra
applied to scenario 2 is 2.0 (600?15 grid), and
smallest grid aspect ratio is 0.02 (30?73 grid). It
deduced from the results of the laboratory scale and
lake scale cases that a nonhydrostatic model needs a
sufficiently large grid aspect ratio [?z/?x > O(10
2
to affect internal wave evolution.
Wadzuk & Hodges: Hydrostatic and Nonhydrostatic Internal Wave Models
54
10
0
10
1
10
2
10
2
10
1
10
0
10
1
10
2
10
3
t/T
?
/
?
m
9 x 9 nh
9 x 9 h
27 x 9 nh
27 x 9 h
?
/
?
m
10
0
10
1
10
2
10
2
10
1
10
0
10
1
10
2
10
3
t/T
?
/
?
m
9 x 9 nh
9 x 9 h
9 x 27 nh
9 x 27 h
9 x 81 nh
9 x 81 h
?
/
?
m
?
/
?
m
10
0
10
1
10
2
10
2
10
1
10
0
10
1
10
2
10
3
t/T
?
/
?
m
27 x 9 nh
27 x 9 h
27 x 27 nh
27 x 27 h
27 x 81 nh
27 x 81 h
?
/
?
m
?
/
?
m
10
0
10
1
10
2
10
2
10
1
10
0
10
1
10
2
10
3
t/T
?
/
?
m
9 x 81 nh
9 x 81 h
27 x 81 nh
27 81 h
81 h
x
81 x
?
/
?
m
?
/
?
m
10
0
10
1
10
2
10
2
10
1
10
0
10
1
10
2
10
3
t/T
?
/
?
m
9 x 27 nh
9 x 27 h
27 27 nh
27 h
x
27 x
?
/
?
m
?
/
?
m
a)
d)
Nu a
d
1,
7
l
is
eft
b)
c)
e)
Figure 4.20:
for different gri
9, b) vertical gr
d)horizontal gri
Lines represent
(nh), markers r
results (h). The
hydrostatic cas
81 makes the fi
out.
merical viscosity of lakesc
d resolutions. a)Vertical gri
id = 27, c) vertical grid = 8
d = 9, e) horizontal grid = 2
nonhydrostatic model resu
epresent hydrostatic model
instability in all non
es where the horizontal grid
gures unclear and so were l
le
=
.
ts
Chapter 4. Comparison of models
55
10
0
10
1
10
2
10
1
10
0
10
1
E
B

E
B
(0
)
/
E
D
(0
)
t/T
9 x 9 nh
9 x 9 h
27 x 9 nh
27 x 9 h
81 x 9 nh
E
B

E
B
(0
)
/
E
D
(0
)
E
B

E
B
(0
)
/
E
D
(0
)
10
0
10
1
10
2
10
2
10
1
10
0
10
1
E
B

E
B
(0
)
/
E
D
(0
)
t/T
9 x 27 nh
9 x 27 h
27 x 27 nh
27 x 27 h
81 x 27 nh
E
B

E
B
(0
)
/
E
D
(0
)
E
B

E
B
(0
)
/
E
D
(0
)
10
0
10
1
10
2
10
2
10
1
10
0
10
1
E
B

E
B
(0
) /
E
D
(0
)
t/T
27 x 9 nh
27 x 9 h
27 x 27 nh
27 x 27 h
27 x 81 nh
27 x 81 h
E
B

E
B
(0
) /
E
D
(0
)
10
0
10
1
10
2
10
2
10
1
10
0
10
1
E
B

E
B
(0
)
/
E
D
(0
)
t/T
9 x 9 nh
9 x 9 h
9 x 27 nh
9 x 27 h
9 x 81 nh
9 x 81 h
E
B

E
B
(0
)
/
E
D
(0
)
E
B

E
B
(0
)
/
E
D
(0
)
a)
b)
d)
e)
Figure 4.21: Background potential energy of
lakescale for different grid resolutions.
rtical grid = 9, b) vertical grid = 27, c)
cal grid = 81, d)horizontal grid = 9, e)
zontal grid = 27. Lines represent non
ostatic model results (nh), markers
esent hydrostatic model results (h).
a)Ve
verti
hori
hydr
repr
10
0
10
1
10
2
1.5
1
0.5
0
0.5
1
1.5
2
E
B

E
B
(0
) /
E
D
(0
)
t/T
9 x 81 nh
9 x 81 h
27 x 81 nh
27 x 81 h
81 x 81 nh
81 x 81 h
E
B

E
B
(0
) /
E
D
(0
)
E
B

E
B
(0
) /
E
D
(0
)
c)
Wadzuk & Hodges: Hydrostatic and Nonhydrostatic Internal Wave Models
56
10
0
10
1
10
2
8
6
4
2
0
2
4
6
8
x 10
4
?
/
?
s
9 x 81 nh
9 x 81 h
27 x 81 nh
27 x 81 h
81 x 81 nh
81 x 81 h
t/T
?
/
?
s
?
/
?
s
10
0
10
1
10
2
10
0
10
1
10
2
10
3
10
4
10
5
10
6
?
/
?
s
t/T
27 x 9 nh
27 x 9 h
27 x 27 nh
27 x 27 h
27 x 81 nh
27 x 81 h
?
/
?
s
?
/
?
s
10
0
10
1
10
2
10
0
10
1
10
2
10
3
10
4
10
5
10
6
?
/
?
s
t/T
9 x 27 nh
9 x 27 h
27 x 27 nh
27 x 27 h
81 x 27 nh
?
/
?
s
?
/
?
s
10
0
10
1
10
2
10
0
10
1
10
2
10
3
10
4
10
5
10
6
?
/
?
s
t/T
9 x 9 nh
9 x 9 h
9 x 27 nh
9 x 27 h
9 x 81 nh
9 x 81 h
?
/
?
s
?
/
?
s
10
0
10
1
10
2
10
0
10
1
10
2
10
3
10
4
10
5
10
6
?
/
?
s
t/T
9 x 9 nh
9 x 9 h
27 x 9 nh
27 x 9 h
81 x 9 nh
?
/
?
s
?
/
?
s
a)
b)
c)
e)
d)
Figure 4.22: Numerical diffusivity of lake
scale for different grid resolutions. a)Vertical
grid = 9, b) vertical grid = 27, c) vertical
81, d)horizontal grid = 9, e) horizontal gr
27. Lines represent nonhydrostatic mode
results (nh), markers represent hydros
model results (h).
grid =
id =
l
tatic
Chapter 4. Comparison of models
57
)
vertical grid = 27, e) vertical grid = 81. Lines represent nonhydrostatic model results (nh), markers represent ults (h).
Figure 4.23: Wavelengths with peak PSD of lakescale for different grid resolutions. a)Horizontal grid = 9, b) horizontal grid = 27m c)vertical grid = 9, d
hydrostatic model res
10
2
10
1
10
0
10
1
10
0
10
1
10
2
t/
T
?/?
o
t/
T
t/
T
9x81 nh
9x81 h
27x81 nh
27x81 h
81x81 h
10
2
10
1
10
0
10
1
10
0
10
1
10
2
t/T
?/?
o
t/T t/T
9x27 nh
9x27 h
27x27 nh
27x27 h
10
2
10
1
10
0
10
1
10
0
10
1
10
2
t
/
T
?/?
o
t
/
T
t
/
T
9x9 nh
9x9 h
27x9 nh
27x9 h
10
2
10
1
10
0
10
1
10
0
10
1
10
2
t/T
?/?
o
t/T t/T
27x9 nh
27x9 h
27x27 nh
27x27 h
27x81 nh
27x81 h
10
2
10
1
10
0
10
1
10
0
10
1
10
2
t/T
?/?
o
t/T t/T
9x9 nh
9x9 h
9x27 nh
9x27 h
9x81 nh
9x81 h
a)
b)
c) d)
e)
Wadzuk & Hodges: Hydrostatic and Nonhydrostatic Internal Wave Models
58
10
0
10
1
10
2
10
2
10
1
10
0
10
1
E
B

E
B
(1
) /
E
D
(1
)
t/T
0.1 nh
0.1 h
1.0 nh
1.0 h
10.0 nh
10.0 h
E
B

E
B
(1
) /
E
D
(1
)
E
B

E
B
(1
) /
E
D
(1
)
10
1
10
0
10
1
10
2
10
3
10
2
10
1
10
0
10
1
E
D
/ E
D
(1
)
t/T
0.1 nh
0.1 h
1.0 nh
1.0 h
10.0 nh
10.0 h
E
D
/ E
D
(1
)
E
D
/ E
D
(1
)
10
0
10
1
10
2
10
1
10
2
10
3
10
4
10
5
?
/
?
s
t/T
0.1 nh
0.1 h
1.0 nh
1.0 h
10.0 nh
10.0 h
?
/
?
s
?
/
?
s
0 10 20 30 40 50 60 70 80 90 100
1000
500
0
500
1000
1500
2000
2500
t/T
0.1 nh
0.1 h
1.0 nh
1.0 h
10.0 nh
10.0 h
?
/
?
m
?
/
?
m
?
/
?
m
a)
b)
c)
d)
Figure 4.24: Varying the density gradient of the lakescale simulation. a)Background potential energy, b) numerical diffusivity, c)
dynamic energy, d) numerical viscosity. Lines represent nonhydrostatic model results (nh), markers represent hydrostatic model
results (h).
Chapter 5. Isolation of nonhydrostatic effects
59
effe
In en
oceans, est
wave
issue a
non
neither
nor pr
press
fou
non
simu
beha
identif
press
advan
hydr
using
could
hydr
simp
me
applicatio
devel
rema
prese
of fou
wave
(Fig
The st
beha
nonhy
is used to identify local regions with nonhydrostatic
behavior. The new nonhydrostatic pressure
isolation approach is based on the separation of
hydrostatic parameter obtained from
hydrostatic pressure isolation method terizes a
cell as hydrostatic or nonhydrostatic. chapter
provides the numerical theory used to p the
nonhydrostatic pressure isolation me d an
application of the method for two diff aves.
Further, the nonhydrostatic pressure n
method is applied to one wave on thre rent grid
resolutions to assess the effect of grid ion on
the nonhydrostatic pressure isolation .
5.1 Numerical Theory
The governing equations (see ?2 eld a
vertical momentum equation for invis :
Chapter 5. Isolation of nonhydrostatic
cts
vironments with hydrostatic aspect ratios (e.g.
uaries, lakes), nonhydrostatic pressure
gradients are only locally significant, such as at steep
fronts. Stansby and Zhou (1998) discussed this
nd suggested isolating regions of significant
hydrostatic pressure gradients; however, they
identified regions of nonhydrostatic behavior
ovided a local solution for the nonhydrostatic
ure. The present research builds the
ndations for linking hydrostatic models with local
hydrostatic solutions to yield practical
lations of nonhydrostatic internal wave
vior. This chapter presents a new method for
ying local regions where the nonhydrostatic
ure is expected to be significant. The key
ce in the stateoftheart is that the non
ostatic region is demonstrated to be identifiable
only a hydrostatic solution. Thus, this method
be used to allow the application of a local non
hydrostatic solver where necessary, while
?hydrostatic? areas are more simply computed with a
ostatic model. The ability to ?bootstrap? from
ler hydrostatic methods to nonhydrostatic
thods should allow faster solution than global
n of nonhydrostatic methods. The
opment of a local nonhydrostatic solver
ins an issue for further research (?6 ff.), so the
nt work can be seen as a necessary development
ndational methods.
As a wave propagates, the fluid around the
front is pushed away from the wave front
ure 5.1); this is where the strongest vertical
velocities are located (Daily and Imberger, 2003).
rong vertical velocities cause the wave to
ve more nonlinearly, and subsequently, more
drostatically. Thus, the vertical acceleration
hydrostatic and nonhydrostatic pressure solutions in
the fractionalstep method, so that the regions with
nonhydrostatic behavior are a priori identified from
the hydrostatic solution. The hydrostatic solution
provides the vertical acceleration used in the non
hydrostatic pressure isolation method. A non
the non
charac
This
develo
thod an
erent w
isolatio
e diffe
resolut
method
.1.1) yi
cid flow
nh
o
PDW 1
Dt z
?
=?
? ?
(5.1)
This equation is a relationship betwee al
acceleration and the nonhydrostatic p vertical
gradient, which ordinarily exists only on
hydrostatic model. The vertical veloc is
representative of a flow field that accounts for non
hydrostatic pressure. In a hydrostatic the
nonhydrostatic pressure is neglected uation
(5.1) is not modeled. The hydrostatic l
velocity is computed diagnostically b ing the
continuity equation (see ?2.1.1). The
hydrostatic pressure isolation method e
hydrostatic vertical velocity and the relationship in
Equation (5.1) to approximate the non static
pressure.
n vertic
ressure
for a n
ity (W)
model,
and Eq
vertica
y enforc
non
uses th
hydro
*
nh
o
PDW 1
Dt z
?
=?
? ?
(5.2)
W* is the hydrostatic vertical velocity; this is the
same vertical velocity used in Equation (2.7).
Therefore, the lefthand side of Equat ) is
known after the hydrostatic model is e
righthand side is unknown and is app with
the nonhydrostatic pressure isolation
comparison of the hydrostatic vertical acceleration
(DW
*
/Dt) and the modeled nonhydrostatic pressure
vertical gradient is seen in Figure 5.2. The
hydrostatic vertical acceleration is the approximated
ion (5.2
applied; th
roximated
method. A
Wave Direction
Figure 5.1: Direction of vertical velocities for a
horizontally propagating wave front.
Wadzuk & Hodges: Hydrostatic and Nonhydrostatic Internal Wave Models
60
,
c
e
ion
ak
dients.
The approximated and modeled nonhydrostatic
pressure vertical gradients are symmetrical around
the wave front, although the modeled non
hydrostatic pressure vertical gradient shows a
slightly wider range of influence around the steep
wave front. Aside from the principle wave front,
there is some steepening associated with a higher
mode wave at the right end of the basin in Figure 5.2.
Both the approximated and modeled nonhydrostatic
pressure vertical gradients identify this location.
.
nonhydrostatic pressure vertical gradient. Note
?modeled? refers to anywhere the nonhydrostati
model is used and ?approximated? refers to when th
hydrostatic model is used to make an approximat
of the nonhydrostatic pressure. The location of pe
values for the approximated nonhydrostatic pressure
vertical gradient is a reasonable match for the peak
modeled nonhydrostatic pressure vertical gra
Horizontal flow dominates internal wave
propagation and is prognostically determined in
hydrostatic models, thus it is practical to use the
horizontal momentum equation in the development
of the nonhydrostatic pressure isolation method
The nonhydrostatic horizontal momentum equation
contains the nonhydrostatic pressure horizontal
gradient:
nh
x
?
?
oo
z'
PDU 1 1
gdz
Dt x x
?
??
??? ?
?
=? + ?
? ??
??
?
ut an
f the
nh
accel
h
? ?
(5.3)
The nonhydrostatic pressure horizontal
gradient is not known in a hydrostatic model, b
approximation of the nonhydrostatic pressure
horizontal gradient effect is required to identify
regions of significant nonhydrostatic behavior. The
approximation for the nonhydrostatic vertical
gradient can be used to estimate the nonhydrostatic
pressure horizontal gradient. Vertically integrating
Equation (5.2) produces an approximation o
nonhydrostatic pressure:
Figure 5.2: Comparison of hydrostatic/no
Density profile; b) hydrostatic vertical
vertical gradient. The main nonhydrostatic
also some nonhydrostatic effect on the rig
ydrostatic vertical momentum terms. a)
eration; c) modeled nonhydrostatic pressure
influence is at the steep wave front. There is
t side of b and c.
1 2 3 4 5 6 7 8 9 10
2
4
6
Length (m)
He
i
g
h
t
(
m
)
a) Density
b) DW
*
/Dt
c) dP
nh
/dz
He
i
g
h
t
(
m
)
He
i
g
h
t
(
m
)
He
i
g
h
t
(
m
)
Chapter 5. Isolation of nonhydrostatic effects
61
a) Density
b) Approximate dP
nh
/dx
1 2 3 4 5 6 7 8 9 10
2
4
6
c) Model dP
nh
/dx
He
i
g
h
t
(
m
)
*
nh o
DW
Pdz
Dt
=??
?
(5.4)
The horizontal spatial derivative of the
approximated nonhydrostatic pressure [Equation
(5.4)] is:
*
nh
o
PDW 1
dz
xDt x
??
=?
??
?
(5.5)
Equation (5.5) approximates the nonhydrostatic
pressure horizontal gradient, as seen in Figure 5.3.
Similar to the nonhydrostatic pressure vertical
gradient, the approximated nonhydrostatic pressure
horizontal gradient matches the occurrence of the
modeled nonhydrostatic pressure horizontal gradien
t
at the steep wa
e
tatic effect is large. The screening
parameter relates the nonhydrostatic pressure
orizontal gradient to other terms within th
orizontal momentum equation that influe
internal wave evolution. Where the screening
parameter is large, that location is considered to have
signi
the
Fi
pr t; c) model non
hydrostatic pressure horizontal gradient.
ve front.
A screening parameter is developed to quantify
the relevance of the approximated nonhydrostatic
pressure with regards to its influence on internal
wave evolution, thereby identifying regions wher
the nonhydros
h e
h nce
ficant nonhydrostatic behavior.
Nonhydrostatic pressure effects are linked to
the nonlinearity of a wave and are propagated at
speed of the wave (?6 ff.). As the nonhydrostatic
pressure gradients are strongest at the wave front,
their relative importance on internal wave evolution
is dependent on how the nonhydrostatic pressure
gradient scales with the horizontal acceleration.
Therefore, the screening parameter that isolates
gure 5.3: Comparison of nonhydrostatic pressure horizontal gradient. a) Density
ofile; b) approximate nonhydrostatic pressure horizontal gradien
Length (m)
b) Approxi ate dP
nh
/dx
1 2 3 4 5 6 7 8 9 10
c) odel dP
nh
/dx
He
i
g
h
t
(
m
)
Length (m)
b) Approximate dP
nh
/dx
1 2 3 4 5 6 7 8 9 10
2
4
6
c) Model dP
nh
/dx
He
i
g
h
t
(
m
)
Length (m)
b) Approxi ate dP
nh
/dx
1 2 3 4 5 6 7 8 9 10
c) odel dP
nh
/dx
He
i
g
h
t
(
m
)
Length (m)
b) Approximate dP
nh
/dx
1 2 3 4 5 6 7 8 9 10
2
4
6
c) Model dP
nh
/dx
He
i
g
h
t
(
m
)
Length (m)
b) Approxi ate dP
nh
/dx
1 2 3 4 5 6 7 8 9 10
c) odel dP
nh
/dx
He
i
g
h
t
(
m
)
Length (m)
b) Approximate dP
nh
/dx
1 2 3 4 5 6 7 8 9 10
2
4
6
c) Model dP
nh
/dx
He
i
g
h
t
(
m
)
Length (m)
b) Approxi ate dP
nh
/dx
1 2 3 4 5 6 7 8 9 10
c) Model dP
nh
/dx
He
i
g
h
t
(
m
)
Length (m)
Wadzuk & Hodges: Hydrostatic and Nonhydrostatic Internal Wave Models
62
regions of nonhydrostatic behavior compares the
nonhydrostatic horizontal gradient to the hydrostatic
horizontal acceleration (DU
*
/Dt):
nh
o
**
max,n
x
DU DU
Dt Dt
?
??
??
??
??
??
???
??
(5.6)
The difference between the maximum
horizontal acceleration at a time step (n) and the
horizontal acceleration at each cell is used to
tigate the effect of cells with very small horizontal
accelerations that may amplify the screening
eter erroneously. The total horizontal
acceleration is used to include the horizontal wave
agation and the nonlinear steepening. The total
horizontal acceleration can be decomposed into the
horizontal wave propagation, which is characterized
ocal horizontal acceleration:
P1???
mi
param
prop
by the l
*
U
t
?
?
(5.7)
nonlinear steepening, which is characterized
dvection acceleration:
and the
by the a
**
**
UU
UW
xz??
s of the approximated nonhydrostatic
dient, the screening parameter is:
??
+ (5.8)
In term
pressure gra
*
**
max,n
Dt Dt
???
??
DW
dz
xDt
DU DU
???
??
?
??
??
??
?
(5.9)
Figure 5.4 shows the comparison between the
approximated and modeled screening parameter.
The approximated screening parameter slightly over
predicts nonhydrostatic regions through the water
column, whereas the modeled screening parameter is
located only at the steep wave front.
5.2 Application of the Nonhydrostatic
Pressure Isolation Method
The nonhydrostatic pressure isolation method is
applied to two different internal waves to
demonstrate the method?s ability to isolate regions
with significant nonhydrostatic pressure effects.
5.2.1 Setup
The test basin is a 2D rectangular basin with an
initial internal wave setup as specified in Table 5.1.
The two different internal waves (internal wave 1
and internal wave 2) were initialized as a cosine
wave with a threelayer hyperbolictangent density
profile. The amplitude of internal wave 1 is 1.125 m,
yielding an amplitude ratio (a/h) of 0.5. The
amplitude of internal wave 2 is 0.3 m, yielding a/h =
0.1. Internal wave 1 is steep (a/L = 0.11), while
internal wave 2 is much less steep (a/L = 0.04).
Internal wave 1 was simulated with three
different square grids (20?15, 40?30 and 80?60,
where each grid cell has dimensions of 0.5 m ? 0.5
125 m,
sen to avoid
damping the nonhydrostatic effect over a grid with a
small aspect ratio (?4.3). Internal wave 2 was
simulated
ternal wav
7
2.
0
ess h
py
Density gradient ??
m, 0.25 m ? 0.25 m and 0.125 m ? 0.
respectively). Square grids were cho
with the 40?30 grid.
5.2.2 Results
The nonhydrostatic pressure isolation method is an
effective way to screen for and target areas that
require a nonhydrostatic solution. However, the
s of test basin.
Value Units
Table 5.1: Dimension
In
Length L
Height H
Upper layer depth h
Depth ratio h/H
Pycnocline thickn
Characteristic Symbol
e 1 Internal wave 2
10 10 m
.5 7.5 m
25 3.0 m
.3 0.4 
55
44kg/m
3
Chapter 5. Isolation of nonhydrostatic effects
63
Figure 5.4: Comparison of screening parameter, ?. a) Density profile; b,d) ? calculated with the
approximate nonhydrostatic pressure horizontal gradient; c,e) ? calculated with the model non
hydrostatic pressure horizontal gradient. In b and c, as cells scale from dark blue to dark red,
the cell is considered more nonhydrostatic.
1 2 3 4 5 6 7 8 9 10
2
4
6
Length (m)
He
i
g
h
t
(
m
)
d) Approximate ?
e) Model ?
He
i
g
h
t
(
m
)
He
i
g
h
t
(
m
)
He
i
g
h
t
(
m
)
He
i
g
h
t
(
m
)
He
i
g
h
t
(
m
)
He
i
g
h
t
(
m
)
i
He
i
g
h
t
(
m
)
t
t
He
i
g
h
t
(
m
)
l
( )
He
i
g
h
t
(
m
)
?
2
4
6
He
i
g
h
t
(
m
)
a) Density
c) Model ?
Length (m)
He
i
g
h
t
(
m
)
1 2 3 4 5 6 7 8 9 10
0
2
4
6
b) Approximate ?
He
i
g
h
t
(
m
)
He
i
g
h
t
(
m
)
He
i
g
h
t
(
m
)
He
i
g
h
t
(
m
)
Wadzuk & Hodges: Hydrostatic and Nonhydrostatic Internal Wave Models
64
reening parameter requires a criterion for
scriminating hydrostatic and nonhydrostatic
gions. Where the screening parameter exceeds the
riterion, that cell is considered nonhydrostatic,
hile a cell that has a screening parameter below the
riterion is considered hydrostatic. Quantifying the
riterion as a function of density gradients, wave
haracteristics and grid resolution is beyond the
ope of this work. However, to demonstrate the
avior of the nonhydrostatic pressure isolation
method, three different criteria were chosen (? = 1, 5
and 10) and applied a posteriori to model results.
Figures 5.5 and 5.6 show that as the nonhydrostatic
screening parameter is increased, the region that is
considered nonhydrostatic decreases. Figure 5.6
shows spikes in the number of cells that exceed the
criterion about every 400 timesteps. This time
coincides with when the steep wave fronts reflect off
the basin wall, where the wave has large vertical
velocities and therefore the nonhydrostatic effect is
strong.
Three grid resolutions were applied to internal
wave 1. The finest grid resolution (80?60) identified
the most cells as nonhydrostatic, while the coarsest
grid (20?15) identified the least (Figure 5.7). This
does not imply that there are more nonhydrostatic
regions with the finer grid, but rather there are more
cells representing this region. Figures 5.4, 5.8 and
5.9 show that all grids have similar ability to identify
the region of nonhydrostatic behavior.
Internal wave 1 is strongly nonlinear and so the
nonhydrostatic pressure effect is strong at the wave
front. Internal wave 2 is weakly nonlinear,
nonetheless the nonhydrostatic pressure isolation
method identifies the steepest part of the wave, as
seen in Figure 5.10. Figure 5.11 shows a comparison
between the number of cells with a screening
parameter greater than one for internal wave 1 and
internal wave 2. Internal wave 2 shows large rises in
the number of cells considered nonhydrostatic when
the wave is traveling along the basin and has the
least number of nonhydrostatic cells when the wave
is at the basin walls. When the wave is reflecting off
the basin walls, the vertical and horizontal
accelerations are greater than anywhere else in the
basin, so the screening parameter is able to properly
detect these regions with nonhydrostatic behavior.
When internal wave 2 is traveling along the basin, it
should not behave nonhydrostatically, as it is a near
linear wave. However, the horizontal accelerations
are very small so the screening parameter is
amplified. The filter used in the screening parameter
to eliminate small accelerations does not work in the
limit of very small horizontal accelerations. Thus,
the screening parameter delineates regions as non
hydrostatic erroneously.
5.3 Summary
The new nonhydrostatic pressure isolation method
identifies regions where the nonhydrostatic pressure
effect is significant. For a sufficiently steep wave,
the method is able to predict these regions quite well.
For weakly nonlinear waves, the nonhydrostatic
pressure isolation method tends to overpredict non
hydrostatic regions, due to the nature of a weakly
nonlinear system. The method of the screening
parameter takes the difference between the local
horizontal acceleration and the maximum horizontal
acceleration in the field to ameliorate the problem of
small accelerations. Despite this construction, a
weakly nonlinear system has such small
accelerations that the nonhydrostatic screening
parameter is still amplified. Thus, the present
method is limited by the need to identify a
characteristic horizontal acceleration scale, which
herein is taken as the maximum value of the domain.
If a wave is weakly nonlinear, and therefore
hydrostatic, the nonhydrostatic pressure isolation
method is not applicable. However, if a wave has
significant steepening, then the nonhydrostatic
pressure isolation method identifies nonhydrostatic
regions. The results from the nonhydrostatic
pressure isolation method must not be viewed alone,
but with a general understanding of the nature of the
modeled wave.
sc
di
re
c
w
c
c
c
sc
beh
Chapter 5. Isolation of nonhydrostatic effects
65
0 200 400 600 800 1000 1200 1400 1600
0
200
400
600
800
1000
# Timesteps
# Ce
l
l
s
>
?
? = 1
? = 5
? = 10
Figure 5.6: Number of cells considered nonhydrostatic for internal wave 1; ? > criteria for a
40?30 grid
hydrostatic and
Figure 5.5: Internal wave 1 domains with dif
white are regions considered non
ferent nonhydrostatic screening parameters. The
the black are regions considered hydrostatic.
Leng
De
p
t
h
(m
th (m)
)
1 2 3 4 5 6 7 8 9 10
0
2
4
6
? = 5
? = 10
= 5
? = 1
De
p
t
h
(m
)
Wadzuk & Hodges: Hydrostatic and Nonhydrostatic Internal Wave Models
66
0 200 400 600 800 1000 1200 1400 1600
0
1000
2000
3000
4000
# timesteps
# c
e
l
l
s
?
> 1
40x30
20x15
80x60
# c
e
l
l
s
?
> 1
# c
e
l
l
s
?
> 1
0 200 400 600 800 1000 1200 1400 1600
0
200
400
600
800
1000
# Timesteps
# C
e
l
l
s
?
> 1
0
40x30
20x15
80x60
# C
e
l
l
s
?
> 1
0
# C
e
l
l
s
?
> 1
0
0 200 400 600 800 1000 1200 1400 1600
0
500
1000
1500
# timesteps
#
C
e
lls
?
> 5
40x30
20x20
80x60
#
C
e
lls
?
> 5
#
C
e
lls
?
> 5
#
C
e
lls
?
> 5
#
C
e
lls
?
> 5
#
C
e
lls
?
> 5
#
C
e
lls
?
> 5
Figure 5.7: Number of cells considered nonhydrostatic for internal wave 1. Each panel has three
different grids (blue line is 40?30 grid, red line is 20?15 grid and green line is 80?60 grid). The top
panel is for a screening parameter, ? > 1; the middle panel has ? > 5 and the bottom panel has ? > 10.
Chapter 5. Isolation of nonhydrostatic effects
Figure 5.8: Comparison of nonhydrostatic screening parameter, gamma, on the
20?15 grid for internal wave 1. a) Density profile; b) ? calculated with the
approximate nonhydrostatic pressure horizontal gradient; c) ? calculated with
the model nonhydrostatic pressure horizontal gradient.
1 2 3 4 5 6 7 8 9 10
2
4
6
Length (m)
He
i
g
h
t
(m)
b) Approximate ?
e) Model ?
a) Density
He
i
g
h
t
(m)
67
Wadzuk & Hodges: Hydrostatic and Nonhydrostatic Internal Wave Models
68
th
1 2 3 4 5 6 7 8 9 10
2
4
6
Length (m)
He
i
g
h
t
(m)
c) Model ?
b) Approximate ?
a) Density
He
i
g
h
t
(m)
He
i
g
h
t
(m)
Figure 5.9: Comparison of nonhydrostatic screening parameter, gamma, on the
80?60 grid for internal wave 1. a) Density profile; b) ? calculated with the
approximate nonhydrostatic pressure horizontal gradient; c) ? calculated with
e model nonhydrostatic pressure horizontal gradient.
Chapter 5. Isolation of nonhydros
tatic effects
69
1 2 3 4 5 6 7 8 9 10
2
4
6
Length (m)
He
i
g
h
t
(
m
)
He
i
g
h
t
(
m
)
He
i
g
h
t
(
m
)
He
i
g
h
t
(
m
)
He
i
g
h
t
(
m
)
a) Density
b) Approximate dP
nh
/dx
c) Model dP
nh
/dx
Figure 5.10: Comparison of nonhydrostatic pressure horizontal gradient for internal wave 2. a) Density
profile; b) approximate nonhydrostatic pressure horizontal gradient; c) model nonhydrostatic pressure
horizontal
Figure 5.11: Number of cells considered nonhydrostatic for internal wave 1 (blue line) and internal
wave 2 (red line). ? > 1.
gradient.
0 0.5 1 1.5 2 2.5 3 3.5 4
0
500
1000
1500
t/T
# C
e
l
l
s
?
> 1
# C
e
l
l
s
?
> 1
# C
e
l
l
s
?
> 1
# C
e
l
l
s
?
> 1
internal wave 1
internal wave 2
Wadzuk & Hodges: Hydrostatic and Nonhydrostatic Internal Wave Models
70
Chapter 6. Conclusions and Recommendations
71
Chapter 6. Conclusions and
Recommendations
The two main objectives of the present research are:
? quantify the differences between hydrostatic
and nonhydrostatic simulations of internal
wave evolution, and
? develop a method to a priori determine
regions with nonhydrostatic behavior.
Both of these objectives have been accomplished,
and the results provide insight into the mechanisms
that control internal wave modeling. This chapter: 1)
summarizes the work and conclusions, 2) discusses a
numerical issue that limits the rapidity of the
pressure Poisson solution?s convergence and 3)
provides recommendations for future work.
6.1 Summary Discussion
Differences between the hydrostatic and non
hydrostatic representation of internal waves are
considerable (?4.1). The hydrostatic model
simulates half of the physics necessary to capture
internal wave evolution (i.e. nonlinear acceleration).
This allows the wave to nonlinearly steepen, but
there is no physical term to model dispersion. Thus,
the hydrostatic model is limited in capability to
modeling hydrostatic shallow waves that remain
essentially linear damped waves (?4.1, scenarios 1
and 6). However, the hydrostatic model is
incapable of physically capturing the dispersion of a
steepened wave into a train of solitons, as the
hydrostatic model neglects the physics that control
this process. The solitonlike formations observed in
the hydrostatic model (in this work (?4) and Hodges
and Delavan, 2004) change when the grid resolution
for a hydrostaticallymodeled wave is varied. This
indicates that the size and shape of the solitonlike
formations are dependent on the grid (?4.1.3) and
therefore are a model fabrication. This finding is
further confirmed by the power spectral density
distributed over wavelengths; the peak power is not
shifted to the wavelengths of the solitonlike
formations, and any power that is seen in these
waves is a result of harmonics (?4.1.3). Numerical
error within the hydrostatic model accumulates to
reduce the energy in the wave and eventually damp
out the wave?s steepness. Numerical diffusion can
be reduced by vertically refining the grid, but the
hydrostatic system is numerically dissipative,
irrespective of the grid resolution.
Unlike the hydrostatic model, the non
hydrostatic model includes the nonhydrostatic
pressure term, which acts as a dispersive force.
Thus, as a wave steepens nonlinearly, the non
hydrostatic pressure applies an effect that disperses
the wave into a train of solitons. Grid resolution
affects the occurrence of the train of solitons (i.e.
finer horizontal grids simulate more emerging soliton
trains than a coarser grid, Figure 4.8). However, the
characteristic wavelength of the solitons is relatively
unaffected by grid resolution, compared to the
hydrostatic model. The peak power in a wave
system shifts to the characteristic soliton wavelength
after solitons have emerged; the shift in power is
indicative of an energy transfer from basinscale
wavelengths to soliton wavelengths. It is thus
concluded that the nonhydrostatic model is indeed
modeling the physics of internal wave evolution.
Nonhydrostatic models do not eliminate numerical
error, however, the characteristics of numerical error
are different in a nonhydrostatic model than in a
hydrostatic model. Like the hydrostatic model,
vertical grid refinement reduces numerical diffusion.
In contrast to the hydrostatic model, horizontal grid
refinement greatly reduces numerical dissipation in
the laboratoryscale case. The laboratoryscale case
(?4.1.3) decreases dynamic energy by only 8% over
nine wave periods for the finest grid resolution
compared to the coarsest grid resolution which
decreases dynamic energy by 85%. It appears that
the neglect of nonhydrostatic pressure effectively
acts as a dissipative term in the hydrostatic
equations. Perhaps an alternative viewpoint is that
the nonhydrostatic pressure redistribution of
momentum reduces the numerical dissipation that is
associated with shallow water conservation of
momentum.
Grid resolution and the subsequent grid aspect
ratio are the two model parameters that control the
nonhydrostatic model?s ability to capture internal
wave evolution. Most of the lakescale cases had
very small grid aspect ratios [O(10
3
) ? O(10
4
)], and
showed no significant difference between the
hydrostatic and nonhydrostatic models. The one
lakescale case nonhydrostatic model that differed
(in diffusion only) from the hydrostatic model had a
grid aspect ratio of O(10
2
); the laboratory scale case,
where results were significantly different between
the hydrostatic and nonhydrostatic models, had grid
aspect ratios of at least O(10
2
). Due to the
successful use of the nonhydrostatic model for
larger grid aspect ratios, the present research
suggests that the grid aspect ratio must be at least of
O(10
2
) to include the effect of nonhydrostatic
Wadzuk & Hodges: Hydrostatic and Nonhydrostatic Internal Wave Models
72
ial resolution
error non
poral resolution must also
condition.
he
and is pr
gradients ve. The
anal ca son
equation
nce
ves
d
For example, flow
.g. Su, et al., 1999)
e step
tatic
the
me
e
id needed to
effectively apply the nonhydrostatic model renders
the s ture
w
c
ler
e
n
th
n
he
a
n

er
nt
r
2
)
pressure. In addition to the fine spat
needed for a physical, low numerical
hydrostatic model, the tem
be fine enough to maintain a small vertical CFL
T pressure Poisson solution is inherently stiff
oblematic for the motion of stiff pressure
at the front of an internal wa
yti l solution (?3.1) provides a pressure Pois
that is independent of a flow field.
However, the pressure Poisson solution for
hydrodynamic problems (e.g. internal wave
evolution) depends upon the gradient of velocity.
The gradient of velocity generally does not match the
speed of propagation of the nonhydrostatic pressure
gradient. Indeed, this appears to be a key differe
between nonhydrostatic solutions of internal wa
and flows that have been more commonly studie
using nonhydrostatic models.
over a backwards facing step (e
has an unsteady recirculation region behind th
due to nonhydrostatic pressure. The nonhydros
pressure gradients in the recirculation oscillate at
same timescale as the motion of the reattachment
point, which is much longer than the timescale of
advection. Thus, the advection velocity limits the
model time step, from which it follows that local
changes in the nonhydrostatic pressure gradient
from one time step to another are small. As the ti
?n+1? pressure field is a good initial guess at the tim
?n? pressure field in a backwardsfacing step, a
Poisson solver should be able to rapidly converge.
In contrast, when modeling an internal wave, for
common subcritical flows the wave is faster than the
fluid velocity, so the baroclinic wave speed limits the
model time step. As the internal wave propagates
rapidly through the model grid, nonhydrostatic
pressure gradients at the wave front are propagated at
the wave speed. Therefore, the nonhydrostatic
pressure field may change quite dramatically
between time steps, such that the time ?n? pressure
field is a poor approximation of the time ?n+1?
pressure field. This effect is demonstrated by
simulating a monochromatic wave (?3.2) under two
different conditions for the Poisson solution: 1) with
the time ?n? pressure is used as the starting point for
the time ?n+1? pressure solution; and 2) with zero
pressure used as the starting point for all pressure
solutions. Figure 3.9 shows that using the time ?n?
pressure always requires a larger number of
iterations to meet the convergence criterion.
The fine timestep and spatial gr
olution of realworld scale problems a fu
capability. The present research investigated a ne
method to identify regions with significant non
hydrostatic behavior. Using the nonhydrostati
pressure isolation method, a model may be
developed which would ?bootstrap? from a simp
hydrostatic model to a nonhydrostatic model in th
regions identified as ?nonhydrostatic.? Such a
hydrostatic/nonhydrostatic hybrid model should
allow faster solution of internal wave evolution.
However, there is a limit to the ability of the no
hydrostatic pressure isolation method. The method
overpredicts nonhydrostatic regions in shallow
waves. These shallow waves are ably modeled wi
the hydrostatic model, since steepening and solito
formation is not the dominant phenomenon in t
wave?s evolution. Thus, the nonhydrostatic
pressure isolation method should only be applied to
waves which are expected to steepen and form
solitons (i.e. regime 2, Figure 3.10).
It is crucial to water quality modeling to have
hydrodynamic model to skillfully simulate the
internal wave evolution and the flow field in a basin
(Gross, et al, 1999). For damped linear internal
waves, a hydrostatic model may adequately model
the behavior as long as viscous damping is modeled.
For internal waves where nonlinear steepening and
soliton formation are the prevalent features, a no
hydrostatic model provides a physical simulation of
internal wave evolution on a sufficiently fine grid,
whereas the hydrostatic model does not. Thus, non
hydrostatic models should be used when nonlinear
internal waves and solitons are present in the flow
field; this will in turn create a more realistic
representation of all the system processes that
depend on the flow field, such as those in a wat
quality model.
6.2 Conclusions
The conclusions that may be drawn from the prese
work are:
? The present hydrostatic model of the Eule
equations shows solitonlike formations for
nonlinear/nonhydrostatic waves, which are
an artifact of the numerical model and grid
resolution.
? There is some evidence that the grid aspect
ratio (?z/?x) must be greater than O(10
to capture the effect of nonhydrostatic
pressure in the present nonhydrostatic
model.
Chapter 6. Conclusions and Recommendations
73
y
mendations
st be
ulation
n
d
g
eed,
ce
c
hydr
the
gion to
r
e
on
e handled by
buffe
to
t
e
c
o finer grids. Additionally, the
essure solver should be
conditioned conjugate gradient
Ack o
This rese
Naval R
2001 of
011057 tre for Water Research at the
Uni
ELC
research t
of Texas
? A hydrostatic model can be used to identif
internalwave regions with nonhydrostatic
behavior.
6.3 Recom
The research presented here investigates emerging
areas within nonhydrostatic modeling. This work
began to chip away at some of the issues challenging
nonhydrostatic modeling, but much remains to be
investigated. The most notable issue that mu
examined in nonhydrostatic modeling is sim
time. To diminish simulation times, this author
believes there are two main approaches: 1) address
the issue of stiffness, and 2) construct a hybrid model
that uses the nonhydrostatic pressure isolation
method.
The pressure Poisson equation is inherently
stiff, but the stiffness is exacerbated by the advectio
methods used. Presently, the velocity at the previous
timestep is used in the source term of the pressure
Poisson equation; this velocity moves much slower
than the wave and the nonhydrostatic pressure
gradients. One suggestion to improve the source
term of the pressure Poisson equation would be to
use a predicted advection velocity field from the
previous two timesteps. For instance, use a secon
order AdamsBashforth which predicts the ?n+1?
velocity from the ?n1? and ?n? velocities. Using the
past two timesteps may better estimate the strong
velocity gradients near the propagating wave front.
A second approach to improve the source term
in the pressure Poisson equation is to use the non
hydrostatic pressure isolation method to identify the
point at an instant in time with the peak non
hydrostatic pressure gradient. This point should
coincide with the propagating wave front and the
strongest velocity gradients. Identifying the location
with peak nonhydrostatic pressure gradients, alon
with the known timestep, will yield the velocity that
of the peak nonhydrostatic pressure gradient. This
velocity should be characteristic of the wave sp
which may be used to provide a better estimation of
the nonhydrostatic pressure field.
The nonhydrostatic pressure isolation method
developed here only provides the relative importan
of the nonhydrostatic pressure effect in the internal
wave?s evolution. In application, the nonhydrostati
screening parameter needs some criterion to
delineate a region as ?hydrostatic? or ?non
ostatic.? For example, where the screening
parameter is above a preset criterion, a cell is
designated as ?nonhydrostatic.? At this cell, the
nonhydrostatic solver may be called to calculate
nonhydrostatic pressure and the updated velocity
field and freesurface. The criterion for a re
be considered nonhydrostatic must be established.
The author does not know if there will be a global
value for the nonhydrostatic screening paramete
criterion, or if it will be dependent on each wave
characteristics.
Another issue associated with the
implementation of the nonhydrostatic pressure
isolation method is that any change in the velocity
field propagates instantaneously over the entire
domain. Thus, changes within the interior of the
domain may cause perturbations on the freesurfac
that could feed back into the flow in the form of non
conservation in mass and volume. There must be
some transition between the ?hydrostatic? and n
hydrostatic? cells; the transition may b
r cells. The buffer cells will account for any
freesurface perturbation and make a correction
the system to ensure conservation of mass and
volume. The size of the buffer region must also be
investigated to determine if it is a function of the
magnitude of the nonhydrostatic pressure gradients
or the size of the nonhydrostatic region.
Finally, the nonhydrostatic model must be
examined on more realworld scale cases with
different grid resolutions to determine if the
instability found in the lakescale case (?4.2) is a
model error or an issue with the fractionalstep
method. To verify the present research?s assessmen
that a grid aspect ratio of at least O(10
2
) must b
used to observe differences between the hydrostati
and nonhydrostatic model, the lakescale case
should be applied t
present nonhydrostatic pr
adapted to use a pre
method or a multigrid method to improve
convergence time.
n wledgements
arch was supported under the Office of
esearch Young Investigator Program Award
Ben R. Hodges, ONR Contract No. N00014
4. The Cen
versity of Western Australia provided the
OM numerical model used in this work under a
partnership agreement with the Universi y
at Austin.
Wadzuk & Hodges: Hydrostatic and Nonhydrostatic Internal Wave Models
74
Bibliography
75
Bibliography
Apel, J.R., (1981). Satellite sensing of ocean
surface dynamics. Ann. Rev, Earth & Planetary Sci.,
8:303342.
Armfield, S. and R. Street, (1999). The
fractionalstep method for the NavierStokes
equations: the accuracy of three variations. J. Comp.
Phys., 153:2:660665.
Batchelor, G.K., (1967). An Introduction to
Fluid Dynamics. Cambridge Univeristy Press, New
York, 615 pgs.
Blumberg, A.F., & G.L. Mellor, (1987). A
description of a threedimensional coastal ocean
circulation model, in ThreeDimensional Coastal
Ocean Models, N. S. Heaps (ed.), American
Geophysical Union, Washington D.C. pp. 116.
Boegman, L., J. Imberger, J.N. Ivey, & J.P.
Antenucci, (2003). High frequency internal waves in
large stratified lakes. Limnol. & Oceanogr.,
48:12:895919.
Casulli, V., (1999). A semiimplicit finite
difference method for nonhydrostatic, freesurface
flows. Int. J. for Numer. Methods in Fluids, 30:425
440.
Casulli, V. & G.S. Stelling, (1998). Numerical
simulation of 3D quasihydrostatic freesurface
flows. J. Hydraul. Engng. 124:678686.
Chen, X., (2003). A fully hydrodynamic model
for threedimensional freesurface flows. Int. J. for
Numer. Methods in Fluids, 42:929952.
CushmanRoisin, B., (1994). Introduction to
Geophysical Fluid Dynamics. Prentice Hall, New
Jersey, 320 pgs.
Daily, C. & J. Imberger, (2003). Modeling
solitons under the hydrostatic and Boussinesq
approximations. Int. J. for Numer. Methods in
Fluids, 43:231252
Delavan, S.K., (2003). Evolution Accumulation
of Numerical Errors in Hydrostatic Models of
Internal Waves, M.S. Thesis, Department of Civil
Engineering, University of Texas, Austin, Aug,
2003.
De Silva, I.P.D., J. Imberger, & G.N. Ivey,
(1997). Localized mixing due to a breaking internal
wave ray at a sloping bed. J. Fluid Mech. 350:127.
Ezer, T., H, Arango, A.F. Shchepetkin (2002).
Developments in terrainfollowing ocean models:
intercomparisons of numerical aspects. Ocean
Modelling, 4:249267
Farmer, D.M., (1978). Observations of long
nonlinear internal waves in a lake. J. Phys.
Oceanogr., 8:6373.
Ferziger, J.H. & M. Peri?, (2002).
Computational Methods for Fluid Mechanics.
Springer, Berlin, Germany, 423 pgs.
Fringer, O. B., & R. L. Street. 2003. The
dynamics of breaking progressive interfacial waves.
J. Fluid Mech. 494: 319353.
Garrett, C. and W. Munk, (1979). Internal
waves in the ocean. Ann. Review Fluid Mech.,
11:339369.
Gill, A.E., (1982). AtmosphereOcean
Dynamics, Academic Press: New York, 1982.
Gross, E.S., J.R. Koseff, & S.G. Monismith,
(1999). Threedimensional salinity simulations of
South San Francisco Bay, J. Hyd. Engr.,
125:11:11991209.
Haidvogel, D.B., H.G. Arango, K. Hedstrom,
A. Beckmann, P. MalanotteRizzoli, A.F.
Shchepetkin, (2000). Model evaluation experiments
in the North Atlantic Basin: simulations in nonlinear
terrainfollowing coordinates. Dynamics of
Atmospheres and Oceans. 32:239?281.
Hamrick, J.M., (1992). A threedimensional
environmental fluid dynamics computer code:
Theoretical and computational aspects. The College
of William and Mary, Virginia Institute of Marine
Science, Special Report, 317, 63 pgs.
He, P., M. Salcudean, I.S. Gartshore, and P.
Nowak, (1996). Multigrid calculation of fluid flows
in complex 3D geometries using curvilinear grids.
Computers and Fluids, 25:4:395419.
Hirsch, C., (1988). Numerical Computation of
Internal and External Flows, Volume 1:
Fundamentals of Numerical Discretization. John
Wiley & Sons, Inc., New York, 515 pgs.
Hodges, B.R., (2000). Numerical Techniques in
CWRELCOM (code release v.1). CWR Manuscript
WP 1422 BH, Centre Water Res., Nedlands, Western
Australia, Australia.
Hodges, B.R., (2004). Accuracy of semi
implicit shallowwater CrankNicolson
discretization. ASCE J. Eng. Mech., 130:8:904910.
Hodges, B.R. & S.K. Delavan, (2004).
Numerical diffusion and dissipation in hydrostatic
Wadzuk & Hodges: Hydrostatic and Nonhydrostatic Internal Wave Models
76
waves. ASCE Engineering
ence Proceedings, June 1316,
2004
00). Modeling basinscale internal
wave
gescale interfacial gravity
wave
2). A weakly nonlinear model of
long h.,
Imbe
d.],
paradigm of planetary problems.
Elsev
77491.
ion
396:183
201.
c Process,
Acad
Geophysical Fluid Flows,
Acad
).
g and breaking over a slope.
J. Flu
d
ress, New York, 730 pgs.
ts with a
ydraul. Engrgr., 129:215224.
Dover
ges, and R.
Stock s:
ic
rpolation. Computer Meth. in Appl.
Mech
cons
8), Waves in Fluids,
Cam
dal flow
ng,
. Street,
(199 l.
(199
Adcroft,
(199
Z. (2001). Dynamics of
Inter
onto
models of internal
Mechanics Confer
.
Hodges, B.R., J. Imberger, A. Saggio and K.
Winters, (20
s in a stratified lake, Limnology and
Oceanography, 45:7:16031620.
Horn, D.A., J. Imberger, & G.N. Ivey, (2001).
The degeneration of lar
s in lakes. J. Fluid Mech., 434:181207.
Horn, D. A., J. Imberger, G. N. Ivey, and L. G.
Redekopp. (200
internal waves in closed basins. J. Fluid Mec
467:269287.
Hutter, K., G. Bauer, Y. Wang & P. Guting,
(1998). Forced motion response in enclosed lakes. In
Physical Processes in Lakes and Oceans (ed. J.
rger), AGU.
Imberger, J., (1994). Transport processes in
lakes: A review, p. 79193. In R. Margalef [e
Limnology now: A
ier Science.
Imberger, J. & G.N. Ivey, (1993). Boundary
mixing in stratified reservoirs. J. Fluid Mech.,
248:4
Javam A., J. Imberger & S.W. Armfield,
(1999). Numerical study of internal wave reflect
from sloping boundaries. J. Fluid Mech.,
Kantha, L.H. & C.A. Clayson, (2000a).
Numerical Models of Oceans and Oceani
emic Press, New York, 940 pgs.
Kantha, L.H. & C.A. Clayson, (2000b). Small
Scale Processes in
emic Press, New York, 888 pgs.
Kao, T.W., F. Pan, & D. Renouard, (1985
Internal solitons on the pycnocline: generation,
propagation, and shoalin
id Mech., 159:1953.
Kim, J. and P. Moin, (1985) Application of a
fractionalstep method to incompressibleNavier
Stokes equations. J. Comp. Phys., 59:308323.
Kinsman, B., (1965). Wind Waves: Their
Generation and Propagation on the Ocean Surface.
Dover Publications, Inc., New York, 676 pgs.
Kreyszig, E, (1999). Advanced Engineering
Mathematics. John Wiley & Sons, Inc., New York,
1156 pgs.
Kundu, P.K. and I.M. Cohen, (2002). Flui
Mechanics, Academic P
Laval, B., B.R. Hodges, & J. Imberger,
(2003a). Reducing numerical diffusion effec
pycnocline filter. J. H
Lamb, H., (1932). Hydrodynamics.
Publications, Inc., New York, 738 pgs.
Laval, B., J. Imberger, B.R. Hod
er. (2003b). Modeling Circulation in Lake
Spatial and Temporal Variations. Limnology and
Oceanography 48:3:983994.
Leonard, B.P. (1979). Stable and accurate
convective modeling procedure based on quadrat
upstream inte
. & Engr., 19:1:5998.
Leonard, B.P. (1991). The ultimate
ervative difference scheme applied to unsteady
onedimensional advection. Computer Meth. in
Appl. Mech. & Engr., 88:1774.
Lighthill, J. (197
bridge University Press, Cambridge, 504 pgs.
Lin, B. and R.A. Falconer, (1997). Ti
and transport modeling using ULTIMATE
QUICKEST scheme. J. Hyd. Engr., 123:4:303314.
Long, R.R., (1972). The steepening of lo
internal waves. Tellus, 24:8899.
Lorenz, E.N., (1955). Available potential
energy and the maintenance of the general
circulation. Tellus, 7:157167.
Mahadevan, A., J. Oliger, & R. L
6a). A nonhydrostatic mesoscale ocean mode
Part I: Wellposedness and scaling. J. Phys.
Oceanogr., 26:18681880.
Mahadevan, A., J. Oliger, & R. L. Street,
6b). A nonhydrostatic mesoscale ocean model.
Part II: Numerical implementation. J. Phys.
Oceanogr., 26:18811900.
Marshall, J., C. Hill, L. Perlman, & A.
7). Hydrostatic, quasihydrostatic, and
nonhydrostatic ocean modeling, J. Geophys. Res.,
102:57335752.
Miropol?sky, Yu.
nal Gravity Waves in the Ocean, Kluwer,
Dordrecht, 406 pgs.
Moore, G. (1965). Cramming more components
integrated circuits. Electronics, 38:8.
Bibliography
77
and the respective biogeochemical processes
in th
. Internal
nce, 208:4443:451
460.
ng.
124:1
autics, 41:4:595604.
sing Toolbox User?s Guide,
(199
lems. Int. J. Numer. Methods in
Fluid
).
simulation by the gaskinetic
schem
). A coastal ocean model intercomparison study
for a
fluids.
nonli

ng
18, 2003.
astal Modeling:
Proc rence,
ar
.
86). The
lake.
y & E.A.
mixi
ith
nt. J.
5.
Nishri, A., J. Imberger, W. Eckert, I.
Ostrovsky, & Y. Geifman, (2000). The physical
regime
e lower water mass of Lake Kinneret, Limnol.
Oceanogr., 45:4:972981.
Osborne A.R. and T.L. Burch, (1980)
solitons in the Andaman sea. Scie
Roache, P. J., (2002). Code verification by the
method of manufactured solutions. J. Fluids E
:410.
Roy, C.J., (2003). Grid convergence error
analysis for mixedorder numerical schemes. Amer.
Inst. Aeronautics Astron
Segur, H. & J.L. Hammack, (1982). Soliton
models of long internal waves. J. Fluid Mech.,
118:285304.
Signal Proces
8). The MathWorks, Inc.
Stansby, P.K. & J.G. Zhou, (1998). Shallow
water flow solver with nonhydrostatic pressure: 2D
vertical plane prob
s, 28:542563.
Su, M. D., K. Xu, and M. S. Ghidaoui. (1999
Lowspeed flow
e. J. Comput. Phys., 150:1:1739.
Tartinville, B., E. Deleersnijder, P. Lazure, R.
Proctor, K.G. Ruddick, and R.E. Uittenbogaard,
(1998
threedimensional idealized test case. App.
Math. Modelling., 22:3:165182.
Turner, J. S. (1973). Buoyancy effects in
Cambridge Univ. Press, NY USA.
van Haren, H., (2004). Some observations of
nearly modified internal wave spectra. J.
Geophys. Res., 109(C03045).
Wadzuk, B.M. & B.R. Hodges, (2003).
Comparing hydrostatic and nonhydrostatic Navier
Stokes models of internal waves. ASCE Engineeri
Mechanics Conference, July 16
Weilbeer, H. and J.A. Jankowski, (2000). A
threedimensional nonhydrostatic model for free
surface flows ? development, verification and
limitations. Estuarine and Co
eedings of the 6th International Confe
Nov. 35, 1999.
Whitham, G.B., (1974). Linear and nonline
waves, Wiley, NY, USA, 636 pgs
Wiegand, R.C. and E.C. Carmack, (19
climatology of internal waves in a deep temperate
J. Geophys. Res., 91:C3:39513958.
Winters, K.B., P.N. Lombard, J.J. Rile
D?Asaro, (1995). Available potential energy and
ng in densitystratified flows. J. Fluid Mech.,
289:115128.
Yuan, H. & C.H. Wu, (2004). A two
dimensional vertical nonhydrostatic ? model w
an implicit method for freesurface flows. I
Numer. Meth. Fluids, 44:81183