Modeling the influence of river discharge on salt intrusion and residual circulation in Danshuei River estuary, Taiwan

(1)

Continental Shelf Research 27 (2007) 900–921

Modeling the inﬂuence of river discharge on salt intrusion and

residual circulation in Danshuei River estuary, Taiwan

Wen-Cheng Liu

a,

, Wei-Bo Chen

b

, Ralph T. Cheng

c

,

Ming-Hsi Hsu

d

, Albert Y. Kuo

e

a

Department of Civil and Disaster Prevention Engineering, National United University, Miao-Li 36003, Taiwan b_{Department of Civil Engineering, 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} e_{National Center for Ocean Research, National Taiwan University, Taipei 10617, Taiwan} Received 3 May 2006; received in revised form 23 November 2006; accepted 4 December 2006

Available online 27 December 2006

Abstract

A 3-D, time-dependent, baroclinic, hydrodynamic and salinity model was implemented and applied to the Danshuei River estuarine system and the adjacent coastal sea in 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 verification in model simulations of barotropic and baroclinic flows. The turbulent diffusivities were ascertained through comparison 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 freshwater discharge on residual current and salinity intrusion under different freshwater inflow condition in the Danshuei River estuarine system. The model results reveal that the characteristic two-layered estuarine circulation prevails most of the time at Kuan-Du station near the river mouth. Comparing the estuarine circulation under low- and mean flow conditions, the circulation strengthens during low-flow period and its strength decreases at moderate river discharge. The river discharge is a dominating factor affecting the salinity intrusion in the estuarine system. A correlation between the distance of salt intrusion and freshwater discharge has been established allowing prediction of salt intrusion for different inflow conditions.

Keywords: 3-D model; Residual circulation; Salt intrusion; Danshuei River estuarine system and coastal sea; Freshwater discharge

1. Introduction

Estuaries, which act as the transition zone between upland rivers and the coastal ocean, are

important nursery and feeding grounds for a very large number of marine species. However, due to increasing urban development, they are being used more and more as repositories for direct point source discharges of contaminants, indirect pollu-tant input through non-point land sources and atmospheric pollutant deposition. Neither the estuaries nor the coastal ocean is capable of

www.elsevier.com/locate/csr

Corresponding author. Tel.: +886 37 381674; fax: +886 37 326567.

(2)

assimilating pollutants indefinitely and the environ-mental concerns now require that any pollutant released in the coastal zone should be carefully regulated and properly managed (Horrevoets et al., 2004; Bilgili et al., 2005). The long-term, large-scale transport and distribution of salt, sediment and pollutants are primarily determined by the residual net circulation and the dispersion processes under the influence of tidal flow. In the Danshuei River, the discharge of runoff from its drainage basin indicates a net movement of water toward the ocean (Liu et al., 2001). However, the intrusion of seawater from the ocean often results in a net estuarine circulation (or gravitational circulation) of upriver flow in the lower layer and downriver flow in the upper layer of the water column (Pritchard, 1956; Dyer, 1973). Due to the usually complex basin geometry, numerical models are often used in studies of water movement, tidal and residual circulation in estuaries, tidal rivers and coastal seas.

Hydrodynamic processes in estuaries and coastal seas are generally 3-D tidal ﬂows in a setting of complex geometry and bottom topography. Tidal circulation and solute transport in estuaries are strongly affected by complex geometry and bathy-metry. Hence, to accurately simulate tidal and residual circulation in coastal seas and estuaries driven by tides, winds and density gradients, the numerical model must be able to accurately and efﬁciently resolve the dynamics of various vertical boundary layers and the complex geometry and bottom topography.

A number of 3-D hydrodynamic models are available to simulate the tidal current, salinity and temperature in estuaries and coastal seas (e.g. Blumberg, 1986; Casulli and Cheng, 1992; Casulli and Walters, 2000; Ji et al., 2001; Casulli and Zanolli, 2002; Chen and Liu, 2003; Chen, 2003; Sankaranarayanan and McCay, 2003; Sankara-narayanan, 2005; Zhong and Li, 2006; Oliveira et al., 2006). In the present study, a 3-D high-resolution hydrodynamic model, UnTRIM (Cheng and Casulli, 2001; Casulli and Zanolli, 2002), was implemented and applied to simulate the water surface elevation, tidal current, residual circulation and salinity distribution in the Danshuei River estuarine system and its adjacent coastal sea in northern Taiwan. Model calibration and veriﬁca-tion were conducted for tidal ranges, water surface elevation, current and salinity by compa-ring them to available ﬁeld data. The model was

then applied to investigate the residual circula-tion and salt water intrusion under different freshwater discharge conditions from upstream boundaries.

2. Study site

The Danshuei River estuary (Fig. 1) is the largest estuarine system in Taiwan, with its drainage basin including the capital city of Taipei. The tidal influence 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 Stream, the Hsintien Stream and the Keelung River. 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. 1). 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 M2 tide 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-sectional contraction and wave reflection, the mean tidal range may reach a maximum of 2.39 m within the system. The phase relationship between tidal elevation and tidal flow is close to standing wave characteristics (Hsu et al., 1999).

Seawater intrudes upriver as a result of tidal dispersion and residual circulation (Liu et al., 2004). Salinity varies on the intra-tidal time scale in response to the ebb and flood of tidal flows as well as in response to changing freshwater inflows. The limit of the 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 (Liu et al., 2001; Hsu et al., 2003). The estuarine circulation strengthens as the river flow decreases (Kuo and Liu, 2001).

A 3-D model is used to cover the Danshuei River estuarine system and adjacent coastal sea where the bottom topography gradually deepens away from the coastline (Fig. 1), and the model uses an unstructured model grid (Fig. 2(a)). Wang (2002)

(3)

used a towed-ADCP to measure the tidal current in 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 on 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.

3. Model description

3.1. Governing equation

The governing equations describing 3-D free-surface ﬂows are the well-known shallow water equations. Such equations express the physical conservation principles of volume, mass and mo-mentum.

The volume conservation is expressed by the following incompressibility condition:

qu qxþ qv qyþ qw qz ¼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 (Casulli and Cheng, 1992): qZ qtþ q qx Z Z h u dz þ q qy Z Z h v dz ¼0, (2)

where t is the time, hðx; yÞ is the prescribed mean water depth measured from the undisturbed water surface and Zðx; y; tÞ is the free-surface elevation. Thus, Hðx; y; tÞ ¼ hðx; yÞ þ Zðx; y; tÞ is the total water depth.

The mass conservation of salt (solutes) is ex-pressed by the following equation:

qS qt þu qS qxþv qS qyþw qS qz ¼ q qx K h qS qx þ q qy K h qS qy þ q qz K v qS qz , ð3Þ Zhu-Wei Mangrove Zone

Gauging Station Tu-Ti-Kung-Pi Kuan-Du Bridge Kuan-Du Wetland Taiwan Strait Pa-Ling Bridge Keelung River Ta-Chih Bridge Chiang-Pei Bridge Wu-Tu Chung-Cheng Bridge Bao Bridge Taipei City Husi-Lang Bridge Hsintien Stream Heng Stream Ru-KouWeir Sonhaia Stream Tahan Stream Taipei Bridge Sanying Bridge

Cheng-Ling Bridge Hsin-Hai Bridge

Flood Diversion Canal

Pacific Ocean TAIWAN Coastal Sea 121° E 24° N Fu-Gi Harbor River Mouth Danshuei Coastal Sea-Pile

Coastal Sea 7 KM 6 5 4 3 2 1 0 Danshuei River

(4)

where S is the salinity, Kh _{and K}v _{are prescribed} non-negative horizontal and vertical diffusivities, 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Þ ¼ p_aðx; y; tÞ þ g½Zðx; y; tÞ z þg Z Z z r r₀ r₀ dz, ð4Þ East Boundary ( Taipower 1 ) Northeast Northwest Boundary Boundary West Boundary ( Zhu-Wei) Danshuei River Keelung River Tahan Stream Hsintien Stream

Fig. 2. (a) An unstructured model grid representing the modeling domain and (b) contour of bottom topography in the Danshuei River estuarine system.

(5)

where p_a is the atmospheric pressure, g is the gravitational acceleration, r is the water density, r₀ is the freshwater density.

The horizontal momentum equations can also be written as Du Dtfv ¼ qpa qx g qZ qxg q qx Z Z z r r0 r dB þvh q 2_u qx2þ q2u qy2 þ q qz v v qu qz , ð5Þ Dv Dtþfu ¼ qp_a qy g qZ qyg q qy Z Z z r r₀ r dB þvh q 2 v qx2þ q2v qy2 þ q qz v v qv qz , ð6Þ

where D( )/Dt is the substantive derivative, f is the Coriolis parameter; vh _{and v}v _{are, respectively, the} coefﬁcients of horizontal and vertical eddy viscos-ities.

The vertical momentum is expressed as Dw Dt ¼v h q2w qx2þ q2w qy2 þ q qz v vqw qz . (7)

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

r ¼ r₀ð1 þ kSÞ, (8)

where k is constant ð¼ 7:8 104ppt1_Þ. 3.2. Turbulence closure model

The vertical turbulence diffusion coefficients (vv and KvÞ strongly affect the flow field, salinity distribution and vertical stratification. Their values can vary over several orders of magnitude at a fixed point in an estuary or coastal sea within a tidal cycle (Odd and Roger, 1978).

Here the mixing length concept is used for vv_and Kv. For ﬂow with two parallel plane boundaries in a wide channel of depth h,Rossby and Montgomery (1935) proposed the following form for vv and Kv based on mixing length theory:

vv¼aZ2 1 Z h 2 qu qz fMðRiÞ, (9) Kv¼aZ2 1 Z h 2 qu qz fsðRiÞ, (10)

where Z is the distance from the surface, a is a constant to be determined empirically and f_M and f_s are the stability functions for momentum and mass, respectively, 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 r qr qz _qu qz 2 . (11)

The stability functions in Eqs. (9) and (10) account for the inhibition of the vertical exchange of momentum and mass (salt) by a stable density

struc-ture. Many forms of f_M and f_s have been

proposed: for example, see Munk and Anderson (1948). Although it is clear that the effect of the stratiﬁcation on the vertical turbulent exchange is a function of Ri, there is no theoretical derivation for these functions (Perrells and Karelse, 1981).

The various formulations for f_M and f_s suggest the following general forms (Blumberg, 1986):

f_M ¼ ð1 þ b_MRiÞqM, (12)

fs¼ ð1 þ bsRiÞqs, (13)

where bM; bs, qM and qs are constants to be determined empirically. In the numerical modeling,

these constants can be determined through

model calibration or taken from previous studies. In this investigation, qM ¼ 12 and qs¼ 32, which indicates that a stable density gradient reduces the vertical turbulent exchange of mass more than that of momentum. The constants b_M and b_s are evaluated through comparison of model results with ﬁeld measurements of a conservative substance such as salt. Furthermore, there is no particular reason that b_M is different from b_s. Since the difference in the turbulence mixing characters for momentum and salt is already included in the constants, q_M and q_s, it is justiﬁable to assume that bM equals bs.

The resulting formations for vvand Kv are vv¼aZ2 1 Z h 2 qu qz ð1 þ bRiÞ1=2, (14) Kv¼aZ2 1 Z h 2 qu qz ð1 þ bRiÞ3=2. (15) The horizontal mixing coefﬁcients (vh _{and K}h_Þ range from 1 to 100 m2_{=s (}_{Bowden and Hamilton,}

(6)

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 values for vh_{and K}h_{are used and they are adjusted} within the range of 1–100 m2_{=s through model} calibration.

3.3. Numerical approximation

A semi-implicit finite-difference scheme intro-duced in TRIM3D by Casulli and Cheng (1992) is used to obtain an efficient numerical algorithm whose stability is independent of free-surface gravity waves, vertical viscosity and bottom friction. The numerical algorithm of UnTRIM is fundamen-tally the same as TRIM3D (Casulli and Cheng, 1992; Casulli and Cattani, 1994), except that the finite-difference treatment of the governing partial differential equations is performed over an unstruc-tured grid mesh. The horizontal model domain is covered by a set of non-overlapping convex poly-gons. Each side of a polygon is either a boundary line or a side of 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 (see Fig. 3). One such grid is called an unstructured orthogonal grid (Casulli and Zanolli, 2002; Casulli and Walters, 2000). The center of a polygon does not necessarily coincide with its geometrical center. The special cases of unstructured orthogonal grids include, of

course, the rectangular ﬁnite-difference grids as well as a grid of uniform equilateral triangles. In these particular cases, the center of each polygon can be identiﬁed with its geometrical center.

When the numerical scheme is extended to an unstructured computational grid, (shown inFig. 3), the momentum equations written in the normal direction to each face are ﬁnite-differenced in the normal direction of each vertical face along oa-, ob-and oc-directions. The momentum equations relate the gradient of water surface elevation between adjoining polygons and the face velocity (velocity normal to the face) on the common face between these polygons. The vertical mixing and bottom friction are discretized implicitly for numerical stability (Cheng and Casulli, 2001). An explicit ﬁnite-difference operator is used to account for the contributions from the discretization of the advec-tion and horizontal dispersion terms. A particular form for this operator can be given in several ways, such as an Eulerian–Lagrangian scheme (Casulli and Cheng, 1992). For stability, the implicitness factor y has to be chosen in the range 1

2pyp1

(Casulli and Cattani, 1994). Along the vertical direction, the model domain is discretized by fixed levels, not necessarily 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 computational cells are accomplished.

The free-surface equation, Eq. (2), is discretized implicitly by the y-method (Casulli and Cattani, 1994; Casulli and Walters, 2000), and only the face-normal velocities are needed to complete the finite-volume integration for the balance of total finite-volume in the polygon. For each polygon, by substituting the finite-differenced momentum equations on all faces into the continuity equation (finite-volume method), the resultant matrix equation governs the water surface elevation distribution for the entire domain. This matrix equation is strongly diagonally dominant, symmetric and positive definite; thus its unique solution can be efficiently determined by preconditioned conjugate gradient iterations until the residual norm becomes smaller than a given tolerance (Golub and van Loan, 1996). Once the free surface for the next time level has been

Fig. 3. Orthogonal unstructured grid (oa, ob and oc present the intersection in the normal direction of each vertical face).

(7)

calculated, the normal velocities on the faces of the prisms are calculated by back substitution into the momentum equations. The details of the finite-difference equations can be found in Casulli and Walters (2000)and Casulli and Zanolli (2002) and thus are not presented here. In summary, by judiciously choosing certain terms in the continuity and momentum equations and by substituting the finite-differenced momentum equations into the continuity equation, a large and sparse matrix equation arises for the water surface elevation for each of the polygons. In this treatment, the baroclinic (salinity/density) terms are treated ex-plicitly, and the resulting matrix is diagonally dominant, positive definite. Once the water surface elevations are solved, they are substituted into the momentum equations to compute the velocity values and the loop of one time-step solution is completed.

For 3-D barotropic flow, the salt transport equation is uncoupled from the momentum equa-tions. A non-negative bottom friction coefficient is specified by the Manning–Chezy formula. If the advection terms are treated by an Eulerian–Lagrangian scheme with positive bottom friction, then the numerical scheme is unconditionally stable. For baroclinic flows, the density gradient terms are expressed explicitly in the momentum equations, and the solutions for the transport variables lagged by one time step are solved. In this case, the numerical scheme is subject to a weak stability condition due to the explicit treatment of the density gradient terms in the momentum 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. In addition, the numerical scheme is also subject to a weak stability condition due to the explicit treat-ment of the horizontal diffusion in the motreat-mentum equations.

3.4. Model implementation

The bottom topography is probably the most important factor that affects the ﬂow properties in environmental modeling. Hence, an accurate representation of bottom topography by the model grid is the most important and fundamental requirement in a successful modeling study (Cheng et al., 1991, 1993). This is particularly true for the Danshuei River estuarine system and adjacent coastal sea, where the bottom topographic

varia-tions are very complex. The model grid must represent accurately the characteristics of the Danshuei River estuarine system that constitutes the model domain.

Because the UnTRIM model uses an unstruc-tured grid, some aspects of the grid definition file are similar to those used in finite-element applications; therefore, the literature on grid generation for finite elements applies in UnTRIM applications. For the purpose of generating an unstructured model grid for UnTRIM model application, a commercial software for mesh generation, Argus software (Bradbury, 1997), has been adopted.

The original bottom topography data in the coastal sea and Danshuei River estuarine system (Fig. 2(b)) used in this study were obtained from the National Center for Ocean Research and Water Resources Agency, Taiwan, respectively. The dee-pest 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 estuarine system and its adjacent coastal sea consists of 8181 3- and 4-sided polygons (Fig. 2(a)). Higher resolution grids are used in the Danshuei River estuary and coarser grids are used in the coastal sea. Sixty vertical layers are speciﬁed with the layer thickness varying from 1 m in the top 10 layers to 2 m for the other layers. For this model grid, a 360 s time step was used in simulations without any sign of numerical stability.

4. Model calibration and veriﬁcation

The numerical model was first calibrated and verified to ascertain the model’s validity before practical applications. A large set of accurate field data are required to calibrate the numerical model and to verify its capability to predict the water surface elevation, current and salinity.

4.1. Preliminary calibration

Friction dissipates energy as the tidal wave propagates landward from the open sea. The tidal range variation along the Danshuei River estuarine system is determined by geometry of the Danshuei River estuarine system and coastal sea, tidal range speciﬁed at the open sea boundaries and the values of bottom roughness. Therefore, the properties of tidal range as a function of distance from the Danshuei River mouth were used to calibrate the bottom friction for the model.

(8)

A single tidal constituent was used as the forcing function at the open sea boundaries for preliminary model calibration. The use of a single constituent expedites the extraction of tidal ranges from the model outputs for comparison with measured data. Because an M2 tide is dominant in the coastal sea and river system, it was used for model calibration. The measured mean tidal ranges at the Taipower 1 and Zhu-Wei station (Fig. 2(a)) are 0.735 and 1.69 m, respectively. Halves of these values were used as the amplitudes of the M2 tide to set the boundary conditions at the two shoreline corners of the open sea boundary. According to the report by Jan et al. (2001), the amplitudes of the M2 tide gradually increase from the shoreline seaward along the west and east boundaries. The amplitudes of the M2tide were interpolated from the east (0.368 m) to northeast boundaries (0.58 m) and from the west (0.845 m) to northwest boundaries (1.95 m) (Fig. 2(a)). The amplitude along the northern boundary of the sea was interpolated between the two end points. The upstream boundary conditions at three tributaries (i.e. Tahan Stream, Hsintien Stream and Keelung River) were given as the mean annual discharges of year 2000. The model was run for 50 tidal cycles to reach a dynamic equilibrium.

The difference between the maximum and minimum surface elevation in the last tidal cycle was considered as the mean tidal range for each grid cell along the axis of Danshuei River estuary. The measured and modeled tidal range distributions along the axes of the estuary and its tributaries are compared (Fig. 4). The tidal range is affected by the river geometry as well as bottom friction. The friction dissipates tidal energy and reduces tidal range as the tide propagates upriver. At the upper-most reach of all three branches the tidal range decreases sharply when the river depth decreases rapidly. Finally the tidal wave reaches its limit when the river bottom rises above high tidal level. The tidal range distribution was satisfactorily predicted by the numerical model (Fig. 4), the average absolute value of differences (mean absolute errors) between computed and measured tidal ranges is 4.0 cm, while the root-mean-square error is 2.0 cm.

Noteworthy is that half of mean tidal range, not the amplitude of the M2 tide from harmonic analysis, was used as the amplitude of the M2 tide to force the model calibration. In the case of the Danshuei River estuarine system and adjacent coastal sea, the other tidal constituents account for a signiﬁcant amount of tidal energy. If the M2

0 10

Distance from Danshuei River mouth (km) 0 0.5 1 1.5 2 2.5 3 Tidal Range (m) D Simulation H Simulation K Simulation D Measured H Measured K Measured D H K D and K confluence D and H confluence 40 20 30

(9)

tide amplitude obtained from harmonic analysis was used for model calibration, it would under-estimate the tidal energy and result in

overestima-tion of bottom fricoverestima-tion coefﬁcient. This result would cause problems when a multiple constituent tide was used to ﬁne-tune the model.

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 (degree) Amplitude (m) Phase (degree) Amplitude (m) Phase (degree) Amplitude (m) Phase (degree) 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 Table 2

The comparison of computed and observed amplitudes and phases of tidal constituents Constituent Danshuei coastal sea

pile

Fu-Gi Harbor Danshuei 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 (degree) 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

Constituent Ru-Kou Weir Hsin-Hai Bridge Chung-Cheng Bridge Ta-Chih Bridge AME RMSE Observed Simulated 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 (degree) 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

(10)

08/21 08/22 08/23 08/24 08/25 08/26 08/27 08/28 08/29 08/30 08/31 0 200 400 600 800 1000 1200 1400 Discharge (m 3/s) Tahan Stream Hsintien Stream Keelung R iver

a

b

c

d

e

f

Simulation Measurement River Mouth Simulation

Measurement SimulationMeasurement

Simulation

Measurement SimulationMeasurement

Tu-Ti-Kung-Pi Taipei Bridge

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) Ru-Kou Weir 8/21 8/31 Hsin-Hai Bridge 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) 8/21 8/31 0 24 48 72 96 120 144 168 192 216 240 Time (hr) Time (day) -2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5 Elevation (m) 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) 8/21 8/31 0 24 48 72 96 120 144 168 192 216 240 Time (hr) -2.5 -2 -3 -1.5 -1 -0.5 0 0.5 1 1.5 2 3 2.5 Elevation (m) 8/21 8/31

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

(11)

4.2. Fine-tune calibration

The bottom roughness was determined by fitting the mean tidal range. Even though a single M2tidal constituent can simulate the tidal energy propaga-tion up the estuary, the time varying water surface variation it represents is far from the measured data. To focus on a particular aspect of astronomical tidal forcing and simplify the calibration process, it was proposed that the model be fine-tuned with the astronomical tide only (Hsu et al., 1999). The approach also allows the use of field data acquired in different time periods for model verification.

Five tidal constituents (i.e. M2, S2, N2, K1 and O1) are sufficient to represent the tidal dynamics in the Taiwan Strait (Jan et al., 2001). The five tidal constituents 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 in Table 1. Amplitudes and phases of these five tidal constitu-ents were used to generate time series of water surface elevation along the open boundaries for a 2-month model simulation. Harmonic analysis was performed on the time series of the model simulated water surface elevation at various locations. After a slight local adjustment of the bottom roughness, the amplitudes and phases of harmonic constants between the computed and observed tides at nine locations are compared (Table 2); the absolute mean difference and root-mean-square errors are also presented. The absolute mean difference and

root-mean-square errors of tidal constituents for ampli-tudes are in the range of 0.01 to 0.02 m and for phases are in the range of 2.4 –6:7. The differences in amplitudes and phases are quite small, and the water surface elevations are reasonably validated. 4.3. Verification 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 constituents (given inTable 1). The model was run for a 1-year simulation. The model results of time-series surface elevation, current and salinity were verified with field data in the same periods when field data were available. Because of the typhoon events and storm runoff, the peak fresh-water discharges occurred on August 23 and 30. The freshwater discharge inputs from three tributaries (Tahan Stream, Hsintien Stream and Keelung River) during the period of August 21–31, 2000 (Fig. 5(a)), and the computed surface elevations and field data at several stations (Fig. 5(b)–(h)) are compared. In general, the modeling results repro-duce the water level variations very well. The water surface elevations at the upriver stations (Fig. 5(f)–(h)) have much more conspicuous re-sponse to pulses of high freshwater discharge, the modeled water levels follow the low water variations closely throughout the storm event. The comparison

Ta-Chih Bridge Simulation Measurement Simulation Measurement Chung-Cheng Bridge 0 24 48 72 96 120 144 168 192 216 240 Time (hr) -2.5 -2 -3 -1.5 -1 -0.5 0 0.5 1 1.5 2 3 2.5 Elevation (m) 8/21 8/31 0 24 48 72 96 120 144 168 192 216 240 Time (hr) -2.5 -2 -3 -1.5 -1 -0.5 0 0.5 1 1.5 2 3 2.5 Elevation (m) 8/21 8/31 Fig. 5. (Continued)

(12)

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 mea-surements 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. The velocity data were recorded by hand and ‘‘ebb or flood’’ directions were also noted based on visual observa-tions. During slack tides, there were some difficul-ties (uncertainty) assigning current directions. Except 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 five stations:

0 4 8 12 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)

Simulation (Top layer) Simulation (Bottom layer) Measurement (Top layer) Measurement (Bottom layer) Kuan-Du Bridge ebb flood Chung-Cheng Bridge ebb flood -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 4 8 12 16 20 24 Time (hr) 0 4 8 12 16 20 24 Time (hr) 0 4 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) -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

a

b

c

d

16 20 24

Fig. 6. 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.

(13)

(a) Kuan-Du Bridge, (b) Taipei Bridge, (c) Hsin-Hai Bridge, (d) Chung-Cheng Bridge and (e) Pa-Ling Bridge,Fig. 6(a)–(e), for May 5, 2000. While there is a certain amount of uncertainties in these measure-ments (due to instrument limitations), these com-parisons 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 circulation and modifies mixing processes. Only limited long-term salinity time series exists at the Zhu-Wei mangrove zone collected by the Taiwan Industrial Technology Institute; this time series was used for model calibration and verification. In the numerical model, the salinity values at open boundaries at the coastal sea were set to 35 ppt. The salinity boundary conditions at the heads of the three tributaries were set to 0 ppt along with the specification of daily freshwater discharges. The wind condition was considered in the model simulation (Fig. 7(b) and (e)). Two 1-week time periods were examined.

During November 15–November 22, 2000, fresh-water discharge was decreasing from 300 to about 150 m3_{=s. The simulated salinity time series} com-pared very favorably with the discrete salinity measurements at the Zhu-Wei mangrove zone, Fig. 7(c). The simulated salinities reproduce the pattern of observed salinity variations in a dynamic variation of salinity between 0.0 and 33 ppt over a tidal cycle. The dynamic variations of salinity covered a range of about 30 ppt, and the mean salinity increased in response to a decrease in freshwater discharge. The salinity was measured at 1.5 m below surface and the simulated salinity was taken from the surface layer, because salinity is generally well mixed in the vertical at this site. In these comparisons, the constants, a ¼ 0:0115 and b ¼ 0:75, gave best results in the turbulence closure model. In the week between December 15 and 22, 2000, freshwater discharge was in an increasing trend from 150 to 300 m3_{=s. While the measured} salinity was as high as 37 ppt suggesting that substantial evaporation might be present but the evaporation mechanism was not built into the model. Nonetheless, the UnTRIM model reﬂected the large dynamic variations of salinity over a tidal cycle, and decreasing trend of mean salinity as freshwater increased near the end of this week (Fig. 7(e)). Fig. 7(c) and (f) also presents the simulated salinity in the absence of wind effect. The comparison of simulated salinity with and without wind conditions shows that there is little difference in salinity at the Zhu-Wei mangrove zone.

5. Model investigations

5.1. Residual current and tidally averaged salinity In most estuaries, residual currents which are masked over by the visible oscillatory tidal current exist due to nonlinear interactions of tidal currents. Gravitational (density) currents are formed if a significant longitudinal salinity gradient exists. The typical estuarine circulation (gravitational current) is characterized by the upriver movement of more saline water in the lower layer and the downriver movement of less saline water in the upper layer (Pritchard, 1956). To remove the tidal and higher frequency fluctuations, a low-pass filter with a cut-off frequency of 36 h1 (Thompson, 1983) was applied to the simulated velocity and salinity time series at three depth levels at Kuan-Du Bridge,

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

Pa-Ling Bridge ebb flood 0 4 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) Fig. 6. (Continued)

(14)

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

Zhu-Wei mangrove zone

Simulation ( wind effect ) Simulation ( no wind ) Measurement

Zhu-Wei mangrove zone

Simulation ( wind effect ) Simulation ( no wind ) Measurement

a

_b

c

d

e

_f

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

(15)

Taipei Bridge and Pa-Ling Bridge (only top and bottom levels) for June 1–September 30, 2000.Fig. 8 depicts the freshwater hydrograph, the residual velocities and tidally averaged salinities. The re-sidual current structure at the Kuan-Du station shows the typical estuarine circulation with the downriver pointing velocity near the surface and upriver pointing velocity near the bottom for most of the time. The strength of the estuarine circulation

was variable and mainly influenced by freshwater discharge. When the freshwater discharge was below 100 m3_{=s, the two-layered flow is persistent} and the water column was stratified with 2–5 ppt difference from top to bottom, Fig. 8(b). In this period, a few runoff events took place, the river discharges increased to 500 m3_{=s or above (for} example, June 13, June 20, July 10, August 23, August 30 and September 10), the salinity was

0 500 1000 1500 2000 2500 Time (hr) -0.4 0 0.4 0.8 1.2 1.6 Residual velocity(m/s) Top layer Middle layer Bottom layer Discharge 0 500 1000 1500 2000 2500 Di s c har ge ( m 3 /s )

Kuan-Du Bridge Kuan-Du Bridge

6/16/13 6/20 7/10 8/23 8/30 9/10 9/30 0 500 1000 1500 2000 2500 Time (hr) 0 6 12 18 24 30 36 42 48 54 60 Salinity (ppt) 0 500 1000 1500 2000 2500Di s c har ge ( m 3 /s ) Top layer Middle layer Bottom layer Discharge 6/16/136/20 7/10 8/238/309/10 9/30 0 500 1000 1500 2000 2500 Time (hr) -0.7 0 0.7 1.4 2.1 2.8 3.5 Residual velocity(m/s) Top layer Middle layer Bottom layer Discharge 0 400 800 1200 1600 2000 Dis char ge ( m 3 /s ) 6/16/13 6/20 7/10 8/23 8/30 9/10 9/30 0 500 1000 1500 2000 2500 Time (hr) 0 5 10 15 20 25 30 Salinity (ppt) 0 400 800 1200 1600 2000 D ischar ge (m 3 /s ) Top layer Middle layer Bottom layer Discharge Taipei Bridge Taipei Bridge 6/16/13 6/20 7/10 8/23 8/30 9/10 9/30

a

b

c

d

Fig. 8. Temporal variation of (a) residual velocity and (b) tidally averaged salinity at Kuan-Du Bridge, (c) residual velocity and (d) tidally averaged salinity at Taipei Bridge and (e) residual velocity and (f) tidally averaged salinity at Pa-Ling Bridge.

(16)

severely depressed,Fig. 8(b), and consequently the estuarine circulation disappeared changing to downstream pointing flow for the entire water column, Fig. 8(a). It is particularly clear, between August 22 and 30, 2000, that the salinity values were near zero from top to bottom, the estuarine circulation disappeared. After August 30, the river discharge decreased to about 250 m3_{/s for a few} days, and consequently the estuarine circulation developed that led to the restoring of stratification. The time required to change from stratified two-layered flow to freshwater in the entire water column is very short (ca. 2 days), which implies that the residual time of the Danshuei River is very short.

The characteristic two-layered circulation was not developed at Taipei Bridge and Pa-Ling Bridge, the net residual currents are all downriver because the water depths at these sites were relatively shallow and the baroclinic forcing (longitudinal salinity gradient) fails to overcome the barotropic forcing (surface slope). Consequently, the full estuarine circulation did not develop at these sites. The vertical stratiﬁcation did not develop because of shallow water depth and turbulent mixing, Fig. 8(c)–(f). When the river discharge was low, weaker downstream velocities are seen near the surface, the velocities in the bottom layer decreased to near zero attempting to move up the river.

5.2. Effect of freshwater discharge

The calibrated and verified 3-D hydrodynamic model of Danshuei River estuary was used to investigate the residual current and salinity re-sponses to variations of freshwater discharge con-ditions. Five tidal constituents, M2, S2, N2, K1 and O1, were used to generate a time series of water surface elevations as open boundary conditions and a constant salinity of 35 ppt was specified at the coastal sea for a 2-month simulation. For low-flow conditions, the model was run with freshwater discharges Q₇₅ (a flow with an exceedence prob-ability of 0.75) of 4.79, 25.32 and 5:0 m3_{=s at the} upstream boundaries of the Tahan Stream, Hsintien Stream and Keelung River, respectively. The averaged residual velocities over two spring–neap cycles, i.e. 29 days, indicate a general upriver residual flow along the deep channel in the lower portions of the Danshuei River estuary,Fig. 9. The average salinity distribution shows that the limits of salt intrusion (i.e. 1 ppt isohaline) are located at the distance of about 25, 4 and 9 km from the Danshuei River mouth, Hsintien River mouth and Keelung River mouth, respectively,Fig. 10.

A moderate discharge condition was simulated using historical mean freshwater discharges of 59.2, 87.1 and 32:8 m3_{=s at the upstream boundaries of} the Tahan Stream, Hsintien Stream and Keelung River, respectively. All other conditions for tides

0 500 1000 1500 2000 2500 Time (hr) Time (hr) -0.2 0 0.2 0.4 0.6 0.8 1 Residual velocity (m/s) 0 40 80 120 160 Di s c har ge ( m 3/s ) Top layer Bottom layer Discharge

Pa-Ling Bridge Pa-Ling Bridge

6/16/13 6/20 7/10 8/23 8/30 9/10 9/30 0 500 1000 1500 2000 2500 0 3 6 9 12 15 18 21 24 27 30 Salinity (ppt) 0 40 80 120 160 D is c h a rg e (m 3/s ) Top layer Bottom layer Discharge 6/16/13 6/20 7/10 8/23 8/30 9/10 9/30

e

f

Fig. 8. (Continued)

(17)

and salinity at open boundaries were kept the same as the low-ﬂow condition. The simulated residual velocities and salinities are depicted (Figs. 11 and12, respectively). Increased freshwater discharge pushes the limits of salt intrusion downriver toward river mouth. In moderate discharge conditions, the limits of salt intrusion are located 11 and 3 km in the Danshuei River and Keelung River, respectively. There is no salt water intrusion in the Hsintien Stream. The higher freshwater discharge prevents

salinity transport upriver to Hsintien Stream. A two-layered circulation exists near Kuan-Du Bridge of the Danshuei River–Tahan Stream and the lower portions of the tributaries under the low-flow condition (Q₇₅Þ, Fig. 9. Under the mean discharge condition, only at the deeper, lower river station at Kuan-Du shows a two-layered circulation, with a weaker upriver flow than that under the low-flow condition, Fig. 11. Comparing between Figs. 10(a) and 12(a), the salinity is high

0 5 10 15 20 25 30

Distance from Danshuei Rivermouth (km)

Danshuei River-Tahan Stream

-10 -5 0 Depth (m) 25 Q75 flow River mouth Kuan-Du Bri. Taipei Bri. Hsin-Hai Bri. cm/s 0 2 4 6 8 10 12

Distance from HsintienStream mouth(km) Hsintien Stream -10 -5 0 Depth (m) 30 Q75 flow Chung-Cheng Bri. cm/s 0 5 10 15 20 25

Distance from KeelungRiver mouth (km) Keelung River -10 -5 0 Depth (m) 15 Q75 flow Pa-Ling Bri. Ta-Chih Bri. cm/s

Fig. 9. Simulation of residual velocity averaged over two spring–neap cycles with Q₇₅freshwater discharge: (a) Danshuei River–Tahan Stream, (b) Hsintien Stream and (c) Keelung River.

(18)

and water depth is large near the Kuan-Du station. The baroclinic effects increase with fresh-water discharge at lower ﬂow, while decrease at higher ﬂow.

5.3. Relationship between freshwater discharge and limit of salt intrusion

The Danshuei River (UnTRIM) numerical

model was used to simulate several freshwater discharge scenarios to establish a relationship

between freshwater discharge and the limit of salt water intrusion. The open boundary condi-tions were kept the same as in the previous cases. Based on the predicted limit of salt water intrusion as a function of freshwater discharge from three tributaries (Fig. 13), an exponential least-square regression ﬁt is developed for each tributary, by assuming the limit of salt intru-sion is exponentially (negatively) correlated with the freshwater discharge. The regression correlations yielded the maximum salinity intrusion

0 5 10 15 20 25 30 -10 -5 0 Salinity : ppt Q75 flow River mouth Kuan-Du Bri. Taipei Bri. Hsin-Hai Bri. 0 2 4 6 8 10 12 -10 -5 0 Salinity : ppt Chung-Cheng Bri. 0 5 10 15 20 25 -10 -5 0 Salinity : ppt Pa-Ling Bri. Ta-Chih Bri.

Distance from Danshuei Rivermouth (km) Danshuei River-Tahan Stream

Depth (m)

Distance from HsintienStream mouth(km) Hsintien Stream

Depth (m)

Distance from KeelungRiver mouth (km) Keelung River Depth (m)

a

b

c

Q75 flow Q75 flow

Fig. 10. Simulation of salinity distribution averaged over two spring–neap cycles with Q₇₅ freshwater discharge: (a) Danshuei River–Tahan Stream, (b) Hsintien Stream and (c) Keelung River.

(19)

distance ðY Þ for a given freshwater discharge ðX Þ at each tributary at

Danshuei River– Tahan Stream:

Y ¼ 23:69 eð0:011X Þ; R2 ¼0:97. (16) Hsintien Stream:

Y ¼ 11:42 eð0:054X Þ; R2 ¼0:90. (17)

Keelung River:

Y ¼ 11:98 eð0:046X Þ; R2¼0:90. (18) The high correlation coefﬁcients indicate that freshwater discharge plays dominant role in salt water intrusion in the Danshuei River estuarine system. The regression equations can be used as

-10 -5 0 25 River mouth Mean flow Kuan-Du Bri. Taipei Bri. Hsin-Hai Bri. cm/s 0 2 4 6 8 10 12 -10 -5 0 30 Mean flow Chung-Cheng Bri. cm/s 0 5 10 15 20 25 0 5 10 15 20 25 30 -10 -5 0 15 Mean flow Pa-Ling Bri. Ta-Chih Bri. cm/s

Distance from Danshuei River mouth (km)

Depth (m)

Distance from Hsintien Stream mouth(km) Hsintien Stream

Depth (m)

Distance from KeelungRiver mouth (km) Keelung River

Depth (m)

Fig. 11. Simulation of residual velocity averaged over two spring–neap cycles with mean freshwater discharge: (a) Danshuei River–Tahan Stream, (b) Hsintien Stream and (c) Keelung River.

(20)

simple tools for prediction of the inﬂuence of freshwater discharge on salt intrusion.

6. Conclusions

A 3-D numerical hydrodynamics model, Un-TRIM, was implemented and applied to the Dan-shuei River estuarine system and its adjacent coastal sea in northern Taiwan. The preliminary calibration

of bottom friction uses a single constituent tide M2 to reproduce the variation of mean tidal range along the Danshuei River estuarine system. The model was then calibrated against five major astronomical tides by adjusting the bottom friction coefficient. 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 comparing directly with

0 5 10 15 20 25 30

-10 -5 0 Salinity : ppt Mean flow River mouth Kuan-Du Bri. Taipei Bri. Hsin-Hai Bri. 0 2 4 6 8 10 12 -10 -5 0 Chung-Cheng Bri. Salinity : ppt Mean flow No salt intrusion 0 5 10 15 20 25 -10 -5 0 Salinity : ppt Mean flow Pa-Ling Bri. Ta-Chih Bri.

Distance from Danshuei Rivermouth (km)

Depth (m)

Distance from HsintienStream mouth(km) Hsintien Stream

Depth (m)

Distance from KeelungRiver mouth (km) Keelung River

Depth (m)

a

b

c

Fig. 12. Simulation of salinity distribution averaged over two spring–neap cycles with mean freshwater discharge: (a) Danshuei River–Tahan Stream, (b) Hsintien Stream and (c) Keelung River.

(21)

measured water surface elevation, current and salinity under variable daily freshwater discharges from upriver of the Tahan Stream, Hsintien Stream and Keelung River in 2000.

The verified model was applied to investigate the influence of freshwater discharge on residual circula-tion and salt water intrusion. The numerical model was run for 2 months and the simulated model results were averaged to generate the subtidal velocity and tidally averaged salinity distributions. A characteristic two-layered circulation prevails at the Kuan-Du station under low-flow conditions. The downriver directed barotropic forcing is smaller than the upriver directed baroclinic forcing in deeper portion of the river, resulting in upriver flow along the river bottom. As freshwater discharge increases, the barotropic forcing increases throughout the estuary, while the baroclinic forcing responds differently at various stations. At the downriver station, particularly at Kuan-Du, where the salinity is high and water depth is large, baroclinic effects increase with freshwater discharge at lower flow, while decrease at higher flow. However, the increase in baroclinic forcing with freshwater discharge is still not enough to overcome the increase of barotropic forcing, thus, resulting in weaker gravitational circulation as flow increases. The magnitude of river discharge is shown to be a

dominant factor affecting salinity intrusion in the estuarine system. Under mean ﬂow conditions, no salt water intrudes into the Hsintien Stream. Reduction of the river discharge to Q₇₅ increased the extent of salt intrusion by 14, 4 and 6 km in the Tahan Stream, Hsintien Stream and Keelung River, respectively.

Using the UnTRIM Danshuei River model, a correlation between the freshwater discharges and the distance of the salt intrusion was established. The distance of the salt intrusion decreases expo-nentially as the freshwater inflow increases. The correlation coefficients of the salt intrusion and freshwater inflow regressions at the Danshuei River–Tahan Stream, Hsintien Stream and Keelung River were greater than 0.9, these regression equations are thus suitable to be used as a simple tool to quickly predict the distance of salt intrusion under the influence of freshwater outflows from major tributaries of Danshuei River.

Acknowledgments

The project, under which this study was con-ducted, was supported by the National Science Council, Taiwan, under Grant numbers NSC 93-2211-E-239-006 and 94-2211-E-239-011. The ﬁnan-cial support is highly appreciated. We also thank the

0 20 40 60 80 100 120 140 160 180 200 Discharge (m /s)3 0 5 10 15 20 25 30

Distance of limit of salt intrusion (km)

D : model result s H : model result s K : model result s D : Y = 23 .6 9*e (- 0. 01 1 * X ) R2 = 0.97 H : Y = 11 .4 2*e (- 0. 05 4 * X ) R2 = 0.90 K : Y = 11 .9 8*e (- 0. 04 6 * X ) R2 = 0.90

Fig. 13. The relationship between the distance of limit of salt intrusion and freshwater discharge input (D: Danshuei River–Tahan Stream, H: Hsintien Stream, K: Keelung River).

(22)

Taiwan Water Resources Agency, Industrial Tech-nology Institute and Central Weather Bureau for providing the measured data.

References

Bilgili, A., Proehl, J.A., Lynch, D.R., Smith, K.W., Swift, M.R., 2005. Estuary/ocean exchange tidal mixing in a Gulf of Maine estuary: a Lagrangian modeling study. Estuarine, Coastal and Shelf Science 65 (4), 607–624.

Blumberg, A.F., 1986. Turbulent mixing processes in lakes, reservoirs and impoundments. In: Gray, W.G. (Ed.), Physi-cal-Based Modeling of Lakes, Reservoirs, and Impound-ments. ASCE, New York, pp. 79–104.

Bowden, K.F., Hamilton, P., 1975. Some experiments with a numerical model of circulation and mixing in tidal estuary. Estuarine and Coastal Marine Science 3 (3), 281–301. Bradbury, R., 1997. User’s Guide for Argus One. Argus

Interware, Inc.

Casulli, V., Cattani, E., 1994. Stability, accuracy and efficiency of a semi-implicit method for three-dimensional shallow water flow. Computers and Mathematics with Applications 27, 99–112. Casulli, V., Cheng, R.T., 1992. Semi-implicit finite difference

methods for three-dimensional shallow water ﬂow. Interna-tional Journal for Numerical Methods in Fluids 15, 629–648. Casulli, V., Walters, R.V., 2000. An unstructured grid, three-dimensional model based on the shallow water equations. International Journal for Numerical Methods in Fluids 32, 331–348.

Casulli, V., Zanolli, P., 2002. Semi-implicit modeling of nonhydrostatic free-surface ﬂows for environmental pro-blems. Mathematical and Computer Modelling 36 (9–10), 1131–1149.

Chen, C., Liu, H., 2003. An unstructured grid, ﬁnite-volume, three-dimensional, primitive equations ocean model: applica-tion to coastal and estuaries. Journal of Atmospheric and Oceanic Technology 20, 159–186.

Chen, X., 2003. An efﬁcient ﬁnite difference scheme for simulating hydrodynamic in narrow rivers and estuaries. International Journal for Numerical Methods in Fluids 42 (3), 233–247.

Cheng, R.T., Casulli, V., 2001. Evaluation of the UnTRIM model for 3-D tidal circulation. In: Proceedings of the Seventh International Conference on Estuarine and Coastal Modeling, Florida, pp. 628–642.

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

Cheng, R.T., Garnter, J.W., Casulli, V., 1993. Tidal, residual, intertidal mudﬂat (TRIM) model and its applications to San Francisco Bay, California. Estuarine, Coastal, and Shelf Science 36 (3), 235–280.

Dyer, K.R., 1973. Estuaries: a Physical Introduction. Wiley, New York, 140pp.

Golub, G.H., van Loan, C.F., 1996. Matrix Computations, third ed. J. Hopkins, London.

Horrevoets, A.C., Savenije, H.H.G., Schuurman, J.N., Graas, S., 2004. The inﬂuence of river discharge on tidal damping in alluvial estuaries. Journal of Hydrology 294, 213–228.

Hsu, M.H., Kuo, A.Y., Kuo, J.T., Liu, W.C., 1999. Procedure to calibrate and verify numerical models of estuarine hydro-dynamics. ASCE Journal of Hydraulic Engineering 125 (2), 166–182.

Hsu, M.H., Fu, J.C., Liu, W.C., 2003. Flood routing with real-time stage correction method for ﬂash ﬂood forecasting in the Tanshui River, Taiwan. Journal of Hydrology 283, 267–280. Jan, S., Wang, Y.H., Chao, S.Y., Wang, D.P., 2001. Develop-ment of a nowcast system for the Taiwan Strait (TSNOW): numerical simulation of barotropic tides. Ocean and Polar Research 23 (3), 195–203.

Ji, Z.G., Morton, M.R., Hamrick, J.M., 2001. Wetting and drying simulation of estuarine processes. Estuarine, Coastal and Shelf Science 53 (5), 683–700.

Kuo, A.Y., Liu, W.C., 2001. Investigation of hydrodynamic characteristics in the Tanshui River estuarine system. In: The 2001 Agricultural Engineering Annual Conference, Depart-ment of Agricultural Engineering, National Taiwan Univer-sity, Taipei, Taiwan, pp. 1–13 (in Chinese).

Liu, W.C., Hsu, M.H., Kuo, A.Y., 2001. Investigation of long-term transport in Tanshui River estuary, Taiwan. ASCE Journal of Waterways, Port, Coastal, and Ocean Engineering 127 (2), 61–71.

Liu, W.C., Hsu, M.H., Wu, C.R., Wang, C.F., Kuo, A.Y., 2004. Modeling salt water intrusion in Tanshui River estuarine system—case-study contrasting now and then. ASCE Journal of Hydraulic Engineering 130 (9), 849–859.

Munk, W.H., Anderson, E.R., 1948. Notes on a theory of the thermoclinic. Journal of Marine Research 7 (3), 276–295. Odd, N.V.M., Roger, J.G., 1978. Vertical mixing in stratiﬁed

tidal ﬂows. ASCE Journal of Hydraulics Division 104 (HY3), 337–351.

Oliveira, A., Fortunato, A.B., Rego, J.R.L., 2006. Effect of morphological changes on the hydrodynamics and ﬂushing properties of the Obidos lagoon (Portugal). Continental Shelf Research 26 (8), 917–942.

Perrells, P.A.J., Karelse, M., 1981. A two-dimensional laterally averaged model for salt intrusion in estuaries. In: Fisher, H.B. (Ed.), Transport Models for Inland and Coastal Waters. Academic Press, San Diego, pp. 483–535.

Pritchard, D.W., 1956. The dynamic structure of a coastal plain estuary. Journal of Marine Research 15, 33–42.

Rossby, C.G., Montgomery, R.B., 1935. The layer of frictional inﬂuence in wind and ocean currents. Physical Oceanography and Meteorology 3 (3), 1–10.

Sankaranarayanan, S., McCay, D.F., 2003. Three-dimensional modeling of tidal circulation in Bay of Fundy. ASCE Journal of Waterways, Port, Coastal, and Ocean Engineering 129 (3), 114–123.

Sankaranarayanan, S., 2005. A 3D boundary-ﬁtted barotropic hydrodynamic model for the New York Harbor region. Continental Shelf Research 25, 2233–2260.

Thompson, R.O.R.Y., 1983. Low-pass ﬁlters to suppress inertial and tidal frequencies. Journal of Physical Oceanography 13, 1077–1083.

Wang, Y.H., 2002. Mapping ﬂow using towed-ADCP in coastal water of Taiwan. In: Proceedings of the 24th Ocean Engineering Conference in Taiwan, pp. 485–490 (in Chinese). Zhong, L., Li, M., 2006. Tidal energy ﬂuxes and dissipation in the Chesapeake Bay. Continental Shelf Research 26 (6), 752–770.