Citation
Three-dimensional conservative eulerian-lagrangian model for coastal and estuarine circulation

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
Creator:
Lee, Jun
Place of Publication:
Gainesville, Fla.
Publisher:
Coastal & Oceanographic Engineering Dept. of Civil & Coastal Engineering, University of Florida
Language:
English

Record Information

Source Institution:
University of Florida
Holding Location:
University of Florida
Rights Management:
All applicable rights reserved by the source institution and holding location.

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 Kij in 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 f ind 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

. . . . . . . ii

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

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 oacre
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 I-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 I-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 I-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 I-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 nonconservative 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 nonconservative 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
O= 1.
Line : Analytical Solution
Circle : Numerical Solution 90
Figure 5.6 Two-dimensional simulation of water velocity
of tidal propagation at three locations with O= 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
O= 1.
Line : Analytical Solution
Circle : Numerical Solution 92
Figure 5.8 Three-dimensional simulation of water velocity
of tidal propagation at three locations with O= 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
O= 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 = 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
O= 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 =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 =1I.
Line Analytical Solution




Circle :Numerical Solution14

Figure 5.20 Three-dimensional model results of water velocity
in X direction at three locations with Coriolis effect when

06=1 .
Line :Analytical Solution

Cicl Nue.a 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

06=1I.
Line :Analytical Solution Circle :Numerical Solution12

Figure 5.23 Two-dimensional model results of water velocity
in X direction at three locations with bottom friction when

0 =1I.
Line :Analytical Solution Circle :Numerical Solution

Figure 5.24 Three-dimensional model results of water surface
elevation at three locations with bottom friction when

06=1I.
Line Analytical Solution
Circle Numerical Solution

. . 122

. . . 123

Figure 5.25 Three-dimensional model results of water velocity
in X direction at three locations with bottom friction when

06=1I.
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




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 )T/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 0 =1.0
Line CH3D
Circle Present model ...... ................ 141
Figure 6.3 Water surface elevation of propagation test with
CH3D when 0=0.501
Line CH3D
Circle Present model ...... ................ 142
Figure 6.4 Depth-averaged water velocity of propagation test
with CH3D when 0 =0.501
Line CH3D
Circle Present model ...... ................ 143
Figure 6.5 Water surface elevation of non-linear test
with CH3D when 6 =1.0
Line : CH3D
Circle : Present model ...... ................ 144
Figure 6.6 Water surface elevation of Nonlinear test
with CH3D when 6 =0.501
Line : CH3D

xiii




Circle :Present model14

Figure 6.7 Depth-averaged water velocity of non-linear test
with CH3D when 06=0. 501
Line :CH3D
Circle :Present model.................146
Figure 6.8 The test of Coriolis effect with CH3D when 06=1.0
Line :CH3D
Circle :Present model.................147

Figure 6.9 Comparison of water surface elevation with
Coriolis ef fect with CH3D when 6 =0. 501
Line CH3D
Circle Present model .............

Figure 6.10 comparison of depth-averaged velocity with
Coriolis ef fect 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 e-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 quasiconservative 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 onedimensional 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(a -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 threedimensional 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 U -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 Obiectives
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,
9 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,
0 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 THREEDIMENSIONAL 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
+ +- 0 (2.1)
dx dy dz
X-momentum equation,
du d u du du 1 dp d 2u d 2u
X- -A(T 2)
dt dx dy dz po dx dx + y
d du
+ -(AV )+ (2.2)
Y-momentum equation,
dv v + v dv 1 dp + 2v 2v
-y+ u- + v- + w-z ---+ AH( + )
dt dx dy dz po dy dx2 2
d dv
+ z(Av-)-fu (2.3)
Z-momentum equation,
-p= -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




Vph = gp(xy,,t)V (xyt)+ gV p(xy,,t)d7
VPh = gp(X,Yliqt)Vq1(Xlylt)+ gJVp(x,,,t)1

(2.5)

Finally the momentum equations in x and y directions are

du du v-+ wdy dz

g d +
fo d + po dx

d 2 2 2H 2T
A (- + 2-- )

dv
dt

dv Ux+
dx

dv vd+ dy

g pd 2v d 2V Pd +AH(x2 +2 po dy x y

Salinity equation,

dS aS +u +

V S +wS = 2S d 2S v +w-z =1(T + 2) dy dz dx2 y

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, q(x,y,t) is the free surface elevation, g is the

du
dt St

du U-+
dx

8 du + -(A, -)+ fv
dz vdz

dv 89 w g y

(2.6)

d dv
+ (AV )- fu dz dz

(2.7)

d dS
+-(D -) dz v dz

(2.8)

p = p(T,S)

(2.9)

dq
-gx




8
gravitational acceleration and A,, DH 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, p0 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.1OT)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 = po(1+ aS) (2.11)
where a is a function of temperature, and for a temperature of 15C, a is 7.57xi0-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
w w
by the prescribed wind stresses Tx and Tw
2 2
pATz T x pCda Uw +Vw (2.12)
dv w 122
pAvl PCdaVw w +WVw (2.13)
Wz w
where Txw and Ty are the wind stresses in x and y directions
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, FT :
Cda = 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




Dij d) diq d) l)(215 w-= +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 :
+ vdz 0 (2.16)
dt dx dy h
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)+ ql(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 b
pAVz = x yu (2.17)
Iv
d pAz Y (2.18)




ii
g U2 + V2
where u and v are velocities at the bottom
layer, and C2 is the Chezy friction coefficient which can be formulated as
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 (rx, ,)= PCdb(112 + V2 (1,,V) (2.20)
where Cdb is the bottom friction coefficient, and tj and V1 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)
Cdb = k In (2.21)
L
where k is the von Karman constant.




12
The open boundary condition for the surface elevation is prescribed by
q(x,y,t)= -An cosr + 0 (2.22)
where An, 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)
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 yz, 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 l-H 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 G is introduced such that for 0 = 1, this method is fully implicit and for 0 = 1/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 8scheme.
3.2 Finite difference form of momentum ecquations
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 e-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 q 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.
n+ j A t n+1 n+l in ,j in
g At L Az (p",j,,- pi, ) Az Pi, ,m- Pi ,m
poAX m= k i-y2, j i t+ ,j 2
n+1 n+1 Un+1 n+1
A i+ '2,j,k+1 i+yjk i+ j,k i+yj,k
Vk+J4 A k-)E AZ+k2
A V k + 2 i + j ,k + Y 2A 2 A i + j ,k Y
+At
(+ ,j,k

vn
+ f,+ j,k

(3.1)




n+1 = nAt [ n ++l )+qn ]
i,+Y,k i,= Fvj+Y2,k I ,+ 0- )X ,";+ _,")

A t n pn.
M- '
- g t Az"2 (p"2i p" ,,) -Az
o 7 m= t~gm +1, ij~r A ;gm

n+1 n+1
v-v
i,j+ Y2,k+1 i,j+ 2,k
Ak+Y Az ij+ ,k+Y -_

+At

- Un+k
i, j,+ ,k

n n
ij+,m Pi,jm
2

n+1 n+1
v -v
i,j+ 2,k i,j+ 2,k-1
Si,j,+ 2k-12

Azl,+>k

(3.2)

n+1 n t M n+1 M n+ 1
- kn ++ ,n k +1 ,k k= z Un j'-i
7, = rl Ax 0.L m A ~jk i+Ak ,- i-), ,Ii-2,Ikj
AX [k=n i k-iM 2
- t(-6 Az. u"kV -/ Az u"- ki~-/,
At L M n+1
A I Z Az+)~12,k tij+ 12k Zijy, A11 -Y2I
Ay t[k=,,I ,I/2 2, k=m
A t (1 Az "
A = Az ij+ y2,k i+y,jk k=m ij- 2,jk i-,jk
--(1-0 k=in
At n
---(1 0)[ A z. ,ki+ Ui,k -I Az, -, it" ]
Ay /___ki 2'k i'- 2
(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 Xmomentum 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 M+ A i.%.+Y2 X

A V .'j ,-M

n+l un+l
ui+Y2,j,M+1- u+ 2,j,M
Az.
i+ 2,j,M+12
n+l n+l
V ,J +12M Y2ij= y.
AzijM+ 1

n+1 n+1
SJm u ym-1 n+ / n+1l
Az i+ 2,j+,,m 2 i+2,,m i+ 2,,m

(3.4) (3.5)

(3.6)

(3.7)

n+1 n+1
iY +2,m ij+ = n+2, n+1
Avidj+gm- ij+ 2,m i,j+ 2,m
i,1+ 2,m- 2

n+2 = (n+ 12 n+ ( 12

n+2 n At n
i+ 2,j,m i+ 2,j,n 2Ax qi+l,i
n+g2 Vn At ( .
i,j+ /2,m i,j+ 2,m 2Ay "ij

and

where

and

(3.8)

(3.9.)
(3.10)

qij




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
At r ,+
A" U"+' = G" g--O(LJ -+ ")AZ.+ + AtfV (3
i+ ,j i+ ,j ,i+ ,j AXj i ,j i+/,j f +2/,j
A "+ = G"n- A ( n+I )AZ fUn+ (3.12)
.".= G+ [ (Z U AZ" ) U+
.1j J A i+/, i+Y2,j i- ,Y
At
- AZ" V. 1 A (jy V" (3.13)
we U Ay id i aj jwhere U, V, AZ, G, S AndA are defined as follows




Un+I ,n+I Un+I ui+,/nl

Sn+1 i,j+Y2, M
Vn+1 iQ+Y2,M-1

SAZ=

n+1
, n+1 i,+Yr ,m

AzM AzM-1

.Azm )

A ,t 0 ) ,, ]
AzM Fut" +M A) + C" + BC" + At
F At 1
Az Fu," g x(1- )( i I' ,)+ C" +BC".M
[ I At ./ MI]
AZM_1 [g.,+ ...,~~- g-~ o IV27+.-, ,J)+c ,M..,+ l + .1.+Y g/

Az..i[ FU" Az,,, Fu"'
I Az 2,jm
Az,, [Fit",4

At
At 1 ") + C" + BC, 1 Ax +Y
At
Ax g )(," -,) + C+B

At
AzMI FVI, g- (1- O)(q,"j+, qi,j) + BC M+ AIt t
A, 1,j+Y2- m
At
Az,,f1 Fvn g(1- ) i) + BC"
[ At
Az,+, Fvn -+ g -0 (- 0)(qfi+1 G i)+ BC m+
[~ At
Azm Fv",; g-0(1-0)(q";+1 i) +BC', + m

where C's are Coriolis terms in the X-momentum equations

(3.14)

(3.15)

(3.16)

vn+
i,j+ Y2=

G"Y2 =

Gn
I, =-




Cn1, =

AtfV+Y,j,M AtfV ,M-

A v ,", +,.,, l Ay te+ jm+
and BC's are the baroclinic terms

Az (p' p,, ) Az
i+A,j,k i+1,j,k ,k i+ ,jM

At
- Z Y2 (PinIjk P1,jk) j
At z (pil P,,J)- Az" 2
- o k= M1 i2 ,jk i+I,j,k j,k)- i+)j, M-1
PA k= +1 + ,2~ j, k n~ +2, jlm+1

AtFm
gT
PoXk=.

nAz+ ki k Pnjk AZn
Az+y2,j,k (Pi,+l,,,k ij,k) AZ+y2,j,,.

n1
Pi+I,j,M Pi,j,M
2
n n
Plj,- Pj,M-I
2
n n 1
Pi+l,j,m+l Pi,j,m+l
2 n
,l n.
Pi+l,j,m Pi,],n
2
(3.18)

(3.17)

AtF
-g[ k=MA

BC"
i+Y2,jlk




At [
- o k
p0AY kA

AZij+ ,k Pi,j+l,k i,j,k) ij+_,nM

-g P--'At [k=Mm-+ A i,j+2,k(Pi,j+l,kAt m
-g 1m Az.J"k (Pj+lk
poAy L k=M- i j+ ,k inj+1,k
At[ M
oA L k i+1 j+ ,k j+1,k gmA j+ k( i j+1,k

ij,k)- A+2.P,,j,k) Azj+,+ i j,k j+,m

n np
Pij+1,M Pi,,M
2
n pm
PiJ+IM-1 Pi,j,M2

,n ii
9ij+1,m+1 -- Pi.j,n+1
2
n n
Piylm- PijmI
2
(3.19)

)[ Z T+ l2,j TU ,
itj i+ 2,j

At F
--(1- ) Az )
Ayij+

- (Az U"j T ]
a- ,j,j

T
V" Az,
i,j+ 1,1- 2

SAzAI + a,_ a Y
- a f AzA,- + amY + aMY

- a My

- am+y Az, + am+y + rAt)

(3.21)

k At
with ak = Azk

BC =
i, j+Y2,k

At
3". = (11,j by A

and

(3.20)

V"
1.1Y




Formal substitution of Equations (3.1) and (3.2)into Equation (3.3) yields :
qn+1 (n
,fj ,j
At [AAZTG g A 8(i,7 I 71+ )]A-'Az} ]
S AzTY AG} g (1 ,n+l -A Zn
A ~ ~ g__ yn A [ O(77i,.+,_ _,,,I -A 1
- Iz A -g A (n J ]) A 2
- At [[A {A 'G}t g 0( ,+, -,,l) {A Az}i+ ,j
+ L j+{A2 AG} g- [ ( +",1 4,j+ )]A A'z 2
+ AI { G}- g [(I +l-j- q+l) A AAz
+ i G- g ( "j )]A i,j, { AzY2
At rO ,.+1 ++I -n+1Ii
- jzy 1 Al-G g- ( ) Az
- [A i {A G}n g ( i +1 )]A A i2,j+
+ AAt [AG g (n 1+ -1 -A'Az1
I+ A7 {A G}i2 g)[e(l, -l, A jJ)] A i+'Az ,
+ A A'G g [( n+.1 +1 1 -'AAz }
(3.22)




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+1 g 2 A[ [(AZ)A Z]i (iq, n+l [(A ) ),_,j~
- (Az A-'Az]i (n,"-i
- y At 2 2TA- IAz (,j ', ',j )
Ayj + 1(3.23) [(Az) A Az.. Qnl, + ;,)]
-- Ii Y2 i ( lf qli,j .i,.=" At 0 [(A Z) AG] [(Az)A-'GA ]
A x Li+ ,j i ,j
At 0~ [[(Az)T A-'G.y (Az)TA-1G] ]
Equation (3.22) can be re-written into a compact form :
d ,"f+'- S n+1 S n+1 S .n+1 S qn+l n
i j i Y2,j i+,j iY2,ji-lj iJ+ 2i,J+ iJ- 2iJ-1 = qi,j
(3.24)
where
S" At2 n T 1 n
",= g Ax [(Az)TA 'AzYj
12n
Sn = g 9 2[(Az)T AzA+n
i,- 2
d, = I+ S"n +S +S"n +n S(3.25)
j +Y ,j -Y ,j i+ +Y (3 .25
q". = 8.n. At [[(Az)T AG] (Az) A-'G ]
Ax ,i+ ,J i-2,J
At 0[(Az) AG y [(Az)TA-'G],
Ay Q+ i j 2




26
If Coriolis terms are considered, Equation (3.24) becomes more complicated :
C +n 1 1j1n l j n l in1n+1 n n l
G in-1,I q q,-Ij -q -i- 7nl +l,~ j(3.18)
_-, ~ C ",7 ~ _ n q n + l - ( n n + 1 l? ~ S 1 ni+l,j-1 8 i+l,j +,j+1 = RHS"
where
C" = FIPn P"n
Q i,-2 i-Y2,j-1
C" = -CM" + FIP" P" FIP" P1
2 i-Yj i,j+ i- Y2,j ij-2 i-2,j
C=-FIpn +pn"
3 i, j+ 2 i- 2,j+)
C= HMn + FIP" P" -FIPn" 2Pn
4 i, j-Y i, j-Y i+2, j-1 ij-Y i-2, j-1/
C = 1+ CMn + CM" +HM" +HM"
C5 L i+Y2,j i-Y2,j+H idj+Y2+ i,j-Y
- FIP Pn FIPn pn FIP Pn FIP P1
i,j+Y2 i+2,j i,j+Y2 i-2,j ij-2 i+ Yj ij-2 i-2,j]
6 1,j+ i,j+Y2 i+y2, j+1 ij+ y i-2,j+1
C7 = FIP" P"
7 i,j- 2 i+ y2,J -1
C"(= CM" + FIP P -FIP" P"
C" = FIP +P" .
\i, j+)/ it y, j+1

(3.27)




with
At2 At2 2 gA2
C g-x2 H= g 2 2
SAy AxAy (3.28)
P= (AzA-'AZ), M= (AzA-'AZ), F-= At
4
and
RHS" = 6," 6I 0(Az) A-'G]". [(AZ) AG]i
lAX k 2J "2
- At [(Az) A-'G (Az)T A_G]n]
Ay ij+ i,
+ {A A{:zTA} j+j -'G}>;+, + +{AzA} j+ { A-G}1,+ ]
+ Atf At O[{AzTA}1 y {A'G} + {AzTA}' 2{A_'G},J
4 Ay
+{AzTAI'A A G,_ + {Az A} { A-'G} I
- ~Q 8AzA)_A 2G AzA 1, A-Y2(
+A zAZ _- A 'Gn +AzT A'G-1
(3.29)
Since the coefficient matrix A is symmetric and positive definite, A-' is also positive definite, and therefore
(Az)TA-IAz 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 ninediagonal 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 tridiagonal 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
= Ax
,n+1 n+1 1+ ,j,k i+ xjk i- j,k i- ,j,k
i,j,k+y ij,k-Y2A
Azn v -AZn Vn
i,j+ ,k i, j+ y,k i,j-(,k i,j-3,k
Ay (3.27)

k= m,m+ 1,...,M




3.4 Conlugate 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 papular 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
I
f(x)= xT Ax-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>O, 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




a
d f (x)
d
f'(x) Ax2f (x) (3.29)
d f (x)
With equation (3.29), the gradient of equation(3.28) is
1 1
f'(x) = ATX+ I 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)
](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 Sn
n+I i+12,] an n+1 S 2,J n n+l
" d d+i,j
S"1 Sn a
Sd d
Ili 1+111 "'J 2, I n n+ _' .+l S 2 d n n+1 qi
Id:d:+ aw'J' + d -,ad,,:_, ,j~
,Jn nliJ-I -n j n+ 1
which, by changing the variable 1 to e,j is equivalent to
ij a+ e+l, -a e,_,j a ye,,+l a,,j e,,j_ = b, (3.34) where
Sn Sn
'-2,,J ai,j+-2 qi,,
a i+ij /iji+l ia,j+ n d~d~+ bi'j n,
ityd d" d Vd d"jd
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 (0)
i,J
(2) Set




(0) (0) e(0) (0) (0)
ij ij ij i+ 2,j i+ ,j 2,j -1,j
(0) (0) b
- a i,+Y ,j+1 a1,j_ i,j- ij

(3.35)

(3) Then for k=0,1,2, .. and until (r(k),r(k)) < calculate

()(k) (k)
where a (k) = ((),Mp(k))
(kt), Mp(k))

(k+1) (k) (k) Mp(k)ij
r.. ri. (k (li ,()..

p(k+.) (k+1) + p(k)() where (k) ij =j + f ij h

((k+1) (k+1))
(k),r(k)

In equation (3.37) each element of the vector Mp which is pre-conditioner is simply given by
Mp k.) = p(k) ak)a (k)
1,j = "i, ai+,l2,]l"Pi+l,] i 'l2,j pi-l, ..

3.3 9)

_ -(k) ( a () i,j+ 21i,j+1 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)
i,yj = i,j 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 DSijk+1 i,j,k ijk i,,k-1
i,j,k AZ vij,k2+ AZ i,j,k-2 Azn k
Ij,k i,j,k+ ijk=FS,"jk + DHAzik S-1,,k- 2S,jk + Si+ljk (3.40)
D+ t Azi k ( -. )
-DHt \i,I-1,kn 2S~ik + S Ilk
+2 i j-I,k 2 i,j,k i, j+I,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 S, 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-Lagrancian 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 EulerianLagrangian 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-La rangian Method (ELM)
Consider the diffusion-free non-conservative advection equation
d+ 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 onedimensional 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 dx8 Dt dt dt dx

(4.3)

and the definition of the velocity,

dx
- U

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

-#

Figure 4.2 An Eulerian Lagrangian Mesh

A finite difference scheme for Equation (4.2) is simply

I-a i-n

uAt




39
n+1= cfl (4.5)
where C~a = 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 Cina" The stability, numerical diffusion and unphysical oscillations of (4.5) depend on the interpolation formula chosen. If a linear interpolation between (i-i) and
(i) is used to estimate Ci-a' 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 C7a 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,

0 p<1. In this case Equation (4.5) becomes




40
k+1 k k k
C. C p( c__) (4.6)
or, equivalently
Cik+1 I- p)in + PCik 1 (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 0p 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 tn+. 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-1 = xs_ ruk (xs), XN = Xi
s NN-, (4.8)
s= N,N- 1,N- 2,...,2,1




where uk(xs) is interpolated with any interpolation formula.
x. x
Then, at x, a can defined by aAx
The finite difference equation (4.6) possesses the minmax 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 k A e (4.9)
where A f is the amplification factor, I is the grid point index and i=". 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
IAk12 = 1- 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-i), 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 82C
S- p(p- 1)-T+ HOT (4.12)
Dt 2Az'
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 this case Equation (4.5) becomes (Casulli, 1987)
k+1 i(p + p)c,+1 + (1- p2) 1 2
Ci 2 ( +A '- f2




Note that, if jai: 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-Larangian 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+ aVu= 0 (4.14)
To obtain a higher order monotone numerical scheme, the FluxCorrected-Transport (FCT) method (Boris and Book 1973, Zalesak 1979) is used. Following the FCT method, the solution is u n+I Tn+I pn+1 4 5
L,k + F k(4.15)
where U+I 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 L k
is the low order approximation at the departure point at t, and Fn+1 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',,
l'n+l
2) Compute a higher order approximate solution UH,k by a standard ELM of order higher than 1, 3) At any instance t"+l set

Fk nl n+1 (U U )n+ -' k (t"H,k --"L,k

(4.16)

where (k 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(U ',, Uk2, U3, Uk4)

(4.17)

and U as the minimum U value of the neighboring points
i.e.

U- = min(Ul, Uk2, U;3, U;4 )

(4.18)




4.2) Set
+ = U+ Un+
Q L,k
Q-= U--Un+ (4.19)
P = U n+ U n+
"H,k "L,k
4.3) If P>O,
C = min(1, )
If P<0,
k = min(1, (4.20)
'P
If P= 0,
U n+1 T]n+l qn+
U = Ui ,i. e, = 0
Hk UL,k 'k 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
U = akUHk + (1- ak)ULk (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 ak. The procedure to define ak is as follows (Priestley, 1993




and 1994) 1) Define

f U(x)dx = Uodx = C

, = (UH,k UL,k)Sk (4.22)
where Sk is the area associated with node k 2) Compute
( -u s U -1usc (4.23)
Zak kHUk)Sk -UL,kSk=C (4.23)
Assume that
kmax (4.24)
If this is not the case, then the definitions are changed to Ak A (4.25)
C* C*
3) determine 'k
Step 1
If <0, then a nax. Otherwise k=
If A:5 0 then a,= a i .Otherwise a k=-- 0

Step 2 : Define surplus
SURPLUS = C* Zakk

Step 3 : Define the average value of a as :




SURPL US
aa vg
Step 4
max max
If aavg< akm then ak :avg Otherwise ak =ak
The lower order interpolation functions include linear(1D), bi-linear(2D), and tri-linear(3D) interpolations. The higher order interpolation functions include quadratic interpolation function(lD), bi-quadratic(2D), and triquadratic(3D) interpolation functions. The lower and higher order interpolation function for I-D, 2-D, and 3D cases are : 1) One dimensional case Linear interpolation : 2 points used
U'= (1- p)U_ + pU1_,,_ (4.26)
Quadratic Interpolation : 3 points used
U'+ = 0.5(p2 + p)Ui_,_z + (1- p2)U1_, + 0.5(p2- p)U_, (4.27) 2) Two dimensional Case Bi-linear Interpolation : 4 points used
U!.' = (1- p)[(1- q)Utn,j-m + qUtn,jm_1]
+ p[(1-q)Ul-n-l,j-m + qU-n-1,_-m-] (4 28)

Bi-quadratic Interpolation : 9 points used




U t+U
I'
(1+ p)[( (1+ q) -n-,j-m- + (1- q )Ui-,-l (1- q)Ui -n- ,j-"z+l]
+(1- p2)[ (1+ q)U-n, + (1- q)U,' -,j-m q) i-n,j-m+ ]
2_-inj-m 2
-_ (1- p)[(2(1 + q)U_ i-n+,-,n-I + (1-q 2) -,,1+,-,,, ( q) U,-+,, -,,, +l]
(4.29)
3) Three dimensional case Tri-Linear Interpolation : 8 points used
U" =
i,j,k
(1- r){(1 p)[(l q)Utl,j-m,k-n + qUt-1,j-m-1,k-n
+ p[(1- q) Ull1,j-m,k-n + q U-1-1,j-m-1,k-n }
+r{(1- p)[(1- q)Ul,j-m,k-n-1 + qUl-1,j-m-1,k-n-1
+ p[(1 q)U,--1,j-m,k-n-1+ qUllj-m-1l,k-n-11l]} (4.30)
Tri-Quadratic Interpolation : 27 points used
UI = (p2 + p)(q2 + q)@
11
,j,k- 4
2 2U~~~~~~~_ 1 2I
2 r 2 r U,'1-1, m-1,-n-1" Ui-1-1, j-m-1,k-n 2 r 2F r Ui'-1-1,j-m-1,k-n+1}
1
+ 2 +-(p +p)(1-q
1 1
(r2 + + U(1-1j-m,k-n-1 r2)U l -1,j-m,k-n + 2 r) U'--I,j-m,k-n+I
2 2




1
+( + p)(q q).
4
1 1
2) _I,j-m+I,,_n + 2( U ,jj ,n+,k-+I}
{2(2 r)Ui'--,j-m+,k-n- + (1- r2)Uil1m+kn r)U
1
+ -(1- p2)(q + q)
2
1 1
2 (r2 + r)U-11,k + ( 2)U, 1,k-,,_ + (r2 r)Ul,j-,,,,k-,+
2 2
+ (1- p2)(1- q2).
1 1
2 (r-2 + r)Uil,j-m,k-n-1 2)Uil,jm,kn + I(2 r)Uitl,j-m,k-n+1I
+ r)Ugl,,,,kfl + (1- 2U
2 i--n,k-n 2~ )~~lklI
1
+ -(1- p2)(q2 q).
2
1 1
{--(r2 + r)Uijm+knl_ + (1- r2 U_ + 2(r + r)Ui-,j-m+,k-n-1 2 U,'_1,j-m+1,k-n 2 2 r)Uil,j-m+I,k-n+1
1
+ -(p2 p)(q + q).
4
lr 1
2(r2 +rU-l+1,j-m-1,k-n-1+(I r 2)U-l+1,j-m-1,k-n + 2 r)U l+1,j-m-1,k-n+1
1
+ (p2
+ -(p2 p)(1- q2)
2
1 1
(2 2+ r)U-+,j-,k-- 2)Ui-l+1,j-m-I,k-n + I(2 r)U -l+,j-m,k-n+1
2+ 2(p2(2
2
2 rUl+ +,k1 + -2)Ul+, j-m+1,k-n 2- )U-l+,j-m+,k-n+
(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
Figure 4.4 Analytical solution of linear advection
eqato atT 25 e
2o7
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 .0 .6
0.2 :
0 7
0 003
x
Figure 4.6 Results of the 1-D linear advection equation
obtained with upwind method at T = 50 Sec
.6 0 i i a p i, 7I.
1:-4
0 .0 0 .2 0 I
00 000 0O
Figure 4.7 Results of the 1-D linear advection equation
obtained with upwind method at T 250 Sec
..... .m ...... .a u t 0 .... ....
0 9
0 a 0,7 0 .0
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.
00
0.2 0 7 0
Figure 4.9 Results of the 1-D linear advection equation
obtained with the Leith's method at T = 50 Sec
0 0
0 o0
00 0 .00
02 0
00 10 (0 00
Figure 4.10 Results of the I-D linear advection equation
obtained with the Leith's method at T = 250 Sec
0 .7 0 .6 220 5 .0
0 .5
.4
0
C2 250 500 760 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
A " 4-....- I 14 Io Io
o 5
Figure 4.13 Results of the 1-D linear advection equation
obtained with the TVD method at T =250 Sec
A.
0 .s 00
.. . '77'-, B . .
Figure 4.14 Results of the l-D linear advection equation
obtained with the TVD method at T = 250 Sec
A
7j
o .
o .5
. 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.
* LS .II II.
t,a
I:
10. .
0.s
2. 0 7 01 0
.tx

Figure 4.15 Results of the 1-D linear
obtained with the ELM with linear at T = 55 sec

advection equation interpolation

.. ......... .... ..
o .
3
Figure 4.16 Results of the I-D linear advection equation
obtained with the ELM with linear interpolation at T = 253 sec
E .. P d:9
0 0.0
0 .4
S070.00
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, ...
~o .s o
.0.1
250 10 750 000

Figure 4.18 Results of the I-D
obtained with the ELM with at T = 55

linear advection equation linear interpolation Sec

1 2LU I I I I d. 0
v o< .8% ,. : .
Io
0 4
o .5
o. 7 A8 1
Figure 4.19 Results of the I-D linear advection equation
obtained with the ELM with linear interpolation at T 253 Sec
0 .4
o I
42
o ,
Figure 4.20 Results of the I-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. o ..
o .7
1 0 .0
Figure 4.21 Results of the I-D linear advection equation
obtained with the Conservative QMSL at T = 55 Sec
0 ...... .n ll.I ...0
U0
0
0 0
0 4
o20 000 000 100
Figure 4.22 Results of the I-D linear advection equation
obtained with the Conservative QMSL at T = 253 Sec
U0o
0
0) 0.0
~~ .5 o'i
0 I 0 000I i 00 ? 0 1000
Figure 4.23 Results of the I-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.

o,+A- -x=L -+1 dxdt = o

(4.35)

Thus :

J x=L(+ c j) +][(uc)Rc(c .d
fx=O _C c-)dx+ (uc)Lldt = 0

(4.36)

where subscript R, 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

S0Lcq




60
If the convective velocity is constant, the convective equation is the same in either conservative or nonconservative form. Hence, to check the conservation of
numerical scheme, the convective velocity should be nonconstant.
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
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 1i-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
-+ uk- 0
whereas the conservative form of governing equation is
1
S u2
+ -0

(4.38)

(4.39)

o
5E5512
-1E-11. T
Elo 2 d0 300 400
E psd Tim.




The initial and boundary conditions are assume to be :
u(x,0) = (1- x) 0 x< 1
u(x,O) = 0 1< x 3 (4.40)
u(0,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 15 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 t_ 1 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 0:5x 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
Figure 4.28 Results of the non-conservative
upwind method at T 1.0 sec
'0A
o 7
0 6
0
a 3
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 3
0 2
or
Figure 4.30 Results of the conservative
upwind method at T 0.5 sec
0 .6 S 0.;
0 4
0 .2
Figure 4.31 Results of the conservative
upwind method at T 1.0 sec
0 3
0 .2 0 .1
0

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.8
.6
W 0
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 .2
i o ..
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 ....< -. .. .. .
D.3
0.2
.1
0o
x
Figure 4.36 Results of the Godunov method for non-linear advection equation at T = 0.5 sec
0 .4 04 3
.3
0.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
0 .2 0 .1 =On 2
oX

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 0 5.. ,. .
UO 4
Figure 4.39 Results of the non-linear advection equation
of ELM with linear interpolation at T =0.5 sec
.0 E
0 0
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.
o .7
Figure 4.43 Results of the non-linear advection equation
of ELM with quadratic interpolation at T 0 .0 sec
0 .2
0 .0
Figure 4.44 Results of the non-linear advection equation
of ELM with quadratic interpolation at T = 3.5 sec
o .7 v T ,c 1L s.
a 5 0.4 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.
07 0.0
o
o 4 1 I Ia
Figure 4.45 Results of the non-linear advection equation of conservative QMSL method at T : 0.5 sec
I 4 00
.o
.sx

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

.4
.o
o .2 o I
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.
I+AI 01l 1 X =L 011t
adt + -l= ---dx=0 (4.41)
Flux difference in time should be Zero.
( Icells+ I n+ ( Icells+1 )n+1 Icells+1 n Icells+1 n i=0 U i+1) i=ii= 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.
Conservation ofCons.rvaIv, ELM
1.0005
1 .0 0 3 i .00 02 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(OOt) = 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 10
Totol concentration

(4.45)




z
cF~nnevtv-101 ELM -Y
[ X

0.0
0.6 0.4

X(Krn)

Figure 4.49 Numerical solution of 2D non-conservative ELM

2 E 1 2 1 E -1 2
0
-1 E -1 2
-2 E -1 2

Figure 4.50 Conservation error of 2D non-conservative ELM

2(0 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
x100
75 50 25
0o

2D Non-Conservative
iF!
2 days
1 lday ]Inital Condition 1:11 I.II

5 days after

I I I I I

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.

SELM




I
0.8 o 0.6 0.4 0.2
0

z
200

X(Krn)

Figure 4.52 Numerical solution of 2D conservative ELM
2 E -12 S1E-12
o
0 0
0
-1 E -1 2
-2 E -12 ......

Figure 4.53 Conservation error of 2D conservative ELM




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.

)fiHlI

175 150
125
EIO
.Sd 100
75 50 25

U 50 100 150 200
X(km)

Figure 4.54 Contour map of numerical solution of conservative ELM

2D Conservative ELM
5 days after
4 days
S3 days
:0 2 days
0 1 day
Initial Condition I I 1 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., 100xlOOxlO 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 nonconservative 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
6yx

Figure 4.55 Numerical solution of 3D non-conservative ELM
0
*.0- .5
1i 5
.2
-2.-5 7
-3 -,. . . . . .

Figure 4.56 Conservation error of 3D non-conservative ELM




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

1
4 days a fte r
2 day after
In a Itia. C.ndltipn 7
IIi ii. ___l_____.....'JHH '....."....i~/ "... , ,
0 S0 100 150 200
X (km)

X
50 25
00
N

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
6y-x

Figure 4.59 Numerical solution of 3D conservative ELM

2 E .1 2 1 .8 E -1 2 1 .6 E -1 2 1 .4 E .1 20

1 2 3 4
T im e (D a ys )

Figure 4.60 Conservation error of 3D conservative ELM

- - - - - - - - - - - - - - -




Figure 4.61 and 4.62 show the numerical solution of the nonconservative ELM in X-Z and Y-Z plane, respectively. Figures show conservative ELM has less damping than non-conservative ELM.

3D Conservative ELM

150 125
100
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 N f

Y (kin)

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 quasimonotone 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.