UFL/COEL2000/014
A THREEDIMENSIONAL CONSERVATIVE
EULERIANLAGRANGIAN MODEL FOR COASTAL
AND ESTUARINE CIRCULATION
by
Jun Lee
Thesis
2000
Copyright 2000
by
Jun Lee
ACKNOWLEDGMENTS
First of all, I would like to thank my advisor, Dr. Sheng,
for his guidance and financial support throughout my graduate
research. He encouraged me when I was frustrated doing my
research. In addition, I would like to thank my committee
members, Dr. Dean and Dr. Thieke, for their careful review of
this thesis.
I would like to thank all of colleagues in my research
group and in particular : Justin and Sun, for their help with
discussing, developing, debugging and verifying the model; Dave
and Du, for doing field work together; Chenxia, Vadim and Jack,
for studying and working together; and my Korean friends Kijin
and Yeonshik, for studying and relaxing together.
I would like to thank everyone in the administrative
offices, Becky, Lucy, and Sandra. Special thanks are sent to
Helen who helped me to find many references. Grateful thanks go
to everyone in the Coastal Engineering Laboratory and
particularly Sidney and Vic who are my diving buddies.
Finally, I would like to thank my parents, my wife and my
son, and God, for without their love and endless support, I
would never be where I am.
TABLE OF CONTENTS
. . . . . . . i i
ACKNOWLEDGMENTS
LIST OF TABLES
LIST OF FIGURES
ABSTRACT . .
1 INTRODUCTION
1.1 Objectives . . . . . . . . .
1.2 Organization of This Study . . . . .
2 GOVERNING EQUATIONS AND BOUNDARY CONDITIONS IN THREE
DIMENSIONAL CARTESIAN COORDINATES . . . . .
2.1 Introduction . . . . . . . . .
2.2 Governing Equations in Cartesian Coordinates
2.3 Boundary Conditions in Cartesian Coordinates
3 A THREEDIMENSIONAL NUMERICAL MODEL . . . . .
3.1 Introduction . . . . . . . . .
3.2 Finite Difference Form of Momentum Equations
3.3 Solution Algorithm . . . . . . .
3.4 Conjugate Gradient Method . . . . .
3.5 Salinity Scheme . . . . . . . .
. . V
. vi
. xvi
5
. . 3
. . 3
. . 5
4 EULERIANLAGRANGIAN (SEMILAGRANGIAN) METHOD . . .. 34
4.1 Introduction . . . . . . . . . 34
4.2 Review of EulerianLagrangian Method (ELM) . 34
4.3 Derivation of EulerianLagrangian Method (ELM) 37
4.4 Conservative EulerianLagrangian Method (ELM) 43
4.5 Numerical Experiments . . . ... . . . 51
4.5.1 Linear Advection Case . . . . .. 51
4.5.2 Conservation Test . . . . . . 59
4.6 Momentum Advection Equation . . . . .. 61
4.6.1 NonLinear Advection Case . . . .. 61
4.6.2 Momentum Conservation Check . . . .. 70
4.7 TwoDimensional Conservation Test . . . .. 71
4.8 ThreeDimensional Conservation Test . . .. 76
4.9 Summary . . . . . . . . . . 81
5 MODEL VERIFICATION WITH ANALYTICAL SOLUTION . . .. 82
5.1 Introduction . . . . . . . . . 82
5.2 Wind Setup . . . . . . . . . . 83
5.3 Tidal Propagation . . . . ... . . . 86
5.4 Tidal Propagation with NonLinear Advection . 98
5.5 Tidal propagation with Coriolis Effect . . .. 104
5.5.1 Kelvin Wave Propagation . . . . .. 104
5.5.2 The Governing Equation for Kelvin Waves . 105
5.5.3 Comparison between Numerical and Theoretical
Solutions . . . . . . . . 109
5.6 Tidal Propagation with Bottom Friction Effect . 117
5.7 DensityDriven Circulation . . . . . . 125
5.7.1 Horizontal Diffusion Test .. . ... 125
5.7.2 Density and Tidally Driven Circulation 129
5.8 Wetting and Drying Test over Tidal Flats . .. 132
5.8.1 Determining Flooding and Drying Cells . 133
5.8.2 Comparison with Analytical Solution . .. 135
6 COMPARISON WITH CH3D . . . . . . . .. . 139
6.1 Introduction . . . . . . . . . 139
6.2 Comparison with CH3D in a Rectangular Basin . 140
7 APPLICATION OF MODEL TO LAKE OKEECHOBEE . . . .. 150
8 CONCLUSIONS AND FUTURE STUDY . . . . . ... 161
8.1 Conclusions . . . . . . . . ... 161
8.2 Future Study . . . . . . . . . 162
APPENDIX A NUMERICAL SCHEMES FOR LINEAR ADVECTION EQUATION 165
APPENDIX B DESCRIPTION OF PROCEDURES OF THE PRESENT MODEL 171
LIST OF REFERENCES . . . . . . . . . . 175
BIOGRAPHICAL SKETCH . . ... . . . . . . 182
LIST OF TABLES
Table page
Table 5.1 Summary of test for analytical solutions and model
simulations . . . . . . . . . . 83
Table 5.2 Comparison between analytical and simulated wind
setup. . . . ... . . . . . . .. 84
LIST OF FIGURES
Figure page
Figure 3.1 Schematic diagram of stencil . . . .. 15
Figure 4.1 Characteristic lines and their feet . .. 35
Figure 4.2 An Eulerian Lagrangian Mesh . . . .. 38
Figure 4.3 Analytical solution of linear advection equation
at T = 50 Sec . . . . . . . . . . 52
Figure 4.4 Analytical solution of linear advection equation
at T = 250 Sec . . . . . . . . . 52
Figure 4.5 Analytical solution of linear advection equation
at T = 500 Sec . . . . . . . . . 52
Figure 4.6 Results of the 1D linear advection equation
obtained with upwind method at T = 50 Sec . . .. 53
Figure 4.7 Results of the 1D linear advection equation
obtained with upwind method at T = 250 Sec . .. 53
Figure 4.8 Results of the lD linear advection equation
obtained with upwind method at T = 500 Sec . .. 53
Figure 4.9 Results of the 1D linear advection equation
obtained with the Leith's method at T = 50 Sec . 54
Figure 4.10 Results of the 1D linear advection equation
obtained with the Leith's method at T = 250 Sec . 54
Figure 4.11 Results of the 1D linear advection equation
obtained with the Leith's method at T = 500 Sec . 54
Figure 4.12 Results of the 1D linear advection equation
obtained with the TVD method at T = 50 Sec . . .. 55
Figure 4.13 Results of the 1D linear advection equation
obtained with the TVD method at T = 250 Sec. . . 55
Figure 4.14 Results of the 1D linear advection equation
obtained with the TVD method at T = 500 Sec. . . 55
Figure 4.15 Results of the 1D linear advection equation
obtained with the ELM with linear interpolation at
T = 55 sec . . . . . . . . . . 56
Figure 4.16 Results of the 1D linear advection equation
obtained with the ELM with linear interpolation at
T = 253 sec . . . . . . . . . . 56
Figure 4.17 Results of the lD linear advection equation
obtained with the ELM with linear interpolation at
T = 506 Sec . . . . . . . . . . 56
Figure 4.18 Results of the 1D linear advection equation
obtained with the ELM with linear interpolation at
T = 55 Sec . . . . . . . . . . 57
Figure 4.19 Results of the 1D linear advection equation
obtained with the ELM with linear interpolation at
T = 253 Sec . . . . . . . . . . 57
Figure 4.20 Results of the lD linear advection equation
obtained with the ELM with linear interpolation at
T = 506 Sec . . . . . . . . . . 57
Figure 4.21 Results of the 1D linear advection equation
obtained with the Conservative QMSL at T = 55 Sec . 58
Figure 4.22 Results of the 1D linear advection equation
obtained with the Conservative QMSL at T = 253 Sec 58
Figure 4.23 Results of the lD linear advection equation
obtained with the Conservative QMSL at T = 506 Sec 58
Figure 4.24 Computational cell used in conservation test 59
Figure 4.25 Conservation error of the conventional ELM 61
Figure 4.26 Conservation error of the QMSL . . .. 61
Figure 4.27 Results of the nonconservative upwind method
at T = 0.5 sec . . . . . . . . . 63
Figure 4.28 Results of the nonconservative upwind method
at T = 1.0 sec . . . . . . . . . 63
vii
Figure 4.29 Results of the nonconservative upwind method
at T = 1.5 sec . . . . . . . . . 63
Figure 4.30 Results of the conservative upwind method
at T = 0.5 sec . . . . . . . . . 64
Figure 4.31 Results of the conservative upwind method
at T = 1.0 sec . . . . . . . . . 64
Figure 4.32 Results of the conservative upwind method
at T = 1.5 sec . . . . . . . . . 64
Figure 4.33 Results of the Leith's method for nonlinear
advection equation at T = 0.5 sec . . . . .. 65
Figure 4.34 Results of the Leith's method for nonlinear
advection equation at T = 1.0 sec . . . . .. 65
Figure 4.35 Results of the Leith's method for nonlinear
advection equation at T = 1.5 sec . . . . .. 65
Figure 4.36 Results of the Godunov method for nonlinear
advection equation at T = 0.5 sec . . . . .. 66
Figure 4.37 Results of the Godunov method for nonlinear
advection equation at T = 1.0 sec . . . . .. 66
Figure 4.38 Results of the Godunov method for nonlinear
advection equation at T = 1.5 sec . . . . .. 66
Figure 4.39 Results of the nonlinear advection equation of
ELM with linear interpolation at T = 0.5 sec . .. 67
Figure 4.40 Results of the nonlinear advection equation of
ELM with linear interpolation at T = 1.0 sec . .. 67
Figure 4.41 Results of the nonlinear advection equation of
ELM with linear interpolation at T = 3.25 sec . 67
Figure 4.42 Results of the nonlinear advection equation of
ELM with quadratic interpolation at T = 0.5 sec . 68
Figure 4.43 Results of the nonlinear advection equation of
ELM with quadratic interpolation at T = 1.0 sec . 68
Figure 4.44 Results of the nonlinear advection equation of
ELM with quadratic interpolation at T = 3.25 sec . 68
viii
Figure 4.45 Results of the nonlinear advection equation of
conservative QMSL method at T = 0.5 sec . . .. 69
Figure 4.46 Results of the nonlinear advection equation of
conservative QMSL method at T = 1.0 sec . . .. 69
Figure 4.47 Results of the nonlinear advection equation of
conservative QMSL method at T = 3.25 sec . . .. 69
Figure 4.48 Conservation error of the conservative QMSL 70
Figure 4.49 Numerical solution of 2D nonconservative ELM 72
Figure 4.50 Conservation error of 2D nonconservative ELM 72
Figure 4.51 Contour map of numerical solution of non
conservative ELM . . . . . . . . . 73
Figure 4.52 Numerical solution of 2D conservative ELM 74
Figure 4.53 Conservation error of 2D conservative ELM 74
Figure 4.54 Contour map of numerical solution of non
conservative ELM . . . . . . . . . 75
Figure 4.55 Numerical solution of 3D nonconservative ELM 77
Figure 4.56 Conservation error of 3D nonconservative ELM 77
Figure 4.57 Numerical solution of 3D nonconservative ELM in
XZ plane . . . . . . . . . . 78
Figure 4.58 Numerical solution of 3D nonconservative ELM in
YZ plane . . . . . . . . . . 78
Figure 4.59 Numerical solution of 3D conservative ELM . 79
Figure 4.60 Conservation error of 3D conservative ELM . 79
Figure 4.61 Numerical solution of 3D conservative ELM in XZ
plane . . . . . . . . . . ... 80
Figure 4.62 Numerical solution of 3D conservative ELM in YZ
plane . . . . . . . . . . . . 80
Figure 5.1 Comparison between analytical and numerical
solutions of water level for wind setup at three locations
within the domain(2D case) . . . . . .. 85
Figure 5.2 Comparison between analytical and numerical
solutions of water level for wind setup at three locations
within the domain(3D case) . . . . . .. 85
Figure 5.3 Schematic diagram of a rectangular basin . 87
Figure 5.4 Computational grid system for tidal propagation
test . . . . . . . . . . ... 88
Figure 5.5 Twodimensional simulation of water surface
elevation of tidal propagation at three locations with
8= 1.
Line : Analytical Solution
Circle : Numerical Solution . . . . . .. 90
Figure 5.6 Twodimensional simulation of water velocity
of tidal propagation at three locations with 8= 1.
Line : Analytical Solution
Circle : Numerical Solution . . . . . .. 91
Figure 5.7 Threedimensional simulation of water surface
elevation of tidal propagation at three locations with
8= 1.
Line : Analytical Solution
Circle : Numerical Solution . . . . . .. 92
Figure 5.8 Threedimensional simulation of water velocity
of tidal propagation at three locations with 0= 1.
Line : Analytical Solution
Circle : Numerical Solution . . . . . .. 93
Figure 5.9 Twodimensional simulation of water surface
elevation of tidal propagation at three locations with
8= 0.501.
Line : Analytical Solution
Circle : Numerical Solution . . . . . .. 94
Figure 5.10 Twodimensional simulation of water velocity
of tidal propagation at three locations with = 0.501.
Line : Analytical Solution
Circle : Numerical Solution . . . . . .. 95
Figure 5.11 Threedimensional simulation of water surface
elevation of tidal propagation three locations with
8= 0.501.
Line : Analytical Solution
Circle : Numerical Solution
Figure 5.12 Threedimensional simulation of water velocity
of tidal propagation at three locations with = 0.501.
Line : Analytical Solution
Circle : Numerical Solution . . . . . .. 97
Figure 5.13 Twodimensional simulation of water surface
elevation of nonlinear advection at three locations with
0= 1 .
Line : Analytical Solution
Circle : Numerical Solution . . . . . .. 100
Figure 5.14 Twodimensional simulation of water velocity
of nonlinear advection at three locations with 0 =1
Line : Analytical Solution
Circle : Numerical Solution . . . . . .. 101
Figure 5.15 Threedimensional simulation of water surface
elevation of nonlinear advection at three locations with
0= 1 .
Line : Analytical Solution
Circle : Numerical Solution . . . . . .. 102
Figure 5.16 Twodimensional simulation of water velocity
of nonlinear advection at three locations with 0 =1
Line : Analytical Solution
Circle : Numerical Solution . . . . . .. 103
Figure 5.17 Twodimensional model results of water surface
elevation at three locations with Coriolis effect when
0= 1 .
Line : Analytical Solution
Circle : Numerical Solution . . . . . .. 112
Figure 5.18 Twodimensional model results of water velocity in
X direction at three locations with Coriolis effect when
0=1.
Line : Analytical Solution
Circle : Numerical Solution . . . . . .. 113
Figure 5.19 Threedimensional model results of water surface
elevation at three locations with Coriolis effect when
0=1.
Line : Analytical Solution
Circle : Numerical Solution
Figure 5.20 Threedimensional model results of water velocity
in X direction at three locations with Coriolis effect when
0=1.
Line : Analytical Solution
Circle : Numerical Solution
. . . . . . . 115
Figure 5.21(a) water level(contours) and velocity field vector
after 48 hours of simulation with Coriolis effect when
0=1.
Contour lines represent water surface elevation . 116
Figure 5.21(b) water level(contours) and velocity field vector
after 60 hours of simulation with Coriolis effect when
0=1.
Contour lines represent water surface elevation . 116
Figure 5.22 Twodimensional model results of water surface
elevation at three locations with bottom friction when
S= 1.
Line : Analytical Solution
Circle : Numerical Solution
Figure 5.23 Twodimensional model results of water velocity
in X direction at three locations with bottom friction when
0=1.
Line : Analytical Solution
Circle : Numerical Solution . . . . . .. 122
Figure 5.24 Threedimensional model results of water surface
elevation at three locations with bottom friction when
S=1.
Line : Analytical Solution
Circle : Numerical Solution
Figure 5.25 Threedimensional model results of water velocity
in X direction at three locations with bottom friction when
0=1.
Line : Analytical Solution
Circle : Numerical Solution
Figure 5.26 Comparison between analytical and numerical
solutions of diffusion equation in X direction
Figure 5.27 Comparison between analytical and numerical
. . . . . . 124
127
xii
. 114
. . . 121
. . . . . . 123
solutions of diffusion equation in Y direction . 127
Figure 5.28 Numerical result of temperature distribution in
rectangular plate . . . . . . . . . 128
Figure 5.29 Comparison between the numerical and analytical
solution of the density driven currents.
Line : Analytical solution
Circle : Numerical solution . . . . . .. 131
Figure 5.30 Definition of total computational depth . 134
Figure 5.31 Wave propagation on a linearly sloping beach 135
Figure 5.32 Nondimensional comparison between wave profiles
as predicted by theory and numerical model of Wetting and
drying(time=0 1 /2) . . . . . . . . 137
Figure 5.33 Nondimensional comparison between wave profiles
as predicted by theory and numerical model of Wetting and
drying(time = 2z/3 g ) . . . . . . . 138
Figure 6.1 Comparison to CH3D with wind forcing . . 140
Figure 6.2 Water surface elevation of propagation test with
CH3D when 6 =1.0
Line : CH3D
Circle : Present model . . . . . . . . 141
Figure 6.3 Water surface elevation of propagation test with
CH3D when =0.501
Line : CH3D
Circle : Present model . . . . . . . . 142
Figure 6.4 Depthaveraged water velocity of propagation test
with CH3D when =0.501
Line : CH3D
Circle : Present model . . . . . . . . 143
Figure 6.5 Water surface elevation of nonlinear test
with CH3D when 0 =1.0
Line : CH3D
Circle : Present model . . . . . . . . 144
Figure 6.6 Water surface elevation of Nonlinear test
with CH3D when 0 =0.501
Line : CH3D
xiii
Circle : Present model
Figure 6.7 Depthaveraged water velocity of nonlinear test
with CH3D when 0=0.501
Line : CH3D
Circle : Present model . . . . . . . . 146
Figure 6.8 The test of Coriolis effect with CH3D when 0=1.0
Line : CH3D
Circle : Present model . . . . . . . 147
Figure 6.9 Comparison of water surface elevation with
Coriolis effect with CH3D when =0.501
Line : CH3D
Circle : Present model . . . . . . .
Figure 6.10 comparison of depthaveraged velocity with
Coriolis effect with CH3D when 0=0.501
Line : CH3D
Circle : Present model . . . . . . . . 149
Figure 7.1 The Bathymetry of Lake Okeechobee
(From the 1987 contour map of Lake Okeechobee) . . 151
. 148
Figure 7.2 Computational grid
. . . . . 151
Figure 7.3 Numerical solution of total water depth as a
function of time at location A, B, C, D, E, and F in Lake
Okeechobee subject to an impulsive wind . .. . . 153
Figure 7.4 Numerical solution of water surface elevation
with time at location A, B, C, D, E, and F . . 154
Figure 7.5 Velocity fields and dry regions at time 1 and 3
Hours. Blank cells indicate the dry cells. . . . 155
Figure 7.6 Velocity fields and dry regions at time 6 and 9
Hours. Blank cells indicate the dry cells . . . 156
Figure 7.7 Velocity fields and dry regions at time 12 and 15
Hours. Blank cells indicate the dry cells . . . 157
Figure 7.8 Velocity fields and dry regions at time 18 and 24
Hours. Blank cells indicate the dry cells. .. . . 158
Figure 7.9 Velocity fields and dry regions at time 30 and 36
Hours. Blank cells indicate the dry cells. . . . 159
xiv
. 145
Figure 7.10 Mass conservation error . ..
. . . . 160
Abstract of Thesis presented to the Graduate School of the
University of Florida in Partial Fulfillment of the
Requirements for the Degree of Master of Science
A THREEDIMENSIONAL CONSERVATIVE EULERIANLAGRANGIAN MODEL
FOR COASTAL AND ESTUARINE CIRCULATION
By
Jun Lee
August 2000
Chairperson: Y. Peter Sheng
Major Department: Civil and Coastal Engineering
This study presents a threedimensional conservative
Eulerian Lagrangian hydrodynamic model with wetting and drying
which is significantly enhanced over the nonconservative
EulerianLagrangian model developed by Casulli and Cheng.
Semiimplicit finite difference method has been used for the
numerical solution of threedimensional shallow water
equations. The advection terms and horizontal diffusion terms
in the momentum equations are discretized with the explicit
conservative EulerianLagrangian Method (ELM). The barotropic
pressure gradient terms in the momentum equations and the
velocities in the vertically integrated continuity equation
are discretized with the 8scheme, while the baroclinic
xvi
pressure gradient terms in the momentum equations are
discretized with the explicit scheme. The vertical diffusion
terms are discretized with the implicit numerical scheme.
Conservation is an important property when studying the
transport of salt, heat, and other conservative species in
coastal and estuarine numerical modeling. If a numerical model
does not have conservative property, the model is unable to
capture a front correctly and the total amount of conservative
substances will not remain unchanged. The conventional
EulerianLagrangian numerical schemes do not have the
conservative property and hence cannot be used for practical
longterm simulations. The present model possesses the global
conservative property because of the use of the quasi
conservative EulerianLagrangian method.
The Coriolis terms in the momentum equations are
discretized explicitly in the X direction and implicitly in
the Y direction, as opposed to the fully explicit treatment in
Casulli and Cheng. This method is conditionally stable but the
time increment is much bigger (about 3 hours) than the other
time step limits. If the Coriolis terms are not considered,
the present model is composed of a large fivediagonal system
of linear equations. The matrix is symmetric and positive
definite. This fivediagonal system of equations is solved
efficiently by the preconditioned conjugate gradient method.
xvii
If the Coriolis terms are considered, the ninediagonal
matrix is positive definite but not symmetric, and is solved
by the successive overrelaxation (SOR) method.
In addition this study compares several advective schemes
versus analytical solutions. "Wetting and drying" features,
which are important for coastal modeling, have been
implemented in the model and compared with analytical solution
as well. Several additional analytical tests were conducted to
demonstrate the accuracy of the model.
After successful comparisons with analytical solutions,
the model is applied to the winddriven circulation of Lake
Okeechobee, Florida. The results demonstrate that the present
model represents the wetting and drying dynamics well and
possesses the conservative property for water mass in Lake
Okeechobee.
xviii
CHAPTER 1
INTRODUCTION
A large number of numerical models of circulation in
coastal, ocean, estuaries, and lakes have been developed over
the last few decades. The earlier models are generally one
dimensional or twodimensional. With the advancement in
computer technology, more recent models are threedimensional
and timedependent.
One pioneering threedimensional coastal and estuarine
model was developed by Leendertse (1967, 1975). The model uses
a timeexplicit method to solve the threedimensional
hydrostatic equation of motion in Cartesian coordinates. Sheng
(1983) developed a threedimensional model of coastal and
estuarine circulation in horizontally Cartesian and vertically
stretched(G grid) coordinates. They solved the external
mode(water level and vertically integrated velocities which
govern the surface gravity waves) using a factorized implicit
scheme and the internal mode(vertical structure of three
dimensional velocities and temperature and salinity) using a
verticallyimplicit scheme.
Blumberg and Mellor(1983, 1987) developed a three
2
dimensional estuarine and coastal Ocean circulation model
using both Cartesian and orthogonal curvilinear coordinates
in the horizontal plane and U grid in the vertical direction.
They solved the external mode using an explicit method and the
internal mode with a verticallyimplicit method. More recently
they used an implicit scheme developed by Casulli(1992) to
solve the external mode.
CH3D (Curvilinear Hydrodynamic Three Dimensional) model
developed by Sheng(1986) is a stateofart finite difference
numerical model which uses boundaryfitted curvilinear
coordinates in the horizontal directions and 7 coordinate
system in the vertical direction. The CH3D model also includes
a sediment transport model and water quality model.
TRIM (Tidal Residual Intertidal Mudflat Model) was
developed by Casulli and Cheng(1992). A Cartesian coordinate
system is used in both horizontal and vertical directions. The
model does not separate the external mode from the internal
mode but solves them together with a timeimplicit method.
However it uses a nonconservative EulerianLagrangian
method(ELM) for the advection terms in the momentum equations
and the linear advectiondiffusion equations. Hence for long
term simulation, there will not be conservation of mass,
momentum, heat, and salinity. In addition their model uses an
explicit scheme for the Coriolis terms which could render the
3
solution unstable.... The model allows wetting and drying in
shallow waters.
1.1 Objectives
The objectives of this study include the following :
* Produce a robust threedimensional model by improving the
Casulli and Cheng(1992) model by implementing a
conservative EulerianLagrangian Method(ELM) and adding
baroclinic terms and salinity equation,
* Improve the EulerianLagrangian method by implementing an
explicit scheme in Xmomentum equation and an implicit
scheme in Ymomentum equation for the Coriolis terms,
* Test the model with analytical solutions, and
* Compare the model results with those of CH3D model.
1.2 Organization of This Study
In Chapter 2, the governing threedimensional shallow
equations which are momentum, continuity and salinity
transport equations are written for a threedimensional
Cartesian coordinates with boundary and initial conditions.
In Chapter 3, the governing shallow equations and the
transport equation will be written by finite difference
equation with boundary conditions. And a preconditioned
conjugate gradient method to solve system of equations
In Chapter 4, the EulerianLagrangian method (ELM), or
4
semiLagrangian numerical method to solve linear and nonlinear
advection equations will be discussed in detail. Model results
are compared to analytical solution.
In Chapter 5, several analytical problems are developed
and compared to model results. Wind test, tidal propagation,
nonlinear term, Coriolis effect, bottom friction effect,
diffusion test, baroclinic effect and wetting and drying test
over tidal flats for both two and three dimensional cases.
In Chapter 6, the present model results will be compared
to those of CH3D model results.
In Chapter 7, the present model will be applied to
simulate winddriven circulation in Lake Okeechobee, Florida.
Chapter 8 will present a summary of the model and
conclusions.
CHAPTER 2
GOVERNING EQUATIONS AND BOUNDARY CONDITIONS IN THREE
DIMENSIONAL CARTESIAN COORDINATES
2.1 Introduction
In this chapter, the basic threedimensional equations
and boundary conditions for estuarine and coastal circulation
are described with respect to Cartesian (x, y, z) coordinate
system.
2.2 Governing Equations in Cartesian Coordinates
The four basic assumptions used for deriving the
governing equations are : (1) the flow is incompressible so
that the flow is nondivergent, (2) the vertical accelerations
are negligible compared tc gravity, thus, the vertical
pressure distributions satisfies hydrostatic approximation,
(3) density is approximated by its mean value except when
multiplied by gravity, i.e., the Boussinesq approximation is
valid, and (4) the eddy viscosity concept is used with
constant turbulent eddy viscosity and diffusivity.
Using a righthanded Cartesian coordinate system (x, y,
z) with the origin at the water surface and z measured upward,
the equations are written as following nonconservative form
(Lamb,1945) :
Continuity equation,
du dv dw
+ + o (2.1)
dx dy dz
Xmomentum equation,
du du du + u 1 dp 2u 2u
+ +v + w + A( +
dt dx dy dz po dx dx dy
d du
+ (AV )+ (2.2)
Ymomentum equation,
dv v dv dv 1 dp d2v d2
+ u + v + w= + AH( + )
dt dx dy az po dy dz px dy2
d dv
+ (Av )fu (2.3)
Zmomentum equation,
9p
pg (2.4)
dz
Additionally, by using the Leibnitz integration rule, the
hydrostatic pressure gradient can be split into the barotropic
and the baroclinic components as follows
h = gp(x,y,,t)V(x,y,t)+ gV px,y,,t)d
VPh = gp(X,y,q, t)V q(Xx y t)+ gJVp(xjy,t)d1
(2.5)
Finally the momentum equations in x and y directions are
du du
v + w
dy dz
g pd +
Pd +
po Z dx
d2 d2 2
AH( +2
9x 2\
dv
dt
dv
u +
dx
dv
v +
dy
g p d d2v d 2V
P d+ AH(dx2 +d2
po dy
Salinity equation,
dS dS
+ ux +
9t dx
dS dS d2S d 2S
v +w d ,( + 2 )
dy dz x dy2
Equation of state,
Where u(x,y,z,t), v(x,y,z,t) and w(x,y,z,t) are
the velocity
components in the horizontal x, y and vertical z direction, t
is the time, 7(x,y,t) is the free surface elevation, g is the
du
dt
St
du
U+
dx
d du
+ (A, )+ fv
dz dz
dv dq
Wd 
(2.6)
d dv
+(A, ) fu
dz dz
(2.7)
d dS
dz dz
(2.8)
p = p(T,S)
(2.9)
dq
gdx
8
gravitational acceleration and A,, D, and A,, D, are the
horizontal and vertical turbulent eddy coefficients, f is the
Coriolis parameter defined as 2Q sino where Q is the
rotational speed of the earth and 0 is the latitude, T is
temperature, S is salinity, p, is the density of freshwater, p
is the density, p is pressure, respectively. Various forms of
the equation of state can be used. Eckart (1958) used :
p = P/(a + 0.698P)
P= 5890+38T0.375T2+3S (2.10)
a = 1779.5 + 11.25T 0.0745T2 (3.8 + 0.10T)S
where T is in C, S is in psu and p is in g/cm3. CH3D uses
equation (2.10), while Blumberg(1978) used
p = p,(1+aS) (2.11)
where a is a function of temperature, and for a temperature
of 15C, a is 7.57x104 The value of a can be derived from
U, tables (Bowden, 1967). In the present model, Blumberg's
equation is used because the present model doesn't have a
temperature equation.
2.3 Boundary Conditions in Cartesian Coordinates
The boundary conditions at the free surface are specified
by the prescribed wind stresses Tx and T
du 2
pA, Tz = r= pCdauwU + vw (2.12)
pA, y = pCdaw u + v (2.13)
dz
where Tx and T" are the wind stresses in x and y directions
x y
at the free surface, uw and vw are the components of wind
speed measured at some distance above the free surface and C,
is the drag coefficient. The drag coefficient is normally a
function of the roughness of the sea surface and the wind
speed at some height above the water surface. For this study,
the empirical relationship developed by Garratt (1977) is
used. Garratt defined the drag coefficient as a linear
function of wind speed measured at 10 meters above the water
surface, W,
Cd = 0.001(0.75 + 0.067W,) (2.14)
where W is the wind speed in meters per second.
At the free surface, the kinematic condition states that
Dr dq d] dc)rl d(r
= +u +v (2.15)
Dt dt dx dy dt
Integrating the continuity equation over the depth and using
the kinematic condition at the free surface leads to the
following free surface equation :
+ udz + vdz = 0 (2.16)
dt dx dy
where h(x,y) is the water depth measured from the undisturbed
water surface and H(x,y,t) is the total water depth, given by
H(x,y,t) = h(x,y)+ rl(x,y,t).
The boundary conditions at the sedimentwater interface
are given by specifying the bottom stress in a form of the
ChezyManning formula
du
pAz = = yu (2.17)
dz
dv
pAd z b 7Y (2.18)
p4 =y
11
g u2 + v2
where = u and v are velocities at the bottom
Cz
layer, and C2 is the Chezy friction coefficient which can be
formulated as
RY
C = 4.64 (2.19)
n
where R is the hydraulic radius given in centimeters and n
is Manning's n. In shallow estuaries, the hydraulic radius
can be approximated by the total depth. In CH3D model, the
bottom stresses are calculated differently :
(rbb)= PCdb(11 + v2)1/2 (1, ) (2.20)
where Cdb is the bottom friction coefficient, and tj and v,
represent the velocity components at the first grid point
above the bottom. Taking z, as half of the bottom layer
thickness (which starts at the bottom roughness height, z0),
Cdb, for a hydraulic rough flow, is given by (Sheng, 1983)
2 ( Z (2.21)2
Cdb = k2 In (2.21)
L z0J _
where k is the von Karman constant.
12
The open boundary condition for the surface elevation is
prescribed by
N ( 2nt /
q(x,y,t)= Z Ancos + (2.22)
n=1 Tn
where A Tn and On are the amplitude, period and phase angle
of the tidal constituents. When open boundary conditions are
given interms of the surface elevation, the normal velocity
component is assumed to be of zero slope while the tangential
velocity component may be either (1) zero, or (2) zero slope
or (3) computed from the momentum equations. In the present
study, it is assume that the velocity gradients are zero at
the open boundary. At the fixed boundary, no flux through the
boundary is allowed.
Salinity is prescribed as a function of time and depth
at the open boundary or salinity is calculated from the one
dimensional advection equation. In this study, the salinity at
the open boundary is prescribed. At the water surface and the
bottom the normal salinity flux vanishes, i.e.
aS
D, 0 (2.23)
dz
If the estuary is connected to freshwater river flow then
the boundary conditions at the interface are given as
S(x, y,z,t) = 0
(2.24)
13
u(x, y,z, t) (2.25)
(h + q)b
where q is the river flow, usually a constant, (m3/s) and b
is width of the corresponding river cells. At the open
boundary the salinity is prescribed as
S(x,y,z,t)= So (2.26)
where So is the prescribed salinity value in psu.
CHAPTER 3
A THREEDIMENSIONAL NUMERICAL MODEL
3.1 Introduction
The governing equations and boundary conditions shown in
Chapter 2 are discretized using finite difference method
similar to that described in Casulli and Cheng (1992). A
simple stability analysis of the twodimensional, vertically
averaged shallow water equations showed that the celerity term
g7i in the equation affects the stability of the finite
difference equations. Results of this analysis led to the
development of an implicit method (Casulli and Cheng, 1992)
and a semiimplicit method (Casulli and Cattani, 1994) which
have been proven to be unconditionally stable.
When the implicit method is used, there is undesirable
numerical damping in the numerical solution. To improve the
accuracy of the implicit method, an implicitness parameter 6
is introduced such that for 0 = 1, this method is fully
implicit and for 8 = /2, this method is semiimplicit and
remains unconditionally stable and doesn't have undesirable
numerical wave damping.
15
In this chapter, the threedimensional semiimplicit
circulation model is developed by using the so called 8
scheme.
3.2 Finite difference form of momentum equations
The horizontal gradients of surface elevation in the
momentum equations (2.6), and (2.7) and the velocity gradients
in the free surface equation (2.16) will be discretized by the
8scheme, and the vertical mixing terms will be discretized
implicitly because the vertical length scale is much smaller
than the horizontal one.
Figure 3.1 Schematic diagram of stencil
As shown in Figure 3.1, the spatial mesh consists of
rectangular boxes of length Ax width Ay and height Azk
16
i.e., a cellcentered grid system is used. Each box is
numbered at its center with indices i, j, and k. The discrete
u velocity is defined at half integer i and integers j and k;
velocity v is defined at integers i, k and half integer j;
vertical velocity w is defined at integers i, j and half
integer k. Water surface elevation ] is defined at the center
of each column i and j. The density of water is defined at the
center of each box i, j and k. The water depth h(x,y) is
specified at the u and v horizontal points.
Then, parametrized semiimplicit discretization for the
momentum equations (2.6) (2.7) and for the free surface
equation (2.15) takes the following form.
= g i j +lFm Pi2j,m
u+i+ 2j,k+,jk y 1 jk k k ,,
AZ (pnj. A A Z il i
Vk+J Az ykj, Azi+ +,2
i+ jk+k i+ ,jk
++
+(+ 3j,k
n ,
i+ fv,+ j,k
(3.1)
n+1 =At [1) + (l ). 1 ]
i,+Y,k= Fvjj+Y,k x 1 ,";j W+ j ;,
=k '
o m k ta^ j1,m i,' n Az m
n+1 n+1
vv
i,j+ Y2,k+1 i,j+ 2,k
A V A
i,j+y ,k+Y
+At
 Un+1:
JU ,s ,
n n
Pi,j+,m Pi,j,m
2
Vn+1 n_ +1
i,j+ 2,k i,j+ 2,k1
Si,j,+ 2k12
Az ,+
i,/j+,k
(3.2)
n+1 n At M n+1 M n+ 1
+ n On A +1 V Az U7'1
j Ax 0 ZAz+k= ij+ j,k ij,k E zi,j,ki.,j,kj
S Z i1 j+12,kV i,+ 1,k Z i,j ,k ij2,k
Ay [k=nI /2 k=m ' 
where m and M denote the limit of index representing the
Ax (1 ) Az+ 12,k i+yj,k M AZ2,1jk ,,k
(10) z ,/ un/ AZ n / ," /
At M M
Ay i2y i2zk j 2zk i/J_ 2' z'
(3.3)
where m and M denote the limit of kindex representing the
bottom and the top finite difference stencil, respectively. In
equation (3.1) and (3.2), Fu and Fv are explicit, nonlinear
18
finite difference operators which include the contribution
arising from the momentum advection terms which are
discretized using EulerianLagrangian method and the
horizontal diffusion terms.
Here, it is important to discuss the treatment of
Coriolis terms in the equations. The Coriolis term in the X
momentum equation is discretized explicitly, while the term in
the Ymomentum equation is discretized implicitly (Sheng,
1983). Although the Coriolis term in the Ymomentum equation
is discretized implicitly, the actual computation of the term
is explicit because the velocity Un+1 is already obtained from
the Xmomentum equation. Stability analysis of this method
shows that the time step of the scheme should be less than
1 which is about 3 hours and much bigger than the other
time step limits, e.g. the horizontal diffusion limit. Since
the numerical grid is staggered, v is not defined at the u
grid point and u is not defined at the v grid point. A simple
way to define the undefined velocity at a grid point is to
take the average of the velocities at four neighboring
points.
The fictitious values of u and v above the free surface
and below the bottom in (3.1) and (3.2) can be eliminated by
using the surface and bottom boundary conditions as shown in
(2.17) and (2.18). The difference forms of (2.17) and (2.18)
are
A
Aii+.M+Y
AV..jA,
zi+ 2,j,M+12
n+l n+l
VI ,J,+ 1 + ,M J= .y
i ,jj+. M+Y wM
^^^^.t~MiX
n+l n+1
U+Jm u+ ,m1 n+ / n+1
i+2, i+ 2,,m
(3.4)
(3.5)
(3.6)
(3.7)
n+1 n+1
+, +2 +,m i, j+ =l n+ n+1
Avi'j+ Azij+ i 2,+ m ij+ l2,m
,1+ 2,m 2
n+2 = F un+ 12 + (i 12
12n+2 n At n
i+ 2,j,m i+ 2,j,n 2Ax i+l,
n+2 Vn At q .
i,j+ 2,m i,j+ 2,m 2Ay j
and
where
and
(3.8)
(3.9.)
(3.10)
qi:j)
By using the boundary conditions (3.4) through (3.7), the
fictitious velocities in equation (3.1) through (3.3) will be
replaced by values of u and v defined within the
computational domain.
3.3 Solution algorithm
To solve equations (3.1), (3.2) and (3.3)
computationally, the equations are first written in matrix
form as
A ,+ 8 + n )AZn/ AU (3.12)
[0 (qi+Ij q _.n(3+l.11)
=." AZ7" U ""+ AZI" i 7(+.l1
At
I AZJ V7. A" (i 13.13)
Atwhere [ And A are defined as follows
Ay ij+1 Q2 I j2 j2
where U, V AZ, G, 8 AndA are defined as follows :
Un+I
,n+I
Un+I
ui+Y,1/nl
V n+
i,j+ Y2, M
ij+ ,M1
SAZ=
n+1
Vn+1
iij+Y/r )
AzM
AzMI1
Az,., )
A t 1,,
AzM Ft" g(1 )(~ ,) + C + BC + Atr,
Az, g t (1 0)( ) ,",) + MC + BC" ,
AZM~~~~~~ ~~~~1 L F., ,_1 IV2 +,. J,M+l+., ,+ < .._
1.+Y g/
Az [ F"'+2,jm
Az,, [Fit" iI
A (1 ) "i) + C" + BC+, .
At + BC"
 g ( 0)(,, ,.n.) + C +,
At
Az., FI ,Y g (1 0)(,,. ,.) + BC M + At
At
Az, FV ,m g(1 )( i," ) + BC" ,
SAz, [,Fv g (1 0)(q~i+ + BC~ + ]
where C's are Coriolis terms in the Xmomentum equations
(3.14)
(3.15)
(3.16)
vn+1
i,j+/2
G"n =
+,A2,
G" =
i.J+/2
Ci ,j =
AtfyV>+I,j,M
AtfV" M
Y,.'
A tfV+",i+ jm+
and BC's are the baroclinic terms
Az (p' p,,k ) Az
i+A,j,k i+I,j,k )k i+ ,jIM
At M
 pZA Y2 A (PinIjk P1,jk) j
At (pl P,,J) Az
 p k=M1 +,jg k i+"I,j,k ,j,k)" i+4 .j 
ik +1 n+ 2 k n, j, k injk n+2, jlm+
SAt m
K ^'
Azi+,j,k (Pi+,j,k Pij,k) AZ+ ,,.
n
Pi+I,j,M Pi,j,M
2
n n 
Pi+I,j,M PiJ,MI
Pi+l,j,m+l Pi,j,m+l
2
,l n
Pi+l,j,m i,],
2
(3.18)
(3.17)
AtL
gpx k=M
BC"
i+Y2,j,k
At [
g p 1y
pAYk=
AZij+ k Pi,j+l,k i,j,k) ij+,M
 P'A y k' =n A z=M\ i,j+ 2,k( i,j+l,k
At "
g Az" 'J (P'lj+l'k
, k=Ml ij+,k ij+,k
At [
poy k+1 A,,j+,ki j+1,k
p0,AY Lk=m
inj,k j+YuM1
Pij,k j+,) A
P,,,) AzIj+~. ,
Pi'j,k) Znij+Ym
n n
Pij+,M Pij,M
2
n _pm
Pi,j+I,M1 ~ Pi,J,M
2
,n ii "
i,j+1,mn+l i.j,m+
2
n n
tily+lm Pij,m
2
(3.19)
A) [ Z+2,j) T /U" ,
i i+ 2,j
t (1) Az,
Ay ij+
 (Az 2,j T
A ,2 ,j
V" Az
i,j+ 2 i 
SAzA, + aMY am_
a y AzA + a _y+ aM_
 a M
 am+ Az,. + aml + /rAt
m+^ "'*/!2
(3.21)
A, At
with ak Az
BC =
i,j+ ,k
At
4,. = ", (I
i,j 'y A
and
(3.20)
Vn I]
jY2_
Formal substitution of Equations (3.1) and (3.2)into
Equation (3.3) yields :
n+1 (5n
,j = ,j
At [ AZ A'G n g (,I 7 '+ )]AAz' ]
S AzT AG g [( n+ 
A 6 AzAG g y q _, ]A Az
At [[AG Ag [0( ,)]. {AAz}+ ,.,
S ['j+ AG} [ g[ (7 ] AJAx, 'z
+ [AI { G}i g (2, [ ,i )] A' 'A z
2  g (  j)] "i,j+, { A AzY2 ,
AAt r ,n+l ,++ I
zJZ 1AG g ^,) A.
[Ai,+{A G} g x[8( l+)]A i+ {A z}+ y
+ [AAt A G [0(un 1)] A 1Az}
+' 9' j ]y
+ Id 2 {dIGi+Y2,j A A A ) A J A }i+i2',
+ A A G g [ n+ ( +1 1)]A 'Az}
(3 .22
[A,'_y G ij 77i+il)] 2 A 2,
25
Neglecting the Coriolis terms and rearranging the
Equation (3.22) by moving "n+l" terms to the lefthandside,
Equation (3.22) becomes
n", 2 2[ [(AZ)TA Z] (iij )i,
[(Az)A Az]i/ (i," 
[(Az)TAAz]. (,+:; +;)]
= o [(A) A, n(AzAG1n
= [(A 'A (+,ji,:;j0 ]
At [([A(Az)A A'G y (Az)TA1G]I.
Equation (3.22) can be rewritten into a compact form :
d, "+1 S n+l S n+1 S .n+1 S qn+l n
i,ij iY2,ji+l,j i2 ,j i,J+ i Si/J1 = q=,j
(3.24)
where
S At2 n T 1 A n
S" g,j x2[(Az A 'Az
'12 n
S/ = g 82[(Az)T AziA]
'ij 2 g' 12
d, = I+ S +SI +S" +Sn (3.25)
Sj i+Yj iYj i+ ,+Y (3.25)
q. = At [(Az)T AG] [(Az)T A'G
Ax 1,+ ,J i2,J
r 
At [(Az) Aol,/ [(Az)TA'G],
Ay Q+ i j 2
26
If Coriolis terms are considered, Equation (3.24) becomes
more complicated :
Cn n+l n. n+l nn+1l 1 n+l
Glni,IjI 2 li,j qilj+l ,/ij q /iq
qn+l i nn+l n(3.18)
C,n.n _+l C",n7i n+l  nqn+l n.n+l j? O"
 6 j+l n +1j1 C8 ]i+l,j '9 li+1,j+l = RHS"
where
C" = FIP"n_ P"n
SiY2 i Y2,j1
C" = CM" + FIP" P" FIP" P",
2 I ,j ij+ iY j ij J ' 2,J
C"~= FIpP" P"
3 i, j+ 2 i 2,j+.)
C = HMn" + FIP" P" FIP" ,P"/ n l
4 i ijY i,jY2 i+ ,jl1 i jY2 ' j1
C"= 1+ CM" + CM" + HM" + HM"
5 < i+ ,j iXj i+ + ij
FIP" P" FIPn P" FIP. P/ FIP P"
i,j+Y2 i+2,j i,j+ 2 i ,j i, 2 i+ j i, 2 i 2,
S= ( HMn" + FIP P" FIP" P"
6 j+j+= i+yj+1 i, j+ y2i j+
C = ( FIP" P"
C7 ji,j y2Y i+ y j +21
C"= CM" FIP" P" FIP" /P"
C" = +FIP".Y i,j
9 ij+ i+ y,j+1]
(3.27)
with
At2 At2 gA2
C= g2 H= g 2 1 2
HJg Ay AxAy (3.28)
(3.28)
A tf
P= (AzAAZ), M= (AzA'AZ), F= 
4
and
RHS" = S" 6 (Az) A G] [(AZ) AG
At [() A'G [(Az) AG]
Ay i[j+ Gi+ ,2
+{ AzzA}j+{ c'G}>+, + {AzTA};:,+{ A'G}1,+ ]
+ Af At O[{AzTA}1 /{AG}, +{AzT'A}1' _{A'G}i
4 Ay L
+{AzT'A} {A 'G}". + {Az A}' ,/{ A'G}" ,
+l A [z '. _ A n +' 1Az 'Gn
(3.29)
Since the coefficient matrix A is symmetric and positive
definite, A' is also positive definite, and therefore
(Az)TA'Az is a nonnegative number. Moreover the coefficient
matrix A is a tridiagonal matrix, and the inversion of the
matrix A is straightforward and efficient by using Thomas
algorithm. Finally Equation (3.23) constitutes a linear five
28
diagonal system of equations for new water surface elevation
"' and the associated matrix is symmetric and strictly
diagonally dominant with positive elements on the main
diagonal and negative ones elsewhere. Thus the system is
positive definite and has a unique solution. If the Coriolis
terms are considered, Equation (3.26) consists of a nine
diagonal system of equation, i.e., a block tridiagonal
matrix. The ninediagonal matrix is positive definite but not
symmetric. The fivediagonal system of equations can be solved
very efficiently by a preconditioned conjugate gradient
method. For the ninediagonal matrix case, the block tri
diagonal algorithm or successive overrelaxation(SOR) method
can be applied. Present model uses SOR method. After obtaining
the new water surface elevation, the new horizontal velocities
u and v can be calculated by using Equation (3.11) and
(3.12).
Finally, from the continuity equation, the new vertical
velocity w can be obtained by setting wn1 =0 and
n A nx
n+1 n+1 AZ +2,j,k i+ ,jk i ,j,k i'Y ,k
i,j,k+2 i,j,k2 AX
Azn v AZn Vn
i,j+ ,k ij+y2,k i,jy,k i,jY,k
Ay (3.27)
k= m,m+ 1,.,M
3.4 Coniugate gradient method
Equation (3.24) constitutes a very large system of
equations, which regime a large amount of computing time for
its solution. Since the system of equations (3.24) has
symmetric and positive definite sparse matrix, the equation
can be solved efficiently by the i.e., conjugate gradient
method.
Conjugate gradient method ( Beckman, 1959. Casulli, 1992)
is the most popular iterative method for solving large system
of equation Ax= b where the matrix A is square, symmetric,
and positive definite (or positive indefinite).
Consider a quadratic function of a vector with the form
1
f(x)= xrAxbTx+ c (3.28)
2
where A is a matrix, x and b are vectors, and c is a scalar
constant. If the matrix A is symmetric and positive definite,
f(x) is minimized by the solution to Ax= b The definition
of positive definite matrix is that if xTAx>0, for every
nonzero vector x, then the matrix A is said to be positive
definite matrix.
The gradient of a quadratic form (J. S. Ahewchuk, 1994)
is defined to be
f(x)
dx
a
f(x)
f(x)= ax2 (3.29)
Sf (x)
With equation (3.29), the gradient of equation(3.28) is
1 1
f'(x) = ATx+ Axb (3.30)
2 2
If A is symmetric, this equation reduces to
f'(x) = Ax b (3.31)
Setting the gradient to zero, we obtain Ax= b which we want
to solve. Therefore, the solution to Ax=b is a critical
point of the function f(x) If A is positive definite as
well as symmetric, this solution is a minimum of f(x), so
Ax = b can be solved by finding an x that minimizes f(x) It
should be noted that if A is not symmetric, the conjugate
gradient method finds the solution to the system (J. S.
Ahewchuk, 1994)
1
(AT+ A)x= b (3.32)
2
31
To solve Equation (3.24) using the conjugate gradient
method, the equation can be written in the normalized form
(Casulli and Cheng, 1992) :
Sn S"
n n+ L+'J in .n+l S /2J n+l
dV iJ d,+ijj
(3.33)
S" S" /
d Sn s 1 j
,J+ pz n+1 Si'l2 n+1 qij
n n nJ1 iJI
which, by changing the variable i1 to e,j is equivalent
to
ei, a ei +lj a e,, e e,,+l  a,,j e,,j_ = b, : (3.34)
where
Sn" S
Sj2,, I ,d" d i,
ai+ j d i j a ,j+ /2 vdijd i j+i biJ
Now the equation (3.34) is applied to the conjugate
gradient method to obtain the unknown value of ei.j The
conjugate gradient algorithm to solve the system of equations
(3.26) takes the following steps.
(1) Guess e(0)
(2) Set
(2) Set
(0) (0) (0) (0) (0)
ij ij i, i+ ,j i+ ,j ,j ,j
(0) (0) b
a,i i,j+l ai,j_ i,j ij
(3.35)
(3) Then for k=0,1,2, ... and until (r(k),r(k)) < calculate
(k) (k) (k))
where a' =(p() M(k))
(k) Mp(k)
(k+1) = (k) (k) M (k) ,
l,J i,J (k ( J()
p(k+l) (k+l) + (k) ) P (), where (k)=
ij = + i,j
(r(k+) (k+1))
(rk),r(k))
r ), r
In equation (3.37) each element of the vector Mp which is
preconditioner is simply given by
Mp() = p(k) (k) (k)
ij "i, i+ 2,Pi+,j i Y2,j Jil
(. 9)
(k) (k)
i,j+ 21i,j+l j 2 i,j1
3.5 Salinity Scheme
The finite difference equation (2.8) is similar to the
e(k+1) (k) (k) (k)
ei,j = ij a Pi,]
(3.36)
(3.37)
(3.38)
33
hydrodynamic finite difference equation. The vertical
diffusion term is treated implicitly, while the horizontal
diffusion and advection terms are treated explicitly. The
advection term will be solved by the conservative ELM method.
The finite difference equation for the salinity equation
(2.8) has the form as
.n+ At D Si,,k+ i,,k ijk Sijk1
S D S+D ' _'
i,j,k Az i,j,k+ 2 Azn i,j,k2 Az" k
Ijk i,j,k+ ij,k
= FS"jk +2 ik S ,jk 2 ,ijk + I+1,j,k (3.40)
DAtAzi,,k ( )
SDHt ijk S 2S" + S +"k
Ay2 i j1,k ijk i,j+1,k
where FSi, k is an explicit discretization of the substantial
derivatives of convective term by ELM. The left hand side of
Equation (3.40) constitutes a tridiagonal matrix and the
solution procedure for the unknown new value of Sj is
exactly the same as that for the momentum equation (2.3).
CHAPTER 4
EULERIANLAGRANGIAN (SEMILAGRANGIAN) METHOD
4.1 Introduction
One major challenge in numerical modeling of free surface
flows is the accurate and efficient treatment of the advection
/convection terms in the governing equations. There are
numerous numerical schemes for solving the advection/
convection terms. In this chapter, numerical treatment of the
advection/convection equations by the EulerianLagrangian
Method (ELM) will be discussed. Other methods such as upwind,
Leith's method and total variation diminishing (TVD) method
are briefly presented in Appendix A.
4.2 Review of EulerianLaqrancian Method (ELM)
The EulerianLagrangian Method uses the Lagrangian form
of governing equations in an Eulerian computational grid
system. The method (Starniforth and Cote, 1991) is generally
referred to as semiLagrangian method in the community of
numerical weather prediction, but is called Eulerian
Lagrangian Method (ELM) by estuarine and coastal modelers,
e.g., Baptista (1987), Casulli and Cheng (1992).
The ELM solves the equation for a variable e.g., a
velocity component or concentration by following a fluid
particle along the characteristic backward in time and then
interpolating at its foot (departure point of the particle) to
obtain the values of the variables.
Figure 4.1 Characteristic lines and their feet
The ELM was also used by Tang and Adams (1995) with
finite difference method and by Baptista (1984) with finite
element method.
The major advantage of the ELM is that the method is
unconditionally stable although it is explicit. In addition,
it is transportive and accurate (Casulli, 1994). Despite these
advantages there are several disadvantages of this method.
Although the method is accurate, its accuracy is dependent on
the interpolation function used. For instance, if a linear
interpolation function (first order accuracy) is used, there
is excessive numerical diffusion in the numerical solution,
with the maximum diffusion on the order of Ax2/(8At) in the
ID case, although the diffusion can be controlled by reducing
the spatial increment or by increasing the time step size. The
second problem is that if higher order (second order of
higher) interpolation is used the numerical solution contains
unphysical numerical oscillations. Godunov (1959) was the
first to show that any approximation methods of order higher
than one is unable to deliver a solution free of oscillations.
This result, which is known as Godunov's theorem, is a major
stumbling block for the development of integration schemes for
hyperbolic problems.
The most significant problem of ELM is that the scheme
doesn't have the conservative property which is that momentum,
mass, and any species should be conserved locally and globally
throughout the numerical simulation period. Roache (1972)
indicated that all methods of characteristics (including ELM)
which solve the nonconservative form of momentum or transport
equations are nonconservative. If the governing equation is
written in conservative form then many numerical schemes can
have conservative property.
There are several ELM methods which yield solutions free
of oscillations and satisfying conservative property
(Priestley 1993 and 1994, Gravel and Staniforth 1994, Rasch
1994, Tang and Adams 1995, Leslie and Purser 1995, Lin and
37
Rood 1996 and 1997). Although the methods have different names
the basic idea is the same. To control numerical oscillations
these methods enforce the solution to be "monotone" i.e. free
of oscillations. To achieve the conservative property, the
conservative ELM uses FCT (Flux Corrected Transport)like
method and checks for global conservation at every time step.
4.3 Derivation of EulerianLaqrangian Method (ELM)
Consider the diffusionfree nonconservative advection
equation :
dc
+ uVc= 0 (4.1)
dt
where the c(x,t) represents any scalar quantities, u is
velocity field, and V is a gradient operator. Most of the
fundamental equations in fluid dynamics can be derived from
first principles in either Lagrangian form or Eulerian form.
Lagrangian equations describe the evolution of the flow that
would be observed following the motion of an individual parcel
of fluid. While Eulerian equations describe the evolution that
would be observed at a fixed point in space. The one
dimensional Eulerian advection equation (4.1) can be expressed
in Lagrangian form as
Dc
= 0 (4.2)
Dt
The mathematical equivalence of Equations (4.1) and (4.2)
follows from the definition of total derivative,
D 8 dx 8
Dt ~t dt dx
(4.3)
and the definition of the velocity,
dx
(4.4)
When a Lagrangian numerical treatment is applied to Equation
(4.2) the computational grid will be continuously deforming in
the general case when u is nonconstant. For operational
advantage, Equation (4.2) will be discretized on a fixed
Eulerian grid system. Figure 4.2 shows the Eulerian Lagrangian
mesh system.
d'
]
/'.
I]
]
/*
]
]
/*
/]
,^ L _______ _
Figure 4.2 An Eulerian Lagrangian Mesh
A finite difference scheme for Equation (4.2) is simply
'a in
uAt
39
Cn+1 = n (4.5)
where c = c[(ia)Ax,nAt] and n is previous time step, i is
the index at a grid point and a is the CFL number. In general
the CFL number is not an integer, therefore (ia ) is not the
index of a grid point and a proper interpolation formula must
be used to define ina The stability, numerical diffusion and
unphysical oscillations of (4.5) depend on the interpolation
formula chosen. If a linear interpolation between (i1) and
(i) is used to estimate Ca_ we obtain the first order upwind
scheme. If a quadratic polynomial fit is used to interpolate
between (i1),(i), and (i+1), one obtains the Leith's
method.
The ELM uses a generalization of the interpolation
concept of cina between two or more grid points which do not
necessarily include the point (i). Consider first the case
that ci,7 is taken to be a linear interpolation using one node
upstream and one downstream. For a given a 0 let n be the
integer part of a and p the decimal part, then a = n+p,
0Op<1. In this case Equation (4.5) becomes
40
k+1 k k k
Ci Ci Pn(i C _) (4n6)
or, equivalently
Cik+1 i PC + PCikn (4.7)
Note that if a<1, then n=0, p=a and finite difference
equation of Equation (4.6) and Equation(4.7) reduce to the
first order upwind method, since 05p<1.
Since the velocity u is generally nonuniform, the
correct value of a can be found from the solution of the
ordinary differential equation (4.4) using any backward
trajectory computation. The velocity u is known only at time
level tn, hence it will be assumed to be constant over one
time step. Then, at each grid point (i), Equation (4.4) will
be integrated numerically from t" to t"n In this study the
backward Euler method is used for computation of the
trajectory. To compute the velocity, the time step At is
divided into N equal increments, T = At/N, and Equation
(4.4) is discretized backward as
xs X_ ruk (xS), N Xi
X X UX X (4.8)
s= N,N N 2,...,2,1
where uk(xs) is interpolated with any interpolation formula.
x. X
Then, at xi, a can defined by a=
Ax
The finite difference equation (4.6) possesses the min
max property (Casulli, 1987), thus this method is entirely
free of oscillations. Moreover the minmax principle also
implies stability, hence the ELM shown in (4.20), though
explicit, is unconditionally stable. Rigorous stability
analysis can be performed using Von Neumann method. Assuming
a solution to equation (4.1) as
Ck A ei(j (4.9)
where A' is the amplification factor, I is the grid point
index and i= v" Substituting Equation (4.9) into (4.1)
Ak = [1 a( i ei(k))ei(knAx) (4.10)
The magnitude of Equation (4.10) is
Ak = 21 2a(1 a)(1 coskAx) (4.11)
If the grid CFL number satisfies 0 a<1, the scheme is
stable. Because the departure point lies within the
interpolation interval (in, in1), the magnitude of (4.11)
is always less than one, consequently the scheme is
unconditionally stable (Starniforth and Cote, 1991).
To examine the artificial diffusion introduced by
Equation(4.6), each term in (4.6) is expanded into Taylor
series to produce the following :
Dc Ax2 d2C
=t p(p 1) + HOT (4.12)
Dt 2A z2
where HOT stands for higher order term.
Since p<1, the least upper bound for the artificial
Ax2
diffusion coefficient is which is the same as the upwind
8At
method. However, since the ELM is unconditionally stable and
the value of At can be arbitrarily large, the maximum
numerical diffusion can be controlled either by reducing the
spatial increment or by increasing the time step.
The ELM with linear interpolation function has no
spurious numerical oscillation and little numerical diffusion,
but it is only first order accurate. If three point
interpolation function is used, higher order ELM can be
obtained.
Let a=n+p, where n is integer, IP <1 and np! 0. In
this case Equation (4.5) becomes (Casulli, 1987)
k+ 1 2 Ck + 2 k 2 +113)
k+l C (p + p) _,,+1 + (1 p2c + (P
Cj \' l
Note that, if an< 1 then (4.13) reduces to the Leith's
method and it is of second order accuracy. However, when p 0
one gets a negative coefficient on the right hand side of
Equation (4.13), consequently this method does not posses the
minmax property and hence it will introduce spurious
oscillation.
4.4 Conservative EulerianLagrancian Method (ELM)
EulerianLagrangian Methods described above do not
possess the conservation property because the finite
difference equations are obtained from the nonconservative
form of the differential equation. To achieve conservation,
some special treatments are needed. Several conservative ELMs
were introduced in the past decade, inducing Priestley
(1993,1994), Gravel and Starniforth (1992,1994), Rasch (1994),
Tang and Adams (1995), Leslie and Purser (1995), and Lin and
Rood (1996,1997). The conservative ELM presented here is based
on the QuasiMonotone SemiLagrangian (QMSL) scheme (Bermejo
and Staniforth, 1992) and QuasiConservative SemiLagrangian
(QCSL) scheme (Priestley 1992, Garcia and Priestley 1994).
The conservative ELM consists of two steps. The first
step obtains a low order solution and a monotone higher order
44
solution, and the second step combines the two solutions. As
seen in the previous section, the higher order solution does
not have numerical diffusion but has numerical oscillation.
The low order solution is a monotone scheme.
Consider the scalar, linear advection equation for u(x,t)
u,+a Vu=0 (4.14)
To obtain a higher order monotone numerical scheme, the Flux
CorrectedTransport (FCT) method (Boris and Book 1973, Zalesak
1979) is used. Following the FCT method, the solution is :
+ = ULk + Fk1 (4.15)
where U"+ is the value of the approximate solution at the
grid point xk at tn+l i.e., foot or departure point of the
particle which is obtained by Euler backward method and ULk
is the low order approximation at the departure point at t",
and Fkn+ is a higher order correction term such that it
approximates the exact solution with an order higher than that
of the low order solution. The higher order approximation is
any solution which is higher than order of one. The higher
order solution used here is a quadratic second order solution.
The overall procedure is :
1) Obtain a monotone low order approximate solution U' ,
2) Compute a higher order approximate solution U", by a
standard ELM of order higher than 1,
3) At any instance t"+ set
Fkn+l Cn+1 (Unl Un+l
k = k (t H,k L,k
(4.16)
where Cc is a grid point coefficient that controls the amount
of correction to be added to the low order approximation in
order for the numerical solution to be monotone. The procedure
of determining the coefficient which is called a limiting
procedure is as follows :
4.1) Define U+ as the maximum U value of the
neighboring points i.e.
U = max(U1,, Uk22, Uk3, Uk4
(4.17)
and U as the minimum U value of the neighboring points
i.e.
U = min(Uf, U"2, U";, U;4
(4.18)
4.2) Set
+ = U U n
Q= U Un (4.19)
T V L,k
p = Un"+ n+1
SH,k L,k
4.3) If P> O,
n+1
Cq = min(l,)
If P
e
C' = min(l, (4.20)
Pk
If P= 0,
UH = ui ,i. e, Fk+ = 0
Hk UL,k Fk' 0 
Once the low order solution and the monotone higher order
solution are obtained, the next step is to recover a
conservative solution. The final quasiconservative solution
can be written as
Uk = LkUHk + (1 ak)U"L (4.21)
where ak (0 ak < 1) is a constant to be determined by the
recovering process.
The procedure to recover conservation is to compare local
flux with the overall flux of the quantity which is related to
aCk. The procedure to define ak is as follows (Priestley, 1993
and 1994)
1) Define
SU(x)dx = Udx =C
k = (UH,k UL,k)Sk (4.22)
where Sk is the area associated with node k
2) Compute
ak(U U)Sk = C C UL,kSk = C* (4.23)
Assume that
kmax (4.24)
If this is not the case, then the definitions are changed to
k k A(4.25)
C*  C*
3) determine (k :
Step 1
If /i < 0 then a, = ax Otherwise ak = 0
Step 2 : Define surplus
SURPLUS = C akk
Step 3 : Define the average value of a as :
SURPLUS
aavg
Step 4
max a=max
If avg < ak then kk avg Otherwise ak = a
The lower order interpolation functions include
linear(lD), bilinear(2D), and trilinear(3D) interpolations.
The higher order interpolation functions include quadratic
interpolation function(lD), biquadratic(2D), and tri
quadratic(3D) interpolation functions. The lower and higher
order interpolation function for lD, 2D, and 3D cases are :
1) One dimensional case
Linear interpolation : 2 points used
U+ = (1 p)Uf + pU, (4.26)
Quadratic Interpolation : 3 points used
U:+ = 0.5(p2 + p)Ui_,_ + (1 p2)U1_, + 0.5(p2 p)U0_,2 (4.27)
2) Two dimensional Case
Bilinear Interpolation : 4 points used
U!' = (1 p)[(l q)Utn,j + qU'_.. jm]
+ p[(1 q)U_nl,jm + qUln1,_m, (4.28)
Biquadratic Interpolation : 9 points used
U t+U 
I,'
P(1+ p)[( (1+ q)UI  + (1 q2 )Uil (1 q)Uin2 ,jz
+(1p2)[ (1+ q)U, + ( )U,'jm q)Uinjm+
(1 p)[( (+ q) n I + (1 q) ,,+ ,,, (1 q) Ufn+,,n+l]
(4.29)
3) Three dimensional case
TriLinear Interpolation : 8 points used
U+=I
(1 r){(l p)[(l q)Ut[,jm,kn + qUtl,jm1,kn
+ p[(1 q) Ugll,j,,kn + q Ul11,jm1,kn }
+r{(1 p)[(1 q)Ul,jm,kn1 + qUil,jm,kn]
+p[(1 q)U _,1,jm,kn1 + qUll jm1l,knl]} (4.30)
TriQuadratic Interpolation : 27 points used
UIt = (P2 + p)(q2 + q)@
11
U ^ k ,2 2
+ + ()Ui1jk_1(1r + 
(r 2 + rU11, jm1,kn1 r2 ) ,',jm,k. + (r 2 r)UiCl,jm.,k_n+l
2 2
+ (p + p)(l q2)
1 1
(r2 + j,kn1 r2)Ui ,jm,kn + (r r)Ui ,j,kn+
2 2
1
+ I(p + p)(4 q).
4
1 1
(r2 + r)Ui_,_ + (1 r)U j1 + (r
2
{ (r2 + r)U,j,k + ( r2)Ulj + I(r 2 r)Uj,kn+_
2 )i__ (nkn 2+ r) ^ijki }
1
+ ( )(q2 q)
2
1 1
S(r2 + r)Uijm,knI + (1 r 2)Ui+1,1,kn + (r2 r)U/+,jm,,kn+}
2 2
1 1
((r2 r+ r)U+,jl,kI + (1 r2)Uiljm + k 2 r)Uil+,jm,kn+
1
2 2
+ (p2 p)(l2 q).
2
1 1
(2 (r2 U+1,jm+1,kn1 2)Ull+1,]jm+I,kn 2 2 r)1Uil+I,jm+1,kn+1
(4.31)
where U can be any quantity such as velocity, salinity or
concentration.
The ELM can be applied to the nonlinear momentum equation
2 2
(4.31)
where U can be any quantity such as velocity, salinity or
concentration.
The ELM can be applied to the nonlinear momentum equation
51
with nonlinear advection. The procedure is exactly the same as
that for the linear advection equation.
4.5 Numerical Experiments
4.5.1 Linear Advection Case
To test the various algorithms which are described above
and Appendix A, various numerical experiments were performed.
The computational domain is 1000 units in length divided into
1000 computational cells. The initial and boundary conditions
are :
Initial condition :
c( 0 < x < 50, 0) = 0.0
c(50 < x < 150, 0) = 1.0 (4.32)
c(x < 150, 0) = 0.0
that is a square wave.
Boundary condition :
c(0,t) = 0.0
c(l,t) = 0.0 (4.33)
The advection velocity u = 1.0. The analytical solution to
Equation (4.1) is given by
c(x,t+ At)= c(x,t) (4.34)
All the numerical experiments were performed by setting the
CFL number to 0.5 because all the schemes discussed in
52
Appendix A except the ELM are unstable when the CFL number is
greater than 1.0. In the ELM test, the CFL number was set to
5.5. The numerical solutions are compared to the analytical
solution at time T = 50, 250, and 500 seconds.
Figures 4.3, 4.4, and 4.5 show the analytical solution to
Equation (4.1) when time T = 50, 250 and 500 sec.
Figure 4.3 Analytical solution of linear advection
equation at T = 50 Sec
Figure 4.4 Analytical solution of linear advection
equation at T = 250 Sec
o [7
Figure 4.5 Analytical solution of linear advection
equation at T = 500 Sec
53
Figures 4.6, 4.7, and 4.8 show the analytical and upwind
solutions of the 1D linear advection equation (4.1) at time
t=50, 250, and 500 sec. It is evident that there is numerical
diffusion in the numerical solution.
0 .
0 .6
0 .4
0.2 
.1 
x
Figure 4.6 Results of the 1D linear advection equation
obtained with upwind method at T = 50 Sec
o .6
0 .7
0 .3 
0 I ,
e 0.1 T0
Figure 4.7 Results of the 1D linear advection equation
obtained with upwind method at T = 250 Sec
07
00 .
0 .2 ,
Figure 4.8 Results of the 1D linear advection equation
obtained with upwind method at T = 500 Sec
Figures 4.9, 4.10, and 4.11 show the results of the
Leith's method. Although there is no numerical diffusion in
the solution,.excessive oscillation is evident.
V .
1o
0 4
A ___
1 .2
0. 2  0' 0
Figure 4.9 Results of the 1D linear advection equation
obtained with the Leith's method at T = 50 Sec
1 2 J t
0 2
V 0.0
:0
0. 
  @. 
obtained with the Leith's method at T = 250 Sec
.2 0 A
S.1  ,
0 1
0 .
o .1 
2050 000 000 1000
Figure 4.11 Results of the 1D linear advection equation
obtained with the Leith's method at T = 500 Sec
Figures 4.12, 4.13, and 4.14 show the results of the TVD
method. The numerical and analytical solutions are in very
good agreement.
Figure 4.12 Results of the 1D linear advection equation
obtained with the TVD method at T = 50 Sec
0 1y1 " 1I 0 o
o 7

0 .2
0 .1
 '  5i't ''  s t'l . ..  7 S't . ..1 F
Figure 4.13 Results of the 1D linear advection equation
obtained with the TVD method at T = 250 Sec
TV o mD th e
o .,
o 5
0.0 
0 .3
0 .2
Figure 4.14 Results of the 1D linear advection equation
obtained with the TVD method at T = 500 Sec
56
Figures 4.15, 4.16, and 4.17 show the numerical results
obtained with the ELM with linear interpolation and a CFL
number of 5.5. The numerical and analytical solutions are
clearly in good agreement.
E l.. i1. . .. ...
0 .9
t .o
S0.S
00 5
2505027 0 00
0.1 
Figure 4.15 Results of the 1D linear
obtained with the ELM with linear
at T = 55 sec
advection equation
interpolation
S .3 
.2 
  t'  s' T' ,,, rI ,
Figure 4.16 Results of the lD linear advection equation
obtained with the ELM with linear interpolation
at T = 253 sec
0 .4
S0.8
0 .1 
Figure 4.17 Results of the 1D linear advection equation
obtained with the ELM with linear interpolation
at T = 506 Sec
57
Figure 4.18, 4.19, and 4.20 show the numerical solution of the
Equation (4.5) obtained with ELM with quadratic interpolation
and a CFL number of 5.5. The numerical solutions show wiggles
in the vicinity of sharp gradients.
... ....
o .s
0 .
.101
250 010 750 000
Figure 4.18 Results of the lD
obtained with the ELM with
at T = 55
linear advection equation
linear interpolation
Sec
I 2 s
Figure 4.19 Results of the 1D linear advection equation
obtained with the ELM with linear interpolation
at T 253 Sec
*^ Unly ti .la.tn
)..9 
S0.0 
o.4 
o 0 000.
Figure 4.20 Results of the 1D linear advection equation
obtained with the ELM with linear interpolation
at T = 506 Sec
00 .7 
Figure 4.20 Results of the 1D linear advection equation
obtained with the ELM with linear interpolation
at T = 506 Sec
58
Figures 4.21, 4.22, and 4.23 show the numerical solution
of the conservative QMSL method. The numerical solution has no
spurious numerical oscillation and agrees well with the
analytical solution.
 o . i . .
5.5 1
.3 *
Figure 4.21 Results of the 1D linear advection equation
obtained with the Conservative QMSL at T = 55 Sec
0. ...... ,ll I,...
0,
0.4 4
2,0 5 0 0 o1000
Figure 4.22 Results of the 1D linear advection equation
obtained with the Conservative QMSL at T = 253 Sec
  t' ^ m '.  A
.0 .
.8i
0.40
!0 
0 20 0 00 0 00 1000
Figure 4.23 Results of the 1D linear advection equation
obtained with the Conservative QMSL at T = 506 Sec
4.5.2 Conservation Test
If one wants to simulate some conservative quantities
numerically, it is necessary to use numerical schemes that
possess the conservative property. Conservative property of a
numerical scheme is that conservative species satisfy the
integrated form of the conservative equation i.e.
+Af Jx=L(01 ~+1Cdxdt=0
(4.35)
Thus :
J x=L ut +rAti
x=L (+ c")dx+ [(uc) (uc)L]dt = 0
ax=O 
(4.36)
where subscript Rc and L, are right and left face of
computational cells as shown in Figure 4.24.
Figure 4.24 Computational cell used in conservation test
t "< I
60
If the convective velocity is constant, the convective
equation is the same in either conservative or non
conservative form. Hence, to check the conservation of
numerical scheme, the convective velocity should be non
constant.
The test problem for the lD linear advection equation is
the same as that discussed in section 4.4.1. However, the
convective velocity is nonconstant
x
u= 1+ (4.37)
L
where L is the length of the computational domain.
Figure 4.25 shows the conservation error of the
conventional ELM with linear interpolation, while Figure 4.26
shows the result of the conservative QMSL method. It is seen
that the conservation error of QMSL fluctuates around zero
with the amplitude on the order of 1013 while that for the
conventional ELM does not converge to zero. In both figures,
conservation error is defined as :
Conservation Error = Total Initial Mass Computed Mass
Figure 4.25 Conservation error of the conventional ELM
Figure 4.26 Conservation error of the QMSL
4.6 Momentum advection equation
4.6.1 NonLinear Advection Case
The EulerianLagrangian Method described in previous
sections have also been applied to strictly nonlinear
problems. The invicid Burgers equation is used for model
tests. The nonconservative form of the governing equation is
du du
+ u= 0
whereas the conservative form of governing equation is
I u2
+ 0
at 9x
(4.38)
(4.39)
5E12
1E 11 
100 200 300 400
Elapsed Tim
The initial and boundary conditions are assume to be
u(x,0)= (1 x) 0 x< 1
u(x,0)= 0 1< x 3 (4.40)
u(O,t)= 1 t> 0
The analytical solution of the problem, for t>0, is given
by :
u(x,t)= 1 for x t
u(x,t)= [(1x)/(1 t)] for t< x<1 (4.41)
u(x,t)= 0 for 1< x
At time t=l a shock develops in x=l and propagates to the
right at speed 0.5. Thus the analytical solution for t2> is
given by :
u(x, t) = 1 for x< 1+0.5(t 1)
u(x,t)= 0 for x> 1+ 0.5(t1) (4.42)
The spatial domain 05x< 3 has been divided into 150
intervals of Ax = 0.02. Nonconservative upwind, conservative
upwind, Leith's method, Godunov method, nonconservative and
conservative ELM have been implemented and compared to the
analytical solution.
The first scheme is the nonconservative upwind scheme.
The stability of the scheme is restricted by the CFL number.
The time step is set to 0.01 sec, with the max CFL number of
0.5.
63
Figures 4.27, 4.28 and 4.29 show comparison of the
numerical and analytical solutions when elapsed time T=0.5, 1,
and 1.5 second. The results show that the scheme is accurate
but there is some numerical diffusion in the numerical
solution. This scheme does not capture the shock, i.e., the
shock does not propagate after the shock develops.
Figure 4.27 Results of the nonconservative
upwind method at T = 0.5 sec
0 3.
0 2
x
Figure 4.28 Results of the nonconservative
upwind method at T = 1.0 sec
o 7.
0 6
V .
Figure 4.29 Results of the nonconservative
upwind method at T = 1.5 sec
64
Figures 4.30, 4.31 and 4.32 show numerical solution of
the conservative upwind numerical scheme. The solution is
accurate and capture the shock, but the propagation speed of
the shock is incorrect. The numerical solution is correct
before the shock is developed.
0. 
0.S
04
0.2
Figure 4.30 Results of the conservative
upwind method at T = 0.5 sec
I 3
00 ..6*
0.5
0 .2
0.
Figure 4.31 Results of the conservative
upwind method at T = 1.0 sec
5 
0X
X
Figure 4.32 Results of the conservative
upwind method at T = 1.5 sec
Figures 4.33, 4.34 and 4.35 show the results of the
Leith's method for the nonlinear advection equation. This
method is accurate until the shock develops but does not
capture the shock. Moreover the solution has unphysical
oscillation.
0 6
0.1
00
Figure 4.33 Results of the Leith's method for
nonlinear advection equation at T = 0.5 sec
Figure 4.34 Results of the Leith's method for
nonlinear advection equation at T = 1.0 sec
o .8
0.6
0 .4
o .3
Figure 4.35 Results of the Leith's method for
nonlinear advection equation at T = 1.5 sec
Figures 4.36, 4.37 and 4.38 show the results of the
Godunov method, which is a kind of TVD. It captures the shock
correctly, and the solution is accurate.
0 .' . . .. .
.3
0.2
S.1
0o
Figure 4.36 Results of the Godunov method for
nonlinear advection equation at T = 0.5 sec
0 4
.3 
S .2
o .I
Figure 4.37 Results of the Godunov method for
nonlinear advection equation at T = 1.0 sec
0 .4
0 .3
O !
Figure 4.38 Results of the Godunov method for
nonlinear advection equation at T = 1.5 sec
67
Figures 4.39, 4.40 and 4.41 show the results of
conventional nonconservative ELM with linear interpolation.
The maximum CFL number is 2.5. It captures the shock
correctly, and solution is accurate.
0 5 .,. ..
.o 4
0.
0.2
Figure 4.39 Results of the nonlinear advection equation
of ELM with linear interpolation at T = 0.5 sec
0 .
0 
I
of ELM with linear interpolation at T = 1.0 sec
0.7
0.4 
0.2 
1 =
Figure 4.41 Results of the nonlinear advection equation
of ELM with linear interpolation at T = 3.25 sec
68
Figures 4.42, 4.43 and 4.44 show the results of the
conventional ELM with quadratic interpolation. It captures the
shock correctly, but there are some numerical oscillation near
the shock.
0 .3
0 .2
x
Figure 4.42 Results of the nonlinear advection equation
of ELM with quadratic interpolation at T = 0.5 sec
0 .8
0 .2
0 .1
ox
x    
Figure 4.43 Results of the nonlinear advection equation
of ELM with quadratic interpolation at T = 1.0 sec
o.7
.4
o.3
o o
Figure 4.44 Results of the nonlinear advection equation
of ELM with quadratic interpolation at T = 3.25 sec
69
Figures 4.45, 4.46 and 4.47 show the results of the
conservative QMSL. It captures the shock correctly, and its
results are accurate.
Figure 4.45 Results of the nonlinear advection equation
I 
of conservative QMSL method at T =0.5 sec
.400
0 
0. 
o ,,
K, ,
Figure 4.46 Results of
of conservative
the nonlinear advection equation
QMSL method at T = 1.0 sec
Figure 4.47 Results of the nonlinear advection equation
of conservative QMSL method at T = 3.25 sec
o I
0 .9
08
0 .7
0.0
0.5
.4
0 .3
0 .2
0 .1
a
4.6.2 Momentum Conservation Check
To test the momentum flux conservation, the following
integral form of equation and conservation form of finite
difference equation are used.
r+A( 011 1 x=L ch11
dt + 2 dx=0 (4.41)
t< 2 '=0 ch
Flux difference in time should be Zero.
( Icells+l n+I Icells+l +l Icel ls+ n Icells+l ) n
i=0 U i+1) i=0 uu) i+1 i=0 i =0
(4.42)
The conservation error is defined as
Conservation Error = Total Flux Computed Flux = 0
Figure 4.48 shows the conservation error of the
conservative QMSL algorithm. The error converges to zero with
time.
C on vatio. I C on. rv II. EL
1.0005
1 .0004 
1 .003
i .' 0 0 2
1 .0 0 0 1
Figure 4.48 Conservation error of the conservative QMSL
4.7 TwoDimensional Conservation Test
The conservative and nonconservative EulerianLagrangian
Method (ELM) were tested in twodimensional space. The
physical domain is 200 km in X direction and 200 km in Y
direction. The computational domain is divided into 1 Km in X
and Y direction, i.e., 200x200 cells. The initial and boundary
conditions are :
Initial conditions
c (5 km < x < 25 km, 0) = 1.0
c (5 km < y < 25 km, 0) = 1.0 (4.43)
C = 0, elsewhere
Boundary conditions :
c(0,0,t) = 0.0
c(200 km, 200 km, t) = 0.0 (4.44)
The advection velocity is uniform, i.e., u = v = 30.5 cm/s.
Figure 4.49 is numerical solution of the nonconservative ELM
with time increment At = 1 hours, i.e, CFL number is 1.098
for 5 days simulation and Figure 4.50 shows the global
conservation error of the nonconservative ELM. The global
conservation error is defined by
Computed concentration Total concentration
Error(%) = X 100
Totol concentration
(4.45)
z
E X
L ,xI
X(Km)
Figure 4.49 Numerical solution of 2D nonconservative ELM
2E12
1E 1 2
0
1 E 1 2
2 E12
Figure 4.50 Conservation error of 2D nonconservative ELM
S20 40 60 80 100 12
Tim e(H ours)
73
Figure 4.51 shows the contour map of numerical solution
of nonconservative ELM. There are some numerical damping in
the numerical solution.
200
175
150
125
100
>
75
50
25
0
0
50 100 150 200
X(km)
Figure 4.51 Contour map of numerical
nonconservative ELM
solution of
Figure 4.52 shows the numerical solution of the
conservative ELM and Figure 4.53 shows the conservation error
of the conservative ELM. Figure 4.52 indicates the numerical
solution of conservative ELM has less damping and more
accurate than nonconservative scheme.
2D NonConservative ELM
,~.
!1] ., '',s
(4 days
 h*<
I 3 days
 2 days
Iday
SInitial
Condition
I I.. I , I
5 days after
' 1 i I
)"
Z
Co;.;1;"."1 [ Y

0
m
X(Km)
Figure 4.52 Numerical solution of 2D conservative ELM
2E12
1E12
0 o
0
1E 12 
2 E 12 ...
Figure 4.53 Conservation error of 2D conservative ELM
o ,
75
Figure 4.54 shows the contour map of numerical solution
of nonconservative ELM. The numerical damping of the
conservative ELM is much less than that of nonconservative
ELM.
)111f
175
150
125
E
.100
75
50
25
0 50 100 150 200
X(km)
Figure 4.54 Contour map of numerical solution of
conservative ELM
2D Conservative ELM
7
S5 days after
4 days
i 3 days
2 days
0:y
1 day
Initial
Condition
SI I I I I t I I
4.8 ThreeDimensional Conservation Test
The conservative and nonconservative EulerianLagrangian
Method (ELM) were tested in threedimensional rectangular
basin. The physical domain is 200 km in X direction,200 km in
Y direction and depth is 150 meters. The computational domain
is divided into 1 Km in X and Y direction and 15 meters in Z
direction, i.e., 100x100x10 cells. The initial and boundary
conditions are :
Initial conditions
c (10 km < x < 50 km, 0) = 1.0
c (10 km < y < 50 km, 0) = 1.0 (4.46)
c(20 m < z < 70 m, 0) = 1.0
c = 0, elsewhere
Boundary conditions :
c (0,0,0,t) = 0.0
c (200km,200km,150m, t) = 0.0 (4.44)
The advection velocity is uniform, i.e., u = v = 30.5 cm/s
and w = 0.012 cm/s. Figure 4.55 shows the numerical solution
of the nonconservative ELM with time increment At = 2 hours,
i.e, CFL numbers in X and Y direction are 1.098 and CFL number
in Z direction is 0.0576 for 5 days simulation and Figure 4.56
shows the global conservation error of the nonconservative
ELM. The numerical solutions of 3D conservative and non
conservative ELM is very close, but the conservation error of
77
the two different schemes is significant. Figure 4.56 shows
the mass of concentration in nonconservative numerical scheme
lose its original concentration with time.
3D NonConservative ELM
z
Figure 4.55 Numerical solution of 3D nonconservative ELM
1
0 .5
1 .5 
S 2 
2
2 .5
3 , I. I. . I
Figure 4.56 Conservation error of 3D nonconservative ELM
78
Figure 4.57 and 4.58 show the numerical solution of the
nonconservative ELM in XZ and YZ plane, respectively.
3D NonConservative ELM
X
150
125
75
100
50
25
X(km)
Figure 4.57 Numerical solution of 3D nonconservative ELM
in XZ plane
3D NonConservative ELM
z
Y
Y(km)
Figure 4.58 Numerical solution of 3D nonconservative ELM
in YZ plane
79
Figure 4.59 shows the numerical solution of the conservative
ELM. Figure 4.60 Shows the conservation error of the
conservative ELM. Conservation error of the conservative ELM
is almost zero.
3D Conservative ELM
z
kx
Y^ 
Figure 4.59 Numerical solution of 3D conservative ELM
2E .12
1 .8 E 1 2
w
1 .6 E 1 2
1 .4 E 1 2
1 2 3 4
T im e(D ays)
Figure 4.60 Conservation error of 3D conservative ELM
  o o  oex tot oo 3 l*
Figure 4.61 and 4.62 show the numerical solution of the non
conservative ELM in XZ and YZ plane, respectively. Figures
show conservative ELM has less damping than nonconservative
ELM.
3D Conservative ELM
2x
150
125
100
75
50
25
X(km)
Figure 4.61 Numerical
solution of 3D nonconservative ELM
in XZ plane
3D Conservative ELM
Y
50
25
D0
5 i
Y(knm)
Figure 4.62 Numerical solution of 3D nonconservative ELM
in YZ plane
4.9 Summary
All the above presented numerical schemes for the 1, 2
and 3dimensional linear and nonlinear advection equations are
relatively accurate. The conventional ELM is accurate but it
does not possess the conservative property. If a higherorder
interpolation function is used, there is spurious unphysical
oscillation in the numerical solutions. The total variation
diminishing (TVD) method produced the best results in terms of
both accuracy and conservation. However, the TVD method is
restricted by the CFL condition. The conservative quasi
monotone semiLagrangian (QMSL) algorithm is relatively
accurate, has the conservative property, and has no numerical
oscillations with the help of FCT algorithm. In addition, the
conservative QMSL is unconditionally stable.
Because of its numerous advantages, the conservative QMSL
method will be used for the present threedimensional model
discussed in chapter 3.
