Effect of channel connection on flow and salinity distribution of Danshuei River estuary

(1)

Eﬀect of channel connection on ﬂow and salinity

distribution of Danshuei River estuary

Wen-Cheng Liu

a,*

, Ming-Hsi Hsu

b

, Albert Y. Kuo

c

, Hsiao-Ying Hung

b

a

Department of Civil and Disaster Prevention Engineering, National United University, 1 Lien Da, Kung-Ching Li, Miao-Li 36003, Taiwan b

Department of Bioenvironmental Systems Engineering, National Taiwan University, Taipei 10617, Taiwan c

National Center for Ocean Research, National Taiwan University, Taipei 10617, Taiwan Received 1 May 2005; received in revised form 1 December 2005; accepted 20 February 2006

Available online 5 June 2006

Abstract

Numerical models are often used to evaluate the potential impact of human alteration of natural water bodies and to help the design the alternation to mitigate its impacts. A vertical (laterally integrated) two-dimensional hydrodynamic model was expanded to include the capability of simulating river loops as well as tributaries. The model was performed and applied to the Danshuei River estuarine system in northern Taiwan which consists of three major tributaries: the Tahan Stream, Hsintien Stream, Keelung River, and one river loop under the Chung-Hsin Bridge. The expanded model was reverified with observational field data of 2000. The verified model was then used to hindcast the river hydrodynamic conditions with a loop connection between the Danshuei River and Keelung River, which existed prior to 1965. It was found that the configuration of river loop connection has significant impacts on the residual transport along the connecting channel and the salinities in the connected river branches. The results show that the model may provide an ideal tool for management decision.

Keywords: Estuary; Channel connection; Hydrodynamic model; Danshuei River system; Salinity distribution; Taiwan

1. Introduction

The estuary, together with the tidal freshwater river upstream, comprises a pathway for exchange of water and materials between a drainage basin and coastal sea. However, various pollution problems seriously dam-age the water resources and ecological environment in a number of estuaries. To remedy this harmful situa-tion, it is necessary to understand the transport processes in the estuary such that the most eﬀective management strategy may be devised.

A numerical model is a powerful tool for the understanding of the hydrodynamic characteristics in estuar-ies. However the majority of past and existing vertical (laterally averaged) two-dimensional numerical models

* _{Corresponding author. Tel.: +886 37 381674; fax: +886 37 326567.} E-mail address:wcliu@nuu.edu.tw(W.-C. Liu).

(2)

of estuarine system are suitable only to a single stem estuaries, or estuaries with tributaries[1–12]. These mod-els do not simulate the channel network in some natural rivers. The purpose of this study is to remedy the deﬁciency by expanding a vertical two-dimensional model to simulate complicated river loops or distributaries as well as tributaries.

The Danshuei River estuarine system (Fig. 1) is the largest estuary in Taiwan, with the capital city of Taipei

in the drainage basin. 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 the three major tributaries: the Tahan Stream, the Hsintien Stream, and the Keelung River, and a loop under the Chung-Hsin Bridge. Except for the period of ﬂood 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. 1). Tidal propagation

is the dominant mechanism controlling the water surface elevation, and ebb and ﬂood ﬂows. The M2tide is the

primary tidal constituent at the river mouth, with mean tidal range of 2.17 m, and up to 3 m at spring tide. Because of the cross-sectional contraction and wave reﬂection, the mean tidal range may reach a maximum of 2.39 m within the system. The phase relationship between tidal elevation and tidal ﬂow is close to producing standing wave.

Sea water intrudes upriver as a result of tidal dispersion and the classical two-layer estuarine circulation. Salinity varies at intra-tidal time scale in response to the ebb and flood of the flows as well as in various longer time scales in response to freshwater inflow. The limit of salt intrusion may reach beyond 25 km in Tahan Stream from the river mouth during low flow. The baroclinic pressure gradient set up by the salinity distribu-tion is large enough to push the denser salt water upriver along the bottom layer of the estuary. This results in the classical two-layer circulation of net upriver flow in the bottom layer and net downriver flow in the upper

layer. The resulting estuarine circulation strengthens as the river ﬂow decreases[13].

The branched vertical two-dimensional estuarine hydrodynamic model developed by Hsu et al. [14,15]

employs a z coordinate in the vertical direction and uses a two-time level, ﬁnite diﬀerence numerical method. The model was applied to study the Danshuei River estuary. The model had previously been calibrated and

veriﬁed for the recent conditions of 1994 and 1995[13]. The model was expanded to include the capability of

Fig. 1. Map of the Danshuei River system (lines across the river are model transects at 0.5-km intervals; numbers are model segment numbers).

(3)

modelling river loop. The expanded model was reverified with the prototype data of 2000 using the bathymet-ric configuration measured in the same year. The verified model was then used to hindcast the river loop con-nection between the Danshuei River and Keelung River.

2. Model formulation

Many numerical models have been developed to simulate the hydrodynamic behaviors of estuaries using one-dimensional, two-dimensional, or three-dimensional framework. Although the three-dimensional models

[16–20]have increasingly been used, a two-dimensional laterally integrated model oﬀers and eﬃcient and prac-tical tool for the narrow and partially mixed estuaries. A branched, laterally integrated, two-dimensional, real-time model of hydrodynamics and salinity was developed for application to the tidal portion of the Danshuei

River estuarine system[14,15]. The model is based on the principles of conservation of water volume, water

momentum, and water mass. The model was expanded to include the simulation of river loop system. The following presents only the additions dealing with the river loop.

2.1. Basic equations

In a right-handed Cartesian coordinate system, with the x-axis directed seaward and the z-axis directed upward, the equations for the hydrodynamic model are as followings:

laterally integrated continuity equation oðuBÞ

ox þ

oðwBÞ

oz ¼ qp; ð1Þ

laterally integrated momentum equation oðuBÞ ot þ oðuBuÞ ox þ oðuBwÞ oz ¼ B q op oxþ o ox AxB ou ox þ o oz AzB ou oz ; ð2Þ hydrostatic equation op ox¼ qg og oxþ g Z g z oq oxdz; ð3Þ

sectionally integrated continuity equation o otðBggÞ þ o ox Z g H ðuBÞdz ¼ q; ð4Þ

laterally integrated mass balance equation for salt oðsBÞ ot þ oðsBuÞ ox þ oðsBwÞ oz ¼ o ox KxB os ox þo oz KzB os oz þ S0; ð5Þ equation of state q¼ q0ð1 þ ksÞ; ð6Þ

where x is distance seaward along river axis; z is distance upward in vertical direction; t is time; qpis lateral

inﬂow per unit lateral area; q is lateral inﬂow per unit river length; g is position of the water surface above mean sea level; s is laterally averaged salinity; u and w is laterally averaged velocity components in the x

and z directions, respectively; B is river width; Bgis width at the water surface including side storage area;

H is total depth below mean sea level; p is pressure; g is gravitational acceleration; Azand Kzare turbulent

viscosity and diﬀusivity, respectively, in z-direction; Ax and Kx are dispersion coeﬃcients for momentum

and mass, respectively, in the x-direction; q and q0are water density and freshwater density, respectively; k

is constant relating density to salinity (=7.5· 104ppt1); S0is source and sink of salt due to exchange with

(4)

Eqs.(1), (2) and (5)are integrated over a layer of ﬁnite thickness and then solved with Eqs.(3), (4) and (6)

by a two-time level, finite difference scheme with spatially staggered grid. The implicit treatment of the vertical diffusion terms results in a tri-diagonal matrix in the vertical direction. To improve numerical stability, the

pressure gradient term in Eq.(3) is evaluated using g at a new time step. The method of solution is detailed

in Park and Kuo[7,21]and Hsu et al.[14,15]. The following presents only the additions dealing with river

loop.

2.2. Treatment of free surface elevation at the loop junction

At the junction of the mainstem and loop (Fig. 2), the ﬁnite diﬀerence equation for free surface elevation,

discretized from of continuity equation, is written as

gi;2¼ gi;1

þ Dt

Bi;1þ STBiþ BEJi;1

Bi;1þ Bi;2 2 þ BEJi;1 wi;1;1 1 DxðTMPiþ1 TMPiþ TMPBmÞ þ Q_i;1 Dx ; ð7Þ

where gi,2and gi,1represent surface elevations at new and old time levels, respectively; Dt is time step; Dx is

segment length; Bi,kis width of the kth layer in the ith segment; STBiis equivalent width of the side storage

area; BEJi,1is additional equivalent width contributed by branching; wi,k,1is vertical velocity at the bottom of

the kth layer at old time level; Qi,kis lateral inﬂow to the kth layer. The advective terms are

TMPi¼ Bi1;1þ Bi;1 2 ðh1þ gi1;1Þ þ ðh1þ gi;1Þ 2 ui;1;1; TMPBm¼ TBm1;1þ TBm;1 2 ðTh1þ T gm1;1Þ þ ðTh1þ T gm;1Þ 2 Tum;1;1;

where h1is the thickness of the ﬁrst layer below mean sea level; ui,k,1is longitudinal velocity at old time step;

TB, Tg, Th and Tu is width, surface elevation, layer thickness and longitudinal velocity, respectively in the

(5)

loop. The subscripts 1, i, and m denote the mainstem segment, loop segment and the top layers, respectively. The side storage area in each segment acts as a source for the main channel on the falling tide and as a sink on the rising tide. The junction segment is considered as the end segment of the loop such that the longitudinal

velocity for the loop Tum,k,2is 0, and

Tg_m;1¼ gi;1;

Tg_m;2¼ gi;2:

ð8Þ 2.3. Treatment of the conservation of water mass at the loop junction

At the junction of mainstem and loop (Fig. 3), the ﬁnite diﬀerence form of the laterally integrated

continu-ity equation, Eq.(1), is written as

wi;k;2¼ 1 Bi;kþBi;kþ1 2 þ BEJi;k Bi;kþ1þ Bi;kþ2 2 þ BEJi;kþ1 w_i;kþ1;2hkþ1 Dx Bi;kþ1þ Biþ1;kþ1 2 uiþ1;kþ1;2 Bi1;kþ1þ Bi;kþ1 2 ui;kþ1;2þ TBm;kþ1þ TBmþ1;kþ1 2 Tum;kþ1;2 þQi;kþ1 Dx ð9Þ

for k = 2 to N (bottom layer) and wi,N,2= 0.

2.4. Treatment of the conservation of longitudinal momentum Tum+1at the loop

Because of the spatially staggered grid used in the model, no representative longitudinal velocity is situated at

the junction point of the model (TuminFig. 2). The momentum balance at the section between the loop and

mainstem is handled by neglecting the horizontal advective term o

oxðukBkukhkÞ and turbulent diﬀusion term

o

ox½ðAxÞkBkhkðou_oxÞk. The magnitudes of these two terms in momentum equation are assumed to be locally

negli-gible in comparison to other terms. Computational tests demonstrate that this is a realistic assumption[14,22].

2.5. Treatment of the conservation of salt at the loop junction

At the junction segment, the additional fluxes of salt from or to the loop are fully accounted for. The QUICKEST scheme is used to express the finite difference form of advective term. It is based on a conservative

Upriver River loop Downriver 2 , 1 − i 1 , 1 , 1 − i u 1 , 1 , i u 1 , 1 , 1 + i u 1 , 1 , 2 + i u 2 , 1 + i

η

2 , 2 , m i Tη η = 2 , 1 + m Tη η 1 , 1 , 2 + m Tu Tum+1,1,1 Junction

(6)

control volume formation and estimates cell wall concentrations with a quadratic interpolation using concen-trations in the two adjacent cells and that at the next upstream cell. The same method is used to treat the salt concentration at the junction of mainstem and loop. Since the ﬁrst segment of the loop is the junction segment in the mainstem, therefore

Tsm;k;1¼ si;k;1; Tsm;k;2¼ si;k;2 ð10Þ

and

Tsm1;k;1¼ si1;k;1;

Ts_m1;k;2¼ si1;k;2:

ð11Þ

3. Model calibration and veriﬁcation

For modelling, the Danshuei River–Tahan Stream is regarded as the mainstem of the river, while the Hsint-ien Stream and the Keelung River are treated as the first and second tributaries, respectively. A river loop exists under the Chung-Hsin Bridge as the result of a small island in the middle of river. A field survey in 2000 by the Taiwan Water Resources Agency measured the cross-sectional profiles approximately every 0.5 km along the tidal portion of the river system. The cross-sectional profiles were used to schematize the estuary in the geometric file for model input. The estuary is divided into 66, 29, 75, and 2 segments (Dx = 0.5 km) for the Danshuei River–Tahan Stream, Hsintien Stream, Keelung River, and the river loop under the Chung-Hsin Bridge, respectively. Because surface elevation at low tide is about 1.5 m below mean sea level at spring tide, the thickness of the top layer is 2.0 m to maintain the surface elevation above the bot-tom of the layer at all times. The thicknesses of the remaining layers are 1.0 m. The measured cross-sections in

0 5 10 15 20 25 30 35

Distance from Danshuei River mouth (km)

0 50 100 150 200 250 300

Mean tidal range (cm)

Simulation

Measurement

Simulation Measurement

Simulation Measurement (a) Danshuei River-Tahan Stream

0 2 4 6 8 10 12 14

Distance from Hsintien Stream mouth (km)

0 50 100 150 200 250 300

(b) Hsintien Stream

0 5 10 15 20 25 30 35 40 45

Distance from Keelung River mouth (km)

0 50 100 150 200 250 300 (c) Keelung River

(7)

each segment were averaged and schematized. First, the bottom elevation was determined at an integral num-ber of meters below the mean sea level. Then the width at each layer was determined while maintaining the cross-sectional area approximately the same as that in the measured proﬁles. The time step increment (Dt) of 25 s, which guaranteed stability, was used for all model simulations.

To avoid the upstream boundary conditions being inﬂuenced by tide, the computational domain is extended beyond the tidal limits in the tributaries as well as the mainstem. A single constituent, the dominant

M2tide, was used as a forcing function at the river mouth for model calibration. The use of a single

constit-uent expedites the extraction of tidal ranges from model outputs for comparison with prototype data. The

downstream boundary condition is an M2tide with amplitude being half of the mean tidal range, 108 cm,

cal-culated from the 365 days of measurement at the river mouth in 2000. The upstream boundary conditions at the three tributaries of Danshuei River–Tahan Stream, Hsintien Stream and Keelung River are 59.2, 87.1, and

32.8 m3/s, respectively, the mean discharges for the year 2000. The model was run for 30 tidal cycles to reach

dynamic equilibrium. The diﬀerence between the maximum and minimum surface elevations in the last tidal cycle was output as the tidal range for each segment. The friction coeﬃcient was adjusted until the model

out-puts agreed satisfactorily with the mean tidal ranges calculated from the 2000 data.Fig. 4presents simulated

and measured variation of the tidal range along the axes of the river and the tributaries. The tidal range is affected by the river geometry as well as bottom friction. The friction dissipates tidal energy and reduces the tidal range as the tide propagates upriver. At the uppermost reaches of all three branches the tidal range decreases sharply as the result of the rapid rise in the river bottom. Finally the tide reaches its limit where the river bottom rises well above high tide level. The average absolute values of the differences (mean absolute errors) between computed and measured tidal ranges are 1.4 cm, while the root-mean-square error is 1.9 cm. The model was verified by simulating the prototype conditions of 2000 without changing the values of the friction coefficient determined by model calibration. Time-varying boundary conditions were specified for the

35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 Number of days since 1/1/2000 -300 -250 -200 -150 -100 -50 0 50 100 150 200 250 300

Water surface elevation (cm) Water surface elevation (cm)

Water surface elevation (cm)

Simulation Measurement Simulation Measurement Simulation Measurement Simulation Measurement (a) Tu-Ti-Kung-Pi 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 Number of days since 1/1/2000

-300 -250 -200 -150 -100 -50 0 50 100 150 200 250 300 (b) Hsin-Hai Bridge 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 Number of days since 1/1/2000

-300 -250 -200 -150 -100 -50 0 50 100 150 200 250 300 (c) Chung-Cheng Bridge 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 Number of days since 1/1/2000

-300 -250 -200 -150 -100 -50 0 50 100 150 200 250 300 (d) Ta-Chih Bridge

Fig. 5. Model reveriﬁcation results: comparisons of measured and computed surface elevation: (a) Tu-Ti-Kung-Pi; (b) Hsin-Hai Bridge; (c) Chung-Cheng Bridge; and (d) Ta-Chih Bridge.

(8)

0 5 10 15 20 Time (hour) -150 -100 -50 0 50 100 150 200 Longitudinal v elocity (cm/s)

Top Layer (Simulation) Bottom Layer (Simulation) Top Layer (Measurement) Bottom Layer (Measurement) (a) 0 5 10 15 20 Time (hour) -150 -100 -50 0 50 100 150 200 Longitudinalv elocity (cm/s)

Top Layer (Simulation) Bottom Layer (Simulation) Top Layer (Measurement) Bottom Layer (Measurement) (b) 0 5 10 15 20 Time (hour) Time (hour) -100 -50 0 50 100 150 Longitudinal v elocity (cm/s)

Top Layer (Simulation) Top Layer (Measurement)

(c) 0 5 10 15 20 Time (hour) -150 -100 -50 0 50 100 150 200 Longitudinal v elocity (cm/s) Longitudinal v elocity (cm/s)

Top Layer (Simulation) Bottom Layer (Simulation) Top Layer (Measurement) Bottom Layer (Measurement) (d) 0 5 10 15 20 -150 -100 -50 0 50 100 150 200

Top Layer (Simulation) Bottom Layer (Simulation) Top Layer (Measurement) Bottom Layer (Measurement)

(e)

Fig. 6. Model reveriﬁcation results: comparison of measured and computed longitudinal velocity on May 5, 2000: (a) Kuan-Du Bridge; (b) Taipei Bridge; (c) Hsin-Hai Bridge; (d) Chung-Cheng Bridge; and (e) Pa-Ling Bridge.

Table 1

Mean absolute diﬀerences and root-mean-square diﬀerences between computed and measured longitudinal velocity

Layer Kuan-Du Bridge Taipei Bridge Hsin-Hai Bridge Chung-Cheng Bridge Pa-Ling Bridge Mean absolute difference (cm/s) Root-mean-square difference (cm/s) Mean absolute difference (cm/s) Root- mean-square difference (cm/s) Mean absolute difference (cm/s) Root- mean-square difference (cm/s) Mean absolute difference (cm/s) Root- mean-square difference (cm/s) Mean absolute difference (cm/s) Root- mean-square difference (cm/s) Surface 9.00 12.53 12.24 18.1 6.82 8.58 6.87 8.76 7.41 9.50 Bottom 8.03 10.09 8.57 12.28 – – 8.70 9.74 9.88 11.52

(9)

model run. Hourly measurements of water surface elevation at the Danshuei River mouth and daily average freshwater discharges upriver of the tidal limits were used for the boundary conditions. These conditions allow

the investigation of the model response to the interaction of tidal forcing and varying river discharges.Fig. 5

shows computed surface elevation, together with ﬁeld data at ﬁve locations. The comparison shows that the model can faithfully reproduce tidal propagation to within ±1.9 cm. An intensive survey was conducted by the Taiwan Water Resources Agency on May 5, 2000. Half-hourly measurements of water surface elevation and

velocity were made continuously for 13 daylight hours.Fig. 6shows the comparison of the time series data of

the longitudinal velocity for May 5, 2000. Table 1 presents the mean absolute diﬀerences and

root-mean-square diﬀerences between the computed and measured time series.Fig. 6shows that there are discrepancies

between model and field data when judging with point by point comparison. However, the model does prop-erly simulate the temporal variation of velocity in terms of current amplitude and phase. The calibrated and verified model has values of the Manning friction coefficient ranging from 0.026 to 0.035 for the Tahan River– Tahan Stream, 0.015 to 0.023 for the Hsintien Stream, 0.016 to 0.023 for the Keelung River and 0.034 for the river loop. The calibrated Manning friction coefficients are slightly different from previous calibration results

[13]due to the signiﬁcant bathymetric changes between year 2000 and 1994.

The estuarine salinity distribution is affected by turbulent diffusion coefficients as well as the freshwater dis-charges. The turbulent mixing term in the vertical direction is the dominant factor that determines stratifica-tion in the water column. The recommended procedure to calibrate the coefficients in the turbulent model is by

matching the observed and computed salinity distributions.Fig. 7shows a comparison of time series salinity

distribution between computed and measured data near the Zhu-Wei station (Fig. 1). The surface and bottom

layers are presented in the ﬁgure while 1.2 below the water surface of ﬁeld data is measured. The model results seem to underestimate the observed data. However the numerical model can favorably mimic the trend of salinity distribution.

4. Model project

There existed a river cutoﬀ connecting the Danshuei River and Keelung River prior to 1965 (Fig. 8). The

island surrounded by the two rivers and the cutoff suffered severe flood damage by the Typhoon Gloria in 1962. To mitigate the flood disaster, Taiwan government filled in the cutoff, thus eliminated the river loop con-figuration. The numerical model was used to hindcast the hydrodynamic conditions of the river cutoff loop to investigate if significant change has resulted from the elimination. Due to the lack of cross-sectional data, the cutoff was assumed to have a uniform rectangular cross-section, 75 m wide and 4 m deep, and divided into 2 segment with Dx = 0.5 km. A constant Manning friction coefficient of 0.015 was used in the cutoff.

To simulate the ﬂow conditions and salinity distributions, the model simulations were conducted using nine-constituent tides. These constituents were extracted from a harmonic analysis using the ﬁeld data at

320 322 324 326 328 330 Number of days since Jan. 1 2000

0 5 10 15 20 25 30 35 40 45 50 Salinity (ppt) Zhu-Wei Station

Simulation (surface layer) Simulation (bottom layer) Observation

(a)

November 16 November 26

350 352 354 356 358 360 Number of days since Jan. 1 2000

0 5 10 15 20 25 30 35 40 45 50 Salinity (ppt) Zhu-Wei Station

Simulation (surface layer) Simulation (bottom layer) Observation

(b)

December 16 December 26

Fig. 7. Model reveriﬁcation results: computed of measured and computed time series salinity distribution at the Zhu-Wei station: (a) November 16–26 and (b) December 16–26.

(10)

the mouth of the Danshuei River in 2000. The nine constituents are M2(12.42 h), S2(12 h), N2(12.9 h), K1

(23.93 h), Sa(8765.32 h), O1(25.82 h), K2 (11.97 h), P1(24.07 h), and M4 (6.21 h). Amplitudes and phases

of tidal constituents were speciﬁed to generate a time series of surface elevation for the downstream boundary condition lasting a one-year (705 tidal cycles). A high tide salinity of 35 ppt (parts per thousand) at the Dans-huei River mouth was used for model simulation. The long-term average river discharges at the tidal limits of the three major tributaries, Tahan Stream, Hsintien Stream, and Keelung River, were speciﬁed as 41.04, 57.46,

and 25.25 m3/s, respectively.

Fig. 9(a) presents the longitudinal velocity at the surface layer in the river cutoff. The figure reveals that the flow is in the ebb direction most of the time, that is, flowing to the Danshuei River. The water surface eleva-Fig. 8. River loop connection between the Danshuei River and Keelung River prior to 1965 (lines across the river are model transects at 0.5-km intervals; numbers are model segment numbers).

18.0 18.5 19.0 19.5 Time (day) -100 -50 0 50 100 Longitudinal velocity (cm/s) (a) Segment 182 ebb flood 18.0 18.5 19.0 19.5 Time (day) -150 -100 -50 0 50 100 150

Water surface elevation (cm)

-10 -5 0 5 10 Difference (cm) WSE (segment 181) WSE (segment 184) Difference (segment 181-184) (b)

(11)

18.0 18.5 19.0 19.5 Time (day) 0 2 4 6 8 10 Salinity (ppt)

(a) Segment 42 (no loop)

Segment 42 (looped estuary)

18.0 18.5 19.0 19.5 Time (day) 0 2 4 6 8 10 12 Salinity (ppt)

(b) Segment 172 (no loop)_{Segment 172 (looped estuary)}

Fig. 10. Time series salinity variations at (a) segment 42 (Danshuei River–Tahan Stream) and (b) segment 172 (Keelung River).

0.0 5.0 10.0 15.0 20.0 25.0 30.0

-12.0 -10.0 -8.0 -6.0 -4.0 -2.0 0.0 Depth (m)

(a) Danshuei River-Tahan Stream

0.0 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0 10.0

-5.0 -4.0 -3.0 -2.0 -1.0 0.0 Depth (m) (b) Hsintien Stream 0.0 5.0 10.0 15.0 20.0 25.0 Distance from Keelung River mouth (km)

-8.0 -6.0 -4.0 -2.0 0.0 Depth (m) (c) Keelung River

Fig. 11. Model hindcast of salinity distribution with looped estuary: (a) Danshuei River–Tahan Stream; (b) Hsintien Stream; and (c) Keelung River (the numbers on the contours refer to the salinity in parts per thousand).

(12)

tions at the junction segments in the Danshuei River and Keelung River are plotted inFig. 9(b). The difference in the water surface elevations between these two segments is consistent with the flow acceleration (or decel-eration) in the cutoff. When the water surface at segment 181 (in the Keelung River) is higher than that at segment 184 (in the Danshuei River), the flow at segment 182 (in the cutoff) is increasing in the ebb direction or decreasing in the flood direction, and vice versa.

Fig. 10compares the salinities in both the Danshuei River (segment 42) and Keelung River (segment 172) for the looped estuary with those for the existing condition (no river loop between). Because the cutoﬀ diverted the ﬂow to the Danshuei River, the salinity with the looped estuary prior to 1965 was higher in the Keelung

River and lower in the Danshuei River than those under existing condition.Figs. 11 and 12show the salinity

distributions averaged over 58 tidal cycles for the looped estuary and for the existing condition, respectively. The salt intrusion of the looped estuary was farther upriver than that of the existing condition in the Keelung River, because the flow was diverted through the cutoff reducing the freshwater discharge in the lower reach of the Keelung River. The difference in the overall salinity distributions in the Danshuei River–Tahan Stream is not significant, and there is little change in the Hsintien Stream.

0.0 5.0 10.0 15.0 20.0 25.0 30.0

-12.0 -10.0 -8.0 -6.0 -4.0 -2.0 0.0 Depth (m)

(a) Danshuei River-Tahan Stream

0.0 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0 10.0

-5.0 -4.0 -3.0 -2.0 -1.0 0.0 Depth (m) (b) Hsintien Stream 0.0 5.0 10.0 15.0 20.0 25.0

Distance from Keelung River mouth (km)

-8.0 -6.0 -4.0 -2.0 0.0 Depth (m) (c) Keelung River

Fig. 12. Model computed salinity distribution with no loop (existing condition): (a) Danshuei River–Tahan Stream; (b) Hsintien Stream; and (c) Keelung River (the numbers on the contours refer to the salinity in parts per thousand).

(13)

5. Conclusions

A vertical (laterally integrated) two-dimensional hydrodynamic numerical model was expanded to include the simulation of river loops and applied to the Danshuei River estuary in Taiwan. The model had previously been calibrated and verified with respect to the 1995 bathymetry using observed time series of water level, velocity, and salinity data of 1994 and 1995. The model was further reverified for the 2000 bathymetry using limited data sets collected in 2000. The results of the reverified model show reasonable agreement between the model predictions and measured data. The reverified model was then used to hindcast the effects of a connect-ing channel between the Danshuei River and Keelung River which existed prior to 1965. Two alternative cases of connection between the Danshuei River and Keelung River were also simulated.

The model simulations were conducted with long-term mean river discharges at the upstream boundaries of the three tributaries. A nine-constituent tide and a high tide salinity of 35 ppt at the river mouth were used as the downstream boundary conditions. The model results reveal that the longitudinal velocity in the pre-1965 river cutoff was directed toward the Danshuei River during most of the tidal cycle. This resulted in a stronger salt intrusion in the lower reach of the Keelung River should the cutoff have not been filled. It suggests that careful planning is required if the pre-1965 river cutoff is to be restored. The expanded model provides a tool to assist with water management through the understanding of the water flow and salinity distribution in the looped estuary.

Acknowledgements

This research was conducted as part of a grant supported by National Science Council, Taiwan, Grant Nos. 91-2211-E-002-038 and 92-2211-E-002-057. The writers also thank the Taiwan Water Resources Agency for providing the prototype data.

References

[1] J.F. Festa, D.V. Hansen, A two-dimensional numerical model of estuarine circulation: the eﬀects of altering depth and river discharge, Estuar. Coast. Mar. Sci. 4 (1976) 309–323.

[2] A.F. Blumberg, The inﬂuence of density variations on estuarine tides and circulations, Estuar. Coast. Mar. Sci. 6 (1978) 209–215. [3] D.P. Wang, D.W. Kravitz, A semi-implicit two-dimensional model of estuarine circulation, J. Phys. Oceanogr. 10 (1980) 441–454. [4] P.A. Parrels, M. Karelse, A two-dimensional laterally averaged model for salt intrusion in estuaries, in: H.B. Fischer (Ed.), Transport

Methods for Inland and Coastal Waters, Academic, San Diego, CA, 1981, pp. 483–535.

[5] M. Ford, J. Wang, R.T. Cheng, Predicting the vertical structure of tidal current and salinity in San Francisco Bay, California, Water Resour. Res. 26 (5) (1990) 1027–1045.

[6] J.W. Lavelle, E.D. Cokelet, G.A. Cannon, A model study of density intrusions into and circulation within a deep, silled estuary: Puget Sound, J. Geophys. Res. 96 (C9) (1991) 16779–16800.

[7] K. Park, A.Y. Kuo, A vertical two-dimensional model of estuarine hydrodynamics and water quality, Speciﬁc Report Applied Marine Science and Ocean Engineering, No. 321, Virginia Institute of Marine Science, Gloucester Point, VA, 1993.

[8] R.W. Cole, E.M. Buchak, CE-QUAL-W2: A two-dimensional, laterally averaged, hydrodynamic and water quality, User’s Manual, Instruction Report EL-95-1, US Army Engineer Waterways Experiment Station, Vicksburg, MS, 1995.

[9] J. Wang, A two-channel laterally averaged estuarine circulation model, J. Geophys. Res. 103 (C9) (1998) 18381–18391. [10] V. Botte, A. Kay, A model of the wind driven circulation in Lake Baikal, Dyn. Atmos. Oceans 35 (2002) 131–152.

[11] J.P. Sierra, A. Sanchez-Arcilla, P.A. Figueras, J. Gonzalez Del Rio, E.K. Rassmussen, C. Mosso, Eﬀects of discharges reductions on salt wedge dynamics of the Ebro river, River Res. Appl. 20 (2004) 61–77.

[12] W.C. Liu, J.T. Kuo, A.Y. Kuo, Modelling hydrodynamics and water quality in the separation waterway of the Yulin oﬀshore industrial park, Taiwan, Environ. Modell. Software 20 (2005) 309–328.

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

[14] M.H. Hsu, A.Y. Kuo, J.T. Kuo, W.C. Liu, Study of tidal characteristics, estuarine circulation and salinity distribution in Tanshui River system (I), Technical Report No. 239, Hydrotech Research Institute, National Taiwan University, Taipei, Taiwan, 1996 (in Chinese).

[15] M.H. Hsu, A.Y. Kuo, J.T. Kuo, W.C. Liu, Study of tidal characteristics, estuarine circulation and salinity distribution in Tanshui River system (II), Technical Report No. 273, Hydrotech Research Institute, National Taiwan University, Taipei, Taiwan, 1997 (in Chinese).

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

(14)

[17] A.F. Blumberg, L.A. Khan, J.P. John St, Thee-dimensional hydrodynamic model of New York harbor region, J. Hydraul. Eng., ASCE 125 (1999) 799–816.

[18] A.H. Al-Rabeh, R.W. Lardner, N. Gunay, GULFHYDRO version 2.0: a software package for hydrodynamics in the Arabian Gulf, Environ. Modell. Software 15 (2000) 23–34.

[19] A. Ikala, K. Myberg, Comparison of hydrodynamic models of the Gulf Finland in 1995: a case study, Environ. Modell. Software 17 (2002) 237–250.

[20] R. Perianez, GISPART: a numerical model to simulate the dispersion of contaminants in the Strait of Gibraltar, Environ. Modell. Software 20 (2005) 797–802.

[21] K. Park, A.Y. Kuo, A numerical model study of hypoxia in the tidal Rappahannock River of Chesapeake Bay, Estuar. Coast. Shelf Sci. 42 (5) (1996) 563–581.

[22] A.J. Elliot, A numerical model of the internal circulation in a branching estuary, Chesapeake Bay Institute Special Report No. 54, The Johns Hopkins University, 1976.