Modelling the impact of wind stress and river discharge on Danshuei River plume

(1)

Modelling the impact of wind stress and river discharge

on Danshuei River plume

Wen-Cheng Liu

a,*

, Wei-Bo Chen

b

, Ralph T. Cheng

c

, Ming-Hsi Hsu

d

a

Department of Civil and Disaster Prevention Engineering, National United University, Miao-Li 36003, Taiwan

b

National Center for Ocean Research, National Taiwan University, Taipei 10617, Taiwan

c

US Geological Survey, Menlo Park, CA 94025, USA

d

Department of Bioenvironmental Systems Engineering, National Taiwan University, Taipei 10617, Taiwan Received 1 July 2006; received in revised form 1 March 2007; accepted 26 March 2007

Available online 22 April 2007

Abstract

A three-dimensional, time-dependent, baroclinic, hydrodynamic and salinity model, UnTRIM, was performed and applied to the Danshuei River estuarine system and adjacent coastal sea in northern Taiwan. The model forcing functions consist of tidal elevations along the open boundaries and freshwater inflows from the main stream and major tributaries in the Danshuei River estuarine system. The bottom friction coefficient was adjusted to achieve model calibration and ver-ification in model simulations of barotropic and baroclinic flows. The turbulent diffusivities were ascertained through com-parison of simulated salinity time series with observations. The model simulation results are in qualitative agreement with the available field data.

The validated model was then used to investigate the influence of wind stress and freshwater discharge on Dasnhuei River plume. As the absence of wind stress, the anticyclonic circulation is prevailed along the north to west coast. The model results reveal when winds are downwelling-favorable, the surface low-salinity waters are flushed out and move to southwest coast. Conversely, large amounts of low-salinity water flushed out the Danshuei River mouth during upwell-ing-favorable winds, as the buoyancy-driven circulation is reversed. Wind stress and freshwater discharge are shown to control the plume structure.

Keywords: River plume; Three-dimensional model; Danshuei River estuarine system and coastal sea; Freshwater discharge; Wind stress; Taiwan

1. Introduction

The costal ocean is the recipient of freshwater and land drained materials that are primarily brought in through river discharge. From the dynamical viewpoint the discharge site (the ‘‘river mouth’’) can be

doi:10.1016/j.apm.2007.03.009

* _{Corresponding author. Tel.: +886 37 381674; fax: +886 37 326567.}

E-mail address:wcliu@nuu.edu.tw(W.-C. Liu).

Applied Mathematical Modelling 32 (2008) 1255–1280

(2)

considered the source both momentum and buoyancy, produced by the release of light ﬂuid into a heavier (more dense) ambient. The resulting dynamical structure is referred to as a (buoyant) river plume[1].

River waters have been regards as means for drinking water supply, maintenance of wild life and transpor-tation, while their natural content in perceptible amounts of nutrients has always influenced the ecology of the receiving waters, which directly affects the food supply. In modern times, population growths of coastal states and technology-related development have added a plethora of man-made waste products on the river intake. These ‘‘materials’’ are carried by rivers (in solution, in suspension or as bed load) to their discharge site. It is becoming increasingly necessary to be able to control the impact of this type of ‘imported’’ pollution in the coastal areas, which continue to be favored as rapidly growing population and recreation sites. Furthermore, the freshening of nearshore waters due to river runoff is one of the mechanisms that control circulation in coastal areas[2]. The present study is motivated by the need to understand and predict the phenomena that control the fate of riverine waters and related materials, taking into account the complex geomorphology that characterizes the Danshuei River system and its coastal sea.

Numerical studies of river plume have mostly concentrated on idealized model domains, so that results could be validated analytically and nondimensional numbers characteristic of plume dynamics could be derived. Such studies concentrate on plume dynamics[1,3–7]. These studies pointed out that, when the buoy-ant discharge is the only forcing mechanism, a ‘‘typical’’ mid-latitude plume develops after (a) the initial rig-orous mixing caused by turbulence and friction and (b) the subsequent balance between Coriolis and the density-induced pressure gradient. The structure of a river plume includes an offshore bulge that extends anti-cyclonically seaward from the river mouth and the coastal current region that develops along a narrow near-shore strip to the right of the river mouth. Wind stress may either enhance or oppose the density-driven coastal current. Kourafalou et al.[1]showed that downwelling favorable winds strengthen the coastal current plume region and confine the bulge nearshore, while upwelling favorable winds may reverse the coastal current and allow a wide offshore bulge. Kourafalou et al.[8]studied, under realistic conditions, the wind and tidal influ-ences on a river plume along the Southeast US continental shelf during spring season. Mikhailova and Shapiro

[9]included the interaction between riverine low-salinity waters and abyssal upwelled waters in the Black Sea. Garvine [10] investigated the dependence of the alongshelf penetration of an unforced buoyant coastal discharge in parameters such as bottom slope, background diffusivity, tidal amplitude, and river discharge. Berdeal et al.[11]used a three-dimensional model (ECOM3d) to simulate the response of a high discharge river plume to an alongshore ambient flow and wind forcing. Arnoux-Chiavassa [12] used a higher order advection scheme to simulate 3D Rhone river plume. Lacroix et al.[13]described the distribution and vari-ability of the salinity in Belgian coastal waters and determined the relative impact of the Scheldt and Rhine/Meuse freshwater plume. Whitney and Garvine [14] simulated the Delaware Bay buoyant outflow and compared results with observations of estuarine and shelf conditions.

The numerical simulation in realistic cases contributed to the understanding of the dynamics and the struc-ture of the plumes by taking into account real bathymetries. A number of numerical models have applied to investigate the river plume in diﬀerent continental shelves and estuaries overall the world. However there is not any report regarding to the Danshuei River plume in Taiwan. A number of three-dimensional hydrodynamic models are available to simulate the tidal current, salinity, and temperature in estuaries and coastal seas

[15–21]. In the present study, a three-dimensional hydrodynamic model, UnTRIM[18,22], was implemented and applied to examine the impact of wind stress and freshwater discharge on river plume development in the Danshuei River coastal sea of northern Taiwan. Before the performance of model application, model valida-tion was conducted with water surface elevavalida-tion, current, and salinity and compared them with available ﬁeld data in the Danshuei River estuarine system. Because of the unstructured gird applied in the model, the grids generation is easily to produce and to ﬁt the topography. Moreover, the UnTRIM model takes the advantage to develop the water quality and sediment models in the future work.

2. Study site

The Danshuei River estuary (Fig. 1a) is the largest estuarine system in Taiwan, with its drainage basin including the capital city of Taipei. The tidal inﬂuence spans a total length of about 82 km, encompassing the entire length of the Danshuei River and the downstream reaches of its three major tributaries: the Tahan

(3)

Stream, the Hsintien Stream, and the Keelung River. The average river discharges at the upstream limits of tide are 41.7 m3/s, 72.4 m3/s, and 25.5 m3/s, respectively, in the Tahan Stream, Hsintien Stream, and Keelung River. The northeast wind prevails during the autumn and winter seasons, while the southwest wind prevails

Fig. 1a. The map of the Danshuei River estuarine system and its adjacent coastal sea.

(4)

during the spring and summer seasons in Taiwan. The moderate to strong wind speed ranges from 6 to 15 m/s

[23]. Except during a flooding event, the astronomical tide may reach as far upriver as Cheng-Ling Bridge in Tahan Stream, the Hsiu-Lang Bridge in Hsintien Stream and the Chiang-Pei Bridge in Keelung River (Fig. 1a). Tidal propagation into the estuary from coastal ocean is the dominant mechanism controlling the variations of water surface elevation, and ebb and flood flows. The M2tide is the primary tidal constituent

at the river mouth, with a mean tidal range of 2.17 m, and up to 3 m at spring tides. Because of the cross-sec-tional contraction and wave reflection, the mean tidal range may reach a maximum of 2.39 m within the sys-tem. The phase relationship between tidal elevation and tidal flow is close to standing wave characteristics[24]. Sea water intrudes upriver as a result of tidal dispersion and residual circulation. Salinity varies in intra-tidal time scale in response to the ebb and flood of intra-tidal flows as well as in response to changing freshwater inflows. The limit of salt intrusion may reach beyond 25 km in the Tahan Stream from river mouth during the period of low flow. The baroclinic pressure gradient due to longitudinal salinity/density distribution is large enough to push the denser salt water upriver along the bottom layer of the estuary forming the classical two-layer estuarine circulation[25–27].

In this study, a baroclinic 3 D model, UnTRIM, is used to cover the Danshuei River estuarine system and adjacent coastal sea where the bathymetry is gradually getting deeper away from the coastline (Fig. 1a), and the model uses an unstructured model grid as shown inFig. 1b. Wang[28]used a towed-ADCP to measure the tidal current at the coastal sea off the Danshuei River. The maximum tidal current near the river mouth can reach 1.25 m/s and the tidal current fluctuates with flood and ebb tides along the coastline. The measured mean flow at the upper 10 m is about 0.5 m/s. The flow decreases with increasing depth; the mean current is about 0.35 m/s at 10–40 m and 0.2 m/s below 40 m of depth.

3. Model description 3.1. Governing equation

The governing equations describing three-dimensional free-surface ﬂows are the well known shallow water equations. Such equations express the physical conservation principles of volume, mass, and momentum.

The volume conservation is expressed by the following incompressibility condition: ou oxþ ov oyþ ow oz ¼ 0; ð1Þ

where u(x, y, z, t), v(x, y, z, t), and w(x, y, z, t) are the velocity components in the horizontal x, y, and vertical z-directions, respectively,

Integrating the continuity equation(1)over depth and using a kinematic condition at the free-surface leads to the following free-surface equation[16]:

og otþ o ox Z g h udz þ o oy Z g h vdz ¼ 0; ð2Þ

where t is the time; h(x, y) is the prescribed bathymetry measured from the undisturbed water surface and g(x, y, t) is the free-surface elevation. Thus, H(x, y, t) = h(x, y) + g(x, y, t) is the total water depth.

The mass conservation of salt (solutes) is expressed by the following diﬀerential equation: oS ot þ u oS oxþ v oS oyþ w oS oz ¼ o ox K hoS ox þ o oy K hoS oy þ o oz K voS oz ; ð3Þ

where S is the salinity; Khand Kvare prescribed non-negative horizontal and vertical diﬀusivities, respectively. Under the hydrostatic approximation, the vertical accelerations in the vertical momentum equation can be neglected.

The hydrostatic pressure p(x, y, z, t) can be expressed as pðx; y; z; tÞ ¼ paðx; y; tÞ þ g½gðx; y; tÞ z þ g Z g z q q0 q0 df; ð4Þ

(5)

where pais the atmospheric pressure; g is the gravitational acceleration; q is the water density; qois the

fresh-water density.

The horizontal momentum equations can also be written as Du Dt fv ¼ opa ox g og ox g o ox Z g z q q0 q d1 þ vh o 2_u ox2þ o2u oy2 þ o oz v vou oz ; ð5Þ Dv Dtþ fu ¼ opa oy g og oy g o oy Z g z q q0 q d1 þ vh o 2_v ox2þ o2v oy2 þ o oz v vov oz ; ð6Þ

where D( )/Dt is the substantive derivative; f is the Coriolis parameter; vhand vvare, respectively, the coef-ﬁcients of horizontal and vertical eddy viscosity.

The vertical momentum is expressed as Dw Dt ¼ v h o 2 w ox2þ o2w oy2 þ o oz v vow oz : ð7Þ

The system is closed by an equation of state which relates the water density to the concentration of salinity (contributions from temperature variations to density are neglected). The equation of state takes the form of

q¼ q0ð1 þ kSÞ; ð8Þ

where k is constant (=7.8· 104ppt1). 3.2. Turbulence closure model

In most estuaries with a significant freshwater discharge, salinity may serve as an idea natural tracer for calibration of mixing processes. Salinity distribution in an estuary is affected by the tidal current, freshwater discharge, density circulation, as well as turbulent mixing processes. Therefore the salinity distribution reflects the combined results of all processes, and in turn it controls density circulation and modifies mixing processes

[24].

In the present model, the mixing processes are modeled with turbulent diffusion terms. The mixing length concept is used to calculate eddy viscosity vvand diffusion coefficient Kvin vertical direction.

The formations for vvand Kvare vv¼ aZ2 ₁Z h 2 ou oz ð1 þ bRiÞ 1=2 ; ð9Þ Kv¼ aZ2 ₁Z h 2 ou oz ð1 þ bRiÞ3=2; ð10Þ

where Z is the distance from the surface, a is a constant to be determined empirically, and a local Richardson number (Ri) is used to characterize stability due to stratiﬁcation, where the Richardson number is deﬁned as

Ri¼ g q oq oz _ou oz 2 : ð11Þ

The horizontal mixing coeﬃcients (vh and Kh) range from 1 to 100 m2/s[29]. Since the horizontal length scales are several orders of magnitude larger than the vertical length scales, the horizontal mixing terms play a relatively insigniﬁcant role in momentum balance; they are retained in the model nevertheless. Constant val-ues for vhand Khare used and they are adjusted, within the range of 1–100 m2/s, through model calibration. 3.3. Numerical approximation

The numerical algorithm of UnTRIM is fundamentally the same as TRIM3D[16,30], except that the finite-difference treatment of the governing partial differential equations is performed over an unstructured grid mesh. Before discretizing the governing equations, the horizontal domain (x, y) is covered by a set of

(6)

non-overlapping convex polygons. Each side of a polygon is either a boundary line or a side an adjacent polygon. Moreover, it is assumed that within each polygon there exists a point (i.e., center) such that the segment joining the centers of two adjacent polygons and the side shared by the two polygons have a non-empty intersection and are orthogonal to each other (shown inFig. 2). One such grid is called an unstructured orthogonal grid[17,31]. The center of a polygon does not necessarily coincide with its geometrical center. The special cases of unstructured orthogonal grids include the rectangular finite-difference grids, as well as a grid of uniform equilateral triangles. In these particular cases the center of each polygon can be identified with its geometrical center. Another example of an unstructured orthogonal grid is a set of Delaunay triangles where the triangulation includes only acute triangles[32].

A semi-implicit scheme define above is used in order to obtain an efficient numerical algorithm whose sta-bility is independent from the free-surface gravity waves, wind stress, vertical viscosity and bottom friction. Consider a typical polygon,Fig. 2, the momentum equation, Eqs.(5) and (6), is finite-differenced in the nor-mal direction of each vertical face along oa, ob, and oc directions. The momentum equation relates the gra-dient of water surface elevation between adjoining polygons and the face velocity on the common face between these polygons. As stated previously, the wind stress, the vertical mixing and the bottom friction are discret-ized implicitly for numerical stability.

An explicit ﬁnite-diﬀerence operator is used to account for the contributions from the discretization of the advection and horizontal dispersion terms. A particular form for this operator can be given in several ways, such as an Eulerian–Lagrangian scheme[16]. For stability, the implicitness factor h has to be chosen in the range1

26h 61[30]. Along the vertical direction, the model domain is discretized by ﬁxed levels, not

neces-sarily uniform, with which a simple finite-difference discretization is used. The vertical space increment is defined as the distance between two consecutive level surfaces. In general, the vertical thickness of the top and bottom layers can vary depending on the spatial location and the thickness of the top layer can also vary with time. The vertical space increment is allowed to vanish. In fact, this is how the wetting and drying of com-putational cells are accomplished.

For three-dimensional barotropic flow, the salt transport equation is un-coupled from the momentum equations. A non-negative bottom friction coefficient is specified by the Manning–Chezy formula [16]. If the advections are treated by an Eulerian–Lagrangian scheme with positive bottom friction, then the numer-ical scheme is unconditionally stable. For baroclinic flows, since the density gradient terms are expressed explicitly in the momentum equations, and the solutions of the transport variables are solved lagged one time-step. In this case, the numerical scheme is subject to a weak stability condition due to the explicit treat-ment of the density gradient terms in the motreat-mentum equations. For stability, the integration time step must be chosen so that the propagation of internal wave must satisfy the Courant–Friedrich–Lewy (CFL) condition. The necessary and sufficient stability condition is formed as[16]

Dt 6 juj Dxþ jvj Dyþ jwj Dzþ 2v h 1 Dx2þ 1 Dy2 1 : ð12Þ

In addition, the numerical scheme is also subject to a weak stability condition due to the explicit treatment of the horizontal diﬀusion in the momentum equations.

3.4. Grid generation

Many factors aﬀect ﬂows in the estuarine system and its adjacent coastal sea, but the bathymetry is the most important factor. An accurate bathymetric representation by the model grid is the most important and fun-damental requirement in a successful modelling study[33,34]. This is particularly true for the Danshuei River estuarine system and adjacent coastal sea, where the bathymetric variations are very complex. The model grid must represent accurately the characteristics of the Danshuei River estuarine system that constitutes the model domain.

The UnTRIM model uses an unstructured grid, some aspects of the grid definition file are similar to those used in finite-element applications, therefore the literatures on grid generation for finite-elements apply in

(7)

UnTRIM applications. For the purpose of generating an unstructured model grid for UnTRIM model appli-cation, a commercial software for mesh generation, Argus software[35], has been adopted.

Fig. 1c. The contour of Danshuei River estuarine system and adjacent coastal sea.

(8)

The original bathymetry data used in this study were obtained from National Center for Ocean Research and Water Resources Agency, Taiwan. The deepest depth within the study area is 110 m (below the mean sea level) near the northeast corner of the model in the coastal sea. The model mesh for the Danshuei River estu-arine system and its adjacent coastal sea consists of 8181 3- and 4-sided polygons (Fig. 1b). Higher resolution grids are used in the Danshuei River estuary, and coarser grids are used in the coastal sea. According to the bathymetry and topography in the Danshuei River estuarine system and adjacent coastal sea (Figs. 1c and 1d), sixty vertical layers are speciﬁed with the layer thickness varying from 1 m in the top 10 layers, 2 m for the other layers. For this model grid, 360 seconds time step was used in simulations without any sign of numerical stability.

4. Model validation

4.1. Calibration with amplitudes and phases

Hsu et al.[24]used a vertical two-dimensional hydrodynamic model to the Danshuei River estuarine sys-tem. They found that the bottom roughness is a significant parameter that affects the tidal amplitudes and phases along the Danshuei River estuarine system. To focus on a particular aspect of the dynamic forcing and simplify the calibration process, it is proposed that the model be tuned with astronomical tide only. The approach also eliminates the uncertainty caused by measured errors in the field data. The satisfactory comparison with field data will serve as verification that the model properly simulates all of the forcing combined.

Jan et al.[36]reported that ﬁve tidal-constituents (i.e., M2, S2, N2, K1, and O1) are suﬃcient to represent the

tidal dynamics in the Taiwan Strait. Therefore the ﬁve tidal-constituents according to Jan et al. [36] were adopted in the model simulation as a forcing function at the coastal sea boundaries. The amplitudes and phases used for the model simulation in the coastal sea are listed inTable 1. Amplitudes and phases of these ﬁve tidal constituents were used to generate time-series water surface elevation along the open boundaries for two-month model simulation. Harmonic analysis was performed on the time series of the model simulated

Table 1

The amplitudes and phases used for the model simulation at the coastal sea boundaries

Constituent East boundary Northeast boundary West boundary Northwest boundary

Amplitude (m) Phase () Amplitude (m) Phase () Amplitude (m) Phase () Amplitude (m) Phase ()

M2 0.47 172.37 0.6711 171.33 0.1236 181.63 1.3716 180.42

S2 0.12 331.75 0.1866 332.25 0.3452 354.5 0.3987 355.2

N2 0.10 262.74 0.1370 262.35 0.2252 281.373 0.2526 279.41

K1 0.21 232.33 0.2223 234.2 0.2121 246.06 0.2275 246.69

O1 0.17 67.54 0.170 67.05 0.1748 75.02 0.1866 74.40

(9)

water surface elevation at various locations. Through the bottom friction coeﬃcient represented by the Manning–Chezy formula in the model was adjusted carefully, the results are presented in Fig. 3. It shows

M2 S2 N2 K1 O1 0 0.2 0.4 0.6 0.8 1 1.2 1.4 Amplitude (m)

Danshuei Coastal Sea-Pile

M2 S2 N2 K1 O1 0 50 100 150 200 250 300 350 400 Phase ( d egr e e )

Danshuei Coastal Sea-Pile

M2 S2 N2 K1 O1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 A m p li tude ( m ) Fu-Gi Harbor M2 S2 N2 K1 O1 0 50 100 150 200 250 300 350 400 P h ase ( d egr e e ) Fu-Gi Harbor M2 S2 N2 K1 O1 0 0.2 0.4 0.6 0.8 1 1.2 A m p li tude ( m )

Danshuei River Mouth

M2 S2 N2 K1 O1 0 50 100 150 200 250 300 350 400 Phase ( d egree)

Danshuei River Mouth

Fig. 3. Comparisons of amplitude and phase of ﬁve major tidal harmonics computed with the 3D model and obtained from tide measurements at nine stations: (j) observation, ( ) simulation.

(10)

the comparison of the amplitudes and phases of harmonic constants between the computed and observed tides at nine locations.Table 2presents that the calculated absolute mean diﬀerence and root-mean-square errors of

M2 S2 N2 K1 O1 0 0.2 0.4 0.6 0.8 1 1.2 Amplitude (m) Tu-Ti-Kung-Pi M2 S2 N2 K1 O1 0 50 100 150 200 250 300 350 Phase (degree) Tu-Ti-Kung-Pi M2 S2 N2 K1 O1 0 0.2 0.4 0.6 0.8 1 1.2 Amplitude (m) Taipei Bridge M2 S2 N2 K1 O1 0 50 100 150 200 250 300 350 Phase (degree) Taipei Bridge M2 S2 N2 K1 O1 0 0.2 0.4 0.6 0.8 1 1.2 Amplitude (m) Ru-Kou Weir M2 S2 N2 K1 O1 0 50 100 150 200 250 300 350 Phase (degree) Ru-Kou Weir Fig. 3 (continued)

(11)

tidal constituents for amplitudes are in the ranges of 0.01–0.02 m and for phases are in the ranges of 2.4–6.7. The diﬀerences in amplitudes and phases are quite small, and the water surface elevations are reasonably validated. M2 S2 N2 K1 O1 0 0.2 0.4 0.6 0.8 1 1.2 Amplitude (m) Hsin-Hai Bridge M2 S2 N2 K1 O1 0 50 100 150 200 250 300 350 P h ase ( d egr e e ) Hsin-Hai Bridge M2 S2 N2 K1 O1 0 0.2 0.4 0.6 0.8 1 1.2 A m p li tude ( m ) Chung-Cheng Bridge M2 S2 N2 K1 O1 0 50 100 150 200 250 300 350 Phase ( d egree) Chung-Cheng Bridge M2 S2 N2 K1 O1 0 0.2 0.4 0.6 0.8 1 1.2 A m p li tude ( m ) Ta-Chih Bridge M2 S2 N2 K1 O1 0 50 100 150 200 250 300 350 Phase (degree) Ta-Chih Bridge Fig. 3 (continued)

(12)

4.2. Veriﬁcation of water surface elevation, current, and salinity

The model verification was conducted with the daily freshwater discharges of the year 2000 as the forcing at the upstream boundaries. The open boundaries at the coastal sea were specified by five tidal-constituent (given inTable 1). The model was run for one-year simulation. The model results of time-series surface ele-vation, current and salinity were verified with field data in the same periods when field data are available. Because of the typhoon events and storm run-off, the peak freshwater discharges occurred on August 23 and 30. The freshwater discharge inputs from three tributaries (Tahan Stream, Hsintein Stream, and Keelung River) during the period of August 21–31, 2000 are shown inFig. 4a, and during these eleven-days the computed surface elevations are compared with field data at several stations (Fig. 4c–i). The wind con-ditions were included in the model simulation (Fig. 4b). In general, the modelling results reproduce the water level variations very well. The water surface elevations at the upriver stations (Fig. 4g–i) have much more conspicuous response to pulses of high freshwater discharge, the modeled water levels follow the low water variations closely throughout the storm event. It demonstrates that the model can accurately repro-duce water surface elevations even under extremely large variations of daily freshwater discharge input from three tributaries.

The Taiwan Water Resources Agency conducted 13 h (during daylight hours) intensified field measurement at five transects on May 5, 2000. Water velocities were measured half-hourly at several stations. Water speed was measured with handheld current meters by personnel on boats, the water directions were not measured.

Table 2

The comparison of computed and observed amplitudes and phases of tidal constituents

Constituent Danshuei Coastal

Sea-Pile

Fu-Gi Harbor Danshui River Mouth Tu-Ti-Kung-Pi Taipei Bridge

Observed Simulated Observed Simulated Observed Simulated Observed Simulated Observed Simulated

(a) Amplitude (m) M2 1.07 1.05 0.69 0.67 1.06 1.07 1.05 1.08 1.03 1.06 S2 0.31 0.30 0.18 0.17 0.29 0.28 0.26 0.27 0.24 0.25 N2 0.22 0.21 0.15 0.14 0.20 0.21 0.18 0.19 0.20 0.20 K1 0.22 0.20 0.23 0.23 0.20 0.20 0.17 0.18 0.14 0.15 O1 0.17 0.17 0.18 0.17 0.17 0.17 0.13 0.15 0.11 0.11 (b) Phase () M2 176.18 174.19 170.18 169.17 179.34 176.96 191.38 196.08 193.80 198.93 S2 351.8 351.42 342.96 340.16 351.90 352.21 4.88 4.76 11.29 9.81 N2 284.74 287.34 257.24 263.74 278.21 284.53 292.02 303.23 296.98 303.38 K1 235.00 239.93 241.00 234.97 242.30 241.87 246.35 250.50 248.86 254.58 O1 68.47 67.78 68.25 67.44 72.60 68.67 76.35 73.93 68.47 61.82

Ru-Kou Weir Hsin-Hai Bridge Chung-Cheng Bridge Ta-Chih Bridge AME RMSE

Observed Simulated Observed Simulated Observed Simulated Observed Simulated

(a) Amplitude (m) M2 1.04 1.04 1.03 1.03 1.06 1.05 0.96 0.98 0.02 0.02 S2 0.26 0.24 0.26 0.23 0.28 0.25 0.22 0.21 0.02 0.02 N2 0.18 0.19 0.19 0.19 0.19 0.19 0.14 0.16 0.01 0.01 K1 0.15 0.14 0.15 0.14 0.16 0.15 0.12 0.12 0.01 0.01 O1 0.11 0.11 0.12 0.11 0.13 0.15 0.11 0.11 0.01 0.01 (b) Phase () M2 210.68 216.10 205.64 211.55 204.82 207.55 216.05 219.00 3.58 3.94 S2 31.32 26.85 27.10 25.81 24.65 19.60 39.43 33.82 2.39 3.14 N2 311.91 319.02 306.34 313.16 306.26 311.42 318.38 321.99 6.27 6.69 K1 255.28 260.19 253.11 259.14 253.44 257.90 266.81 267.66 4.17 4.62 O1 84.08 82.65 80.50 73.83 84.38 77.47 85.83 78.65 4.08 4.86

(13)

As the velocity data were recorded by hand and ‘‘ebb or ﬂood’’ directions were also noted based on visual observations. During slack tides, there were some diﬃculties (uncertainty) assigning current directions. Except

08/21 08/22 08/23 08/24 08/25 08/26 08/27 08/28 08/29 08/30 08/31 Time (day) 0 200 400 600 800 1000 1200 1400 Dischar ge (m 3/s ) Tahan Stream Hsintien Stream Keelung River 0 24 48 72 96 120 144 168 192 216 240 Time (hr) -2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5 Elevation (m) Simulation Measurement Tu-Ti-Kung-Pi 8/21 8/31 0 24 48 72 96 120 144 168 192 216 240 Time (hr) -12 -9 -6 -3 0 3 6 9 12 Wi n d vec to r ( m / s ) 8/21 8/31 0 24 48 72 96 120 144 168 192 216 240 Time (hr) -2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5 Elevation (m) Simulation Measurement Taipei Bridge 8/21 8/31 0 24 48 72 96 120 144 168 192 216 240 Time (hr) -2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5 Elevation (m) Simulation Measurement River Mouth 8/21 8/31 0 24 48 72 96 120 144 168 192 216 240 Time (hr) -2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5 Elevation (m) Simulation Measurement Ru-Kou Weir 8/21 8/31

Fig. 4. (a) Freshwater discharge inputs from three tributaries during the period of August 21–31, 2000, (b) wind conditions, and the comparison of water surface elevation at: (c) Danshuei River mouth, (d) Tu-Ti-Kung-Pi, (e) Taipei Bridge, (f) Ru-Kou-Weir, (g) Hsin-Hai Bridge, (h) Chung-Cheng Bridge, and (i) Ta-Chih Bridge.

(14)

at those upriver transects close to tidal limits, the velocity measurements were made at two depths near bottom and near surface. The measured velocities near top and bottom are compared with model computed velocities at ﬁve stations: (a) Kuan-Du Bridge, (b) Taipei Bridge, (c) Hsin-Hai Bridge, (d) Chung-Cheng Bridge, and (e) Pa-Ling Bridge,Fig. 5a–e, for May 5, 2000. While there are a certain amount of uncertainties in these mea-surements (instrument limitations), these comparisons show that the model satisfactorily predicted the velocity in both the top and bottom layers along the channel.

In most estuaries when significant freshwater discharges and stratification exist, salinity becomes a natural conservative tracer for studying the mixing processes. Salinity distribution in an estuary is affected by the tidal current, freshwater discharge, as well as tidal and turbulent mixing processes. Therefore the resulting salinity distribution reflects the combined result of all processes, and in turn it controls (density) gravitational circu-lation and modifies mixing processes. Only limited long-term salinity time-series exists at the Zhu-Wei man-grove zone collected by the Taiwan Industrial Technology Institute; this time-series was used for model calibration and verification. In the numerical model, the salinities values at open boundaries at the coastal sea were set to 35 ppt. The salinity boundary conditions at heads of three tributaries were set to 0 ppt along with the specification of daily freshwater discharges. The wind conditions were included in the model simula-tion (Fig. 6b and e). Two one-week time periods were examined. During November 15 through November 22, 2000 (Fig. 6a), fresh-water discharge was receding from 300 m3/s to about 150 m3/s. The computed salinity time series compared very favorably with the discrete salinity measurements at the Zhu-Wei mangrove zone,

Fig. 6c. The computed salinities reproduce the pattern of observed salinity variations in a dynamic variation of

0 24 48 72 96 120 144 168 192 216 240 Time (hr) -3 -2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5 3 Elevation (m) Simulation Measurement Hsin-Hai Bridge 8/21 8/31 0 24 48 72 96 120 144 168 192 216 240 Time (hr) -3 -2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5 3 Elevation (m) Simulation Measurement Ta-Chih Bridge 8/21 8/31 0 24 48 72 96 120 144 168 192 216 240 Time (hr) -3 -2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5 3 Elev at io n (m ) Simulation Measurement Chung-Cheng Bridge 8/21 8/31 Fig. 4 (continued)

(15)

0 8 12 16 20 24 Time (hr) -3 -2 -1 0 1 2 3 Longitudinal velocity (m/s)

Simulation (Top layer) Simulation (Bottom layer) Measurement (Top layer) Measurement (Bottom layer)

Kuan-Du Bridge ebb flood 0 12 16 20 24 Time (hr) -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 Longitudinal velocity (m/s)

Simulation (Top layer)

Simulation (Bottom layer) Measurement (Top layer) Measurement (Bottom layer)

Chung-Cheng Bridge ebb flood 0 12 16 20 24 Time (hr) -1.6 -1.2 -0.8 -0.4 0 0.4 0.8 1.2 1.6 Longitudinal velocity (m/s)

Taipei Bridge ebb flood 0 12 16 20 24 Time (hr) -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 Longitudinal velocity (m/s)

Pa-Ling Bridge ebb flood 0 8 12 16 20 24 Time (hr) -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 Longitudinal velocity (m/s)

Hsin-Hai Bridge ebb flood 4 4 8 8 4 4 8 4

Fig. 5. The comparison of computed longitudinal velocity with time series data at: (a) Kuan-Du Bridge, (b) Taipei Bridge, (c) Hsin-Hai Bridge, (d) Chung-Cheng Bridge, and (e) Pa-Ling Bridge.

(16)

salinity between 0.0 and 33 ppt in a tidal cycle. The dynamic variations of salinity covered a range of about 30 ppt, and mean salinity value increased in response to decrease fresh water discharge. The salinity was measured at 1.5 m below surface and the computed salinity was taken from the surface layer, but salinity is generally well mixed in the vertical at this site. In the week between December 15 and 22, 2000, fresh water

11/15 11/16 11/17 11/18 11/19 11/20 11/21 11/22 Time (day) 0 50 100 150 200 250 300 350 400 Discharge ( m 3/s ) 12/15 12/16 12/17 12/18 12/19 12/20 12/21 12/22 Time (day) 0 50 100 150 200 250 300 350 400 Discharge ( m 3 /s ) 0 24 48 72 96 120 144 168 Time (hr) -12 -9 -6 -3 0 3 6 9 12 Wind vector (m /s ) 11/15 11/16 11/17 11/18 11/19 11/20 11/21 11/22 0 24 48 72 96 120 144 168 Time (hr) -12 -9 -6 -3 0 3 6 9 12 Wind vector (m /s ) 12/15 12/16 12/17 12/18 12/19 12/20 12/21 12/22 0 24 48 72 96 120 144 168 Time (hr) 0 5 10 15 20 25 30 35 40 45 Salinity ( ppt ) 11/15 11/16 11/17 11/18 11/19 11/20 11/21 11/22

Zhu-Wei mangrove zone Simulation Measurement 0 24 48 72 96 120 144 168 Time (hr) 0 5 10 15 20 25 30 35 40 45 50 Salinity ( ppt) 12/15 12/16 12/17 12/18 12/19 12/20 12/21 12/22 Zhu-Wei mangrovezone Simulation Measurement

Fig. 6. The comparison of computed salinity with time series data at Zhu-Wei mangrove zone in 2000 (a) discharge, (b) wind, and (c) salinity from 15 to 22 November; (d) discharge, (e) wind, and (f) salinity from 15 to 22 December.

(17)

discharge was in an increasing trend from 150 m3/s to 300 m3/s (Fig. 6d). While the measured salinity was as high as 37 ppt suggesting that substantial evaporation might be present but the evaporation mechanism was not built in the model. Nonetheless, the UnTRIM model reﬂected the large dynamic variations of salinity in a tidal cycle, and decreasing trend of mean salinity as fresh water increased near the end of this week. In these comparisons, the constants, a = 0.0115 and b = 0.75, gave best results in the turbulence closure model.

5. Model investigations

5.1. Freshwater discharge forcing

River discharge constituents the most important forcing since there would be no estuary or buoyant out-ﬂow without freshwater input. The inﬂuence of river discharge on salinity intrusion and estuarine circulation has been investigated in the Danshuei River estuarine system[25,26], however the river plume characteristics have not really studied. The validated three-dimensional hydrodynamic model was then used to investigate the river plume in Danshuei River mouth to variations of freshwater discharge conditions. Five tidal constituents, M2, S2, N2, K1, and O1, were used to generate a time series of water surface elevations as open boundary

con-ditions and a constant salinity 35 ppt were specified at the coastal sea for a three-month model simulation. The amplitudes and phases used for model simulation at the coastal sea boundaries presented inTable 1. No wind stress forcing was considered in the model. The model was run with low flow (Q50, a flow with an exceedence

probability of 50%) and high ﬂow (Q10) conditions. The freshwater discharges at the upstream boundaries of

the Tahan Stream, Hsintien Stream, and Keelung River were listed inTable 3.

Fig. 7presents the tidally averaged salinity distribution under Q10and Q50ﬂow conditions. A band of

low-salinity water is formed along the coastline. When the freshwater discharge increases, the river plume distance is far from the coastline. The salinity (arbitrarily limited by the 34 isohaline) due to the river plume is located at the distance 3 km from the Danshuei River mouth under Q10ﬂow, while salinity is at 1–2 km distance from

the Danshuei River mouth under Q50ﬂow condition. The high freshwater discharge results in the extension of

low-salinity. The residual currents Q10and Q50ﬂow conditions are shown inFig. 8. The anticyclonic

circula-tion is included along the north to west coast. The plumes and the buoyancy-driven coastal currents from the river are linked. The magnitude of the computed residual current is about 0.15–0.3 m/s that is the same order of magnitude measured by Wang [28]. With the pure freshwater discharge and tidal eﬀects including in the model simulation, the results reveal that freshwater discharge exerts a signiﬁcant impact on the plume distance and salinity distributions in Danshuei River coast region.

5.2. Wind stress forcing

In the experiments that follow, wind stress and tide are the forcing function, so the salinity remains constant and the velocity ﬁeld is purely wind-driven. Five tidal constituents, M2, S2, N2, K1, and O1, are forced at the

sea boundaries (Table 1). In Taiwan, the northeast wind prevails during the autumn and winter seasons, while the southwest wind prevails during the spring and summer seasons. The moderate to strong wind speed ranges

Table 3

Freshwater discharges at upstream boundaries

Upstream boundaries Freshwater discharge (m3_/s)

Q50 Q10

Tahan Stream 11.66 91.7

Hsintien Stream 31.7 149.5

(18)

from 6 to 15 m/s [23]. Therefore the wind-driven experiments include two main wind directions (coming from northeast and southwest) and employ a magnitude of 10 m/s. The fields of model computed mean sur-face elevation are shown inFig. 9. It is interesting to note that significant variation in wind direction cause quite different elevation field and, therefore, circulation pattern. Furthermore, the wind-driven velocity field (Fig. 10) is greatly influenced by the fact that the Danshuei River coastal sea is a topographic steep basin. The present of shallow coastline and deep-sea region enhance the spatial variation of the flow. When the wind direction is from northeast, the mean surface elevation is higher than that wind direction is from southeast

(19)

(Fig. 9). When the wind direction is from northeast, the residual current is higher than that wind direction is from southeast (Fig. 10).

Fig. 11presents the tidally averaged salinity distribution with different wind direction. Because no freshwa-ter discharges are considered in the model simulation, the plume of saline wafreshwa-ter is confined at the Danshuei River mouth. However, the wind direction also affects the upwelling/downwelling of the saline water. In the case of model simulation, the wind stress and tidal forcing are included. The results reveal that the wind direc-tion has a dominant effect on the velocity field in the Danshuei River coastal sea.

(20)

5.3. The combined eﬀect of buoyancy and wind stress

The combined eﬀect of wind-driven (barotropic), buoyancy-driven (baroclinic), and tidal forcing is inves-tigated in the model simulation. Two simulations, corresponding to the main wind directions (southwest and northeast) have been performed. The wind intensity is assumed 10 m/s and the Q10freshwater discharge from

the upriver boundaries (i.e. Tahan Stream, Hsintien Stream, and Keelung River) is imposed. Five tidal con-stituents, M2, S2, N2, K1, and O1, are also forced at the sea boundaries (seeTable 1).

Fig. 9. Mean surface elevation ﬁelds (in meter) for pure wind stress forcing (10 m/s); wind direction is from (a) Northeast and (b) Southeast.

(21)

Fig. 12shows the surface plume of tidally averaged salinity distribution under wind directions, southwest and northeast. When winds are downwelling-favorable, the surface low-salinity waters are flushed out and move to southwest coast. Conversely, large amounts of low-salinity water flushed out the Danshuei River mouth during upwelling-favorable winds, as the buoyancy-driven circulation is reversed (Fig. 13). In this study of wind action on estuarine plume, Chao[5]obtained similar circulation and notes that surface currents are dominated by wind-induced Ekman drifts. He also mentioned that when the wind is parallel to the coast and downwelling favorable, it reinforces the coastal jet, pushing the surface light water and enhancing the stratification.

(22)

The velocity fields (Fig. 13) exhibit the strongest southwest and northeast coastal current depended on the wind directions; this current has a wind-driven (barotropic) and buoyancy-driven (baroclinic) component. The flow field is reduced at the shallow coastline. In both cases the circulation is largely modified by wind stresses. As mentions in Chao[5], the surface circulation is affected by Ekman drift. This circulation leads to a new shape for the surface plume. On the other hand, residual bottom water circulations are less affected by wind surface stress (not shown). In the case of model simulation, the freshwater discharge, wind stress, and tidal

(23)

forcing are included. The results reveal that the combination of wind-driven and buoyancy-driven enhances the plume region along the Danshuei River coastal area.

(24)

6. Conclusions

A three-dimensional numerical hydrodynamics model, UnTRIM, was implemented and applied to the Danshuei River estuarine system and its adjacent coastal sea in northern Taiwan. The model was calibrated

(25)

against five major astronomical tides through the adjustment the bottom friction coefficient represented by Manning–Chezy formula. The final model calibration was accomplished when the tidal amplitudes and phases of these astronomical tidal constituents were accurately reproduced. The model was further verified by com-paring directly with field measured water surface elevation, current and salinity under variable daily freshwa-ter discharges from upriver of the Tahan Stream, Hsintien Stream, and Keelung River in year 2000.

The verified model was then applied to investigate the influence of wind stress and freshwater discharge on Danshuei River plume. The numerical experiments were conduced with three cases with buoyancy-driven forc-ing, wind stress forcforc-ing, and combined effect of buoyancy and wind stress. In these three cases, the tidal influ-ences are included in the model simulations. The results reveal that when winds are downwelling-favorable, the surface low-salinity waters are flushed out and move to southwest coast. Conversely, large amounts of low-salinity water flushed out the Danshuei River mouth during upwelling-favorable winds, as the buoy-ancy-driven circulation is reversed. If the pure freshwater discharge and tidal forcing are included, the mag-nitude of buoyancy-driven significantly affects the plume distance from the Danshui River mouth to adjacent coastal sea. In the case of model simulation, the freshwater discharge, wind stress, and tidal forcing are included. The results reveal that the combination of wind-driven and buoyancy-driven enhances the plume region along the Danshuei River coastal area.

Acknowledgements

The project, under which this study was conducted, was supported by the National Science Council, Tai-wan, under grant number 94-2211-E-239-011. The ﬁnancial support is highly appreciated. We also thank the Taiwan Water Resources Agency, Industrial Technology Institute, and Central Weather Bureau for providing the prototype data. The authors also would like to express their appreciation to the manuscript reviewers; through their comments this paper was substantially improved.

References

[1] V.H. Kourafalou, L.Y. Oey, J.D. Wang, T.N. Lee, The fate of river discharge on the continental shelf. 1. Modeling the river plume and the inner shelf coastal current, J. Geophys. Res. 101 (C2) (1996) 3415–3434.

[2] V.H. Kourafalou, River plume development in semi-enclosed Mediterranean regions: North Adriatic Sea and Northwestern Aegean Sea, J. Mar. Syst. 30 (3–4) (2001) 181–205.

[3] T.W. Kao, The dynamics of small scale fronts; Part I. Shelf water structure due to freshwater discharge, J. Phys. Oceanogr. 11 (1981) 1215–1223.

[4] R.W. Garvine, Estuary plumes and fronts in shelf waters: a layer model, J. Oceanogr. 17 (1987) 1877–1896. [5] S.Y. Chao, River-forced estuarine plumes, J. Phys. Oceanogr. 18 (1988) 72–88.

[6] L.Y. Oey, G.L. Mellor, Subtidal variability of estuarine outﬂow, plume and coastal current: a model study, J. Phys. Oceanogr. 23 (1) (1993) 164–171.

[7] A.E. Yankovsky, D.C. Chapman, A simple theory for the fate of buoyant coastal discharge, J. Phys. Oceanogr. 27 (1997) 1386–1401. [8] V.H. Kourafalou, T.N. Lee, L.Y. Oey, J.D. Wang, The fate of river discharge on the continental shelf 2. Transport of coastal

low-salinity waters under realistic wind and tidal forcing, J. Geophys. Res. 101 (C2) (1996) 3435–3455.

[9] E.N. Mikhailova, N.B. Shapiro, Modelling of propagation and transformation of riverine waters on a northwestern shelf and in an abyssal part of the Black Sea, J. Phys. Oceanogr. 8 (3) (1997) 169–177.

[10] R.W. Garvine, Penetration of buoyant coastal discharge onto the continental shelf: A numerical model experiment, J. Phys. Oceanogr. 29 (1999) 1892–1909.

[11] I.G. Berdeal, B.M. Hickey, M. Kawase, Inﬂuence of wind stress and ambient ﬂow on a high discharge river plume, J. Geophys. Res. 107 (C9) (2002) 3130–3153.

[12] S. Arnoux-Chiavassa, V. Rey, P. Fraunie, Modeling 3D Rhone river plume using a higher order advection scheme, Oceanol. Acta 26 (2003) 299–309.

[13] G. Lacroix, K. Ruddick, J. Ozer, C. Lancelot, Modelling the impact of the Scheldt and Rhine/Meuse plumes on the salinity distribution in Belgian water (southern North Sea), J. Sea Res. 52 (2004) 149–163.

[14] M.M. Whitney, R.W. Garvine, Simulating the Delaware Bay buoyant outﬂow: Comparison with observations, J. Phys. Oceanogr. 36 (1) (2006) 3–21.

[15] A.F. Blumberg, Turbulent mixing processes in lakes, reservoirs and impoundments, in: W.G. Gray (Ed.), Physical-based Modeling of Lakes, Reservoirs, and Impoundments, ASCE, 1986, pp. 79–104.

[16] V. Casulli, R.T. Cheng, Semi-implicit finite difference methods for three-dimensional shallow water flow, Int. J. Numer. Methods Fluids 15 (1992) 629–648.

(26)

[17] V. Casulli, R.V. Walters, An unstructured grid, three-dimensional model based on the shallow water equations, Int. J. Numer. Methods Fluids 32 (2000) 331–348.

[18] V. Casulli, P. Zanolli, Semi-implicit modeling of nonhydrostatic free-surface ﬂows for environmental problems, Math. Comput. Modell. 36 (9–10) (2002) 1131–1149.

[19] C. Chen, H. Liu, An unstructured grid, ﬁnite-volume, three-dimensional, primitive equations ocean model: application to coastal and estuaries, J. Atmos. Ocean. Technol. 20 (2003) 159–186.

[20] S. Sankaranarayanan, A 3D boundary-ﬁtted barotropic hydrodynamic model for the New York Harbor region, Contin. Shelf Res. 25 (18) (2005) 2233–2260.

[21] L. Zhong, M. Li, Tidal energy ﬂuxes and dissipation in the Chesapeake Bay, Contin. Shelf Res. 26 (6) (2006) 752–770.

[22] R.T. Cheng, V. Casulli, Evaluation of the UnTRIM model for 3-D tidal circulation, Proc. 7th Int. Conf. Estuarine Coast. Model., Florida (2001) 628–642.

[23] S.F. Su, Variation of Tidal Current Vertical Structure at Tanshui Estuary Region of Freshwater Inﬂuence. M.Sc. thesis, National Sun Yat-Sen University, 2006 (in Chinese).

[24] M.H. Hsu, A.Y. Kuo, J.T. Kuo, W.C. Liu, Procedure to calibrate and verify numerical models of estuarine hydrodynamics, ASCE J. Hydraul. Eng. 125 (2) (1999) 166–182.

[25] W.C. Liu, M.H. Hsu, A.Y. Kuo, Investigation of long-term transport in Tanshui River estuary, Taiwan, ASCE J. Waterway, Port, Coastal, Ocean Eng. 127 (2) (2001) 61–71.

[26] W.C. Liu, M.H. Hsu, A.Y. Kuo, J.T. Kuo, The inﬂuence on river discharge on salinity intrusion in the Tanshui River estuary, Taiwan, J. Coastal Res. 17 (3) (2001) 544–552.

[27] W.C. Liu, M.H. Hsu, C.R. Wu, C.F. Wang, A.Y. Kuo, Modeling salt water intrusion in Tanshui River estuarine system-Case-study contrasting now and then, ASCE J. Hydraul. Eng. 130 (9) (2004) 849–859.

[28] Y.H. Wang, Mapping ﬂow using towed-ADCP in coastal water of Taiwan, Proc. 24th Ocean Eng. Conf. Taiwan (2002) 485–490, in Chinese.

[29] K.F. Bowden, P. Hamilton, Some experiments with a numerical model of circulation and mixing in tidal estuary, Est. Coast. Mar. Sci. 3 (3) (1975) 281–301.

[30] V. Casulli, E. Cattani, Stability accuracy and eﬃciency of a semi-implicit method for three-dimensional shallow water ﬂow, Comput. Math. Appl. 27 (1994) 99–112.

[31] V. Casulli, G.S. Stelling, Numerical simulation of 3D quasi-hydrostatic free-surface ﬂows, ASCE J. Hydraul. Eng. 124 (7) (1998) 678–686.

[32] S. Rebay, Eﬃcient unstructured mesh generation by means of Delaunay triangulation and Bowyer-Watson algorithm, J. Comput. Phys. 106 (1993) 125–138.

[33] R.T. Cheng, J.R. Burau, J.W. Gartner, Interfacing data analysis and numerical modelling for tidal hydrodynamic phenomena, in: B.B. Parker (Ed.), Tidal Hydrodynamics, John Wiley & Sons, New York, 1991, pp. 201–219.

[34] R.T. Cheng, J.W. Garnter, V. Casulli, Tidal, residual, intertidal mudﬂat (TRIM) model and its applications to San Francisco Bay, California, Est., Coast., Shelf Sci. 36 (3) (1993) 235–280.

[35] R. Bradbury, User’s guide for Argus One, Argus Interware, Inc, 1997.

[36] S. Jan, Y.H. Wang, S.Y. Chao, D.P. Wang, Development of a nowcast system for the Taiwan Strait (TSNOW): numerical simulation of barotropic tides, Ocean Polar Res. 23 (3) (2001) 195–203.