• TABLE OF CONTENTS
HIDE
 Title Page
 Copyright
 Acknowledgement
 Table of Contents
 List of Tables
 List of Figures
 Abstract
 Introduction
 Governing equations and boundary...
 A three-dimensional numerical...
 Eulerian-Langrian (semi-Langrian)...
 Model verification with analytical...
 Comparison with CH3D
 Application of model to Lake...
 Conclusions and future study
 Appendix
 Reference
 Biographical sketch
 Signature page














Title: Three-dimensional conservative eulerian-lagrangian model for coastal and estuarine circulation
CITATION THUMBNAILS PAGE IMAGE ZOOMABLE
Full Citation
STANDARD VIEW MARC VIEW
Permanent Link: http://ufdc.ufl.edu/UF00091419/00001
 Material Information
Title: Three-dimensional conservative eulerian-lagrangian model for coastal and estuarine circulation
Series Title: Three-dimensional conservative eulerian-lagrangian model for coastal and estuarine circulation
Physical Description: Book
Language: English
Creator: Lee, Jun
Publisher: Coastal & Oceanographic Engineering Dept. of Civil & Coastal Engineering, University of Florida
Place of Publication: Gainesville, Fla.
 Record Information
Bibliographic ID: UF00091419
Volume ID: VID00001
Source Institution: University of Florida
Holding Location: University of Florida
Rights Management: All rights reserved by the source institution and holding location.

Table of Contents
    Title Page
        Title Page
    Copyright
        Page i
    Acknowledgement
        Page ii
    Table of Contents
        Page iii
        Page iv
    List of Tables
        Page v
    List of Figures
        Page vi
        Page vii
        Page viii
        Page ix
        Page x
        Page xi
        Page xii
        Page xiii
        Page xiv
        Page xv
    Abstract
        Page xvi
        Page xvii
        Page xviii
    Introduction
        Page 1
        Page 2
        Page 3
        Page 4
    Governing equations and boundary conditions in three-dimensional Cartesian coordinates
        Page 5
        Page 6
        Page 7
        Page 8
        Page 9
        Page 10
        Page 11
        Page 12
        Page 13
    A three-dimensional numerical model
        Page 14
        Page 15
        Page 16
        Page 17
        Page 18
        Page 19
        Page 20
        Page 21
        Page 22
        Page 23
        Page 24
        Page 25
        Page 26
        Page 27
        Page 28
        Page 29
        Page 30
        Page 31
        Page 32
        Page 33
    Eulerian-Langrian (semi-Langrian) Method
        Page 34
        Page 35
        Page 36
        Page 37
        Page 38
        Page 39
        Page 40
        Page 41
        Page 42
        Page 43
        Page 44
        Page 45
        Page 46
        Page 47
        Page 48
        Page 49
        Page 50
        Page 51
        Page 52
        Page 53
        Page 54
        Page 55
        Page 56
        Page 57
        Page 58
        Page 59
        Page 60
        Page 61
        Page 62
        Page 63
        Page 64
        Page 65
        Page 66
        Page 67
        Page 68
        Page 69
        Page 70
        Page 71
        Page 72
        Page 73
        Page 74
        Page 75
        Page 76
        Page 77
        Page 78
        Page 79
        Page 80
        Page 81
    Model verification with analytical solution
        Page 82
        Page 83
        Page 84
        Page 85
        Page 86
        Page 87
        Page 88
        Page 89
        Page 90
        Page 91
        Page 92
        Page 93
        Page 94
        Page 95
        Page 96
        Page 97
        Page 98
        Page 99
        Page 100
        Page 101
        Page 102
        Page 103
        Page 104
        Page 105
        Page 106
        Page 107
        Page 108
        Page 109
        Page 110
        Page 111
        Page 112
        Page 113
        Page 114
        Page 115
        Page 116
        Page 117
        Page 118
        Page 119
        Page 120
        Page 121
        Page 122
        Page 123
        Page 124
        Page 125
        Page 126
        Page 127
        Page 128
        Page 129
        Page 130
        Page 131
        Page 132
        Page 133
        Page 134
        Page 135
        Page 136
        Page 137
        Page 138
    Comparison with CH3D
        Page 139
        Page 140
        Page 141
        Page 142
        Page 143
        Page 144
        Page 145
        Page 146
        Page 147
        Page 148
        Page 149
    Application of model to Lake Okeechobee
        Page 150
        Page 151
        Page 152
        Page 153
        Page 154
        Page 155
        Page 156
        Page 157
        Page 158
        Page 159
        Page 160
    Conclusions and future study
        Page 161
        Page 162
        Page 163
        Page 164
    Appendix
        Page 165
        Page 166
        Page 167
        Page 168
        Page 169
        Page 170
        Page 171
        Page 172
        Page 173
        Page 174
    Reference
        Page 175
        Page 176
        Page 177
        Page 178
        Page 179
        Page 180
    Biographical sketch
        Page 181
    Signature page
        Page 182
Full Text




UFL/COEL-2000/014


A THREE-DIMENSIONAL CONSERVATIVE
EULERIAN-LAGRANGIAN 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 THREE-DIMENSIONAL 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 EULERIAN-LAGRANGIAN (SEMI-LAGRANGIAN) METHOD . . .. 34

4.1 Introduction . . . . . . . . . 34
4.2 Review of Eulerian-Lagrangian Method (ELM) . 34
4.3 Derivation of Eulerian-Lagrangian Method (ELM) 37
4.4 Conservative Eulerian-Lagrangian 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 Non-Linear Advection Case . . . .. 61
4.6.2 Momentum Conservation Check . . . .. 70










4.7 Two-Dimensional Conservation Test . . . .. 71
4.8 Three-Dimensional 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 Non-Linear 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 Density-Driven 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 1-D linear advection equation
obtained with upwind method at T = 50 Sec . . .. 53

Figure 4.7 Results of the 1-D linear advection equation
obtained with upwind method at T = 250 Sec . .. 53

Figure 4.8 Results of the l-D linear advection equation
obtained with upwind method at T = 500 Sec . .. 53

Figure 4.9 Results of the 1-D linear advection equation
obtained with the Leith's method at T = 50 Sec . 54

Figure 4.10 Results of the 1-D linear advection equation
obtained with the Leith's method at T = 250 Sec . 54

Figure 4.11 Results of the 1-D linear advection equation
obtained with the Leith's method at T = 500 Sec . 54

Figure 4.12 Results of the 1-D linear advection equation
obtained with the TVD method at T = 50 Sec . . .. 55


Figure 4.13 Results of the 1-D linear advection equation










obtained with the TVD method at T = 250 Sec. . . 55

Figure 4.14 Results of the 1-D linear advection equation
obtained with the TVD method at T = 500 Sec. . . 55

Figure 4.15 Results of the 1-D linear advection equation
obtained with the ELM with linear interpolation at
T = 55 sec . . . . . . . . . . 56

Figure 4.16 Results of the 1-D linear advection equation
obtained with the ELM with linear interpolation at
T = 253 sec . . . . . . . . . . 56

Figure 4.17 Results of the l-D linear advection equation
obtained with the ELM with linear interpolation at
T = 506 Sec . . . . . . . . . . 56

Figure 4.18 Results of the 1-D linear advection equation
obtained with the ELM with linear interpolation at
T = 55 Sec . . . . . . . . . . 57

Figure 4.19 Results of the 1-D linear advection equation
obtained with the ELM with linear interpolation at
T = 253 Sec . . . . . . . . . . 57

Figure 4.20 Results of the l-D linear advection equation
obtained with the ELM with linear interpolation at
T = 506 Sec . . . . . . . . . . 57

Figure 4.21 Results of the 1-D linear advection equation
obtained with the Conservative QMSL at T = 55 Sec . 58

Figure 4.22 Results of the 1-D linear advection equation
obtained with the Conservative QMSL at T = 253 Sec 58

Figure 4.23 Results of the l-D 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 non-conservative upwind method
at T = 0.5 sec . . . . . . . . . 63

Figure 4.28 Results of the non-conservative upwind method
at T = 1.0 sec . . . . . . . . . 63


vii











Figure 4.29 Results of the non-conservative 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 non-linear
advection equation at T = 0.5 sec . . . . .. 65

Figure 4.34 Results of the Leith's method for non-linear
advection equation at T = 1.0 sec . . . . .. 65

Figure 4.35 Results of the Leith's method for non-linear
advection equation at T = 1.5 sec . . . . .. 65

Figure 4.36 Results of the Godunov method for non-linear
advection equation at T = 0.5 sec . . . . .. 66

Figure 4.37 Results of the Godunov method for non-linear
advection equation at T = 1.0 sec . . . . .. 66

Figure 4.38 Results of the Godunov method for non-linear
advection equation at T = 1.5 sec . . . . .. 66

Figure 4.39 Results of the non-linear advection equation of
ELM with linear interpolation at T = 0.5 sec . .. 67

Figure 4.40 Results of the non-linear advection equation of
ELM with linear interpolation at T = 1.0 sec . .. 67

Figure 4.41 Results of the non-linear advection equation of
ELM with linear interpolation at T = 3.25 sec . 67

Figure 4.42 Results of the non-linear advection equation of
ELM with quadratic interpolation at T = 0.5 sec . 68

Figure 4.43 Results of the non-linear advection equation of
ELM with quadratic interpolation at T = 1.0 sec . 68

Figure 4.44 Results of the non-linear advection equation of
ELM with quadratic interpolation at T = 3.25 sec . 68


viii










Figure 4.45 Results of the non-linear advection equation of
conservative QMSL method at T = 0.5 sec . . .. 69

Figure 4.46 Results of the non-linear advection equation of
conservative QMSL method at T = 1.0 sec . . .. 69

Figure 4.47 Results of the non-linear 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 non-conservative ELM 72

Figure 4.50 Conservation error of 2D non-conservative 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 non-conservative ELM 77

Figure 4.56 Conservation error of 3D non-conservative ELM 77

Figure 4.57 Numerical solution of 3D non-conservative ELM in
X-Z plane . . . . . . . . . . 78

Figure 4.58 Numerical solution of 3D non-conservative ELM in
Y-Z 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 X-Z
plane . . . . . . . . . . ... 80

Figure 4.62 Numerical solution of 3D conservative ELM in Y-Z
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 Two-dimensional simulation of water surface
elevation of tidal propagation at three locations with
8= 1.
Line : Analytical Solution
Circle : Numerical Solution . . . . . .. 90

Figure 5.6 Two-dimensional simulation of water velocity
of tidal propagation at three locations with 8= 1.
Line : Analytical Solution
Circle : Numerical Solution . . . . . .. 91

Figure 5.7 Three-dimensional simulation of water surface
elevation of tidal propagation at three locations with
8= 1.
Line : Analytical Solution
Circle : Numerical Solution . . . . . .. 92

Figure 5.8 Three-dimensional simulation of water velocity
of tidal propagation at three locations with 0= 1.
Line : Analytical Solution
Circle : Numerical Solution . . . . . .. 93

Figure 5.9 Two-dimensional 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 Two-dimensional simulation of water velocity
of tidal propagation at three locations with = 0.501.
Line : Analytical Solution
Circle : Numerical Solution . . . . . .. 95

Figure 5.11 Three-dimensional simulation of water surface
elevation of tidal propagation three locations with
8= 0.501.
Line : Analytical Solution










Circle : Numerical Solution


Figure 5.12 Three-dimensional simulation of water velocity
of tidal propagation at three locations with = 0.501.
Line : Analytical Solution
Circle : Numerical Solution . . . . . .. 97

Figure 5.13 Two-dimensional simulation of water surface
elevation of non-linear advection at three locations with
0= 1 .
Line : Analytical Solution
Circle : Numerical Solution . . . . . .. 100

Figure 5.14 Two-dimensional simulation of water velocity
of non-linear advection at three locations with 0 =1
Line : Analytical Solution
Circle : Numerical Solution . . . . . .. 101

Figure 5.15 Three-dimensional simulation of water surface
elevation of non-linear advection at three locations with
0= 1 .
Line : Analytical Solution
Circle : Numerical Solution . . . . . .. 102

Figure 5.16 Two-dimensional simulation of water velocity
of non-linear advection at three locations with 0 =1
Line : Analytical Solution
Circle : Numerical Solution . . . . . .. 103

Figure 5.17 Two-dimensional 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 Two-dimensional 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 Three-dimensional model results of water surface
elevation at three locations with Coriolis effect when
0=1.
Line : Analytical Solution










Circle : Numerical Solution


Figure 5.20 Three-dimensional 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 Two-dimensional model results of water surface
elevation at three locations with bottom friction when


S= 1.
Line : Analytical Solution
Circle : Numerical Solution


Figure 5.23 Two-dimensional 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 Three-dimensional model results of water surface
elevation at three locations with bottom friction when


S=1.
Line : Analytical Solution
Circle : Numerical Solution


Figure 5.25 Three-dimensional 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 Non-dimensional comparison between wave profiles
as predicted by theory and numerical model of Wetting and
drying(time=0 1 /2) . . . . . . . . 137

Figure 5.33 Non-dimensional 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 Depth-averaged water velocity of propagation test
with CH3D when =0.501
Line : CH3D
Circle : Present model . . . . . . . . 143

Figure 6.5 Water surface elevation of non-linear 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 Depth-averaged water velocity of non-linear 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 depth-averaged 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 THREE-DIMENSIONAL CONSERVATIVE EULERIAN-LAGRANGIAN 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 three-dimensional conservative

Eulerian Lagrangian hydrodynamic model with wetting and drying

which is significantly enhanced over the non-conservative

Eulerian-Lagrangian model developed by Casulli and Cheng.

Semi-implicit finite difference method has been used for the

numerical solution of three-dimensional shallow water

equations. The advection terms and horizontal diffusion terms

in the momentum equations are discretized with the explicit

conservative Eulerian-Lagrangian Method (ELM). The barotropic

pressure gradient terms in the momentum equations and the

velocities in the vertically integrated continuity equation

are discretized with the 8-scheme, 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

Eulerian-Lagrangian numerical schemes do not have the

conservative property and hence cannot be used for practical

long-term simulations. The present model possesses the global

conservative property because of the use of the quasi-

conservative Eulerian-Lagrangian 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 five-diagonal system

of linear equations. The matrix is symmetric and positive

definite. This five-diagonal system of equations is solved

efficiently by the pre-conditioned conjugate gradient method.


xvii










If the Coriolis terms are considered, the nine-diagonal

matrix is positive definite but not symmetric, and is solved

by the successive over-relaxation (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 wind-driven 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 two-dimensional. With the advancement in

computer technology, more recent models are three-dimensional

and time-dependent.

One pioneering three-dimensional coastal and estuarine

model was developed by Leendertse (1967, 1975). The model uses

a time-explicit method to solve the three-dimensional

hydrostatic equation of motion in Cartesian coordinates. Sheng

(1983) developed a three-dimensional 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

vertically-implicit 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 vertically-implicit 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 state-of-art finite difference

numerical model which uses boundary-fitted 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 time-implicit method.

However it uses a non-conservative Eulerian-Lagrangian

method(ELM) for the advection terms in the momentum equations

and the linear advection-diffusion 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 three-dimensional model by improving the

Casulli and Cheng(1992) model by implementing a

conservative Eulerian-Lagrangian Method(ELM) and adding

baroclinic terms and salinity equation,

* Improve the Eulerian-Lagrangian method by implementing an

explicit scheme in X-momentum equation and an implicit

scheme in Y-momentum 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 three-dimensional shallow

equations which are momentum, continuity and salinity

transport equations are written for a three-dimensional

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 pre-conditioned

conjugate gradient method to solve system of equations

In Chapter 4, the Eulerian-Lagrangian method (ELM), or










4

semi-Lagrangian 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 wind-driven 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 three-dimensional 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 non-divergent, (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 right-handed Cartesian coordinate system (x, y,

z) with the origin at the water surface and z measured upward,

the equations are written as following non-conservative form

(Lamb,1945) :

Continuity equation,

du dv dw
+ +-- o (2.1)
dx dy dz

X-momentum 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)

Y-momentum 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)

Z-momentum 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 +
P-d +
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+38T-0.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.57x10-4 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 sediment-water interface

are given by specifying the bottom stress in a form of the

Chezy-Manning 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 THREE-DIMENSIONAL 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 two-dimensional, 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 semi-implicit 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 semi-implicit and

remains unconditionally stable and doesn't have undesirable

numerical wave damping.










15

In this chapter, the three-dimensional semi-implicit

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

8-scheme, 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 cell-centered 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 semi-implicit 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 yk-j, 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
v-v
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,k-1
Si,j,+ 2k-12


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 ij-2,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 AZ-2,1jk ,,k
--(1-0) 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 k-index 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 Eulerian-Lagrangian 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 Y-momentum equation is discretized implicitly (Sheng,

1983). Although the Coriolis term in the Y-momentum equation

is discretized implicitly, the actual computation of the term

is explicit because the velocity Un+1 is already obtained from

the X-momentum 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+ ,m-1 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 j-2 j-2


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+ ,M-1


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 IV-2 +,.- 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 X-momentum 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,M-I







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=M-l ij+,k ij+,k



At [
poy k+1 A,,j+,ki j+1,k

p0,AY Lk=m


inj,k j+YuM-1






Pij,k j+,)- A

P,,,) AzIj+~. ,

Pi'j,k)- Znij+Ym


n n
Pij+,M Pij,M
2

n _pm
Pi,j+I,M-1 ~ Pi,J,M-
2


,n ii "
i,j+1,mn+l i-.j,m+-
2
n n
tily+l-m -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, + aM-Y am_
a y AzA- + a _y+ a-M_


- 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]
j-Y2_









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 '+ )]A-Az' ]


S AzT AG g [( -n+ -

A- 6- AzA-G g- y q _, ]A- Az
At [[AG Ag [0( ,-)]. {A-Az}+ ,.,

S ['j+ A-G} [ 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 {d-IGi+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 re-arranging the
Equation (3.22) by moving "n+l" terms to the left-hand-side,
Equation (3.22) becomes

n", 2 2[ [(AZ)TA Z] (iij )i,

[(Az)A- Az]i/ (i," -


[(Az)TA-Az]. (,+:; +;)]

= o [(A) A-, n(AzAG1n
= [(A -'A (+,j-i,:;j-0 ]


At [([A(Az)A A'G y (Az)TA-1G]I.


Equation (3.22) can be re-written into a compact form :
d, "+1- S n+l S n+1 S .n+1 S qn+l n
i,-ij i-Y2,ji+l,j i-2 ,j i,J+ i- Si/J-1 = 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 Az-iA]
'ij 2 g' -12
d, = I+ S +SI +S" +Sn (3.25)
Sj i+Yj i-Yj i+ ,+Y (3.25)

q. = -At [(Az)T AG] [(Az)T A-'G
Ax 1,+ ,J i-2,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-,Ij-I 2 li-,j qi-lj+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 +1j-1 C8 ]i+l,j '9 li+1,j+l = RHS"


where


C" = FIP"n_ P"n
SiY-2 i- Y2,j-1

C" = -CM" + FIP" P" FIP" P",
2 I- ,j ij+ i-Y j i-j- J '- 2,J

C"~= FIpP" P"
3 i, j+ 2 i- 2,j+.)

C = HMn" + FIP" P" FIP" ,P"/ n l
4 i ij-Y i,j-Y2 i+ ,j-l1 i j-Y2 '- j-1

C"= 1+ CM" + CM" + HM" + HM"
5 < i+ ,j i-Xj 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 +2-1-

C"=- CM" FIP" P" FIP" /P"

C" = -+FIP".Y i,j
9 ij+ i+ y,j+1]


(3.27)









with

At2 At2 gA2
C= g-2 H= g 2 1 2
HJg Ay AxAy (3.28)
(3.28)
A tf
P= (AzA-AZ), 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'A-z is a non-negative number. Moreover the coefficient

matrix A is a tri-diagonal 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 tri-diagonal

matrix. The nine-diagonal matrix is positive definite but not

symmetric. The five-diagonal system of equations can be solved

very efficiently by a preconditioned conjugate gradient

method. For the nine-diagonal matrix case, the block tri-

diagonal algorithm or successive over-relaxation(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,k-2 AX
Azn v -AZn Vn
i,j+ ,k ij+y2,k i,j-y,k i,j-Y,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)= -xrAx-bTx+ 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+ -Ax-b (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 iJ-I


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
S-j2,, 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

pre-conditioner is simply given by



Mp() = p(k) (k) (k)
ij "i, i+ 2,Pi+,j i- Y2,j Ji-l


(. 9)


(k) (k)
i,j+ 21i,j+l j 2 i,j-1


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 Sijk-1
S D S+D -' _'
i,j,k Az i,j,k+ 2 Azn i,j,k-2 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 j-1,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 tri-diagonal 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
EULERIAN-LAGRANGIAN (SEMI-LAGRANGIAN) 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 Eulerian-Lagrangian

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 Eulerian-Laqrancian Method (ELM)


The Eulerian-Lagrangian 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 semi-Lagrangian 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 non-conservative form of momentum or transport

equations are non-conservative. 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 Eulerian-Laqrangian Method (ELM)


Consider the diffusion-free non-conservative 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 non-constant. 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 i-n


uAt









39

Cn+1 = n (4.5)



where c = c[(i-a)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 (i-a ) 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 (i-1) and


(i) is used to estimate Ca_ we obtain the first order upwind


scheme. If a quadratic polynomial fit is used to interpolate

between (i-1),(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 + PCik-n (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 non-uniform, 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 min-max 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- e-i(k))e-i(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 (i-n, i-n-1), 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

min-max property and hence it will introduce spurious

oscillation.



4.4 Conservative Eulerian-Lagrancian Method (ELM)



Eulerian-Lagrangian Methods described above do not

possess the conservation property because the finite

difference equations are obtained from the non-conservative

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 Quasi-Monotone Semi-Lagrangian (QMSL) scheme (Bermejo

and Staniforth, 1992) and Quasi-Conservative Semi-Lagrangian

(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-

Corrected-Transport (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 (U-nl 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 quasi-conservative 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), bi-linear(2D), and tri-linear(3D) interpolations.

The higher order interpolation functions include quadratic

interpolation function(lD), bi-quadratic(2D), and tri-

quadratic(3D) interpolation functions. The lower and higher

order interpolation function for l-D, 2-D, 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

Bi-linear Interpolation : 4 points used

U!' = (1- p)[(l- q)Utn,j + qU'_.. jm]


+ p[(1- q)U_-n-l,j-m + qUl-n-1,_-m, (4.28)


Bi-quadratic Interpolation : 9 points used










U t+U -
I,'


P(1+ p)[( (1+ q)UI -- + (1- q2 )Ui--l (1- q)Ui-n-2 -,j-z


+(1-p2)[ (1+ q)U,- + (- )U,'j-m q)Ui-nj-m+


(1- p)[( (+ q) -n -I + (1- q) -,,+ ,,, (1- q) Uf-n+,-,n+l]
(4.29)
3) Three dimensional case
Tri-Linear Interpolation : 8 points used

U+=I

(1- r){(l- p)[(l- q)Ut-[,j-m,k-n + qUtl-,j-m-1,k-n

+ p[(1- q) Ugll,j-,,k-n + q Ul1--1,j-m-1,k-n }

+r{(1- p)[(1- q)U-l,j-m,k-n-1 + qUil-,j-m-,k-n-]

+p[(1- q)U _,--1,j-m,k-n-1 + qUll- j-m-1l,k-n-l]} (4.30)

Tri-Quadratic Interpolation : 27 points used

UIt = (P2 + p)(q2 + q)@

11
U ^ k ,2 2

+ + ()Ui1-jk_1(1r + -
(r 2 + rU-1-1, j-m-1,k-n-1 -r2 ) ,'--,j-m-,k-. + (r 2 r)UiCl--,j-m.-,k-_n+l
2 2
+ -(p + p)(l- q2)-

1 1
(r2 + -j-,k-n-1 r2)Ui ,j-m,k-n + (r r)Ui ,j-,k-n+
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-,k-n+_
2 )i-__ (-nk-n 2+ r) ^ij-k-i }
1



+ -(- )(q2 q)
2
1 1







S(r2 + r)Uij-m,k-nI + (1- r 2)Ui+1,-1,k-n + -(r2 r)U/+,j-m-,,k-n+}
2 2
1 1



((r2 r+ r)U-+,j-l,k-I + (1- r2)Uilj-m + -k- 2 r)Ui-l+,j-m,k-n+
1



2 2
+ -(p2 p)(l-2 q).




2
1 1





(2 (r2 U-+1,-j-m+1,k-n-1 2)Ull+1,]j-m+I,k-n 2 2 r)1Ui-l+I,j-m+1,k-n+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 1-D 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 1-D 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 1-D linear advection equation
obtained with upwind method at T = 250 Sec






07
00 .



0 .2 ,-





Figure 4.8 Results of the 1-D 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 1-D 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 1-D 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 1-D linear advection equation
obtained with the TVD method at T = 50 Sec










0 1y1 " 1-I 0 o
o 7



-

0 .2
0 .1











---- ' -- 5i't -'' ------- s t'l . ..-- -- 7 S't --. ..--1-- F
Figure 4.13 Results of the 1-D 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 1-D 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 1-D 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 l-D 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 1-D 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 l-D
obtained with the ELM with
at T = 55


linear advection equation
linear interpolation
Sec


I 2 -s












Figure 4.19 Results of the 1-D linear advection equation
obtained with the ELM with linear interpolation
at T 253 Sec
*^----- Unly ti- .la-.--t--n
)..9 --










S0.0 -
o.4 -



o 0 000.





















Figure 4.20 Results of the 1-D linear advection equation
obtained with the ELM with linear interpolation
at T = 506 Sec
00 .7 --








Figure 4.20 Results of the 1-D 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 1-D 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 1-D linear advection equation
obtained with the Conservative QMSL at T = 253 Sec









--- ----- t' ^ m '. ---- A
.0 .
.8-i

0.40
!0 -




0 20 0 00 0 00 1000



Figure 4.23 Results of the 1-D 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 l-D linear advection equation is

the same as that discussed in section 4.4.1. However, the

convective velocity is non-constant

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 10-13 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 Non-Linear Advection Case

The Eulerian-Lagrangian Method described in previous

sections have also been applied to strictly nonlinear

problems. The invicid Burgers equation is used for model

tests. The non-conservative 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)


-5E-12
-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)= [(1-x)/(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(t-1) (4.42)


The spatial domain 05x< 3 has been divided into 150

intervals of Ax = 0.02. Non-conservative upwind, conservative

upwind, Leith's method, Godunov method, non-conservative and

conservative ELM have been implemented and compared to the

analytical solution.

The first scheme is the non-conservative 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 non-conservative
upwind method at T = 0.5 sec





0 3.



0 2

x


Figure 4.28 Results of the non-conservative
upwind method at T = 1.0 sec





o 7.

0 6
V .


Figure 4.29 Results of the non-conservative
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
non-linear advection equation at T = 0.5 sec


Figure 4.34 Results of the Leith's method for
non-linear 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
non-linear 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
non-linear advection equation at T = 0.5 sec










0 4
.3 -
S .2
o .I






Figure 4.37 Results of the Godunov method for
non-linear advection equation at T = 1.0 sec












0 .4
0 .3


O !


Figure 4.38 Results of the Godunov method for
non-linear advection equation at T = 1.5 sec













67


Figures 4.39, 4.40 and 4.41 show the results of


conventional non-conservative 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 non-linear 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 non-linear 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 non-linear 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 non-linear 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 non-linear 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 non-linear 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 non-linear advection equation
QMSL method at T = 1.0 sec


Figure 4.47 Results of the non-linear 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 Two-Dimensional Conservation Test

The conservative and non-conservative Eulerian-Lagrangian

Method (ELM) were tested in two-dimensional 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 non-conservative 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 non-conservative 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 non-conservative ELM


2E-12



1E -1 2



0



-1 E -1 2



-2 E-12


Figure 4.50 Conservation error of 2D non-conservative ELM


S20 40 60 80 100 12


Tim e(H ours)











73

Figure 4.51 shows the contour map of numerical solution

of non-conservative 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
non-conservative 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 non-conservative scheme.


2D Non-Conservative 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





2E-12


1E-12


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 non-conservative ELM. The numerical damping of the

conservative ELM is much less than that of non-conservative

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 Three-Dimensional Conservation Test

The conservative and non-conservative Eulerian-Lagrangian

Method (ELM) were tested in three-dimensional 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 non-conservative 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 non-conservative

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 non-conservative numerical scheme


lose its original concentration with time.


3D Non-Conservative ELM


z


Figure 4.55 Numerical solution of 3D non-conservative ELM









-1
-0 .5




-1 .5 -
S 2 -
-2


-2 .5

-3 -, I. I. . I


Figure 4.56 Conservation error of 3D non-conservative ELM













78


Figure 4.57 and 4.58 show the numerical solution of the


non-conservative ELM in X-Z and Y-Z plane, respectively.


3D Non-Conservative ELM


X-
150


125




75
100





50


25


X(km)


Figure 4.57 Numerical solution of 3D non-conservative ELM
in X-Z plane


3D Non-Conservative ELM


z


Y


Y(km)


Figure 4.58 Numerical solution of 3D non-conservative ELM
in Y-Z 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 X-Z and Y-Z plane, respectively. Figures


show conservative ELM has less damping than non-conservative


ELM.


3D Conservative ELM


2-x
150

125

100

75

50

25


X(km)


Figure 4.61 Numerical


solution of 3D non-conservative ELM
in X-Z plane


3D Conservative ELM


Y
50


25

D0

5 i


Y(knm)


Figure 4.62 Numerical solution of 3D non-conservative ELM
in Y-Z plane












4.9 Summary



All the above presented numerical schemes for the 1, 2

and 3-dimensional linear and nonlinear advection equations are

relatively accurate. The conventional ELM is accurate but it

does not possess the conservative property. If a higher-order

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 semi-Lagrangian (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 three-dimensional model

discussed in chapter 3.




University of Florida Home Page
© 2004 - 2010 University of Florida George A. Smathers Libraries.
All rights reserved.

Acceptable Use, Copyright, and Disclaimer Statement
Last updated October 10, 2010 - - mvs