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

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.

MILLISECOND	CLASS.METHOD	MESSAGE
0	sobekcm_page_globals.constructor
0	sobekcm_page_globals.constructor	Application State validated or built
0	sobekcm_database.verify_item_lookup_object
0	sobekcm_page_globals.constructor	Navigation Object created from URI query string
0	sobekcm_database.verify_item_lookup_object
0	sobekcm_page_globals.display_item	Retrieving item or group information
0	sobekcm_page_globals.get_entire_collection_hierarchy	Retrieving hierarchy information
0	sobekcm_assistant.get_entire_collection_hierarchy
0	cached_data_manager.retrieve_item_aggregation
0	cached_data_manager.retrieve_item_aggregation	Found item aggregation on local cache
0	item_aggregation_builder.get_item_aggregation	Found 'all' item aggregation in cache
0	system.web.ui.page.page_load (ufdc.page_load)
0	sobekcm_page_globals.constructor.on_page_load
0	html_echo_mainwriter.add_style_references	Adding style references to HTML
0	html_echo_mainwriter.add_text_to_page	Reading the text from the file and echoing back to the output stream
2	html_echo_mainwriter.add_text_to_page	Finished reading and writing the file