Front Cover
 Front Cover
 Table of Contents
 List of Figures
 Review of numerical schemes
 Numerical model
 Adaptive grid technique
 Numerical experiments

Group Title: UFL/COEL (University of Florida. Coastal and Oceanographic Engineering Laboratory) ; 90/016
Title: Two-dimensional finite-difference model for estuarine and coastal fronts
Full Citation
Permanent Link: http://ufdc.ufl.edu/UF00084101/00001
 Material Information
Title: Two-dimensional finite-difference model for estuarine and coastal fronts
Physical Description: xi, 127 p. : ill. ; 28 cm.
Language: English
Creator: Luo, Yu
University of Florida -- Coastal and Oceanographic Engineering Dept
Publisher: University of Florida
Place of Publication: Gainesville Fla
Publication Date: 1990
Subject: Navier-Stokes equations -- Numerical solutions   ( lcsh )
Gravity -- Measurement   ( lcsh )
Estuarine oceanography   ( lcsh )
Genre: bibliography   ( marcgt )
technical report   ( marcgt )
non-fiction   ( marcgt )
Statement of Responsibility: by Yu Luo.
General Note: "UFL/COEL-90/016."
Funding: This publication is being made available as part of the report series written by the faculty, staff, and students of the Coastal and Oceanographic Program of the Department of Civil and Coastal Engineering.
 Record Information
Bibliographic ID: UF00084101
Volume ID: VID00001
Source Institution: University of Florida
Rights Management: All rights reserved by the source institution and holding location.
Resource Identifier: oclc - 24229241


This item has the following downloads:

UF00084101 ( PDF )

Table of Contents
    Front Cover
        Front Cover
    Front Cover
        Page i
        Page ii
        Page iii
    Table of Contents
        Page iv
        Page v
    List of Figures
        Page vi
        Page vii
        Page viii
        Page ix
        Page x
        Page xi
        Page 1
        Page 2
        Page 3
        Page 4
        Page 5
        Page 6
        Page 7
        Page 8
        Page 9
        Page 10
        Page 11
        Page 12
    Review of numerical schemes
        Page 13
        Page 14
        Page 15
        Page 16
        Page 17
        Page 18
        Page 19
        Page 20
        Page 21
        Page 22
        Page 23
        Page 24
        Page 25
        Page 26
    Numerical model
        Page 27
        Page 28
        Page 29
        Page 30
        Page 31
        Page 32
        Page 33
        Page 34
        Page 35
        Page 36
        Page 37
        Page 38
        Page 39
        Page 40
        Page 41
        Page 42
        Page 43
    Adaptive grid technique
        Page 44
        Page 45
        Page 46
        Page 47
        Page 48
        Page 49
        Page 50
        Page 51
        Page 52
        Page 53
        Page 54
        Page 55
        Page 56
        Page 57
        Page 58
        Page 59
    Numerical experiments
        Page 60
        Page 61
        Page 62
        Page 63
        Page 64
        Page 65
        Page 66
        Page 67
        Page 68
        Page 69
        Page 70
        Page 71
        Page 72
        Page 73
        Page 74
        Page 75
        Page 76
        Page 77
        Page 78
        Page 79
        Page 80
        Page 81
        Page 82
        Page 83
        Page 84
        Page 85
        Page 86
        Page 87
        Page 88
        Page 89
        Page 90
        Page 91
        Page 92
        Page 93
        Page 94
        Page 95
        Page 96
        Page 97
        Page 98
        Page 99
        Page 100
        Page 101
        Page 102
        Page 103
        Page 104
        Page 105
        Page 106
        Page 107
        Page 108
        Page 109
        Page 110
        Page 111
        Page 112
        Page 113
        Page 114
        Page 115
        Page 116
        Page 117
        Page 118
        Page 119
        Page 120
        Page 121
        Page 122
        Page 123
        Page 124
        Page 125
        Page 126
        Page 127
Full Text




Yu Luo










I would like to express my gratitude to my advisor and

committee chairman Professor Y. Peter Sheng for his

initiating, motivating, directing and financially

supporting my work and for his criticism and patience

during the writing process. I would also like to extend my

thanks and appreciation to my thesis committee members,

Professor Donald M. Sheppard, Professor Wei Shyy and

Professor Robert J. Thieke, for their comments about this

work and patience in reviewing this thesis.

I must thank Mr. Hongjun Li, my husband, for his

numerous important and valuable suggestions about this work

and for his general help on computer programming and during

the writing process.

My gratitude must also be extended to Dr. James

O'Donnell for his useful suggestions on this work.

The support of all my fellow graduate assistants and

postdoctoral associate is also deeply appreciated. I would

like to thank Dr. Donald Eliason, Mr. Hyekeun Lee, Mr.

Steve Peene, Mr. Xingjian Chen and Mr. Jei Choi. My

appreciation must also be extended to Mr. Suburna Malakar


for his help in dealing with the VAX computers.

Finally, I would like to thank my family in China for

their encouragement throughout this study.




ACKNOWLEDGMENTS .......................................... ii

LIST OF FIGURES ........................................ vi

ABSTRACT .............................................. x


1 INTRODUCTION ..................... ................. 1

2 NUMERICAL SCHEME REVIEW ........................... 13
2.1 First-Order Upwind Scheme ................... 15
2.2 Flux-Corrected Transport (FCT) Scheme ...... 15
2.2.1 The First-Order Scheme ................ 18
2.2.2 The Second-Order Scheme ............... 18
2.2.3 The Combination-Hybrid Scheme ......... 19
2.3 Central-Difference Scheme .................... 20
2.4 Total Variation Diminishing (TVD) Scheme .... 21
2.5 Quadratic Upstream Interpolation for Convective
Kinematics with Estimated Streaming Terms
(QUICKEST) Method ........................... 22
2.6 Time Marching Method ...................... 24
2.7 Summary ........................................ 26

3 NUMERICAL MODEL ..................................... 27

3.1 Governing Equation .......................... 27
3.2 Numerical Model ................................ 30
3.2.1 Computation Domain ..................... 31
3.2.2 Descretization of Equations .......... 31
3.3 Finite-Difference Equation and Procedure ...... 34
3.3.1 Find i+lby Solving Continuity Equation.. 34
3.3.2 Solve u ?1 and v+1 ...................... 35
3.3.3 Solve S + ............................... 40
3.4 Stability ...................................... 42

4 ADAPTIVE GRID TECHNIQUE .......................... 44

4.1 Grid Adaptive Algorithm ..................... 45
4.2 Weight Function ................................ 47
4.3 Steady-State vs. Transient Problems ........ 50
4.4 Temporal Coupling ........................... 52

4.5 Verification of The Adaptive Grid System .....

5 NUMERICAL EXPERIMENTS ............................



5.1 "Breaking Dam" Problem ........................ 60
5.1.1 Physical Problems ....................... 60
5.1.2 Numerical Results ...................... 62
5.2 Flow with Linear Density Distribution ....... 64
5.2.1 Physical Problems ......................... 64
5.2.2 Numerical Results ....................... 69
5.3 Sudden Release Problem .......................... 74
5.3.1 Physical Problems ......................... 74
5.3.2 Numerical Results ...................... 77
5.4 Water Discharge Problem ......................... 79
5.4.1 Physical Problems ......................... 79
5.4.2 Numerical Results ........................ 91
5.5 The Application of Adaptive Grid Technique ... 100

6 CONCLUSION .......................................

APPEND IX ...............................................

REFERENCES ............................................

BIOGRAPHICAL SKETCH .....................................






Estuarine front with foam line observed just
north of the lightering area inside
Delaware Bay from Bowman and Iverson(1977) .....

Thermal imagery of Point Beach power plant
from Scarpace et al.(1975) .....................

South Pass effluent into ambient water masses

during flooding and ebbing tides
from Wright and Coleman(1974) .......

The comparison of numerical schemes
at CFL=0.8 .........................

The comparison of numerical schemes
at CFL=0.2 .........................

QUICKEST method ....................

Results of QUICKEST method
(a) Linear equation ................
(b) Nonlinear equation ...........

Computation domain and grid system .

Initial mass distribution .........

Momentum balance principle ........

Case (1) on adaptive grid points 20

Case (2) on adaptive grid points 20

Case (3) on adaptive grid points 20 ...........

The comparison of solutions on uniform grid and
adaptive grid
(a) Case (1) on uniform grid 20 and 80
(b) Case (3) on uniform grid 20 and 40 .........

Initial State of Dam ..........................







........... 12

........... 16

........... 17

........... 23













5.2 "Breaking dam" problem
(a) Analytic solutions at unit intervals of time
(b) Numerical solutions corresponding to (a)... 65

5.3 The numerical solutions for the layer depth with
rotation effects in (a) t=10 and (b) t=25 ..... 66

.5.4 Nonratating "breaking dam" problem
at t=l and t=6 ................................. 67

5.5 With rotating effects, the gravitive
adjustment at t=l and t=3 ........................ 68

5.6 Density Distribution ......................... 69

5.7 Numerical time-sequence showing the isopycnals r
with p= .058x10-3g/cm3. a=15xl0-cm-.
at t=8.Os, 12.3s. 16.7s. 20.9s and 24.9s ....... 71

5.8 Corresponding to Figure 5.2-1, a time-sequence
showing the motion of the isopycnals
in laboratory experiments from
Simpson and Linden(1988) ...................... 72

5.9 A log-log plot of the angle 9 of inclination
of the isopycnals to the vertical against the
dimensionless t'==(Iga) /'2t ..................... 73

5.10 Laboratory Experiment Case 1 ..................... 75

5.11 Laboratory Experiment Case 2 .... ................ 75

5.12 Large Scale Bottom Plume ........................ 76

5.13 Large Scale Surface Plume ...................... 76

5.14 Shadowgraphs of a volume of salt water
with ho=7cm, H=49cm, xo=50cm, g'=47cms-2
collapsing into fresh water at 2.2s, 3.8s
and 6.0s after release ...... .. ... .... ........ 80

5.15 Shadowgraphs of a volume of salt water with
ho/H=1, ho=7ctm, x0=50cm, g'=47cms-2
collapsing into fresh water at 4.7s,
7.7s and 10.7s after release ..... ............. 81

5.16 Numerical solutions corresponding
to Figure 5.14 ................................ 82

sta' En"eee-n, .-AIr^ li

3.17 Numerical solutions corresponding
to Figure 5.15 .................. ............. 83

5.18 Numerical solutions of Case 3 at 20hrs, 40hrs,
60hrs and 80hrs ................................. 84

5.19 Numerical solutions of Case 4 at 20hrs. 40hrs.
60hrs and 0Ohrs .................................... 85

5.20 Numerical solution for surface plume at 20hrs.
4Ohrs, 60hrs and SOhrs by Wang(1984) ............ 86

5.21 Numerical solution for bottom plume at 20hrs,
40hrs. 60hrs and 80hrs by Wang(1984) ............. 87

5.22 Numerical solution for surface plume at 20hrs.
40hrs, 60hrs and SO/rs by current model ......... 88

5.23 Numerical Solution for Bottom Plume at 20hrs,
40hrs, 60hrs and 80hrs Coresponding tp 5.3-12 ... 89

5.24 Water Discharge Problem ...................... 90

5.25 Water discharge problem with Fr=0.04
for overflow case at 2hrs, 4hrs. 6hrs and 8hrs
(a) Velocity vector ............................ 92
(b) Salinity contour ........................... 93

5.26 Water discharge problem with Fr=0.22 for
overflow case at 20hrs 40ihrs. 60hrs and 80hrs
(a) Velocity vector ........................... 94
(b) Salinity contour .............................. 95

5.27 Water discharge problem with Fr=0.057 for
underflow case at 4h/rs. 8hrs, 12hrs and 16hrs
(a) Velocity vector ............................. 96
(b) Salinity contour ........................... 97

5.28 Water discharge problem with Fr=0.64 for
underflow case at 4hrs, 8hrs, 12hrs and 16hrs
(a) Velocity vector ............................. 98
(b) Salinity contour ........................... 99

5.29 Two regular grids with grid size 80x10 and
40x10 for sudden release problem Case 4 ....... 102

5.30 Adaptive grid at 20hrs. 40hrs, 60hrs and 80hrs
with grid size 40x10 ......................... 103

5.31 The salinity contour obtained on uniform grid
80x10 at 20hrs. 40hrs. 60/irs and SOhrs ......... 104

5.32 The salinity contour obtained on uniform grid
40x10 at 20hrs. 40hrs. 60hrs and SOhrs ......... 105

5.33 The salinity contour obtained on adaptive grid
40x10 at 20hrs, 40/hs. 60hrs and 80hrs ......... 106

5.34 Two uniform grids with grid size 50x10 and
25x10 for overflow water discharge problem
at Fr=0.22 ...... .... ... ......................... 107

5.35 Adaptive grid at 2hrs. 4-hrs. 6hrs and Shrs
with grid size 25x10 ............................. 108

5.36 The velocity vector plots at 2hrs. 4hrs, 6hrs and
Shrs on uniform grid 25x10 ..................... 109

5.37 The velocity vector plots at 2hrs, 4hrs, 6hrs and
8hrs on adaptive grid 25x10 .................... 110

5.38 The salinity contour obtained on uniform grid
50x10 at 2hrs, 4hrs, 6hrs and Shrs ............. 111

5.39 The salinity contour obtained on uniform grid
25x10 at 2hrs, 4hrs, 6hrs and Shrs ............. 112

5.40 The salinity contour obtained on adaptive grid
25x10 at 2hrs, 4hrs, 6hrs and Shrs .............. 113

5.41 The salinity contour obtained on uniform grid
25x10 at 2hrs. 4hrs. Chrs and Shrs by Upwind scheme. 114

5.42 The velocity vector plots at 2hrs, 4hrs, 6hrs and
8hrs on adaptive grid 25x10 by Upwind scheme ... 115

A.1 Sketch of characteristic method ............... 119

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



Yu Luo

December 1990

Chairman: Professor Y. Peter Sheng
Major Department: Coastal and Oceanographic Engineering

A two-dimensional numerical model is developed to

predict the properties of gravity currents and resolve the

sharp density gradient at the frontal zone in estuarine and

coastal waters. The governing equations, which are

vertically integrated Navier-Stokes equations of fluid

motion are solved by using a finite-difference scheme with

a staggered grid in space and the Euler forward-time scheme

for time integration. In order to capture the sharp

variation density gradient at the front zone accurately and

efficiently, this study employs a high resolution TVD

(Total Variation Diminishing) scheme, which is simplified

by "freezing" the local constants for complicated system

equations, and an adaptive grid method with weight function

correction. A series of analytical solutions and laboratory

experiment results is used to verify the accuracy of the

model. The model is then employed to investigate the

dynamics of fresh water or cool water discharge into a

large water body.


The dynamics of coastal and estuarine are often

locally dominated by shallow layers of buoyant water or

plumes. These layers are often established by the discharge

of fresh water from a river or bay into estuarine or

coastal water, and as another example, the cooling water or

wastewater from a power plant into a river or estuary. As

the discharging and receiving waters with different

densities (due to different salinities or different

temperatures) meet, they are often separated by a sharp

front, i.e. a narrow zone of sharp density gradient.

Fronts occur over a wide range of spatial scales. In

an estuary, a bottom salinity front (salt wedge) separates

the incoming ocean water from the estuarine water. Near the

river mouth, a surface salinity front (river plume)

separates the spreading, brackish water from the ocean

water. On shallow continental shelves, fronts are formed at

the transition between well-mixed and stratified waters.

And along the shelf break, a semi-permanent shelf-slope

front is formed during the winter months, separating cold,

fresh shelf water from warm, saline open-ocean water. Other

fronts are found on the continental shelf, separating a

zone of coastal water from oceanic water. Fronts also occur

on a larger scale in the deep ocean, between water masses

of different properties (Bowden, 1983).

There have been many observations of plumes and

fronts. The Mississippi River's discharge at South Pass

produces a plume into the Gulf of Mexico. The behaviour of

the plume. with its load of suspended sediment, has been

the subject of a number of investigations, including that

of Wright and Coleman(1974). The spreading river water has

a considerable influence on the coastal band of water out

to the edge of the continental shelf, as marked by the 200m

contour which, in the region of the delta, is about 50 km

from the coast.

The Columbia River flows into the Pacific Ocean as the

boundary between the states of Washington and Oregon, USA.

Duxbury(1965) found that its estuary is partially mixed as

a result of a high river discharge averaging 7,300 m3s-,

combined with strong tidal currents. The outflowing upper

layer forms a plume of relatively low salinity water which

can be distinguished over a large area of the adjacent

ocean. The behaviour of the plume depends largely on the

wind and wind-driven currents and shows a seasonal


Another example of river plume is from Bowman(1978).

The Hudson River plume, which is largely determined by

wind conditions, flows into the New York Bight.

Szekielda and Mitchell (1972) found heavy metal

concentrations three orders of magnitude higher than

background values in frontal zones along Delaware Bay for

chromium, copper, lead, mercury, silver, and zinc. Klemas

and Polis(1977) graphically illustrated how convergent

tidal currents in Delaware Bay captured and advected oil

slicks over period as short as one hour.

Scarpace et al.(1975) presented a study of the

structure of the thermal plume of the Point Beach power

plant in Lake Michigan. Again a major feature of the plume

was the sharp frontal boundaries in the near field.

Secondary fronts were also observed in the far field. These

fronts appeared as boundaries between areas of slightly

different temperatures.

Garvine and Monk(1974) reported the results of

observation of the hydrographic fields and the circulation

in the neighborhood of the Connecticut River front. They

found that the front essentially separates the brackish

plume water from the ambient Long Island Sound water and

that mixing between them takes place within 50m of the


Since the gravity current is an interesting and

important problem of concern to scientists,


environmentalists and geographiers. there have been- many

attempts to model plumes and fronts. For example, Kao et

al.(1978) and Valentine and Kao(1983) examined the detailed

structure of the flowfield produced by a buoyant inflow

into a quiescent receiving basin. They reduced the

primitive equations into stream function and vorticity

equations. Which are then solved by using a conventional

over-relaxation method with the stream function for the

point source potential flow solution as an initial guess.

In their discussion of accuracy related to the advection

term, they defined a false diffusion term which resulted

in truncation error. They analyzed the effects of the false

diffusion by comparing the advective term, the physical

diffusion term and the false diffusion term at various

areas in plume region including frontal zone. They

concluded that the advective term is dominant over the

false and physical diffusion in the headwave region.

However, in the main established region behind the

headwave, the false diffusion becomes important. This

computational diffusion is probably large enough to -mask

the physical diffusion and viscous effect. Their results

showed that the currents reached a steady speed of

propagation and the layer behind the front appeared to

reach a steady finite depth. They discussed the influence

of the densimetric Froude number for the stratified gravity

currents in near field and reported that their results are

consistent with the observation that the frontal zone forms

a strong surface convergence. Furthermore, their model

prediction agreed quantitatively with data from the

Connecticut River Plume provided by Garvine and Monk(1974).

Fischer(1976) developed a numerical model for density

currents in estuaries. He pointed out that the main

difficulty for the numerical model lies in the numerical

treatment of the convection terms. The usual methods for

treating the convection terms (e.g., Forward Time/Centered

Space, Upstream Difference, Lax-Wendroff, ADI, Leapfrog or

Galerkin Finite-Element Method) break down when the

convection of fairly sharp 'density jumps is to be

calculated, as is the case for salt wedges. He described a

new method "Hermite-interpolation method," in which the

spatial density distribution is represented by IIermite

interpolation functions and the density value p and its
local gradient are defined as time dependent variables.
In his model, the advection term is neglected in the

momentum equation which is solved by an explicit leap-frog

method with an implicit treatment of the viscosity term. He

concluded that Hermite-interpolation method gave

satisfactory results in every case.

O'Donnell and Garvine(1983) and O'Donnell(1986)

developed an approach based on the assumption that the

plume may be divided into two dynamically distinct domains.

the inviscid region or plume body, and the frontal zone.

The plume body is considered to be the upper layer of a two

layered flow between which turbulent exchange with the

ambient fluid is negligible. The dynamics of the flow in

this region is dominated by gravitational spreading across

which flow properties are related by algebraic jump

conditions. In their two-dimensional model, they used the

"front patch technique." which employed the

bicharacteristic relationship and the characteristic cones

of the systems in the frontal zone. In the interior of the

plume, the two-step Lax-Wendroff scheme was employed. They

discussed and calculated the 'growth of a plume from a

radial source in a steady crossflow. It was found that the

nearshore areas of the plume are most sensitive to the

effects of rotation. The effect of the oscillation tidal

crossflow on the Connecticut River plume indicated that the

plume is diluted by vertical mixing soon after a change in

the direction of the crossflow.

Wang(1984) developed a primitive equation model to

examine the two-dimensional mutual intrusion of a gravity

current in vertical plain and the formation of a density

front. He used a staggered-grid, central-difference scheme

in space and a leapfrog scheme in time, except for the

advection terms in the salinity equation which were treated

with an upwind scheme. A mode-spliting technique (Sheng et

al.. 1978) was used to achieve computational efficiency.

Later Wang(1984) changed the upwind scheme in the salinity

equation to a diffusion-corrected advection scheme

(Smolarkiewicz, 1983. 1984) for solving the advection term.

He pointed out that the upwind scheme can generate large

numerical diffusion which can be significantly reduced by

using a dffusion-corrected scheme. His solution gave a

realistic description of the frontal structure and the mean

southward current by comparing with Garvin and Monk's(1974)

observations from the New England shelf-slope front.

Wang(1984) extended his two-dimensional model to a three-

dimensional one to study gravity currents produced by

instantaneous release of a buoyant fluid in a rectangular

channel with and without the Coriolis effect. His results

were compared with laboratory experiments provided by

Rottman and Simpson(1983).

Chao and Boicourt(1986) developed a three-dimensional,

primitive-equation model to study the onset of estuarine

plumes. Their numerical solution essentially followed the

scheme of Semtner(1974). Equations and boundary conditions

are finite-differenced in space to second-order accuracies

and leap-frogged in time with lagged friction and

diffusion. In case of outflow boundary condition, upwind

finite differencing is employed for the advection term. The


finite-difference scheme for the advection terms conserves

mass. energy, salt and salt variance. At each time step,

the vertically-averaged current and the deviation from the

mean are found separately. The former was obtained by a

five-point, sucessive over-relaxation of the vorticity of

the depth-averaged flow predicted via a vorticity equation.

The latter was found by integrating the momentum equations

in time while removing the barotropic components. Regular

averaging of two consecutive time steps was carried out to [

avoid the splitting of the developing solutions at even and [

odd time steps. If the ocean becomes gravitationally

unstable, instantaneous vertical mixing was carried out to

achieve a stable stratification. Their results showed that

the plume expands in the direction of propagation of the

coastally trapped waves after the freshwater release. The

seaward expansion of the bulge decreases, and the

undercurrent leaks out of the bulge and propagates behind

the nose of the bore. The three-dimensional structures of

the density-driven estuarine circulation and also of the

bore intrusion along the shelf have been identified.

All these models tried to obtain good results by using

a variety of techniques aimed to improve the accuracy and

efficiency of the numerical solution in vicinity of a

front. This thesis will present a study to further improve

the accuracy and efficiency of the numerical solution of

frontal problems. A two-dimensional numerical model is

developed by using a finite-difference method in which a

high resolution TVD (Total Variation Dissipation) scheme

for simulating advection terms is employed in both the

momentum and the salinity equations and an adaptive grid

technique is applied. To verify the accuracy of this model,

numerical solutions are compared with analytical solutions

and some laboratory experiment results. The capabilities of

this model in solving the realistic problems are

demonstrated with further numerical experiments.

This thesis consists of 6 chapters. An introduction is

given in Chapter 1. In Chapter 2. a brief review of several

numerical schemes for solving advection terms are

presented. A model equation (Bergur's Equation) is used to

verify those schemes separately and the comparisons of

solutions are presented. In Chapter 3, the derivation of

the model from the shallow water equations is presented.

The computational domain and the computational procedure

are explained. Chapter 4 describes an adaptive grid

technique. The effect of weight functions is discussed and

different test cases are evaluated to assess the

effectiveness and efficiency of the numerical grid

adaption. In Chapter 5. results of a series of numerical

experiments are presented. Comparisons with analytical

solutions and laboratory experiment results are presented.


Some more realistic problems are also solved by the model.

Chapter 6 gives the conclusion and suggestions.

- 5.* .
I _* 0

Figure 1.1 Estuarine front with foam line observed just north
of the lightering area inside Delaware Bay from
Bowman and Iverson(1977)

Figure 1.2 Thermal imagery of Point Beach power plant
from Scarpace et al(1975)


Gulf water -
S 00 1000 0 00 1000
met*'$ Momrs

Figure 1.3 South Pass effluent into ambient water masses
during flooding and ebbing tides from Wright and
oor Ae r

FJ~ure1.3Souh Pss fflen int Ben w~e mtse
duringI odn n bbn =ds rm rgh n
Co-len (174


The mathematical model of gravity currents is governed

by Navier-Stokes equations. However, the numerical

solutions of Navier-Stokes equations are plagued with

difficulties caused by the numerical treatment of

advection terms. An ideal numerical scheme should have the

following properties:

1. Accuracy: Second order accuracy is usually


2. Capibility to capture the discontinuities.

3. Efficiency.

4. Robustness: Able to handle different problems.

5. Convenience: Easy to program.

Unfortunately, no scheme can meet all of the above

requirements presently. Many of the schemes work well only

for the special application for which they have been

designed. Wang(1984) used an upwind scheme and a diffusion-

corrected scheme for the advection terms. O'Donnell(1986)

used the Lax-Wendroff scheme for the advections terms to

solve the two layer front problem.

In this chapter, a few schemes which are usually

employed to solve hyperbolic equations are disscused based

on a model equation, i.e., a one-dimensional advection-

diffusion equation:

atO+=F x (2.1)


u Dependent variables

F u

y Viscosity coefficient

t Time

x Spacial coordinate

For inviscid flow, p=0. The equation becomes

Sat> (2.2)

In this equation, the characteristic speed is aF=u.

The initial conditions are

1 x=O, t=O

u= (2.3)

0 x>O. t=O

and the boundary conditions are

1 x=O

u=f (2.4)

0 x=l

All numerical computations are completed on a grid of

20 points. The Courant number, CFL=uy, is chosen 0.8 and


2.1 First-Order Upwind Scheme

The first-order upwind scheme is written as

u('+=u,-A(OF',++( -2o0)F'-( 1 -)F'1!) (2.5)

0 ui>O

Aj A= x (2.6)
1 ui

The results in Figure 2.1(a) and 2.2(a) with different

CFL numbers demonstrates that a non-oscillating shock is

captured, but the shock is smeared due to the high

dissipativity of this scheme.

2.2 Flux-Corrected Transport (FCT) Scheme

McDonald(1984) made an improvement to the first-order

upwind scheme by introducing a flux-correction term which

results in a second order FCT scheme. His method is

separated into two parts, the first-order scheme and the

second-order scheme. The combination of these two parts


1.0 O O O Oe9


0.0 --~-e-o0-9-O-e-e-

-0.50.0 0.1 o.2 0.3 0o. 0.5 0.6 0.7-O.1 0.9 1.0


1.0 i ----,--,--C-,.----- -,-----------


0.0 C -e- --e-

.0 .-. -.2 -.3 .4 o.- 0.6 0. 0.8 0.9 1.0


1.5 1.5

1.0 --e---e-e--e-e -e-o-- -e- -e-e-e-eO-e--- e

.S 0.5 -

0.0 -e-e -e- -e- e- o.o -e-o-e-o- oe--e-

-0.5o.0 0. o2 0.3 O. 0.5 0.16 0.7 0.8 0.9 1.0 -0.50.0 0.1 0.2 0.3 0.1 o.S 06 0.7 0.8 0.9 1.0
(c) (cl)

Figiu re 2.1 IThe complariison of numerical schemes at CIF'l=0.8





a, a


S oP

Uj 3


C )





0>x .



In 0 1 0 41
0, 0 0, a


0 i

411 0 In 0 In
- ~ 0





x c
u, Z






makes a hybrid form by the method of flux-correction

(Zalesak, 1979).

2.2.1 The First-Order Scheme

The numerical difference equation is written as


0 is the numerical flux given by

At F
Ax -



LFX i+l

where w is a variable which is used to determine the upwind

direction. It can be defined as

wi+l/2=-(Fi+l-Fi)(ui+i-ui) (2.9)

2.2.2 The Second-Order Scheme

From the relation

F ,,au
t -^t

a finite difference equation can be written as

U =i -i .+l/2+? i--1/2



Fn+ l/=F"- i (2.11)

The one-sided second-order differences in flux conservative

form is defined as

31"n+l/2 1n+1/2
si "2 i-1 W+1/2 0

Ai+t/= (2.12)
3 n+1/2 1- n+1/2
Fi+l i+2 Wi+1/2<0

2.2.3 The Combination-Hybrid Scheme

u,=ut--i+1/2+0i-1/2 (2.13)

Where, uj represents a first-order update. The second-order

flux-correction is difined as

64i +l/2==
Because of non-monotonicity of this flux-correction, a

filter, i.e., a flux-limiter must be introduced (Boris and

Book, 1973).

645i+1/2=-sgn(6i+1/2)max O, min(16i+1/21 sgn(6 ,i+i/2)(ui+2-u +I),

sgn(6i+,/2)(u$--u;_I)] (2.15)

Then the hybrid scheme is completed as follows

u2+'1= ui-6 4i+12+6i-1/2' (2.16)

The numerical results are shown in Figure 2.1(b) and

Figure 2.2(b). As can be seen from the figures, the shock

is sharper than the first order results and it is also free

of oscillation. However, this scheme requires much more CPU

time in each time step. The extension of this scheme to a

system of hyperbolic equations is quite complicated because

the direction of the characteristic velocity of the system

has to be determined at each grid point.

2.3 Central-Difference Scheme

In comparison with the upwind schemes, the central-

difference schemes require less CPU time and are easier to

program. However, classical second order central difference

schemes are not able to give smooth solutions in the

neighbhood of discontinuities. Generally all central

difference schemes can be expressed as the sum of a flux

term and a high order diffusion term.

uj ,=u-1 ( F --1 (2.17)

where F,4i/3=1(F,+ +F,), F,_1/2=.(Fi,_+Fi). Fi+/2-Fi-1/2 is the flux

term. Different ways to construct the diffusion term D will

result in different schemes.

If Di is expressed as

=i=( A-)2[(Ai+/2(F -F1)-Ai-_/2(F -Fit)]
AX 1+ 9 r I-


Ai+l/2=-- 1(ui+ui+1)
Ai -/2-=-(Ui+ui-t)

the well-known Lax-Wendroff scheme (Lax and Wendroff, 1960)

is obtained. It can be seen from Figure 2.1(c) and 2.2(c)

that oscillations in the solutions are presented obviously.

2.4 Total Variation Diminishing (TVD) Scheme

The high resolution oscillation-free TVD scheme (Yee

and Harten, 1985) is constructed. If diffusion term I( works

as a flux-correction term. As treated in FCT scheme, all

oscillations are "limited" by a function called "limiter".

The TVD diffusion term is given by

i -= i+1/2-- i-1/2 (2.19)

The diffusion function 0 can be written as

1i+1/2=---b(ai+l/2)[Aui+/2 -Qi-1/2] (2.20)

where 0 is a correction term to prevent zero value of a.

( a)= 6=10-4 (2.21)
a'+62 <6

a=a=u (2.22)

Q is the limiter function

Q+1/2=minmod(Aui_l/2 ,I i+1/2 Ui+3/2) (2.23)

ui-1/2=ui--ui-1 (2.24)

The minmod function is defined as

minmod(x, y)=sgn(x)max{ O, min[ljx, ysgn(x)]} (2.25)

As shown in Figure 2.1(d) and 2.2(d), TVD scheme gives

a high resolution discontinuity. There is no oscillation in

the solution. TVD scheme also needs more CPU time than

upwind and Lax-Wendroff schemes do, but it is easier to

program than FCT scheme.

2.5 Quadratic Upstream Interpolation for Convective Kinematics
with Estimated Streaming Terms (QUICKEST) Method

The QUICKEST method is a third order scheme based on

the quadratic upstream interpolation (Leonard, 1979). This

method involves three-point upstream-weighted quadratic

polynomial interpolation for estimating control-volume face

values and gradients.

The finite difference equation can be written as

un+ n=+C,1U-Cru
U i = U7+! ut-- rtar


where C=uAt is the local
and r denote i-1 and i+1
2) 2'

Courant number. The subscripts I

i-1 i i+1

Figure 2.3 QUICKEST method

u;= (uC+u)- (1-C)CURV- ~Cr1GRAD
u;2CC 1 2


The curvature term, .which must be upstream-weighted

for stability, is defined as





uRR 2uR+UC Ur

where subscripts RR, R, C, L and LL denote i+2, i+1, i, i-1 and


The gradient term is defined as

GRAD=ui+-u, (2.29)

u" can be found with corresponding formulas (shifted

one index to the left).

Leonard(1979) gave a result for linear model equation

as shown in Figure 2.4(a). For nonlinear equation, when CFL

is larger than 0.4, 'the numerical results present serious

oscillation. Even when the CFL reduces to 0.2, the result

shown in Figure 2.4(b) is still not good as other schemes

discussed above.

2.6 Time Marching Method

An Euler forward explicit scheme is used for time

marching. Even though implicit scheme can guarantee large

time steps due to unconditional stability, it needs much

more CPU time in each time step. .On the other hand, the

explicit method is more accurate because of the smaller

truncation error O(At). Another reason to chose explicit

method is that the explicit method is much easier to be

implemented than the implicit method for any kind of

computational domain.





'0. 0o 0 0.1 0.2 0.3 0.11 0.s 0.16 0.7 0.8 0.9 1.0

(a) Linear equation



-0.So 0.1 0.2 0.3 0. 9 0 .6 o.l? 0.L--. 1.0-

(b) Non- iin;r t (eq(iiLt. i ill

Figure 2.4 Ilesul-ts of QUICKEST meitl.hod



Eo E 6 9 -

_ iC~___;i_; ~ ___


2.7 Summary

From the comparison of these five schemes shown in

Figure 2.1, Figure 2.2 and Figure 2.4, only the solutions

solved by FCT and TVD have the sharper gradient variations.

However, the FCT method is more complicated than TVD

method. The programing of FCT method is also more difficult

than that of TVD method. Therefore, TVD method is chosen as

the numerical scheme for advection terms in the current

model .


3.1 Governing Equations

The two-dimensional dynamics in the vertical plane

across the density front may be mathematically described in

space and time by means of a set of five equations together

with appropriate boundary and initial conditions. Four of

the equations represent the conservation of volume,

momentum and salt. The fifth relates the density of water

to salinity.

In Cartesian coordinates,

Conservation of mass:

auaw( )
=0 (3.1.1)
dx Oz

Conservation of momentum:

au a(uu) 8(uw)
9t Ox Oz

P= P (A u,+ (A, ) (3.1.2)
Poax 57 z aZtz OX OX

Suva vw (AA- )+ (A v) (3.1.3)
at _x Oz Oz Oz x dx


Conservation of salt:

as O(us) +(ws)
-t x Oz

a as 8 O Os

Equation of state





x,z Subscripts denoting partial derivatives in

the longitudinal and vertical directions

u, w Longitudinal and vertical velocity components

v Latitude velocity component. It exists only when

the Corioli effect is considered

p Pressure

p Density

s Salinity

f Coriolis parameter

Ah,AV Horizontal and vertical eddy viscosity

Kh ,KV Horizontal and vertical eddy diffusivity

0 An empirical constant : 7.8 x 10-4

In equations (3.1.1) (3.1.6), the pressure is

7- P g

assumed to be hydrostatic and the flow is assumed to be

incompressible and isothermal. Also, the Boussinesq

approximation is used.

The boundary conditions are

(1) At the sea surface there is no surface stress.



(2) At the ocean bottom there is bottom friction, and no

flow coming through the bottom



where A is a linear drag coefficient.

(3) At the coast, the normal flow and flux vanishes



all av)=

all az

L9 Z

AU,_x V=O

K.J =0

(4) At the open ocean end, the ambient condition are

specified as


A,(u, v, )=O (3.1.10)
(''x ax Dx


where 7 = Water surface level

=0 mean water level

>0 above mean water level

3.2 Numerical Model

A longitudinal and vertical two-dimensional model for

nonrotating water model is developed in this section.

Equations (3.1.1) through (3.1.5) are solved with a

two-time-level and central-difference scheme in a staggered

grid. In order to reduce the numerical oscillation and to

capture the sharp density variation in front zone, the high

resolution TVD scheme and adaptive grid technique are

employed in the horizontal direction.

3.2.1 Computational Domain

Figure 3.1 shows the computational domain and the

staggered grid. Computation begins from the bottom of the

water body to the surface for the continuity equation. The

momentum equations are then computed from the surface to

the bottom. The i-index is counted from the closed

boundary(left) to the open boundary(right). The k-index is

counted from the bottom to the surface.

3.2.2 Discretization of Equations

The first step in this numerical model is to divide

the water body into horizontal layers of thickness Az.

Equation (3.1.1) (3.1.5) are then applied to each layer

and integrated vertically across the layer thickness

resulting in a set of hybrid differential-finite-difference

equations (Cerco, 1982).

After integration, the continuity equation (3.1.1) is

rewriten on surface layer as

SW 8(uh) (3.2.1)
t = x .

Below the free surface, the equation (3.1.1) becomes

= (uh) (3.2.2)
wt = W x

2 3 4 5 '6

Figure 3.1 Computation domain and grid system


u = Vertically averaged value of velocity in each layer

h = Layer thickness

Az+1 surface layer


Az otherwise

t,b = Subscripts denoting parameter evaluated at the

top or bottom of a layer

Rewrite Eq. (3.1.2) and Eq. (3.1.3) in the conservation



at +Fax Oz
Q^t 9x az


- uuh



vh -u,w, +ubwb -P A +( h)+rT,-r

-fuh-vw, +ubWb+ 0(Ah)+ry,-ry,
Ox Ox

v = Vertically averaged value of velocity

in each layer

V= Iuhl

From Eq. (3.1.5)

8(sh) O(u h)
-+ -+SteW-Sbbw
Ot Ox

=-(K(h )+(K) )-(K0 )b


s = Vertical averaged value of salinity in each layer

P = Integrated pressure over the layer thickness

r7,ry= Shear stresses

r,,,r,-= (Avu A ),

(A"z, A', ),' otherwise

T7,b 'yb -I

(-Au, -Av)b

estuary bottom

3.3 Finite-Difference Equation and Procedure

3.3.1 Find 1?+' by Solving Continuity Equation

Computation begins from the bottom and ends at the


From equation (3.2.2) and equation (3.2.1)

z ( 3 3 1)Z( U
,j Ax Zj-l il- '.f



At the bottom. w'=0. i.e. w'1=0


Sn _(u ~2, .n u -1
i.2-- =_.(U i u -" i-1.2)

At the surface, wt= Thus,

h Tl +Atwu


3.3.2 Solve u1' and v+1

The finite-difference equation for equation (3.2.3) is

t' =?+t(FLUX(W +R) (T?) (3.3.4)

where FLUX term can be separated into two terms



where FLUXC(W ) are the net fluxes solved by central

difference schemes and FLUXT(W ) are the TVD type diffusion


QF.-- F i+i,-F i-i/2


where F i+1/2--

F i_- +F
F i-l/2

In the current development, the concept of TVD scheme

is only defined for scalar conservation laws or constant

coefficient hyperbolic systems. It is impossible to fully

extend the theory of nolinear scalar TVD scheme to the

nonlinear systems. What we can do is to extend the scalar

TVD scheme to the system so that the resulting scheme is

TVD for the "locally frozen" constant coefficient system.

(Yee, 1987; Li and Kroll, 1989) For our two-dimensional

problems, this procedure can be greatly simplified.

The TVD type diffusion term FLUXT(Wi) is expressed as

FLUXT(W,)= (Ti;1/2Ji+l/*>-Ti-,/.9-1/2 (3.3.7)

where T is the right-eigenvector matrix of the Jacobian.

From equation (3.2.4)

Suu u\. (3.3.8)
F LuvhJ =iuV14

For a the "frozen" constant coefficient system, the

Jacobian matrix can be written as

A u (3.3.9)

Obviously, the two eigenvalues of A are

At=u A,=u

and the corresponding right-eigenvector matrix and its

invers are

T = 1 T- 1 0
T 0 1 0 1


In this case, the TVD FLUXT can be reduced to


The diffusion function 0 can be expressed as


Si+1/2=- i+1/(Ali/2 i+/2-0 i+1/2)

AW j+i/2 =T/12( i+i -WW)=' i+-W i






lI.+l/21 < 10-4



FLUXT(W i)=. .Zi+/2-0 '-1/2-

6 is a small positive value.

ai+1/2= 2 -" U+1/2

The limiter function Q can be written as follows

Q i+,/2=minmnod(AW l;_1/, AW i+1/2' i+3/2)


In these calculations. no matrix is involved, so the

computation of FLUXT is much faster.

R(W1) consists of many terms. Writing equation (3.2.5)


1 0P 8 8ub
fvh-u,w,+ubwb Po- x + (A+ -')+r, ,- rb

-fuh-vw,+uIbwb+') (Ah Vh)+ ry-Tyb
O N 1N


t= .(Ui,j +ui,j+1)

vt= ~(ij+i,j +l)

Wt= (W i,j+wi+l,j)

U =. (uij-. +ui, )
2+ U

Vb= 1("ij- Vij)

b = (WJi,j-l+ iW+tj-t)

constants, the stress terms can be



Assuming Ah and Av are


written as

r = I U ij+= Ui.j

9v A =v vi.j+ --vi j
rt --A = AZ
r2 = '^ z~' A A' Az-


rb = Av- lb Avuj zj-
baz Az

Vi -V i.j-1
T yb -A b -= Av .Z- 1
j; b \z


S_ 9uh A ui+ i+/2 -2ujhi-/2.j +ui-l,jhi-3/2,
'kh-'-A )=Ah A2
xA ax Ax2


The pressure grident -lP is calculated from the

surface to the bottom. At surface layer, from hydrostatic

equation (3.1.4)


p(z)=pg(7- z)

ap= 2- g- z)
ax 8r9x

Ox jax JI-ax



By Liebnitz's rule,

api = ,0 20P;
=P9g(i+Az) +g(Azd + xx


=pg( + z) I+ i- g(i+ z)2Pi+ j--Pi-l.J
=pg(?+Az ) "2Ax 4 +' Ax

At the subsurfaces,

S P IP .P Aza
ij.i--- ij4+1 gAz i,+t g (Pi+,j-Pi-1,i)

j=JM-1, ..... 1


3.3.3 Solve S'+1

Once the velocity components u, v are obtained, the

salinity equations can be solved.

From equation (3.2.6). the finite-difference equation

for salinity equation with two flux terms can be written as

n+ =( hi FLUXC(sh)+FLUXT(sh)+RHS(sh)3.3.4)
t,j h,,,

where h=Az+q for surface, and h=Az for subsurface.

Again, FLUXC(sh) is the fluxe solved by central

difference schemes. FLUXT(sh) is the TVD type diffusion


FLUXC(sh)=-_ x- sWt-SBwb)

1 (sh)i+l.,+(sh)i.j (sh)i.j+(sh)i-, i ,
= -A(x -2-u -2 u2-1'

(sw,-sbwb) (3.3.25)


st = (s ij l+s;
wl =wi,j

s=.-(si,j +s

b i j-i j
W"b= ;i~j--1

Since the salinity equation is solved independently,

the FLUXT(sh) can be calculate in the same way as solved in

scale equation discussed in Chapter 2.

FLUXT(sh)= (.i+1/2 _i-/2)


where 0 can be solved by substituting A(sh) into Au in

equation (2.22).


and the limiter function can be given as


The diffusion term of salinity equation is solved with

the assumptions of constant K,, and K1,.

Qi+,/2=minmod(A(sh)i-1/2., A(sh)i+1/2, A(sh) i+3/2)

A(sh)i+/2=. ((sh)i+(sh)i+,)

0H 1h sh= sh
RHS(sh) =Kh -+v
ox" 0z2

S2 =K 1 ((sh)i+ -2(sh)i,+(sh)i-1j)

(h -=K ((sh)i.+-,-2(sh)i,j+(sh);ij-

3.4 Stability

Equations (3.2.3) and (3.2.4) are solved explicitly in

the longitudinal direction. The computations have to be

subjected to several stability and convergence constraints.

Cerco(1982) discussed the following inequalities.


(Ah, ,)h

(A uK2a)
(AAKh) 2




Inequality (3.4.1) is restrictive to the selection of

Ah and Kh. uAx/kh is well known as "cell Reynolds number."

If the "cell Reynolds number" is equal to 2, the equations

(3.4.2) and (3.4.3) satisfy the stability criterion CFL





If this number is greater than 2. equation (3.4.2) will not

satisfy CFL condition. In this case, the computation

remains stable but oscillations develop in the solutions.-'

Inequality (3.4.2) arises from the explicit treatment of

the longitudinal dispersion term and is more restrictive

than inequality (3.4.3) which is due to the explicit

treatment of the advection term. If At satisfies the

inequality (3.4.2), the numerical calculations are stable

and convergent. Inequality (3.4.3) is the dead limitation.

The violation of this constraint results in ultimate

instabilities in numerical computation.

Finally, the most restrictive stability criterion

associated with the wave propagation on the free surface

must be satisfied.

At < Ax (3.4.4)

Actually, this is the dominant stability criterion to

limit the time step At.


Adaptive grid technique is a method by which the mesh

points can be self-adjusted as the numerical solution

evolves in time such that grid resolution increases in

regions with large gradients of solutions but decraeases in

regions with little or no gradients. This adjustment can

reduce the truncation error and improve the numerical

accuracy in regions with sharp gradients such as an

estuarine or coastal front. In addition, grid points will

be reduced in regions with little or no gradients, thus

increasing the numerical efficiency.

Due to its potential for improving the accuracy and

efficiency of the numerical solutions, adaptive grid method

has have rapidly emerged as an important tool in numerical

simulations of flow problem. With adaptive grid, each grid

point acts like an internal observer whose purpose is to

monitor the solution. The grid points distribution over the

computational domain is correspondingly adjusted

dynamically to concentrate grid points in regions of high

solution variability. The improved accuracy is achieved by

the redistribution of the grid points to increase the


number of grid points in the region where the solution has

large gradients. This reduces the truncation error in the

numerical solution. The ability of the adaptive grid to

redistribute the grid points into the regions requiring

finer grid resolution on an as needed basis reduces the

number of grid points required to compute the solution.

Thus requiring less storage and computational time for

performing a numerical simulation.

In the front zone, sharp density gradient often exists

over a spatial scale of tens of meters. Since the frontal

dynamics are affected by the large-scale estuarine and

coastal circulations, it is essential to model both the

large-scale circulation and the front-scale dynamics. If a

non-adaptive grid is used. very fine grid resolution in the

entire computational domain is necessary to obtain high

accuracy solution, thus requiring large computer storage

and huge CPU time. In thus cases, the adaptive grid

technique can be employed to refine the grid locally while

keeping the total grid point number unchanged, to obtain

accurate solutions at the front zone. This should lead to

savings in computer storage and CPU time.

4.1 Grid Adaptive Algorithm

The adaptive grid scheme developed here is developed

based on the idea of "minimal moment" of the moment-

balancing principle according to Connett et al.(1988). It

belongs to the category of equidistribution adaptive grid

method (Thompson et al., 1985).

Let a segment L in a one-demensional grid have a mass

V. then the density P of this segment (not the density of a

fluid) is W/L (see Figure 4.1). The movement of a grid

point, i.e. the center point in the segment L is determined

by the two adjacent grid points and the densities assigned

to the two segments (Ax, and Ax,l,). The algorithm is set up

to move a single grid point at one time. Once the local

algorithm is defined, the entire grid is traversed. The

traversals are repeated until a satisfactory adaption of

overall grid is achieved. The local algorithm is discussed

in the following.


Pi-1/2 Pi+1/2

xi_- Ax, x, Axi+i xi+l


Figure 4.1 Initial Mass Distribution

If the density P is uniform over the entire segment L,

the mass center is on the segment center. Then the moment

of the segment L is minimized.

When the density on either side of x, i.e. -x. and

-x.. is changed. the mass center should be chanced. If

P,. on -:<. is larger than P,./\. the mass center should be

moved to the left to a ne' position x:,) ( Figure 4.2).


new mass center (x.) new

old mass center (x,),od

'-1/2 I Wi+l/2

Figure 4.2 Moment Balance Principle

As shown in Figure 4.2, w.i;/ (wi-1/2= Pil/ 2xi)

w i/ (wi+1/2-=Pi+/12) are weight functions on Ax, and -Ax+i .


Let us define

10o=11+1,= (Ax,+Axxi)


According to moment-balance principle,

wi-l/2 11= "i+ /2 12-


SWil!/2 I ',-k/'

Wi+ /2

-Wi+ (' ." .-i) (4. 1- 3)


S- -) (4.1.4)
-1 2 w-/*+

Therefore. the mass center will be shifted from (xi)o, to

(x) new

(xi) neUw=x .i-I+1- + l

=(Xi)o.1d -A'NSi i+2 :>:>

=(x,)old xi + li (4.1.5)

The adaption will be advanced from i=l to i=Imar. The

adjustment of each point is based on the previous point

positions. Then repeating the same procedure to find 14 and.

(x,)ne. until I|(Ax,)new-(.Ax.i)o dIe (here e is a small value).

4.2 Weight Function

The main idea of the adaptive grid method is to have

the grid points move as the physical solution developed. A

sensor S has to be defined on the grid points, and the grid

is driven to respond to the magnitude of the gradients of

S, so as to obtain small cells in the presence of the

larger gradients of S and large cells in the presence of

the small gradients. Usually, there are many quantities to

be solved in a physical problem. The quantity which has

discontinuty is chosen as the sensor S(r). For the front

problem discussed in this paper. salinity is chosen as the

sensor S.

As described in section 4.1, the weight function is

the product of the density and the length of the interval


wi-1/2=Pi,-1/2- xi (4.2.1)

Where, the density function P is some kind of combination

of the first derivative of the sensor S(z), S, and the second

derivative of the sensor S(x). Si, so that P can present the

large variation of the salinity in solutions. The selection

of P for the adaptive grid method is somewhat arbitrary and

subjective. Two types of P are used in this study:

P=(1 +eYts, js.l) (4.2.2)

p l +aIlS I+3|SI (4.2.3)

The arbitrary constants a and 3 are the adaptation

parameters. When a increases, the sensitivity of P to Szz

will increase. When 3 increases, the sensitivity of P. to S,

will increase. There is no theoretical way to chose a and

3. The determination of their values are based on numerical


Another influence on the sensitivity of P is the

formular of P selected. If P is too sensitive with the

variation of ISzzj and ISZi like the exponential form Eq.

(4.2.2). the adaption may not converge and may sometimes be

unstable. Usually, the formula of P is chosen as polynomial

type (Thompson et al., 1985. Shyy, 1987), Eq. (4.2.3), so

that it is not too sensitive to select the adaption

parameters a and 0.

Satisfactory adaption generally can not be achieved in

one step. During the period of repeated adaptions, new sets

of grid system continue to be built up, by using new weight

functions which are interpolated from the old values.

4.3 Steady-State vs. Transient Problems

Equidistribution adaptive method is very effective for

solving steady state problem. For unsteady case, the grid

system always follow the solution. Unfortunately, since

the adaption is performed based on knowledge of the

previous solutions, the region where large gradient of the

solution occurs moves forward with the development of time.

Thus, the refined grid region always lags behind the new

large-gradient region. Therefore the adaptive effects are

not obvious for this unsteady flow problems.

If the adaption can capture the movement of the large

gradient or the refined grid region can catch up with

large-gradient region, this grid system will be able to

simulate the large variation of the solution solved in next

time step. Thus, an effective adaption can be achieved. To

catch the movement of large gradient, the moving grid

method (Eiseman and Bockelie. 1989) can be used. However,

the moving grid method is more complicated and requires

more CPU time than the equidistribution method. In this

work, an alternative method. which is basically a modified

equidistribution method. is developed to refine the grid

points in large-gradients region.

In general, weight function includes the influence of

the first derivative and the second derivative of the

sensor S. However, at the downstream region of the front,

there is no influence from the upstream region, hence the

two derivatives are very small or even zero, i.e. the weight

functions are very small. In this case, the adaptive grid

cells are always coarse and the fine grid cells lag behind

the sharp "head" zone of the front. Therefore, there must

be a weight function correction to transfer the influence

of the two derivative in front "head" zone forward to the

downstream ambient area. The following weight function at

point i includes the influence of weight functions at two

previous points i-i and i-2.

w,= (Aw,_2+Bwi_+Cwi) (4.3.1)

where A, B and C are arbitrarily chosen constants.

4.4 Temporal Coupling

Between every adaption, interpolations have to be

made. If adaption advances through the entire computation

period, the solution should be satisfactory. However,

interpolation error will always exist whenever adaption is

employed. In order to reduce the interpolation error,

adaption will not be performed at every time step. Starting

from a particular moment, the solutions obtained on a fixed

grid are used as previous knowledge and the adaption

begins. Over a certain time interval, adaption is performed

every time step. At every step, a new grid is developed

based on known solutions. The known solutions are then

interpolated onto the new grid for marching into the next

time step. By successively repeating this procedure, the

forward matching front is captured on the adaptive grid

system until the desired final time instant.

For front problems, the layer in which the salinity

has a typical horizontal gradient distribution representing


a front, is chosen as the reference layer. The sensitive

variable S is then taken from this layer. After adaption.

all variables are interpolated onto the newly established

grid system.

The grid cells will become more and more clustered as

the adaption period increases in time. In order to avoid

the development of infinitesimally small grid cells, a

minimum grid size is chosen to be approximately 10% of the

initial fixed srid -size. and 3 can also affect the

smoothness of the grid cells.

4.5 Verification of The Adaptive Grid System

The one-dimensional Burgers equation is used to verify

the adaptive grid method in capturing the discontinuity.

Three cases are provided here to present the steady state

and transient solutions. The transient case has the same

properties as the front problem.

Consider the Burgers equation.

au 4.ILu (4.5 .1)

The initial conditions are

Case (1)

1 x=0

u = linear O
-1 x=1

Case (2) u = cos[(2k+1),rxl O

Case (3) 1 x=O

u = ^(4.5.4)

0 x>O

In case (1) and case (2). a steady-state solution can

be obtained, whereas case (3) is essentially a transient


Figure 4.3 through 4.6 present the comparison between

the solutions on adaptive grid and fixed grid. For cases

(1) and (2), the weight functions are exponential types,

and solutions are given in Figures 4.3 and 4.4.

respectively. The discontinuity captured in the adaptive

grid with 20 points is as sharp as that in the fixed grid

with 80 points (See Figure 4.6(a)). Figure 4.3 through 4.6

also show the evolution of the solution from the transient

to a steady state. The grid points move with the variation

of the solution until the final steady state is arrived.

Figure 4.5 presents the marching discontinuity. A

polynomial type weight function is employed in this

adaptive grid. The discontinuity captured in adaptive grid

with 20 points is as sharper as that in fixed grid with 40

points shown in Figure 4.6(b).































II ME 0.1193

-I__fl.-LLULLLUU LLJ Jlll__-


I I L I). IJ iUJl


l'i gi;re -1.3 (Cas (1) ,in iLda.pt i ve yi icl I itinl.s 20()

TIE- 0.0000

I I 11111I1I1 II 11





, 0.00


TIHE- 0.25011


I IHE- 0. 19'16

1.001 --


















I*'igis' p -i.1 Ca~s. c (2) oni adiaplt-. iv<: grl Id"| iints.. "-0










T IHE-0. 11167

TIME-0. 7631

I I I 1 I I I

- IJ_ _J _J _L J .. _





: 0.00

-I .n


-z.oo LlIIIIII I I





: 0.00





I IHE- 0.372'1


I'igi're '1. 5 Case .(3) oil adaptLive' g i(I point.H 20

T IME-0.0000




















IIHE-0. 1101



iTi E- 0.2216

1 I li l i t IIllIll I I I I




tx 21.


IX 'iI .

(aL) Case (1) on uniform grid 20 land 80

(b) Ca u (3:) on uni orm grid 20 "ni( -10

'a'lT u coijv e riid l ion ( 1 ii i in i ,gri
adaptive grid





o 0.00































NX- 1.

1I1111111111 II II 11 111111 II ill lll III IIIIIIlII1111 1 II Ill1

Figure I4.0(i


Front problems are very complicated. In general, there

are no theoretical solutions for the real physical

problems. Some theoretical solutions and laboratory

experiment results have been obtained by simplifying the

real problems. In this chapter, four kinds of physical

problems will be solved by the numerical model described

earlier in this thesis. Three of the problems contain

simplifying assumptions while the fouth one is more

realistic. Whenever possible, the solutions are compared

with theoretical solutions, laboratory experiment results

and field observations.

5.1 "Breaking Dam" Problems

5.1.1 Physical Problems

Two fluids with different densities can be adjusted

under gravity from initial unbalanced state to final

equilibrium. This has been an old but a very typical topic

discussed by ocean and atmosphere scientists since the

beginning of 1900's.

Water spreading out from a suddenly broken dam becomes a

plume with the water front at its boundary. Due to its

similar physical feature to the suddenly released fluid

flowing into the surrounding ambient fluid with a

vertically varying density structure, the "breaking dam"

problem can be regarded as a simplified front problem.


Figure 5.1 Initial State of Dam

Assumming the "'breaking dam" flow is one-dimensional

i.e. along the horizontal direction, and neglecting viscosity

and diffusion terms. equation (3.2.1) and (3.2.2) become

8uD+D =0
x Oax

a D uuD -fD= 0D
at +x =O

at dx







x>O. t=O

If Coriolis effect

equations are reduced to

is neglected,

the governing

OD, u=uD_
Ot ax

uD+9 uuDl 8D2
Ot Ox 2 .O

For these simplified equations,

solutions are given in the Appendix.

the theoretical

5.1.2 Numerical Results

In the "Breaking Dam" problem. the finite difference

equations are solved with spatial increament of Ax=0.025.

The time steps are limited by the following condition:


At < CFL( Ax )

where c=IgD and the Courant number CFL=0.6.

Figure 5.2(a) presents the analytic solutions at 10

time instants. The numerical solutions. of equation (5.1.5)




and (5.1.6) are shown in Figure 5.2(b). As can be seen from

the figures, the numerical solutions are in excellent

agreement with the theoretical solutions.

Figure 5.2 shows that immediately after the release of

the initial discontinuity, a hydraulic jump forms. The wave

fronts move away from the initial discontinuity at the

constant speed. Due to the hydraulic balance and nonlinear

effects, the discontinuity propagates like a compression

wave towards lower pressure region. and a rarefaction wave

travels in the opposite direction into the high pressure


Figures 5.3(a) and 5.3(b) demonstrate the effect of

Coriolis terms on the numerical solutions at 10 time

instants and 25 time instants respectively. Figures 5.4 and

5.5 corresponding to the Figures 5.2 and 5.3 present the

numerical solutions at 6 time instants and 3 time instants

respectively. As pointed out by Rossby(1938), the Coriolis

force simply balances the pressure gradient. Based on this

geostrophic balance, the gravitational adjustment process

leads ultimately to a steady state.

Rotation effects contribute an additional dispersive

term to the original one. i.e.

W2= J +k'c (5.1.8)

The wave with this dispersion relation is reffered to as

the "Poincare Wave" (Stoker. 1957). The wave front moving

away from the initial step is followed by a "wake" of waves

that trail behind because of the dispersion. The decay of

the rightward propagating jump is caused by a rarefaction

of Poincare waves behind the compression. The leftward

propagating jump forms also by the trailing wake of

Poincare waves which increase the pressure after the

passage uf the front. At t=6. the second jump appears

because of the nonlinear steepening of the Poincare wave.

From the solutions, it can be seen that there is

energy loss between the initial state and the final

geostrophic equilibrium. The wave fronts carry energy, so

for any finite region the energy is lost through the sides

by "radiation" of Poincare waves until the only energy left

is that to hold the final balance. Cahn(1945), Gill(1982)

and O'Donnell(1988) discussed the linear and nonlinear

rotation adjustment problems by theoretical and numerical


5.2 Flow With Linear Density Distribution

5.2.1 Physical Problem

When a horizontal density gradient exists in a fluid,

vertical gravitational field vorticity is generated, and

109 ----- -__T9 .-.r-- ...

0 --- U ....... .....

S------ ------ -

3 U ------------ -

3 -- -- 3 -

0 | .---- _______ -- .--- ---- -- -- -L--- --

-12 -10 -8 -6 -4 -2 0 2 i 6 8 10 12 -12 -t0 8 -6 -4 -2 0 2 1
x x

(u) (l,)

Figure 5.2 "nreakinig dla" problems

(a) Ana.lytlic -olutionii at unit. int.urval o( I' t.iIiie

(1b) N ntiiurical rol utiioni c:o -r,<; pi)",, i ngu t,< (i.

' I I I I

S .

_-.6 I 10 12_ I_
6 0] 10 12


'I I- I S *I 5 V11(~ 1rgyr


-1 I -10 -8 -L -
-12 -10 -8 -6 -4 -2



2 1 6 10. I



20 ------

10 --- -- -


6 __--------'--

I -------- --- --

S _--__...._.. ... . .
0- -- --- -~___ -- -


- 1 -- 1 --I --1- 1 J-I 1.._ __.l .... -. I._ l._.l 1
12 -12 -10 -8 -6 -4 -2 0 2 I 6 0 10 12

( ')

Figure 5.3 'The 1numerical solutions F'or 1.1t,- layt.: (lt:lt).i itl.l
rIo,.i.ta ion clffct;s in (a) 1.:z 10 i(I (I>) t.=2-

Figure 5.4

Nonratating "breaking dam" problem at t=l and t=6

Figure 5.5

With rotating effects, the gravitive adjustment
at t=1 and t=3

flow is induced no matter how small the density gradients.

As an example, an uniform horizontal density gradient

distribution producing fluid motion in a rectangular tank

is discussed.

Figure 5.6 Density Distribution

As shown in Figure 5.2-1, the density distribution

follows p=p(1--ax), where a=!p is a measure of the

horizontal density gradient and ~ is the value of the

largest density in the tank. The length and depth of the

tank are 360cm and 15cm, respectively. Since the flow is

two-dimensional, the width of the tank is eliminated from

the equation. Simpson and Linden(1988) presented their

theoretical and laboratory results under the same


5.2.2 Numerical Results

The numerical experiment is solved with spatial

increment Ax=10cm and Az=1.5cm. The time step is


Figure 5.7 is the numerical contour plots of salinity

obtained at 8.0s, 12.3s. 16.Ts, 20.9s and 24.9s after

initial state with =1.058 x 10-3g/cm-3. The vorticity

generated by the horizontal density produces a flow towards

the right in the lower half of the tank and towards the

left in the upper half. The vertical stratification

increases with the development of time. The isopycnals

remain straight throughout the motion. Corresponding to the

numerical experiment. Figure 5.3 shows a time sequence of

photographs of laboratory experiment provided by Simpson

and Linden(1988). As can be seen, the numerical results and

laboratory experiment have good agreement. In figure 5.3,

the angle of inclination 9 'of the isopycnals to the

vertical for theoretical, numerical and laboratory results

is plotted as a function of time t. In this log-log plot,

the time has been non-dimensionalized by (.ga)2 According

to Simpson and Linden's(1988) theoretical computation, the

initial angle of inclination is given by tanO=,. The

laboratory experiment is realized by suddenly rotating the

tank from vertical position to horizontal position so that

the motion begins with an initial 0. In the numerical

experiments, the numerical results are recorded from a

initial 0 which is the same value as that in laboratory

experiment. Two sets of numerical results are obtained with

a=6.5x10-scm-l, p=I .058x10-g/cm3 and a=15x10-5cm-, p=1 .0

r~ ----=--=---~==-L~=_~-===---,-_-

IFigiire 5.7 'Nimtecicxl t;imI-sec(tionce showing t.Iicisociias
Wit.ii f= 1.058x 103gl/cn 3, a= I5X 1 0-S CI -
atl. t.=8.0s" 12.3s, 16.7s, 20.9s and 2,1.9s.!)s


r (



II. U I1Y*.)S ** &

Figure 5.8

Corresponding to Figure 5.7 a t.ime-scquc: showing the motion of the isopyca: ls in
laboratory experiments from Simpson iULIl
Lini(Ien(1988). Only the right-hIn tl I( l" f1' tIl O i
tank is visible and the vertical I i nes are iLt.
10cm intervals

.- ---. r=15x10-cm-I
o o 0 I.ab 1
X x x Lab 2

-. 0 0

- -- 0



-0.60 -0.20 0.20 0.60

log(' )

F igture 5 .

A log-log plot of the angle 0 of iicn:l ii t.iion (,1
the isopycnals to lhlie ver.tl i:ic l againsitl t.ill:
cl imensionless t =(ly/2) ti




0.8 -

0.4 -



g/cm3 are shown in the Figure 5.2-4. rhe numerical results

are identical to the laboratory results and theoretical

computation. In the numerical experiments, the top layer of

the tank is treated as a free surface whereas it is rigid

boundary in actual laboratory experiment. Therefore, some

error is brought into the numerical solutions. For rigid

boundary, the forms of governing equations will be changed.

Thus, there are must be the changes in the numerical model.

5.3 Sudden Release Problem

5.3.1 Physical Problem

When a fluid is suddenly released into the ambient

fluid which has a vertically varying density structure, a

gravity front is created. In the laboratory experiment, the

gravity front is formed by releasing fluid of a different

density from behind lock-gates. Since the volume of the

fluid behind the gate is usually much smaller than that of

the other fluid, the experiment approximates the release of

a finite volume of fluid into another fluid of infinite

volume. Rottman and Simpson(1983) presented their

laboratory results. In order to compare with the laboratory

results, four test cases are simulated. The first two cases

are in the same conditions as in the laboratory

experiments, which represent bottom plumes.

Figure 5.10 Laborat

Where xo=50 cm

ho=7 cm

Po=1 .048

ory Experiment Case 1

L=150 cm

H=49 cm

g/cm p1,= g/cm3

Figure 5.11 Laboratory Experiment Case 2

Where xo=50 cm L=150 cm

h0=7 cm H=7 cm

po=1.048 g/cm pl=l g/cm3

For realistic large-scale river and ocean flows, there

Case 1:

Case 2:

are not only bottom plumes, but also surface bouyant

plumes. Therefore. the next t-'o cases describe the Large

scale bottom and surface gravity currents.

Case 3:.


H !

Figure 5.12

Large Scale Bottom Plume

Where x0=3 km L=40 km

ho=7 m H=10 m

Po=1.000546 g/crn

p,=l g/cm3

Case 4:


Pi H


Figure 5.13 Large Scale Surface Plume

where xo=3 kmn L=40 km

ho=T m H=10 m

Po=1 g/cm

p =1.000546 g/cmn3

5.3.2 Numerical Results

For the sudden release problem, the numerical

solutions are given in Figure 5.16 to Figure 5.19. From the

figures, it appears that the two-dimensional non-rotating

gravity current produced by an instantaneous release passes

through two distinct phases. There is an initial adjustment

phase during which the initial conditions are important,

and an eventual self-similar phase. In the initial phase

the "nose" velocity, i.e. the velocity of gravity-current

front, is about constant. In the self-similar phase, the

nose velocity decreases. The transition from the first to

the second phase occurs when a disturbance generated at the

end wall overtakes the front.

Figure 5.14 and Figure 5.15 present Rottman and

Simpson's(1983) laboratory experimental results, which

correspond to Figure 5.10 and 5.11. The numerical solutions

in Figures 5.16 and 5.17 show the same density front

formation as obtained in the laboratory experiments. The

bores, surges and hydraulic jumps exist during the front


propagation. The comparison of computed and measured front

velocities are shown in Table 1.

For realistic large-scale coastal or river density

plume. the numerical solutions shown in Figure 5.18 and

5.19 also can give the same results about the front


The solution shown in Figures 5.18 and 5.19 demonstrate

that the hydraulic jumps occur in both bottom and surface

currents and the two phases of fronts are obvious. The

front velocities are given in Table 2 in comparison with

the laboratory experiment and theoretical solutions. The

theoretical solutions which are solved by characteristic

method are from Rottman and' Simpson(1983). The front

velocities follow the relation

u =F(I) (5.3.1)

Theoretically, the bottom front has larger velocity

than surface front due to the reduced acceleration of

gravity g' g. This can be found from table 2 which gives

both bottom and surface front velocities at the same

initial conditions.

Front Velocities

Table 2. Bottom


Surface Front velocity y

Wang(1984) investigated both surface and bottom front

problems by his numerical model. His results are shown in

Figure 5.20 and 5.21. In comparison with his results,

solutions obtained by using current model are shown in

Figure 5.22 and 5.23.

5.4 Water Discharge Problem

5.4.1 Physical Problem

The discharge water of a different density into

Case Theory Lab Numerical

1 11.79 12.24 12.5

2 8.16 5.S 6.9

Case Theory Lab Numerical

3 8.709 9.68 10.68

4 8.707 9.67 8.47

Table 1.



figure 5. 14' Shado(wgralphs of a. vol ei of Isal 1t. wa ;.t:r w i li
Iho=7cm, II=49cm,, xo=50cm, g' =17cms- col li apsi i ntl
freshly wat.t'r at 2.2s, 3.8s anld 6.() a; t.(:r I.':l:;at.-i
The dlLtel v'ert.ical line indicat. v th.lie ipsition
oC t*lhe I ock galte, tlie uIe wal I his al. f;ir left., i;nld
-th th in v;:rt. ical linlies .ar at. I Ol ().n ii n4ti e val'i
(llotlrnani and Simpson (1983))

Figre 5. 15 Shadowgraphs of a voliae of sait. wtt r Wi Li*
. /11 x* .0 17c 'n .
'. .' *

fl I

Figure 5.15 Shadowgraphs of a volume of salt waltr wiitlh
ho/|i/=l ho==7cm, xo=50cm, g'=417cus-" col lnplsing
into fresh water at 4.7s, 7.7s anmd 10.7s a Tt cr
release (Ilottmlan arid Simpson (1983))

1 _/ ii ___i iii~iii~ __ ~_~__C_____~__I_____I _~ ____i







i'igure 5.16 NPumerical solutions correspol tlnding to lF'i gllic 5 .1 1

v I I Ir E 0

.*,,-- ------

L. .. ... .].. . . .

. .

2 .................
.I... .... .

....... .......
.. .. .
r-~i~i: ..P~r........ .... .............. ..

...... .... ..
;ii F'!I"."M--~-~1'1M02-121% M,'----





Figure 5.17 Numerical solutions corresponding to Figure 5.15

i 7-7 1 -T i E
---- --, ------- -

0 40 km
___________________ ________- ______ -- !t I t sii :

-10.m ... .
.. "2 ,-: -------,-,--

0.5540 0.5080 0. i620 0.4160

Figure 5.13 Numerical solutions of Case 3
at 20hrs. 40hrs, 60hrs and 80hrs

- - - - -
b7- .. .... .......

SI I t J~ I I 1 I I t i ii I t L I [ I i -

I- -

::-: 4

I C [ I. It i~I i I. i I- I~ .i.. i L .i i i t. i. i. '-..

l_ i t Fl~~r LrE I F t I I I It I i

...... ............. ....... ......

.- : . : : : : . . .. . . . . . .... .. . . . . ....~-~~~:~. ~~
.......... ...... ...~:~- ~:
. . .
. . . .





Figure 5.19 Numerical solutions of Case 4 at 20hrs,40ihvs, 60hrs
and 80hrs


o 10 20 30

0o o 20 so

to 20 10

Figure 5.20 Numerical Solution for Surface Plume
at 20hrs. 40hrs, 60hrs and 80hrs by Wang(1984)

0 t0 20 i0













2 0 30

a 10 20 30

Figure 5.21 Numerical Solution for Bottom Plume
at 20hrs. 40hrs, 60hrs and 80hrs by Wang(1984)


Figure 5.22 Numerical .Solution for Surface Plume at 20hrs,
40hrs, 60hrs and SOhrs Coresponding to 5.3-11

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

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