Uchiyama, Y., Dec., 2004: Rep. Port and Airport Res. Inst., Yokosuka, Japan, pp.3-21.
Copyright © 2004 Port and Airport Research Institute. All Rights Reserved.

REPORT OF THE PORT AND
AIRPORT RESEARCH INSTITUTE
Vol.43, No.4 (Dec., 2004)

Modeling wetting and drying scheme based on an extended logarithmic law
for a three-dimensional sigma-coordinate coastal ocean model

Yusuke UCHIYAMA*

Synopsis
A wetting and drying scheme (WDS) is developed to incorporate into a three-dimensional sigmacoordinate coastal ocean model (POM; Blumberg and Mellor, 1983) with proposing an extended logarithmic
profile within and above bed boundary layers to accurately simulate 3D hydrodynamics in extremely shallow coastal seas. While the WDS is simplistically formulated, it can selectively retain the computational
robustness and water mass conservation in the domain.

The developed model, WD-POM, is verified tohave a good agreement with the analytical solution for non-linear long-waves propagating on a sloping
beach proposed by Carrier and Greenspan (1958). WD-POM is also applied to a tidal simulation in San
Francisco Bay, California, USA, which is encompassed by huge intertidal areas. The model-generated surface elevations and tidal currents are compared to the observed data-set by NOS-NOAA (2003) and demonstrate that WD-POM is capable of precisely reproducing the estuarine intertidal hydrodynamics in San Francisco Bay.
The influence of intertidal topographies on estuarine hydrodynamics is next examined, by comparing between the simulated result with WDS and that without WDS. The wetting and drying process is represented to significantly alter development of tidal waves and a current field in San Francisco Bay. The intertidal sloping bathymetry appears to enhance refraction and shoaling of tidal waves propagating
anticlockwise as Kelvin waves. The intertidal effects are evidently exhibited to modify the amplitude andphase speed of the tidal waves, and to shift the current direction to be shore-normal.
Key Words: wetting and drying, intertidal flat, estuary, hydrodynamics, bed boundary layer* Senior Researcher, Littoral Drift Division, Marine Environment and Engineering Department
3-1-1 Nagase, Yokosuka 239-0826, Japan
phone: +81-46-844-5045, facsimile: +81-46-841-9812, email: uchiyama@pari.go.jp-3-

Uchiyama, Y., Dec., 2004: Rep. Port and Airport Res. Inst., Yokosuka, Japan, pp.3-21.
Copyright © 2004 Port and Airport Research Institute. All Rights Reserved.

港湾空港技術研究所報告
第43巻第4号（2004.12）

拡張対数則を導入した干出・冠水スキームの開発と3次元σ座標海洋モデルへの適用

内山 雄介*

要旨

3次元σ座標系海洋流動モデル（POM; Blumberg and Mellor, 1983）に適用可能な干出・冠水スキーム（WDS）を開発した．WDSでは，拡張対数則を新たに導入することにより，干出域を含む潮間帯において水深が極めて浅くなった場合でも底面境界層周辺の流速場を正確に再現することが可能となっている．このWDSはシンプルな形で定式化されているが，領域全体の質量保存と計算上のロバストネスを選択的に確保することができるようにモデリングされている．開発された3次元モデル（WD-POM）の再現性を確認するため，計算結果をCarrier and Greenspan（1958）による
非線形長波の斜面遡上問題に対する理論解と比較したところ，両者は良好に一致した．次いで，広大な潮間帯を有する内湾海域である米国サンフランシスコ湾の潮流計算にWD-POMを適用した
ところ，モデルによる潮位変動と潮流速変動がNOS-NOAA（2003）による観測結果とほぼ完全に一致したことから，WD-POMが潮間帯を含む内湾域の流動を忠実に再現することが可能であることが確認された．
潮間帯の存在が内湾流動におよぼす影響を照査するため，サンフランシスコ湾の地形に対してWDSを適用した場合としない場合の計算を実施し，両者の結果を比較した．潮間帯の存在は湾内における潮汐波の伝播と沿岸部における潮汐残差流の形成に対して重大な影響を及ぼしていることが実証された．すなわち，ケルビン波として湾内を反時計回りに伝播する潮汐波は，潮間帯上の斜面を遡上する際に屈折変形を受けて等深線に垂直な流速成分を生じると同時に，浅水変形により潮位変動振幅が増大し，逆に底面摩擦の影響が相対的に大きくなるため潮汐波伝播の位相が遅れることなどが明らかとなった．
キーワード：冠水・干出，潮間帯干潟，内湾海域，流動，底面境界層* 海洋・水工部〒239-0826

漂砂研究室 主任研究官

神奈川県横須賀市長瀬3-1-1

電話：046-844-5045，FAX：046-841-9812，Eメール：uchiyama@pari.go.jp-4-Uchiyama, Y., Dec., 2004: Rep. Port and Airport Res. Inst., Yokosuka, Japan, pp.3-21.
Copyright © 2004 Port and Airport Research Institute. All Rights Reserved.

CONTENTS

Synopsis·············································································································································3
1. INTRODUCTION························································································································7
2. METHODS··································································································································92.1. Brief description of Princeton Ocean Model··············································································92.2. Wetting and drying scheme······································································································92.3. Extended logarithmic law for bed boundary layer respresentation··············································10
3. RESULTS··································································································································13
3.1. Comparison to the analytical solution·······················································································13
3.2. Tidal simulation in San Francisco Bay·····················································································14
4. DISCUSSION·····························································································································15
5. CONCLUSIONS························································································································18
ACKNOWLEDGEMENTS···············································································································18
REFERENCES·································································································································19-5-Uchiyama, Y., Dec., 2004: Rep. Port and Airport Res. Inst., Yokosuka, Japan, pp.3-21.
Copyright © 2004 Port and Airport Research Institute. All Rights Reserved.

-6-

Uchiyama, Y., Dec., 2004: Rep. Port and Airport Res. Inst., Yokosuka, Japan, pp.3-21.
Copyright © 2004 Port and Airport Research Institute. All Rights Reserved.

tions are necessary and therefore time-consuming in

1. INTRODUCTION

terms of computation. In the wet-dry point treatment
Estuarine intertidal mudflats and salt marshes havemethod, the computational domain covers the entire ba-

been recognized to play a significant role in marine envi-

sin of interest including intertidal area.

ronment and aquatic ecosystems.

Efforts have beenboundary between the land and sea, where water flux

made to investigate hydrodynamics and associated sedi-equals zero, is determined for each time step by the wet-ment transport as well as chemical substances attached todry point treatment technique.

the sediment particles on intertidal areas (Kuwae et al.,are defined using the local total depth: when the depth is1998; Austen et al., 1999; Van Der Lee, 2000; Le Hir, etpositive, it is defined as a wet point; otherwise it is a dry

al. 2000; Uchiyama et al., 2001; Talke and Stacey,

point.

2003). A key factor dominating intertidal hydrodynam-

averaged 2D numerical models for estuaries including

ics is a wetting and drying process as it affects suction ofintertidal areas (Leendertse, 1970, 1987; Flather andpore water in bed sediments influencing on sedimentHeaps, 1975; Leclerc et al., 1990 Cheng et al., 1993;

suspension, evaporation related to thermal and salinity

Luettich, and Westerink, 1995a, 1995b, 1995c; Hubbertenvironments of overlaying seawater, rapid movement ofand McInnes, 1998; Ip et al., 1998; Part et al., 2002;

tidal currents resulting in formation of highly turbid

Davis and Sheng, 2003). Applying the wet-dry point“tidal bore” on the bed, and so on (Whitehouse et al.,technique to 3D primitive equation models requires that

2000). These facts imply that it is important to resolve

the equations satisfy a finite value solution when the to-the vertical structure within and above the bed boundary

tal depth approaches zero. In the z-coordinate system,

layer in order to precisely assess intertidal environments

the method used in the 2D model can be applied easily toincluding sediment suspension and deposition.

the 3D case without any technical difficulties (Casulli

The wet and dry points

This method has been widely used in the depth-

and Cheng, 1991; Casulli and Cattani, 1994; Hervouet

Numerical experiments have been frequently per-and Janin, 1994; Lin and Falconer, 1997; Gross et al.,formed to investigate estuarine hydrodynamics and1999).

transport processes, whereas the wetting and drying

However, to ensure numerical stability, thethickness of the layer closest to the surface must beprocess in intertidal areas is often neglected. A proper

greater than the amplitude of the tidal oscillation (Davies

approach to replicate the flooding-draining processes,

et al., 1997; Lin and Falconer, 1997) as shown in Figure

i.e., submergence at high water and emergence at low

1 (a). This constraint results in simulations difficult towater, must thus be developed. In terms of mathemat-

use the wet-dry point method in an extremely shallow es-

ics, however, this seems intricate because the computa-

tuary like San Francisco Bay where the median depth istional domain varies with time if intertidal process must

about 2m whereas the maximum tidal elevation reaches

be taken into account (Cheng et al., 1993). Attempts

about 2.5m since it is impossible to include the verticalhave been made to solve this problem for the last two orhydrodynamic structure that is frequently essential tothree decades and two different approaches have been

predict the bed shear stresses and baroclinic motion on

developed: a moving boundary method and a wet-dry
point treatment method.

The moving

shallow basins in stratified estuaries.

In the moving boundary

system, is reconfigured for each time step according to

The σ-transformation has been broadly used in coastal
ocean models because it is well-suited for varying bot-the boundary defined by a zero total depth and a zero

tom topography as indicated in Figure 1 (b) (Blumberg

normal transport at an interface between land and waterand Mellor, 1983, 1987; Hedivogal et al., 1991;(Lynch and Grey, 1980; Siden and Lynch, 1988; Austria

Shchepetkin and McWilliams, 1998, 2003). This trans-

and Aldama, 1990; Shi, 1995; Hervouet and Van Haren,formation is no longer valid if the total depth becomes

1996; Bates and Hervouet, 1998). Although this ap-

zero.

proach is more precise than others, it has been used only

treatment program in the σ-transformation model because of this limitation. Two methods have been

method, the computational domain, or the model gridfor some idealized estuarine cases because grid regenera--7-

It seems hard to incorporate a wet-dry point

Uchiyama, Y., Dec., 2004: Rep. Port and Airport Res. Inst., Yokosuka, Japan, pp.3-21.
Copyright © 2004 Port and Airport Research Institute. All Rights Reserved.

z=0

z=0

h

draineddewateringhbathymetry
expressed insimulations(a) z-coordinate models
σ=0
h

z=0

z=0
h

dewatering
z

=

H(

σ= 0
drainedx,y)bathymetry
expressed insimulations
σ = -1

σ = -1(b) σ-coordinate models

Figure 1: Schematic diagram of wetting and drying processes expressed by (a) z-coordinate models employed by Lin
and Falconer (1997) and (b) σ-coordinate models by Xie et al. (2003).
proposed to avoid the occurrence of singular points when

EFDC (Hamrick, 1992), they are considered to containthe local water depth approaches zero in the area of dry

problems in the treatment of representing the bed bound-points.

dry points as a result of infinite horizontal and verticalary condition, since σ = −1 is defined at the bed surface,based upon the conventional logarithmic law of the wall.

diffusion (Xie et al., 2003). This approach is simple

For instance, suppose the total depth becomes 10cm andbut it basically does not guarantee the volume conservation. The second method introduces an additional bot-10 vertical σ layers are used in POM simulations which
employ the Arakawa-C grid with uniform grid spacing,

tom boundary layer and redefines the local water depth

then the height of the horizontal velocities at the lowest

as a sum of the real depth and the thickness of the bottom

grid point becomes 0.5cm above the bed while the typi-boundary layer as proposed by Ip et al. (1998) andcal value for the bed roughness height in estuaries may

Zheng et al. (2003; ZH03, hereafter).

usually be 1.0cm.

The first one assumes that the currents vanish at

The second

This means that the velocities closest

method seems accurate in mass conservation as com-

to the bed are defined below the roughness height, andpared to the first method. ZH03 explains that this addi-

consequently the logarithmic law above the bed is no

tional thin layer, ranging from 1cm to 25cm, could be

longer valid. It is thus required to develop a new meth-considered as a viscous sub-layer (ZH03) or porous me-

odology for reproducing bottom boundary layer precisely

dia (Ip et al., 1998) that lays beneath the logarithmic

in extremely shallow estuaries and intertidal areas.

layer in turbulent boundary layer. However in the for-In the present study, a wetting and drying scheme

mulation of ZH03 as well as the model proposed by Ji et(WDS) based on the wet-dry point method is developedal. (2001), Xie et al. (2003) and Oey (2004), which are

to incorporate into a three-dimensional σ-coordinate

based on POM, ECOM (Blumberg et al., 1993) and-8-

Uchiyama, Y., Dec., 2004: Rep. Port and Airport Res. Inst., Yokosuka, Japan, pp.3-21.
Copyright © 2004 Port and Airport Research Institute. All Rights Reserved.

ocean model to simulate hydrodynamics in estuaries in-heat, salinity and turbulence. Several other techniques

cluding intertidal areas such as mudflats and salt marshes

are also used to conserve numerical accuracy and to(WD-POM).

Attentions are paid particularly to theshorten computational times but they are not essential tomass conservation and hydrodynamic processes near thethe present study (for the detail, see Blumberg and Mel-bed by adapting an extended logarithmic law so as tolor, 1983, 1987).

precisely predict vertical profiles of currents in extremely shallow seas. WD-POM is at first applied to a

2.2 Wetting and drying scheme

run-up problem of non-linear long-waves on a sloping

WD-POM is configured for a pre-defined computa-

beach and is compared to the analytical solution proposed by Carrier and Greenspan (1958).

tional domain that contains both wet and dry grid cells.

Next theA land mask function (LMF) is created to mask out allmodel is used to calculate tidal currents in San Francisco

land cells during model simulations. The LMF is as-

Bay, California, USA, which has large intertidal areas ofsigned zero for all land cells and unity for all water cells.

about 200km2 while the total surface area is 1240km2.

The value of LMF at each grid cell is re-evaluated at

The accuracy of WD-POM is then confirmed with theeach external time step based on a set of wetting-drying

harmonic analysis of the measured and the model-

criteria. In the wetting and drying calculation presented

generated tidal elevations and currents. In addition, be-here, the first step is to define three sub-depth scales, dcr,

cause the influence of intertidal bathymetry on estuarinedmin, and δ, prior to the simulations (Figure 2):hydrodynamics is of primary interest, the model is subsequently applied to simulate tidal waves and currents in
San Francisco Bay.

dcr > dmin ≧ δ ≧ 0.

The model-generated result with

WDS is compared with the result without WDS in order(1)

The second step is to scan the total water depth at eachto exhibit the intertidal topographic effects on hydrody-

grid cell, Di,j = Hi,j + ηi,j, where Hi,j is a vertical elevation

namics in the bay.

in the z-coordinate below MSL (mean sea level) and ηi,j
is a water elevation above MSL in the z-coordinate at a2. METHODS

certain external time step at a horizontal girid cell (i, j).
If Di,j ≤ dcr, the grid cell is regarded as potentially dry
and WDS is immediately applied. In WDS three condi-

2.1 Brief description of Princeton Ocean Model
Princeton Ocean Model (POM) is based on a set oftions are introduced to evaluate each of the potentiallyprimitive equations consisting of the continuity equation,

dry cells to be wet or dry by examining the four adjacent

the 3D Reynolds-averaged Navier-Stokes equations, thegrid cells with;heat and salinity transport equations, Mellor-Yamadaa) min [ηi-1, j , ηi+1, j , ηi, j-1 , ηi, j+1 ] ≤ ηi, j ,level 2.5 turbulence closure model (Mellor and Yamada,b) min [Di-1, j , Di+1, j , Di, j-1 , Di, j+1 ] ≤ dcr , or1982; Galperin et al., 1988) for vertical eddy motion, andc) max [LMFi-1, j , LMFi+1, j , LMFi, j-1 , LMFi, j+1 ] = 0.

the Smagorinsky-type formulation for horizontal eddy
viscosity and diffusivity with the Boussinesq approximation and the hydrostatic assumption.

If at least one of these three conditions is satisfied, theThe governingpotentially dry grid is considered effectively dry and isequations are transformed from the 3D Cartesian coordinate into the 2DH orthogonal curvilinear coordinate andsubsequently removed from the computational domain;

the vertical σ coordinate to smoothly resolve complex
terrestrial geometries and marine bathymetries. The

sidered wet and remains in the computational domain to

the LMF is set to be 0. Otherwise the grid cell is conbe scanned again at the next time step; the LMF is heldmode-split technique is used in POM to solve the gov-

1. When the cell becomes effectively dry, the water

erning equations, comprising the external mode to com-elevation ηi,j is stored in the computer memory and restored in the next flooding. This procedure allows con-pute fast-moving 2DH surface gravity waves and the internal mode to estimate 3D vertical structure of currents,

serving the total water volume in the whole computa--9-Uchiyama, Y., Dec., 2004: Rep. Port and Airport Res. Inst., Yokosuka, Japan, pp.3-21.
Copyright © 2004 Port and Airport Research Institute. All Rights Reserved.

WD-POM by adjusting the value of dmin. An example
shoreline
Di -3 Di -2

lustrated also in Figure 2.
Di-3, Di-2, Di-1 ≤ dcr

Di -1
Di

dra
i

how WDS works in a 2D vertical case is conceptually il-

Di +1 Di +2

2.3 Extended logarithmic law for bed boundary layer
representation

Di, Di+1, Di+2 > dcr
ne
d

i nu
nd
at
ed

The last sub-depth scale in Eqn. (1), δ, is the most important one since it is used to retain the logarithmic velocity profile near the bed even though the total depthdcrδapproaches zero. In other words, δ is necessary to reflect the vertical profiles within and above the bottom

dminboundary layer and is indispensable to correctly estimate

Figure 2: Conceptual illustration of the wetting anddrying scheme presented here for a two dimensionalvertical case. Onshore cells i-3 and i-2 are effectively
dry and the cell i-1 is potentially dry. The condition a)is satisfied at the cell i-1 to be effectively dry while thecell i is still wet. Hence the shoreline (interface between a wet cell and a dry cell) can be defined between
the cell i-1 and i.

tional domain during the simulations.

the bed shear stresses. In the original version of POM
as well as the other versions used in the previous works
by ZH03, Xie et al. (2003) and Oey (2004), the boundary
condition for the horizontal components of 3D currents
at the seabed is expresses as:
K M ⎛ ∂U,⎜D ⎝ ∂σ

[

]

∂V ⎞
2
2 12
(U , V ), σ → −1
⎟ = Cd U + V
∂σ ⎠

(2)

WDS requires

only flux-blocking conditions across the land-sea inter⎤⎡κ2
Cd = MAX ⎢, 0.0025⎥2⎦⎣ [ln{(1 + σ kb −1 )H z0 }]faces when wet cells become dry, while dry cells aretemporarily dormant and are dynamically activated(3)

through the mass and momentum conservation.
where KM is vertical eddy diffusivity, U and V are the

The minimum depth, dmin, is used to avoid the zero

horizontal two components of the 3D velocities, Cd is a

depth because the σ-transformed primitive equations do
not accept zero or negative water depths. When the

bed friction coefficient, κ is the von Karman constant,σkb-1 is the lowest σ level for the horizontal current ve-

horizontal grid spacing is sparse or the time step is set

locities, and z0 is the roughness height usually being 1cm

relatively large, the water surface may rapidly descendin estuaries. In the intertidal areas, variations in the wa-

and be on or below the seabed during an external time

ter surface elevation η may be comparable to or greaterthan the water depth H relative to MSL, so that D should

step. In WD-POM, if the water depth Di,j at a potentially dry cell is rapidly dewatered and becomes less thanbe used instead of H in Eqn. (3). In simulations of es-

dmin during one single time step, the cell is considered ef-

tuarine tidal currents, grid spacing may usually be expo-

fectively dry, the local water depth is held dmin, and LMF

nentially decreased downwards resolve the bed boundary

is set zero. Whereas WDS itself is designed to con-layer; the value of −0.991 for σkb-1 is employed in the
San Francisco Bay simulations described later. The

serve the water mass in the whole computational domain,it may no longer be valid if one employs no-zero dmin.

profile of Cd according to the variable depth can be rep-

However, the simulations become markedly robust withresented by the dotted curve in Figure 3 (a). The value

no-zero dmin if an external time step is set longer or veryof Cd is held constant to be 0.0025 if the total depth D issparse horizontal grid spacing is configured. Introduc-less than several hundred meters since POM is originallying dmin is likely equivalent to defining a thin sub-layer

designed to simulate hydrodynamics in the ocean. Thison the bed as used in ZH03. The present method isexpression is considered to be undesirable for the bound-more flexible since users can control the performance ofary condition at the seabed in the intertidal computations.

-10-Uchiyama, Y., Dec., 2004: Rep. Port and Airport Res. Inst., Yokosuka, Japan, pp.3-21.
Copyright © 2004 Port and Airport Research Institute. All Rights Reserved.

(a) vertical profiles of Cd(b) Velocity profiles on a logarithmic axis

Figure 3: Schematic diagrams of the original (dotted lines) and extended logarithmic (solid lines) profiles near the bed.

becomes infinity) if the total depth is lowered to be extremely shallow. This procedure is explained technically to correspond to elevating the grid system verticallyat δ as displayed in Figure 4. In hydrodynamics no discrepancy is expected because the bed shear stress used in
the bottom boundary condition in POM is evaluated between the bed and the lowest grid points defined at z =(1+σkb-1)D/ z0+δ/z0. In practice the smaller value of δis set, the more accurate estimation of the bed shear
stresses can be achieved; the smallest value of δ =z0 isfavorable as expressed in Eqn. (5).
Cd =

κ2

[ln{(1 + σ kb −1 )D z0 + 1}] 2(5)(c) Velocity profiles on a Cartesian axis

Together with Eqn. (5), leaving the determination of l

Figure 3: Continued.

unchanged, where l is turbulent macroscale, the verticaleddy viscosity KM and diffusivities KH and Kq must be

Therefore an extended logarithmic law represented byslightly altered to close the TKE and macroscale equa-

Eqn. (4) is introduced to reproduce the bottom bound-tions so that:

ary layer properly.
Cd =

κ

2

[ln{(1 + σ kb−1 )D z0 + δ z0 }]2

K M = S M q(l + κδ ) = S M q (l + κz0 )⎫⎪,K H = S H q(l + κδ ) = S H q (l + κz0 ) ⎬K q = S q q (l + κδ ) = S q q (l + κz0 ) ⎪⎭(4)

Equation (4) introduces δ/z0 in the denominator to avoid
the logarithm in Eqn. (3) to be zero or negative (Cd thus

(6)

where SM, SH, Sq and are the stability functions and q issquare root of TKE multiplied by 2, used in the Mellor--11-

Uchiyama, Y., Dec., 2004: Rep. Port and Airport Res. Inst., Yokosuka, Japan, pp.3-21.
Copyright © 2004 Port and Airport Research Institute. All Rights Reserved.

i-1

i

i-1

i+1

i

i+1

σ kb-2σ kb-2σ kb-1σ kb-1(1+σkb-1)Droughnessheight, z0

roughnessheight, z0

δ
(1+σkb-1)D+δ

Figure 4: Definition sketch of δ. In the extended logarithmic law, the grid system is moved upwards at a height of δ.

(a) Superposition of the instantaneous water-surface profiles(b) Temporal variation of the computed shoreline positions
Figure 5: Comparison of the results from WD-POM and CG58.

Yamada level 2.5 sub-model in POM.

The final result

trated in Figure 3 (b), the structure of the bottom bound-

seems similar to that in the 1D model for the oscillatory

ary layer (vertical profiles of the horizontal current ve-boundary layer simulations by Mellor (2002). As illus-

locities) based on the original logarithmic law is much

-12-

Uchiyama, Y., Dec., 2004: Rep. Port and Airport Res. Inst., Yokosuka, Japan, pp.3-21.
Copyright © 2004 Port and Airport Research Institute. All Rights Reserved.

wet

wet

dry

shore-most wet
cell position

dry

actual shoreline defined by
the intersection between
water surface and bed

Figure 6: Definition of the shoreline position in Figure 5 (b) in a vertical 2D case. The model-predicted shore-most
wet cell positions are converted with the linear extrapolation to estimate the actual shoreline positions as an interface
between the wet and dry cells.
different from that based on the extended logarithmic

control the numerical robustness in accordance with thelaw if δ = z0 is employed and (1+σkb-1)D/z0+δ/z0
approaches z0. The velocity profiles estimated with thevalue of dmin. WDS is modeled to precisely reproduce
the benthic boundary layer if the computations are car-original logarithmic law shows that return flow is gener-

ried out for estuaries that have extremely shallow and/or

ated below z0 as shown in Figure 3. The extended loga-

intertidal areas. The numerical solutions derived fromrithmic law yields a little departure in its velocity pro-

WD-POM are next compared to the analytical solution

files from that estimated by the original one when thefor swashing nonlinear long-waves on a linearly sloping

depth decreases.

Note that the extended logarithmic

beach (Carrier and Greenspan, 1958; CG58 hereafter).

law is consistent with hydrodynamics near the seabedbecause the velocity profile estimated with the extended

In the simulation, dcr=1mm, dmin=0 and δ=0 are respectively employed for the three sub-depth scales, and theone coincides with that calculated with the original one if

bed friction is thus neglected in the WD-POM simulation

the profile with the extended one is moved upwards at a

as is in CG58. Figure 5 (a) indicates superposition ofheight of z0 as indicated in Figure 3 (c).

This meanswater surfaces at several phases reproduced by the modelthat the procedure expressed by Eqns. (5) and (6) is consistent with uplifting the grid system at a height of z0.

and at the phases of 0 and 2π by CG58. The modelgenerated water surface seems to have a reasonable

In addition, because δ is not used in the σ-

agreement with that calculated on the basis of CG58 al-transformation, δ is independent of D and can be set atzero if the simulations do not include the intertidal inun-

though minor noises, which may be caused by reflected

dation and drainage.

shoreline, suggesting that WD-POM is capable of simu-

waves from the moving boundary, are produced near thelating the wetting and drying process in a laboratory

3.

scale.

RESULTS

The temporal variations in the shoreline posi-tions calculated by the model and CG58 are shown in
3.1 Comparison to the analytical solution

Figure 5 (b).

Although the model-generated shoreline

is defined at the center of shore-most wet cells displayed

Whereas WDS proposed here is formulated simplisti-with “plus” marks, manipulation should be done to esti-cally, it guarantees the water mass conservation and canmate the actual shoreline defined at the intersection be--13-

Uchiyama, Y., Dec., 2004: Rep. Port and Airport Res. Inst., Yokosuka, Japan, pp.3-21.
Copyright © 2004 Port and Airport Research Institute. All Rights Reserved.

Napa River
Sacramento
River

Carquinez
StraitPort Chicago
San Joaquin River
Richmond

G.G

San
Francisco

Berkeley
Oakland
AlamedaPacific
Ocean

Redwood City

Figure 7: The bathymetric map of San Francisco Bay, California, USA (left) and the computational grid alignment using the horizontal orthogonal curvilinear transformation (right). The caption “G.G.” shows Golden Gate. The
bathymetry is relative to the MLLW level at San Francisco tide observatory (NOS-NOAA, 2003).
tween the water surface and the beach face as schemati-

fluvial freshwater fluxes nor wind forcing are consid-

cally illustrated in Figure 6. Shoreline positions calcu-

ered. Ten vertical sigma layers, in which the downward

lated with the linear extrapolation are also indicated withgrid spacing decreases exponentially, are assigned whilst

circle marks in Figure 5 (b).

45 × 160 horizontal curvilinear cells are used in thesimulation. The time steps for the external mode and

The model-predicted

shoreline movement reasonably agrees with the theory.

the internal mode are set at 5s and 50s. The three subdepth scales are set as dcr=20cm, dmin= 5cm, and δ = z0 =1cm, respectively.

3.2 Tidal simulation in San Francisco Bay
The model is next applied to simulate tidal currents in
The

The model-generated surface water elevations and 3D

bay has the surface area of about 1240km2 and is encom-

horizontal velocities in the whole computational domain

passed by intertidal mudflats of about 200km2.

South

are decomposed into the 17 dominant harmonic constitu-

Bay, the southernmost part of San Francisco Bay, is of-ents and compared to the observed data as listed in Ta-

ten described as "a tidally oscillating lagoon with den-bles 1 and 2. The model result is verified to success-sity-driven exchanges with the northern reach", andfully reproduce propagation of the tidal waves andNorth Bay including San Pablo Bay and Suisun Bay has

variations in the tidal currents observed in San Francisco

huge influx from the Sacramento-San Joaquin drainage

Bay by NOAA-NOS (2003).

basin (Conomos et al., 1985). The simulation is per-is next conducted as represented in Figures 8 and 9.

formed for a total of 40 days including the first 3 days

The volume influx at the open boundary is almost equal

for the simulation spin-up (Cheng et al., 1993) by solely

to the volumetric change integrated for the whole com-

imposing the astronomical tidal elevations on the open

putational domain, indicating that the volume is basicallyboundary situated off Golden Gate (the bay mouth).

conserved except during the spring tide.

The density is held constant during the simulation and nosome spiky but minor errors are found during this period,

San Francisco Bay, California, U.S.A. (Figure 7).

-14-

Mass conservation check

AlthoughUchiyama, Y., Dec., 2004: Rep. Port and Airport Res. Inst., Yokosuka, Japan, pp.3-21.
Copyright © 2004 Port and Airport Research Institute. All Rights Reserved.

Table 1: Four dominant components in the harmonic constituents of the tidal elevations(OBS: observed, MDL: model-predicted, DIF: difference between OBS and MDL)

M2

S2

K1

O1

Pt. Reyes
San Francisco
AlamedaRedwood City
Richmond
Port Chicago
Pt. Reyes
San Francisco
AlamedaRedwood City
Richmond
Port Chicago
Pt. Reyes
San Francisco
AlamedaRedwood City
Richmond
Port Chicago
Pt. Reyes
San Francisco
AlamedaRedwood City
Richmond
Port Chicago

Amplitude (cm)
OBS
MDL
DIF
54.7
55.4
0.7
58.1
57.8 −0.3
68.6
66.3 −2.3
89.3
82.9 −6.4
61.8
59.0 −2.8
51.4
49.2 −2.2
16.8
16.6 −0.2
17.2
17.2
0.0
19.6
18.8 −0.8
24.8
23.1 −1.7
17.4
16.8 −0.6
14.1
12.3 −1.8
32.6
32.5 −0.1
32.0
32.8
0.8
32.7
33.7
1.0
35.3
35.4
0.1
32.1
32.5
0.4
27.0
28.6
1.6
24.0
24.2
0.2
23.8
24.6
0.8
23.7
24.5
0.8
24.3
24.0 −0.3
23.3
23.9
0.6
18.3
19.3
1.0

OBS
190.5
208.0
222.9
238.0
222.1
283.4
188.4
206.6
223.4
243.0
223.0
283.1
207.7
211.5
217.8
224.5
218.4
255.3
203.9
207.4
212,9
219.3
213.4
245.8

Phase (°)
MDL
191.8
196.0
214.5
233.4
216.2
299.0
188.5
193.5
215.2
239.8
214.9
301.1
208.2
209.4
217.9
228.8
220.0
268.7
203.9
204.4
211.2
219.8
213.1
256.5

DIF
1.3
−12.0
−8.4
−4.6
−5.9
15.6
0.1
−13.1
−8.2
−3.2
−8.1
18.0
0.5
−2.1
0.1
4.3
1.6
13.4
0.0
−3.0
−1.7
0.5
−0.3
11.3

Table 2: Four dominant components in the harmonic constituents of the tidal currents(OBS: observed, MDL: model-predicted, DIF: difference between OBS and MDL)

Locations
Richmond

Oakland

M2
S2
K1
O1
M2
S2
K1
O1

Major Axis (cm/s)
OBS MDL DIF
75.8
68.6 −7.2
21.6
19.3 −2.3
23.3
21.7 −1.6
16.7
14.0 −2.7
52.5
48.6 −3.9
14.6
16.7
1.9
10.4
11.8
1.4
4.5
6.0
1.5

Orientation (°)
OBS MDL DIF
118.8 105.2 −13.6
116.0 105.4 −10.6
121.6 105.5 −16.1
120.2 105.4 −14.8
99.8 103.7
3.9
102.0 103.4
1.4
100.0 105.1
5.1
110.7 106.8−3.9

Phase Epoch (°)
OBS MDL DIF
176.4 170.8 −5.6
172.8 168.1 −4.7
161.8 157.9 −3.9
155.0 148.7 −6.3
326.6 325.6
1.0
318.1 318.2
0.1
298.6 309.7 11.1
277.2 297.2 20.0

it returns to have a value of about zero (i.e., no error isand agree with the results by Cheng et al. (1993); the wa-

measured) immediately after the spring tide.

ter surface rises and descends to inundate and drain out
intertidal areas formed on the fringe of the bay and theInstantaneous distributions of 2DH current vectors,

currents are strong along the deeper channel while weakwater surface elevations, and dry cells at two phases dur-on the shallow areas.

ing the spring tide, at ebbing and high slack phases, aredisplayed in Figure 10. The current fields are realistic-15-

Uchiyama, Y., Dec., 2004: Rep. Port and Airport Res. Inst., Yokosuka, Japan, pp.3-21.
Copyright © 2004 Port and Airport Research Institute. All Rights Reserved.

Figure 8: Volume influx at the open boundary, Vin, and volumetric change, ∆Vc.

Figure 9: Errors in volume flux per unit area of dry cells, (Vin-∆Vc)/Adry: Adry is drained area.

Figure 10: Distributions of the model-generated 2DH current velocities, water elevations, and dry cells at an ebbing
phase (left) and at a high slack phase (right) during the spring tide.

uses a modified bathymetry in which the depths in the4shallow areas less than 2m deep are replace by the depth

DISCUSSION

of 2m (i.e., H ≥ 2m). In other words, the first one takes
account of WDS (Case 1: w/ WDS, hereafter), and the

Effects of intertidal topography on estuarine hydrodynamics are of interest. Two computations are thus car-latter does not consider WDS (Case 2: w/o WDS).

ried out in order to examine a role of shallow fringes.

simplify comparison of tidally-averaged features be-

One computation is done with the actual bathymetry of

tween the two simulations, only a monochromatic sinu-

San Francisco Bay shown in Figure 7, whereas the othersoidal tidal wave at a period of 12 h with amplitude of

-16-

To

Uchiyama, Y., Dec., 2004: Rep. Port and Airport Res. Inst., Yokosuka, Japan, pp.3-21.
Copyright © 2004 Port and Airport Research Institute. All Rights Reserved.

Figure 11: (Left) iso-tidal amplitude and (right) iso-phase charts estimated with the harmonic analysis for the simulated
results w/ WDS (Case 1) and w/o WDS (Case 2).

0.8m is imposed at the open boundary off Golden Gate.

tors in San Pablo Bay calculated with the results of two

This tidal wave approximately coincides with the M2 tidecomputations seem to be similar in almost the whole

which is the predominant tidal component in San Fran-

domain as shown in Figure 12, but slightly different par-cisco Bay (Conomos et al., 1985).

ticularly in the shallow regions at the depth of less than

Simulated water

elevations are decomposed with the harmonic analysis to2m.

estimate co-tidal amplitude and co-phase charts as dis-WDS) show that the principal directions of the tidal el-played in Figure 11.

The both panels in Figure 11 ex-lipses are nearly parallel to the shore and that the tidal

hibit significant difference between two cases can bewaves are progressive because the tidal ellipses are more

found particularly in the southern and northern reaches

circular.

of the bay remote from the bay mouth; the tidal ampli-

waves are traveling and being refracted on the shallow

tudes calculated in Case 1 (w/ WDS) are 5 to 10cm lar-basins to generate the shore-normal components.

ger than those in Case 2 (w/o WDS), while the phases

shapes of the tidal ellipses turn into more linear to dem-

are slightly delayed at 5 to 20 degrees (i.e., 10 to 40

onstrate that oscillatory standing waves are predominant

minutes) in Case 1 behind those in Case 2.

in Case 1. In addition, it is interpreted that the phase of

The tidal currents near the shoreline in Case 2 (w/o

By contrast in Case 1 (w/ WDS), the tidal
The

In San Pablo Bay located in the southwest part of San

the tidal waves in Case 1 significantly retards probably

Francisco Bay, the tidal waves are propagating anti-

because of the bed friction on the shallows, while theclockwise as Kelvin waves attenuated by bottom-

tidal amplitude is amplified owing to the reflected waves

frictional drag as presented by Hendershott and Speranza

on the shallow slopes to form partially standing waves.

(1971). In the present simulations, possible processes

The resultant residual currents are also different in thewould be refraction and shoaling of the tidal waves due

shallow fringes; the residual vectors align normal to theto the sloping bathymetry in the shallow fringes, andisobaths in Case 1 (w/ WDS) whereas parallel to thecentrifugal forces induced by horizontally curved geome-shoreline in Case 2 (w/o WDS). As a result, a circulat-

try of the bay (Lacy and Monismith, 2001). Spatial dis-ing current pattern in San Pablo Bay seems to be moretributions of tidal ellipses and tidal residual current vec-

attenuated in Case 1 (w/ WDS) than that in Case 2 (w/o-17-

Uchiyama, Y., Dec., 2004: Rep. Port and Airport Res. Inst., Yokosuka, Japan, pp.3-21.
Copyright © 2004 Port and Airport Research Institute. All Rights Reserved.

Figure 12: (Left) tidal ellipses and (Right) tidal residual current vectors based on the harmonic analysis with the resultof Case 1 (w/ WDS) and Case 2 (w/o WDS).

WDS). These results suggest that the directions of thenificant role in altering the propagation of the tidal

tidal residual fluxes can also be significantly different inwaves and formation of the tidal residual currents. The

response to the intertidal topographies.

tidal waves propagating as anticlockwise Kelvin waves
are found to generate the shore-normal velocity compo-

5nent, and thus become partially standing waves, during

CONCLUSIONS

flood onto the intertidal sloping bathymetry while the

A wetting and drying scheme for a three-dimensional

shore-parallel component is predominant in the tidal cur-terrain-following coastal ocean model (POM, Blumberg

rents without the wetting and drying. The intertidal ef-

and Mellor, 1983, 1987) is developed with consideration

fects are also observed from the simulations so that theof the mass conservation and computational robustness.

tidal amplitude is also enhanced on the intertidal areas

An extended logarithmic law for the intertidal simula-because of shoaling of the tidal waves.

tions is proposed to precisely reproduce the bed shear
stresses and the velocity profiles within and above the
bed boundary layer.

The phases, incontrast, are visibly retarded as the bed frictional dragacts on the intertidal currents more on the intertidal ar-

The model results are verified byeas.

comparing to the analytical solution by Carrier andGreenspan (1958) theory and the observed data in San(Accepted on August 10, 2004)

Francisco Bay, California, USA, by NOS- NOAA (2003)
to show a reasonable agreement in the variations in water
surface elevations, 3D current velocities, and shoreline

ACKNOWLEDGEMENTS

positions.

This research project was partially supported by Japan

Effects of the intertidal topography on estuarine hy-

Society for Promotion of Science and is a part of the au-

drodynamics are examined by performing two simula-thor’s work during his stay at University of California,

tions with and without the wetting and drying process.

Berkeley, USA, since March 2002 till March 2004.

The intertidal topography is demonstrated to play a sig--18-

Uchiyama, Y., Dec., 2004: Rep. Port and Airport Res. Inst., Yokosuka, Japan, pp.3-21.
Copyright © 2004 Port and Airport Research Institute. All Rights Reserved.

The author would like to express sincere appreciation to

Casulli, V. and Cattani, E., 1994. Stability, accuracy andProfessor Mark T. Stacey of UC Berkeley, and Professor

efficiency of a semi-implicit method for three-

Rodney J. Sobey of Imperial College of London, UK, fordimensional shallow water flow, Computers andtheir helpful comments on the manuscript. Thanks are

Mathematics with Application, 27: 99-112.

also due to Drs. Hirofumi Hinata and Keisuke Nakayama

Cheng, R.T., Casulli, V. and Gartner, J.W., 1993. Tidal,of National Institute for Land and Infrastructure Man-

residual, intertidal mudflat (TRIM) model and its ap-

agement, Japan, for their comments on the numerical

plication to San Francisco Bay, California, Estuarine,

modeling.

Coastal and Shelf Sci., 36: 235-280.
Conomos, T.J., Smith, R.E. and Gartner, J.W., 1985. Environmental setting of San Francisco Bay, Hydrobi-

REFERENCES

ologia, 129: 1-12.
Austria, P.M. and Aldama, A.A., 1990. Adaptive mesh

Davies, A.M., Jones, J.E. and Xing, J., 1997. Review ofscheme for free surface flows with moving bounda-

recent developments in tidal hydrodynamic modeling,ries, In Gambolati, G., Rinaldo, A., Brebbiam, C.A.,

I: Spectral models, J. Hydraulic Eng., 123, 278-292.

Fray, W.G. and Pinder, G.F. (eds.), Computational

Davis, J.R. and Sheng, Y.P., 2003. Development of par-Methods in Surface Hydrology, Springer-Verlag, New

allel storm surge model, Int. J. Numerical Methods inYork, USA, 456-460.

Fluids, 42: 549-580.

Austen, I., Andersen, T.J. and Edelvang, K., 1999. The

Flather, R.A. and Heaps., N.S., 1975. Tidal computations

influence of benthic diatoms and invertebrates on thefor Morecambe Bay, Geophys. J. of the Royal Astro-erodibility of an intertidal mudflat, the Danish Wad-nomical Society, 42, 489-517.

den Sea, Estuarine, Coastal and Shelf Sci., 49: 99 –Galperin, B., Kantha, L.H., Hassid, S. and Rosati, A.,

111.

1988. A quasi-equilibrium turbulent energy model for
geophysical flows, J. Atmos. Sci., 45: 55-62.

Bates, P.D. and Hervouet, J.-M., 1998. A new method

Gross, E.S., Koseff, J.R. and Monismith, S.G., 1999.

for moving-boundary hydrodynamic problems in shal-Three-dimensional salinity simulations of South San

low water, Proc. R. Soc. Lond. A, 455, 3107-3128
Blumberg, A.F. and Mellor, G.L., 1983. Diagnostic andFrancisco Bay, J. Hydraulic Eng., 125, 1199-1209.

prognostic numerical circulation studies of the South

Hamrick, J.M., 1992. A three-dimensional environmental

Atlantic Bight, J. Geophys. Res., 88: 4579-4593.

fluid dynamics computer code: theoritical and compu-

Blumberg, A.F. and Mellor, G.L., 1987. A description oftational aspects, Special Report 317, The College ofa three-dimensional coastal ocean circulation model,

William and Mary, Virginia Institute of Marine Sci-

In Three-dimensional Coastal Ocean Models, Vol.4,

ence, Williamsburg, VA, 63p.
Haidvogel, D. B., Wilkin, J.L. and Young, R., 1991. A

Heaps, N. (ed.), Ameri. Geophys. Union, Washingtonsemi-spectral primitive equation ocean circulation

D.C., pp.208.

model using vertical sigma and orthogonal curvilinear

Blumberg, A.F., Signell, R.P. and Jenter, H.L., 1993.

horizontal coordinates, J. Comp. Phys., 94: 151-184.

Modeling transport processes in the coastal ocean, J.

Hendershott, M.C. and Speranza, A., 1971. Co-

Mar. Envir. Eng., 1: 31-52.

oscillating tides in long, narrow bays; the Taylor prob-

Carrier, G.F. and Greenspan, H.P., 1958. Water waves oflem revisited, Deep-Sea Res., Vol.18, pp.959-980.

finite amplitude of a sloping beach, J. Fluid Mech., 4:

Hervouet, J.M. and Janin, J.M., 1994. Finite-element al-

97-109.
Casulli, V. and Cheng, R.T., 1991. A semi-implicit fi-

gorithms for modeling flood propagation, In Molinaro,

nite-difference model for three-dimensional tidal cir-P. and Ntale, L. (eds.), Modeling of Flood Propaga-

culation, In Spaulding, M., Swanson, C., Cheng, R.,tion Over Initially Dry Areas, Ameri. Soc. Civil. Eng.,

Blumberg, A. and Bedford, K. (eds), Proc. 2nd Intl.

New York, NY, 102-113.
Hervouet, J.-M. & Van Haren, L. 1996. Recent advances

Conf. Estuarine and Coastal Modeling, Ameri. Soc.

in numerical methods for fluid flows, In: Anderson,

Civil. Eng., Tampa, FL, USA, 620-631.

-19-

Uchiyama, Y., Dec., 2004: Rep. Port and Airport Res. Inst., Yokosuka, Japan, pp.3-21.
Copyright © 2004 Port and Airport Research Institute. All Rights Reserved.

M.G., Walling, D.E. and Bates, P.D. (eds), Floodplain

CIRC hydrodynamic model: implementation andprocesses, Wiley, Chichester, UK, 183-214.

performance in one-dimensional flows, Report for the

Hubbert, G.D. and McInnes, K.L., 1999. A storm surge

Department of the Army, DACW39-94-M-5869.

model for coastal planning and impact studies, J.

Luettich, R.A. and Westerink, J.J., 1995b. Implementa-

Coastal Res., 15 (1), 168-185.

tion and testing of elemental flooding and drying in

Ip, J.T.C., Lynch, D.R. and Friedrichs, C.T., 1998. Simu-

the ADCIRC hydrodynamic model, Final Report forlation of estuarine flooding and dewatering with appli-

the Department of the Army, DACW39-94-M-5869.

cation to Great Bay, New Hampshire, Estuarine,

Luettich, R.A. and Westerink, J.J., 1995c. Continental

Coastal and Shelf Sci., 47, 119-141.

shelf scale convergence studies with a barotropic tidal

Ji, Z-G., Morton, M.R. and Hamrick J.M., 2001. Wetting

model, Coastal Estuarine Stud., 47: 349-371.

and drying simulation of estuarine processes, Estua-

Lynch, D.R. and Gray, W.G., 1980. Finite element simu-

rine, Coastal and Shelf Sci., 53: 683-700.

lation of flow in deforming regions, J. Computational

Johns, B., Dube, S.K., Sinha, P.C., Mohanty, U.C. andPhysics, 36, 135-153.

Rao, A.D, 1982. The simulation of a continuously de-

Mellor, G.L., 2002. Oscillatory Boundary Layers, J.

forming lateral boundary in problems involving thePhys. Oceanogr., 32: 3075-3088.

shallow water equations, Computer and Fluids, 10,

Mellor, G.L. and Yamada, T., 1982. Development of a105-116.

turbulence closure model for geophysical fluid problems, Rev. Geophys. Space Phys., 20: 851-875.

Kuwae, T., Hosokawa, Y., and Eguchi, N., 1998. Dissolved inorganic nitrogen cycling in Banzu intertidal

National Ocean Service, National Oceanic and Atmos-

sand-flat, Japan, Mangroves and Salt Marshes, 2: 167-

pheric Administration (NOS-NOAA), 2003. http://

175.

www. co-ops.nos.noaa.gov/.
Oey, L-Y., 2004. A wetting and drying scheme for POM,

Lacy, J. and Monismith, S.G., 2001. Secondary currents

Ocean Modeling, 9. (in press)

in a curved, stratified channel, J. Geophys. Res., 106

Park, K., Oh, J.-H., Kim, H.-S. and Im, H.-H., 2002.

(C12): 31,283-31,302.
Leclerc, M., Bellemare, J-F., Dumas, G. and Dhatt, G.,

Case study: Mass transport mechanism in Kyunggi

1990. A finite element model of estuarian and river

Bay around Han River mouth, Korea, J. Hydr. Eng.,
128: 257-267.

flows with moving boundaries, Advances in Water Re-

Shchepetkin, A.F., McWilliams, J.C., 1998. Quasi-

sources, 13: 158-168.
Leendertse, J.J., 1970. A water-quality simulation model

monotone advection schemes based on explicit locallyfor well mixed estuaries and coastal seas: Principles ofadaptive dissipation. Mon. Wea. Rev. 126, 1541–1580.

computation, Report RM-6230-RC, Vol. I, Rand Cor-

Shchepetkin, A.F., McWilliams, J.C., 2003. The regional
ocean modeling system: a split-explicit, free-surface,poration, Santa Monica, CA.

topography-following coordinate ocean model. J.

Leendertse, J.J., 1987. Aspects of SIMSYS2D, a system

Geophys. Res., 108: 3501-3534.

for two-dimensional flow computation, Report R3572-USGS, Rand Corporation, Santa Monica, CA.

Shi, F., 1995. On moving boundary numerical models of

Le Hir, P., Roberts, W., Cazaillet, O., Christie, M., Bas-coastal sea dynamics, Ph.D. Dissertation, Ocean Univ.
of Qingdao, Qingdao, China.

soullet, P. and Bacher, C., 2000. Characterization ofSiden, G.L.D. and Lynch, D.R., 1988. Wave equationintertidal flat hydrodynamics, Continental Shelf Res.,

hydrodynamics o n deforming elements, Intl. J. Nu-

20: 1433 – 1459.

merical Methods in Fluids, 8, 1071-1093.

Lin, B. and Falconer, R.A, 1997. Three-dimensional
layer-integrated modeling of estuarine flows with

Talke, S.A. and Stacey, M.T., 2003. The influence offlooding and drying, Estuarine, Coastal and Shelf Sci.,

oceanic swell on flows over an estuarine intertidal

44, 737-751.

mudflat in San Francisco Bay, Estuarine, Coastal and
Shelf Sci., 58: 541-554.

Luettich, R.A. and Westerink, J.J, 1995a. An assessmentof flooding and drying techniques for use in the AD--20-

Uchiyama, Y., Dec., 2004: Rep. Port and Airport Res. Inst., Yokosuka, Japan, pp.3-21.
Copyright © 2004 Port and Airport Research Institute. All Rights Reserved.

Uchiyama, Y., Kuriyama, Y. and Katoh, K., 2001. Sus-

Xie, L., Pietrafesa, L.J. and Peng, M., 2003. Incorpora-pended sediment and morphological response ontion of a mass-conserving inundation scheme into aBanzu tidal flat, Japan, Proc.4 Int’l Conf. Coastal

three dimensional storm surge model, J. Coastal Res.,th

Dynamics, Ameri. Soc. Civil Eng., Lund, Sweden:

1-17.

1038-1047.

Zheng, L., Chen, C. and Liu, H., 2003. A modeling study

Van Der Lee, W.T.B., 2000. Temporal variation of flocof the Satilla River Estuary, Georgia. I: Flooding-size and settling velocity in the Dollard estuary, Con-drying process and water exchange over the salttinental Shelf Res., 20: 1495 – 1511.

marsh-estuary-shelf complex, Estuaries, 26 (3): 651-

Whitehouse, R.J.S., Soulsby, R.L., Roberts, W. and669.

Mitchener,H.J., 2000. Intertidal Processes, In: Dynamics of estuarine muds, Thomas Telford, London, UK:
163-168.

-21-