- Email: [email protected]
3-D fluid-mud dynamics in the Jiaojiang Estuary, China
Estuarine, Coastal and Shelf Science 65 (2005) 747e762 www.elsevier.com/locate/ecss
3-D fluid-mud dynamics in the Jiaojiang Estuary, China W.B. Guan a,b,*, S.C. Kot b, E. Wolanski c a
Laboratory of Ocean Dynamic Processes and Satellite Oceanography, Second Institute of Oceanography, State Oceanic Administration, Hangzhou 310012, China b Department of Mechanical Engineering, Hong Kong University of Science and Technology, Clear Water Bay, Hong Kong, PR China c Australian Institute of Marine Science, P.M.B. No. 3, Townsville M.C., Queensland 4810, Australia Received 10 March 2004; accepted 27 May 2005 Available online 9 November 2005
Abstract A 3-D model has been developed for the muddy Jiaojiang Estuary and adjoining coastal waters, and verified against field observations. To simulate fluid-mud formation, the model uses a fine resolution grid near the bottom and involves coupling processes between hydrodynamics and fluid mud such as the sediment-induced buoyancy, increasing turbulent kinetic energy sink and kinematic viscosity, mixing by internal waves riding on the lutocline, and non-Newtonian properties of fluid mud. The effective hydrodynamic drag was reduced in the presence of fluid mud. It is shown that the estuary is infilled by tidal pumping and that longitudinal and transversal gradients of suspended sediment concentration, salinity, and currents control the formation of mud banks. Thus a 3-D model is necessary to estimate the fate of mud, although the model results are very sensitive to details of the parameterization of the hydrodynamics-mud feedback processes. Ó 2005 Elsevier Ltd. All rights reserved. Keywords: cohesive sediment; model; turbidity maximum; lutocline; fluid mud; China
1. Introduction The dynamics of fine-grained sediment is of major importance in many coastal areas and estuaries, because this sediment generates siltation and degrades the water quality. Fine-grained sediment (mud) is composed mainly of clay (particle size ! 2 mm) and silt (particle size ! 63 mm). Fluid mud is formed at suspended sediment concentrations greater than 5 g l1 and is the major mechanism of the transport of cohesive sediments in coastal zones and estuaries (Faas, 1984; Odd et al., 1992). Estuarine models that include fluid mud flows seldom deal with the internal structure of fluid mud layers; * Corresponding author. Laboratory of Ocean Dynamic Processes and Satellite Oceanography, Second Institute of Oceanography, State Oceanic Administration, Hangzhou 310012, China. E-mail address: [email protected] (W.B. Guan). 0272-7714/$ - see front matter Ó 2005 Elsevier Ltd. All rights reserved. doi:10.1016/j.ecss.2005.05.017
instead they assume a sharp interface between the mud layer and the water above, with semi-empirical exchange formulations between them for entrainment and settling. For example, Odd and Cooper (1989) developed a depth-averaged, 2-D fluid mud model to describe the formation, movement, re-entrainment and dewatering of fluid mud layers, and the fluid mud was assumed to be a Bingham plastic fluid with a uniform density. Roberts (1993) replaced the Bingham plastic assumption with a high viscosity Newtonian fluid model. Le Normant et al. (1998) coupled a horizontal 2-D fluid mud model to a 3-D model for simulating the turbidity maximum in the Loire Estuary. This model assumed a sharp discontinuity at water/fluid mud interface and no vertical velocity gradient within the fluid mud layer, and it failed to reproduce some observations (Le Hir et al., 2001). In fluid mud flows the suspended sediment concentration, C, is so high that interactions between particles
748
W.B. Guan et al. / Estuarine, Coastal and Shelf Science 65 (2005) 747e762
and water become strong, and a description of its mechanics needs a multiphase approach. This was achieved by Teisson et al. (1992) who modeled the fluid mud layer by using a two-phase flow model. This approach also enables a combination of hydrodynamics and soil mechanics, for which pore-water movement relative to settling particles is considered. However, the computations are very time-consuming; also, unsteady configurations cannot be readily simulated. In many cases of natural fine sediment transport, particles and water have similar horizontal velocity components, so that a unique set of horizontal momentum equations can be used for the water/ sediment mixture. For example, Le Hir (1997) presented an integrated modeling approach, simulating water- and mud dynamics by solving for the mass conservation and momentum equations over the entire column. In their vertical 2-D modeling study of the fluid-mud dynamics and the turbidity maximum in the Gironde Estuary, Le Hir et al. (2001) replaced the integrated formulations by a continuous one, which is probably more representative of a progressive transition between the water column and bottom sediment. Le Hir and Cayocca (2002) adopted the concept of continuous modeling to simulate 3-D turbidity currents on a slope. In their model, a zlevel coordinate system with a stepped bottom topography was chosen, which is perhaps suitable for their simplified problem but not for an actual case, especially a macro-tidal estuary. Based on simple analysis and historical measurements in the Amazon Estuary, Vinzon and Mehta (2001) found that it is essential to account for the strong interaction between flow and sediment dynamics within the lower boundary layer for modeling sediment discharge. In this paper these concepts were applied to a finescale, 3-D model of the topographically-complex, muddy, macro-tidal, Jiaojiang Estuary in China for which an extensive data set on the dynamics of water and suspended sediment was available. It is shown that the model can fairly reliably predict cohesive sediment dynamics, but that the results are very sensitive to details of the parameterization of these various processes. The model thus points to the need for processbased studies of these processes. The extremely turbid Jiaojiang Estuary (Fig. 1) is located about 200 km south of the Yangtze River (Changjiang River) mouth; it drains the Lingjiang and Yongning rivers. The Yongning River is dammed and has minimal runoff. The Lingjiang River drains a catchment area of 6520 km2. Its annual runoff is 6.66 ! 109 m3, corresponding to a mean discharge of 211 m3 s1 (Li et al., 1993). The Jiaojiang Estuary is about 35 km long. In the estuary, the mean width of the channel is about 1.2 km with a maximum of 1.8 km at the mouth. The estuary faces Taizhou Bay, which is shallow (the depth varies between 1 and 3 m at low
spring tide). Semi-diurnal macro-tides prevail with a mean tidal range of about 4 m and a maximum tidal range of 6.3 m. The vertically-averaged tidal current peaks at 2.0 m s1. At Haimen, the flood tide lasts 5.1 h while the ebb tide lasts 7.3 h (Dong et al., 1997). Off Haimen, the maximum flood and ebb currents are, respectively, 2.1 and 1.8 m s1 (Zhu, 1986). The ratio of the freshwater volume to the tidal volume is about 0.04, but more than 0.1 during periods of very high river discharge (Bi and Sun, 1984). A turbidity maximum zone is present where C peaks at 40 g l1 (Bi and Sun, 1984; Fu and Bi, 1989). These authors estimated that the fluvial sediment inflow is about 1.2 ! 106 tons year1, corresponding to a mean sediment concentration of 0.18 g l1. Most fine-grained sediment found in the estuary comes from coastal waters, and originates from the Yangtze River (Bi and Sun, 1984; Zhu, 1986). Siltation in the estuary is about 0.2 m year1, only 5% of which is due to fluvial material (Bi and Sun, 1984; Fu and Bi, 1989; Li et al., 1993; Dong et al., 1997; Guan et al., 1998).
2. Methods 2.1. Field observations An Inter-Ocean S4 electromagnetic current meter (Courtesy, Mark Geneau, Inter-Ocean) was deployed 1.75 m above the bed at mooring station M1 (see Fig. 1 for a location map) during April 12e18, 1991. At station C1 also three backscatterance nephelometers were deployed at elevations of 0.35, 0.55 and 1.67 m above the bed. These three sensors recorded data at an interval of 5 min. The nephelometers were calibrated in situ to yield C (Li et al., 1993; Dong et al., 1997; Guan et al., 1998). At stations C1 and C2, 25-h-long observational stations were maintained at spring tides and neap tides, during which hourly, vertical profiles of salinity, C and currents were obtained using the ship-born CTD cumnephelometer of Wolanski et al. (1988) and a rotor-type current meter (model SLC9-1, made by the Institute of Marine Instrument, Ocean University of China, Qingdao). Longitudinal and transversal transects of temperature, salinity and C were also obtained using the shipborn CTD. 2.2. 3-D model The generalized s-coordinate Princeton Ocean Model (POM) was used (Blumberg and Mellor, 1987; Mellor et al., 2001) to simulate the hydrodynamics in the estuary, and further a sediment transport module was developed to integrate with the hydrodynamics model. The incompressible, hydrostatic and Boussinesq approximations were used. The wateresediment mixture
W.B. Guan et al. / Estuarine, Coastal and Shelf Science 65 (2005) 747e762
749
(a) 28.85 N 28.8 Yongquan
0
2
4
6 km
28.75
Latitude (°N)
Lutou
Taizhou Bay
Shixianfu
M1, C2
28.7 Luocun
C1 LaoshuyuSongpuzha
Shanpu Haimen
Jiaojiang Estuary 28.65
T3
Korea
ina
28.6
Ch
Japan Study site
28.55
121.25
121.3
121.35
121.4
121.45
121.5
121.55
121.6
121.65
Longitude (°E)
(b)
Yongquan
Luocun
Songpuzha
Haimen
0
z (m)
-2
-4
-6
-8 0
5
10
15
20
25
30
35
40
45
Distance along transect (km)
(c)
C2
M1
Fig. 1. (a) Orthogonal curvilinear grid used in the simulations. M1 is the mooring stations; C1 and C2 are anchor stations; T3 is the tide stations. The transect connecting all thick dots will be used to show the model results below. The triangles denote locations of towns. (b) Vertical coordinates for the longitudinal transect indicated in Fig. 2. The triangles denote locations of towns. (c) Bathymetry (in m) for the estuarine part from Songpuzha to Yongquan. A dike built on the central shoal of the estuary has been indicated.
750
W.B. Guan et al. / Estuarine, Coastal and Shelf Science 65 (2005) 747e762
was treated as a single-phase fluid in which all particles follow water except for their settling velocity (Winterwerp et al., 2001). Some non-Newtonian effects of highconcentrated mud suspensions were included. The equations were written in a right-handed Cartesian coordinate system, (x, y, z, t), with the z-axis pointing vertically upwards, and t denoting time. u, v and w are the velocity components in the x-, y- and z-directions, respectively, r is the fluid density, r0 is the reference density, g is the gravitational acceleration, AM is the horizontal turbulent eddy viscosity, KM is the vertical turbulent eddy viscosity, KMB is the background vertical viscosity representing the mixing process due to the internal wave and the increment due to the suspended sediment, and f is the Coriolis parameter. The continuity and momentum equations were:
sediment concentration, C, and therefore a background turbulence value was imposed. The equations for salinity, S, and suspended sediment concentration, C, were:
vu vv vw C C Z0; vx vy vz
where ws is the settling velocity of suspended sediment particles (positive downward), AH and KH are the horizontal and vertical turbulent mass diffusivities, respectively, and KHB is the background value of the vertical diffusivity. The leapfrog scheme used in the original POM was replaced by a total-variation-diminishing (TVD) scheme with the superbee limiter as a weighting function between the upwind scheme and either the LaxeWendroff scheme in the horizontal or the central scheme in the vertical for the advective terms of scalar variables (James, 1996; Luyten et al., 1999). The drying and flooding algorithm of Flather and Heaps (1975) was used. Temperature differences are small in the Jiaojiang Estuary and were thus ignored (Guan et al., 1998). The density r (g l1) of the sedimentewater mixture was computed as (Odd, 1988):
vu vu2 vuv vuw C C C fv vt vx vy vz 1 vp v vu CFx C ðKM CKMB Þ Z ; r0 vx vz vz vv vvu vv2 vvw C C C Cfu vt vx vy vz 1 vp v vv CFy C ðKM CKMB Þ ; Z r0 vy vz vz
0Z
1 vp r g; r0 vz r0
ð1Þ
ð2Þ
ð3Þ
ð4Þ
with v vu v vu vv Fx Z C 2AM C AM ; vx vx vy vy vx v vv v vv vu C Fy Z 2AM C AM : vy vy vx vx vy In these equations, a differentiation is made between background turbulence and flow-generated turbulence. This was necessary because strong sediment-induced stratification results in zero or negligible vertical mixing, when the MelloreYamada level 2.5 scheme is used. However, field data reveal that background level mixing exists due to internal waves riding on the lutocline (Jianhua and Wolanski, 1998; Jiang and Mehta, 2000, 2002). This background turbulence affects also the conservation equations for salinity, S, and suspended
vS vuS vvS vwS C C C vt vx vy vz v vS v vS v vS Z AH C AH C ðKH CKHB Þ ; ð5Þ vx vx vy vy vz vz
vC vuC vvC vðw ws ÞC C C C vt vx vy vz v vC v vC v vC Z AH C AH C ðKH CKHB Þ ; ð6Þ vx vx vy vy vz vz
r rw rZrw CC s Z1000C0:78SC0:62C; rs
ð7Þ
where C is in g l1, rw Z 1000 C 0.78S is the water density (g l1), and rs Z 2650 (g l1) is the sediment density. Horizontal turbulent mixing, which describes the effect of sub-grid scale motion, was represented by a Smagorinsky-type formulation (Smagorinsky et al., 1965). Vertical turbulent mixing was calculated using the Mellor and Yamada (1982) turbulence closure model. This model calculates q2, the turbulence kinetic energy, and q2l, where q is the turbulence velocity scale and l is the length scale. The wall proximity function of Blumberg et al. (1992) was used, to account for the presence of free surface and to reproduce data that suggest that l is larger near the free surface than near a solid wall. Mellor (2001) introduced a correction function in the turbulence closure model based on experiments, and thus a critical constant GHc must be
751
W.B. Guan et al. / Estuarine, Coastal and Shelf Science 65 (2005) 747e762
tuned so that model calculations agree more closely with field data. At the bottom, ws CCðKH CKHB Þ
vC ZFD FE ; at zZ Hðx; yÞ; vz
ð12Þ
m
with ð8Þ
where FD is the deposition rate, and FE is the erosion rate of sediment. To calculate FE, the bed was discretized into a number (NLAY Z 10) of elements. The erosion rate equation developed by Ariathurai (1974) from the flume data of Partheniades (1965) was used: k E½t =t 1 for tb OtE FE Z e b E ; ð9Þ 0 for tb %tE where tb is the shear stress at the bottom, tE is the threshold shear stress, E is the erosion constant, and ke (ke Z 1 at the uppermost bed layer) is the remaining ratio of the capability of flow for erosion of the newly exposed layer. The threshold bed shear strength tE was: tE Zas Cbs s ;
for C%Chind max k1 Cm1 ; wsf m ; ws Z ð1C=Cgel Þ 2 ð1C=rs Þ wsr for Chind !C!Cgel 1Ca Cbm (
ð10Þ
in which Cs is the dry sediment density, as Z 6.85 ! 106 N m2 and bs Z 2.55 (Owen, 1975). The sediment budget calculations start from the top layer. If the dry sediment mass per unit bed area eroded during one time step (FEDt) did not exceed the total amount Fs contained in the top layer, the total erosion rate for this time step was set to be FE. Otherwise, the erosion rate from this layer was set to be Fs/Dt and the value of ke was decreased by Fs/FEDt. The same procedure was applied for the underlying bed layers until ke Z 0 or the bottom bed layer was reached. The total erosion rate for the present time step was the sum of the erosion rates for all 10 (NLAY) layers. The dry sediment densities for these 10 layers increased smoothly downward from 80 g l1 in the top layer to 152 g l1 in the bottom layer, while the maximum loads all were 5 kg m2. The deposition law followed Krone (1962): Cb ws b ½1 tb =tD for tb !tD FD Z ; ð11Þ 0 for tb RtD where FD is the mass deposited per unit area per time, Cb is the value C in the bottommost (‘b’) water layer, and ws b the corresponding settling velocity, and tD is the critical shear stress required for deposition. C was assumed to be smaller than a maximum value Cmax. At Cb O Cmax, FD was decreased so that the resulting Cb Z Cmax. The mud was assumed stationary when Cb O Cmax Z 90 g l1 (Odd et al., 1992). In water with salinity larger than a critical value Sfloc, the floc settling velocity is strongly dependent on C. Accordingly, for S O Sfloc, the settling velocity was
b
1 wsr Zk1 Cm hind
m 1Cam Chind m2 ; 1 Chind =Cgel ð1 Chind =rs Þ
where Chind is the value of C at which hindered settling starts, Cgel is the sediment concentration at which a space-filling network is formed, wsf is the free settling velocity, wsr is a reference settling velocity for hindered settling, k1, m1, am and bm are empirical constants, and m2 is an empirical exponent accounting for non-linear effects of return flow and wake due to the settling of sediment particle (Winterwerp, 2002). In the water column, C should be smaller than Cgel. The settling velocity ws is set to be the free settling velocity wsf for salinity !Sfloc and C ! Chind. Suspended sediment acts as a sink for turbulent kinetic energy. This decreases the values of the eddy viscosity and the sediment diffusivity. Thus, for a given shear stress the velocity gradient will be steeper in the stratified flow case and the drag will be reduced (Hill and McCave, 2001). The MonineObhukhov length was used to account for the effects of sediment stratification: Lz Z
U3b r0 ; kgr0 w0
ð13Þ
where U*b is the friction velocity at the bottom, k Z 0.4 is the von Karman constant and r0 w0 is the turbulent density flux: r0 w0 Z KH
vr r Z 0 PB ; vz g
ð14Þ
where PB is the buoyancy production (or sink), g vr PB Z KH : r0 vz Following Businger et al. (1971), the velocity profile in a stratified flow is given by: Ub HCz aðHCzÞ UZ C ; ð15Þ ln z0 Lz k where a is an empirical constant between 4.7 and 5.2. A quadratic friction law at the bed was assumed, tb =r0 ZU2b ZCD U2b ;
ð16Þ
where Ub is the velocity at bottommost grid point zb. Substituting Eq. (16) into Eq. (15) yields a relationship for CD,
752
W.B. Guan et al. / Estuarine, Coastal and Shelf Science 65 (2005) 747e762
1ZAðCD Þ1=2 B=CD ;
ð17Þ
with 1 HCzb AZ ln ; k z0
BZ
aðHCzb ÞPB : U3b
Thus A and B are known and independent of CD. The value of CD was calculated from Eq. (17). For the non-Newtonian fluid-mud mixture, the dependence of the apparent viscosity on the shear rate was calculated following Le Hir et al. (2001), m0 1Cam Cbm r0 1 aB CbB qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiCKIW ; C r0 2 2 Y0 C ðvu=vzÞ Cðvv=vzÞ
KMB Z
ð23Þ ð18Þ
where aB, bB, and Y0 are empirical coefficients, m0 is the dynamic viscosity of clear water, and in the fluid mud layer h 8 i3 < KMB0 1 Rig =Rigc 2 for 0!Rig !Rigc KIW Z 0 ; for Rig RRigc : KMB0 for Rig !0 ð19Þ where KMB0 is a empirical constant, Rig is the gradient Richardson number, 2 2 g vr vu vv Rig Z þ ; ð20Þ r0 vz vz vz and Rigc is the critical gradient Richardson number (Kantha and Clayson, 2000). Without any guidance from the literature, the background diffusivity coefficient KHB was assumed to be constant. The boundary conditions at the closed lateral boundaries are straightforward. The conditions of zero mass, momentum, salt, sediment, and turbulence fluxes must be satisfied at a solid lateral boundary. Along the shoreline where river inflow or outflow may occur, the depth-averaged velocities can be calculated. Assuming that the vertical distributions of velocities u and v are parabolic, they are specified as follows (Chen, 1994): uZ3u
ðHCzÞ2 ; D2
ð21Þ
2
vZ3v
ðHCzÞ ; D2
vertical mean velocity components used in the external (vertically integrated) mode of the dynamic equations. At the head of the river, zero salinity is applied, and C can approximately be specified as vertically uniform. Generally, zero normal gradient conditions are appropriate for q2 and l at the riverine boundary. At the seaward open boundary, when there is an inflow, S and C need to be prescribed. Immediately after low water slack, the salinity of the inflowing water will not be equal to the salinity of seawater. It will take a certain time before the maximum salinity Smax is found at the boundary. This lag effect in the salinity at the boundary was modeled following Thatcher and Harleman (1972) as: bðtÞSmax C½1 bðtÞSðtlws Þ for tlws %t%tlws CTTH SZ ; Smax for tlwsCTTH %t%thws
ð22Þ
where D Z H C h is the total depth of the water column, and u and v denote the depth-averaged velocities, i.e., the
where TTH is the ThatchereHarleman time lag (return time), tlws and thws are the times of the last low water slack and the next high water slack. For time t within tlws and tlws C TTH, the salinity S(t) is obtained by a weighted interpolation of S(tlws) and Smax. The salinity is equal to Smax during the remaining inflow period from t Z tlws C TTH to t Z thws. The time function b(t) is defined by: t tlws 2 bðtÞZsin 0:5p : ð24Þ TTH The following equation was used to calculate the seaward, open boundary conditions for C: ðCL CCH Þ ðCH CL Þ t tlws p C sin 2pmin ;1 CZ ; 2 2 2 TF ð25Þ where TF is the mean duration of flood current, CL is C at low water slack, and CH is C at peak flood current. At the seaward open boundary during an outflow, v v ðS; CÞCUN ðS; CÞZ0; vt vN
ð26Þ
where N denotes the coordinate normal to the boundary, and UN equals u or v according to the normal direction of the boundary. As far as q2 and l are concerned, any advection at the lateral boundaries can be assumed to be zero. For tidal calculations, either the free surface elevation or the tidal current velocity can be formulated. Observations of sea levels are more available than ocean currents. By modifying two classical open boundary conditions proposed by Ippen (1966) and Reid and Bodine (1968), Ly et al. (2000) gave a new formulation.
753
W.B. Guan et al. / Estuarine, Coastal and Shelf Science 65 (2005) 747e762 Table 1 Specification of some model parameters Parameter
Description
Value
Source
DtE
Time step for the external mode
2s
Dt
Time step for the internal mode
6s
r0 CM CH q
Reference density Dimensionless coefficient of horizontal eddy viscosity Dimensionless coefficient of horizontal eddy diffusivity Filter coefficient in the Asselin filter used in the solving the momentum equations Critical value of the Richardson number GH Critical elevation difference chosen so as to prevent the drying and flooding of a grid box at alternate time steps Minimum allowable value of the total water depth to avoid singularities in some relevant terms as a grid box dries Effective roughness length Empirical constant in the velocity profile equation for stratified flow Dynamic viscosity of clear water Empirical coefficient for the apparent viscosity of mudewater mixture Empirical exponent for the apparent viscosity of mudewater mixture Empirical coefficient for the Bingham stress of mudewater mixture Empirical exponent for the Bingham stress of mudewater mixture Empirical constant for the Moore rheological model of mudewater mixture Empirical constant reflecting the effect of internal wave on vertical mixing Critical gradient Richardson number above which the effect of internal wave should be taken into account Background diffusivity coefficient Critical salinity value above which flocculation is enhanced Critical C value above which the hindered settling occurs Free settling velocity Empirical coefficient for the settling velocity of floc Empirical exponent for the settling velocity of floc Empirical exponent accounting for non-linear effects of return flow and wake due to the settling of sediment Sediment concentration at which a space-filling network is formed Critical shear stress for deposition Erosion constant
1012 g l1 0.05 0 0.1
Determined by the CFL condition Determined by the CFL condition Dong et al. (1997) Mellor (2001) Mellor (2001) Mellor (2001)
120 0.02 m
Model calibration Flather and Heaps (1975)
0.2 m
Flather and Heaps (1975)
0.5 mm 4.9 1 ! 103 N s m2 0.2 1 8.7 ! 107 N m2 2.55 1 ! 103 s1 5 ! 104 m2 s1 0.33
Model calibration Businger et al. (1971) Molecular diffusion Ross (1988) Ross (1988) Ross (1988) Ross (1988) Le Hir et al. (2001) Model calibration Geyer and Smith (1987)
5 ! 105 m2 s1 0.01 10 g l1 0.1 mm s1 0.089 mm s1 1.3 1
Model Model Model Model Model Model Model
120 g l1 0.6 N m2 2.5 ! 103 kg m2 s1
Winterwerp (1999) Model calibration Model calibration
GHc e Dcrit z0 a m0 am bm aB bB Y0 KMB0 Rigc KHB Sfloc Chind wsf k1 m1 m2 Cgel tD E
According to Ly et al. (2000), the depth-averaged current normal to the open boundaries has the form: UN ZUNi Gðg=DÞ
1=2
ðhi hs Þ;
ð27Þ
where UNi and hi are the first interior depth-averaged velocity and sea level in the model domain, respectively. If the normal direction of the boundary is the same as the direction of the coordinate axis, Eq. (27) has a positive sign and otherwise a negative sign. The quantity hs is the specified time dependent tidal elevation using a harmonic expansion of the form: hs ðtÞZA0 C
NT X
Am cosðum tC40m 4m Þ;
ð28Þ
m
where A0 represents the residual sea level, Am, um, 40m, and 4m are, respectively, the amplitude, frequency, initial value of the phase umt C 40m and the epoch of the m-th tidal constituent, and NT is the number of tidal constituents.
calibration calibration calibration calibration calibration calibration calibration
The model domain (Fig. 1a) includes the entire Jiaojiang Estuary, the two tributaries, and coastal waters (Taizhou Bay). In the horizontal, an orthogonal curvilinear grid was used, with 102 ! 29 grid elements. In the vertical, 15 unequally spaced sigma-layers were used in the water, and 10 layers in the substrate, with a much higher resolution near the bottom. The bottommost grid was located 0.06 m above the bed. The simulation started at 0 h on April 10, 1991, and ended at 1200 hours on April 22, 1991. The model parameters are given in Table 1. The initial salinity distribution was assumed to be vertically uniform. From Songpuzha (Fig. 1) to the eastern boundary, the salinity was assumed to be constant along all cross-channels and to increase gradually from 0 to 24. For the other regions, the salinity was set to be 0. Initial values of C were 5 g l1. The initial surface elevations and velocities were set to zero everywhere. The stationary sediment in the bed was assumed to be uniformly 30 kg m2 everywhere.
754
W.B. Guan et al. / Estuarine, Coastal and Shelf Science 65 (2005) 747e762
Velocity (m s-1)
(a) 1 0 -1 -2
20
80
100
120
140
160
180
80
100
120
140
160
180
80
100
120
140
160
180
80
100
120
140
160
180
(b)
16 12 8
Suspended sediment concentration (g l-1)
4 0
24
(c)
20 16 12 8 4 0
50
(d)
40 30 20 10 0
Time (h) Fig. 2. Comparison of observed (solid lines) and predicted (dashed lines) (a) longitudinal velocity at 1.75 m above the bed, and the suspended sediment concentration (SSC, in g l1) at (b) 1.67 m, (c) 0.55 m and (d) 0.35 m above the bed, in April 1991 at station M1.
In addition, no wind, evaporation and precipitation are considered at the sea surface. Since the seaward open boundary is parallel to the depth contours, the surface elevation at this boundary is horizontal. The surface elevation was calculated from 11 tidal constituents (Q1, O1, P1, K1, N2, M2, S2, K2, M4, MS4 and M6) and zero residual sea level (A0). The outer salinity Smax was set to be 25 and the return time TTH was set to be 2.5 h all along the seaward open boundary. The time of the last low water slack tlws was recorded by the numerical model for each grid point of the open boundary. Eq. (25) was used to calculate the C boundary conditions during an inflow. At the sea surface, the C at the low water slack CL was set to 0.5 g l1, while the C at the peak flood current CH was set to 1.0 g l1. At the bottom, CL Z 2.0 g l1 and CH Z 4.0 g l1. In April 1991, the discharge of the Lingjiang River ranges from 71 to 954 m3 s1 with a mean value of
225 m3 s1. During the simulations, the instantaneous freshwater discharge was interpolated from these daily gauging data. The discharge in the Yongning River was zero. The salinity of the freshwater was 0, while its C was 0.18 g l1.
3. Results The model was verified against the field data including tidal levels at the T3 station, the longitudinal velocity and suspended sediment concentration (C ) data at station M1 (see Fig. 2 and see the standard experiment in Table 2) and longitudinal velocity, salinity and C data at stations C1 and C2. The currents at 1.75-m elevation, including the tidal asymmetry, were faithfully reproduced by the model. The flood current duration is shorter and the peak current is greater than the
755
W.B. Guan et al. / Estuarine, Coastal and Shelf Science 65 (2005) 747e762 Table 2 Root mean square (RMS) errors and other statistical values at mooring station M1 Elevation
1.75 m
1.67 m
Parameter
RMS error in longitudinal velocity (m s1)
RMS error in C (g l1)
Field data Standard experiment E1 E2 E3 E4 E5 E6
0.00 0.24 0.28 0.24 0.32 0.23 0.24 0.22
0.00 3.58 15.68 3.72 3.50 3.59 3.95 3.71
0.55 m
0.00 4.39 15.54 4.82 4.87 4.41 5.14 4.29
corresponding value at ebb tide. Ebb velocities peak at about 1 m s1, while flood velocity peak at 1.8 m s1. A similar asymmetry is also found, and faithfully reproduced by the model, at stations C1 and C2. Seaward, the tidal current asymmetry is smaller.
0.35 m
0.00 7.23 15.16 7.39 7.71 7.20 7.56 7.00
0.55 m Minimum of C (g l1)
Maximum of C (g l1)
Mean of C (g l1)
Standard deviation of C (g l1)
1.54 2.07 4.84 2.43 2.90 2.08 2.21 2.17
21.52 20.83 37.90 24.44 23.35 24.38 26.34 19.37
7.39 6.02 20.17 6.74 6.99 5.93 7.45 6.10
3.45 3.43 8.67 4.14 4.26 3.39 4.43 3.34
The data at the two lower elevations (0.55 and 0.35 m above the bed) of station M1 (Fig. 2c and d) show five flood tide peaks in C during the first 60 h of the records. These peaks were reproduced successfully by a vertical 2-D model under the assumption of the existence of
(a)
(b)
(c)
Fig. 3. Time series of the predicted (a) longitudinal velocity (m s1), (b) salinity and (c) suspended sediment concentration (SSC, in g l1) at station M1 during a spring tide.
756
W.B. Guan et al. / Estuarine, Coastal and Shelf Science 65 (2005) 747e762
Amount of Sediment (×109 kg)
4
Suspended particulate matter Deposited mud Total sediment mass
3
2
1
0 80
100
120
140
160
180
Volume-averaged SSC Volume-averaged kinetic energy
8
1.6
1.2 6 0.8 4 0.4
2
Kinetic Energy (m2 s-2)
Suspended Sediment Concentration (g l-1)
Time (h)
0 80
100
120
140
160
180
Time (h)
Fig. 4. Time series of sediment amount, volume-averaged SSC and volume-averaged kinetic energy over the estuarine part from Songpuzha to Yongquan.
a limited soft mud layer deposited at the bed (Guan et al., 1998). Since the 3-D model also includes a conceptually similar sediment bed module, it also successfully reproduces these peaks (Fig. 2bed). The comparison between measured and simulated surface elevations for the T3 station at the mouth near Haimen also shows a good agreement between field data and model results, both in amplitude and phase (not shown). The tidal wave at station T3 shows a tidal asymmetry, although a symmetric tidal wave is forced at the nearby seaward open boundary. The model predicts that during a tidal cycle, the isohaline of 1 migrates between Luocun and Songpuza; it also predicts that at station M1 (see Fig. 3) during peak currents the vertically-averaged suspended sediment concentration (C ) commonly reaches 5 g l1 in the turbidity maximum zone, that the near-bottom C peaks at 90 g l1 in the estuary, that at the beginning of the flood tide the values of C are higher than those at peak currents because of low tidal level, that the vertical profiles of velocity near the bed are not logarithmic when C O 20 g l1, and that there is a profound asymmetry in currents. The model also shows that the lutocline rises during accelerating tidal currents and falls during decelerating currents. A lutocline occurs during the whole tidal cycle except around high water slack. Defining fluid mud as the water with C O 20 g l1 (Mehta and Srinivas, 1993), the distribution and
thickness of fluid mud layer are predicted by the model. Accordingly, at the beginning of the flood tide, fluid mud occurs near Songpuzha and in Taizhou Bay, while there is no fluid mud yet in the estuary from Yongquan to Songpuzha. During swift tidal currents, fluid mud occurs also near Songpuzha and generally throughout the estuary. The maximum thickness of the fluid mud layer is about 0.60 m and this occurs near Haimen. All these predictions match well the observations of Guan et al. (1998). Fig. 4 shows that during the accelerating tidal currents the bottom sediment is eroded and entrained into the mobile sediment, to be deposited upstream at flood tide and downstream at ebb tide. There is a significant exchange of sediment between various areas of the estuary. Overall, the total deposited mud increases with time, as result of net advection from coastal waters. The upper limit of the total suspended sediment mass is about 1.2 ! 106 tons, which is comparable to the annual fluvial sediment discharge. During each tide cycle, there are two peaks in volume-averaged kinetic energy (Fig. 4), with the largest peak corresponding to the peak flood current and the smallest peak to the peak ebb current. The time-series plots of predicted, tidally-averaged salinity, currents and C (Figs. 5e7) suggest that there is little net gravitational circulation in the longitudinal direction and that it is the asymmetry in tidal currents
W.B. Guan et al. / Estuarine, Coastal and Shelf Science 65 (2005) 747e762
757
(a)
(b)
(c)
Fig. 5. Horizontal distributions of the predicted, vertically-averaged (a) velocity, (b) salinity and (c) suspended sediment concentration (SSC, in g l1) time-averaged over several tidal cycles.
that drives the outer sediment into the estuary from the coastal waters. The turbidity maximum is found to extend longitudinally about 16 km, and its center is located off Haimen (Figs. 5 and 6). There, cross-channel gradients in currents and C are apparent, with greater values occurring at the southern side than at the opposed side, indicating the existence of cross-channel residual circulation (Fig. 7). For the tidally-averaged mud fluxes the models predict that, off Haimen the northern side of the estuary is flood-dominant, while the southern side is slightly ebb-dominant. Since there is net flood dominance throughout the estuary mouth, mud is pumped into the estuary from coastal waters and is redistributed asymmetrically within the estuary, resulting in the formation of mud banks. This observation is in quantitative agreement with the location of the mud banks and with port authorities dredging logbooks. 4. Discussion Because of the asymmetry of currents, the eddy viscosity and eddy diffusivity are larger at flood than at
ebb tide. Large Richardson numbers occur in the nearbottom layers which have a strong stratification in C (Jianhua and Wolanski, 1998; Jiang and Mehta, 2000, 2002). Mixing occurs as a result of the breaking of internal waves riding on the lutocline. This process is not included in the turbulence closure model due to strong buoyancy damping at the lutocline. Therefore, background mixing coefficients must be added. Without this, the model would predict that the top of the water column would completely clear of sediment at slack water, a prediction at odds with observations (see also Run E6 discussed below). The numerical experiment described in the last section is called a standard experiment. Following these specifications of the standard experiment (Table 1), eight simulations using the 3-D model were completed to explore parameter sensitivity. Table 3 summarizes the objectives and specifically varied parameters in the eight simulations. Discussion on sediment dynamics associated with these simulations including the standard experiment will be expounded sequentially in this section.
758
W.B. Guan et al. / Estuarine, Coastal and Shelf Science 65 (2005) 747e762
(a)
(b)
(c)
Fig. 6. Longitudinal distributions of the time-averaged (a) along-transect velocity, (b) salinity and (c) suspended sediment concentration (SSC, in g l1) over several tidal cycles.
4.1. Buoyancy effects Time series results from Run E1 (no sedimentinduced buoyancy effects) are plotted in Fig. 8, while the corresponding results from Run E2 (no salinityinduced buoyancy effects) are shown in Fig. 9. Without sediment-induced buoyancy effects (Run E2), the predictions are unrealistic. Indeed, the predicted C values are unreasonably high and sometimes are uniformly equal to 30 g l1 in the vertical, which is unrealistic when compared to field observations (Guan et al., 1998). This finding results from the fact that, when there are no sediment-induced buoyancy effects, there is no damping of turbulence in the presence of fluid mud. The differences between Run E2 and the standard experiment are smaller. Nevertheless, without salinityinduced buoyancy effects, the vertical mixing coefficients are also unrealistically high especially in the upper layers, since all weak stratification near the sea surface has been ignored. 4.2. Hydraulic roughness and drag Run E3 was undertaken to demonstrate how sensitive the model is to the specific value of the
effective roughness length. When this length is increased by a factor of 10, the predicted surface elevation at low tide became unrealistically high, and the predicted ebb tidal currents become unrealistically small. The predicted tidal fluctuations of C also became unrealistic. Run E4 investigated the effect of the stratification on the effective hydraulic drag. At peak currents, stratification in C became stronger and thus the overall hydraulic roughness was reduced. However, this impact was moderate.
4.3. Sensitivity of a tuning constant in the turbulence closure model Guided by laboratory data, Mellor (2001) introduced a Richardson number cutoff and thus a tuning closure constant GHc to the turbulence dissipation in the MelloreYamada turbulence closure model. An empirical value of GHc Z 6.0 was adopted. In this case with GHc Z 6.0 (Run E5), the peaks in C during flood tide became unrealistically small. In the standard experiment, GHc Z 120.
759
W.B. Guan et al. / Estuarine, Coastal and Shelf Science 65 (2005) 747e762
(a)
Table 3 List of the additional experiments. The parameters are identical to those in the standard experiment unless otherwise noted Experiment Objective
Parameters different from the standard experiment
E1
rZ1000C0:78S (g l1)
E2
(b) E3
E4
E5
(c)
E6
E7
E8
Fig. 7. Tidally-averaged, cross-channel distributions in a transect past station M1 of the time-averaged (a) velocity (m s1); the longitudinal component is shown as a contour line, and the cross-channel component as a vector, (b) salinity, and (c) suspended sediment concentration (SSC, in g l1).
4.4. Background vertical mixing Generally, the relative, mud-induced increase of the eddy viscosity is much smaller than that for the eddy viscosity, except for flows with a high C and a low shear rate. Although increasing the mud-induced viscosity is necessary to keep the numerical stability of the model when high C values occur, it does not significantly change the overall result. On the other hand, the Richardson-number-dependent coefficient for the background vertical viscosity, which is considered to account for the effects of the internal wave, can influence the model results near the bed. Run E6 demonstrates that the time-averaged C value near the bed is reduced without the background vertical viscosity associated with internal wave. That is to say, internal waves contribute to the vertical transport of momentum and thereby indirectly affect the transport of mass (Uittenbogaard, pers. comm.).
To investigate sediment-induced buoyancy effects To investigate salinity-induced buoyancy effects To demonstrate how sensitive the model is to the effective roughness length To investigate the effect of stratification on the effective hydraulic drag To demonstrate the sensitivity of the tuning constant GHc in the MelloreYamada turbulence closure model To demonstrate the effect of background vertical mixing controlled by internal wave To understand the consequences of a small freshwater discharge of the Lingjiang River. To understand the consequences of a large freshwater discharge of the Lingjiang River.
rZ1000C0:62C (g l1)
z0 Z 5 mm
aZ0
GHc Z 6
KIW Z 0 m2 s1
Freshwater discharge is set to 71 m3 s1 Freshwater discharge is set to 954 m3 s1
4.5. Effects of the river discharge During a low river discharge (Run E7), the total mass of suspended mass decreases, and this results in increased sedimentation. Salty water moves further upstream, and the vertically-averaged and time-mean values of C are reduced. During high river discharge (Run E8), the sediment is gradually carried out of the estuary. The kinetic energy increases during ebb and decreases during flood. The masses of mobile sediment, stationary sediment, and total sediment are reduced, while turbidity maximum moves seaward to Taizhou Bay. 5. Conclusions A 3-D, POM-based model has been developed here for the highly turbid, macro-tidal, Jiaojiang Estuary and adjoining coastal waters. The model was verified against field observations of current and suspended sediment concentration, C. A fine resolution grid near the bottom was necessary to simulate fluid-mud formation. Sediment-induced buoyancy effects are found to be significant. Turbulence is significantly dampened at the high C. Because of
760
W.B. Guan et al. / Estuarine, Coastal and Shelf Science 65 (2005) 747e762
(a)
(b)
(c)
Fig. 8. Same as Fig. 3 for Run E1.
the strong stratification at the lutocline, the Mellore Yamada scheme tends to under-predict vertical mixing, thus requiring background mixing to be added. The MonineObukhov length scale appears successful at calculating the influence of sediment-induced stratification on the effective hydrodynamic drag. The hydraulic roughness is predicted to be measurably decreased in the very turbid estuary. The effect of mud in increasing the viscosity is small but found to be necessary for numerical stability. The TVD scheme avoids all negative values of salinity and C, which would occur frequently otherwise in the presence of large horizontal gradients of salinity and vertical gradients of C. A lutocline occurs almost over the entire tidal cycle, except near high water slack. During accelerating currents, the lutocline rises, and during slack water it falls. The total mass of mobile sediment trapped in the estuary is around 1.2 ! 106 tons, a value comparable to the annual river discharge. Asymmetry in tidal currents generates tidal pumping and estuarine infilling from sediment pumped from coastal waters. During high flow periods, the sediment is
carried to coastal waters; this sediment is advected back in the estuary during low flow periods. The 3-D model appears able to realistically predict the location of fluid mud banks as the areas where nearbottom currents advect mud either towards the deeper regions of the estuary or in coastal waters near the mouth. The simulations are very sensitive to details of the parameterization, suggesting that predictive 3-D modeling of mud dynamics remain an art rather than a hard science, and that improvements in the model requires process-based studies to improve the parameterization of these processes.
Acknowledgments This study was supported by the China Natural Science Foundation under contract No. 40206014, and by the National Key Basic Research Program from the Ministry of Science and Technology of China under contract No. G1999043702. The field data were collected under the Sino-Australia cooperation between the
W.B. Guan et al. / Estuarine, Coastal and Shelf Science 65 (2005) 747e762
761
(a)
(b)
(c)
Fig. 9. Same as Fig. 3 for Run E2.
Second Institute of Oceanography, People’s Republic of China and the Australian Institute of Marine Science. It is a pleasure to thank Dr. A. Mehta for editing this paper. Two anonymous reviewers helped improve this paper.
References Ariathurai, R., 1974. A finite element model for sediment transport in estuaries. Ph.D. thesis, University of California, Davis, 192 pp. Bi, A.H., Sun, Z.L., 1984. A preliminary study on the estuarine processes in the Jiaojiang River, China. Journal of Sediment Research 3, 12e26 (in Chinese). Blumberg, A.F., Galperin, B., O’Connor, D.J., 1992. Modeling vertical structure of open channel flows. Journal of Hydraulic Engineering 118, 1119e1134. Blumberg, A.F., Mellor, G.L., 1987. A description of a three-dimensional coastal ocean circulation model. In: Heaps, N. (Ed.), ThreeDimensional Coastal Ocean Models. Coastal and Estuarine Science, vol. 4. American Geophysical Union, Washington, DC, pp. 1e16. Businger, J.A., Wyngaard, J.C., Izumi, Y., Bradley, E.F., 1971. Fluxe profile relationship in atmospheric surface layer. Journal of Atmospheric Science 28, 190e201.
Chen, X., 1994. Effects of hydrodynamics and sediment transport process on nutrient dynamics in shallow lakes and estuaries. Ph.D. Dissertation, University of Florida. Dong, L.X., Wolanski, E., Li, Y., 1997. Field and modeling studies of fine sediment dynamics in the extremely turbid Jiaojiang River estuary, China. Journal of Coastal Research 13, 995e1003. Faas, R.W., 1984. Time and density-dependent properties of fluid mud suspensions, NE Brazilian Continental Shelf. Geo-Marine Letters 4, 147e152. Flather, R.A., Heaps, N.S., 1975. Tidal computations for Morecambe Bay. Geophysical Journal of the Royal Astronomical Society 42, 489. Fu, N.P., Bi, A.H., 1989. Discussion on the problems in the sediment transport of the Jiaojiang River. Journal of Sediment Research 8, 51e57 (in Chinese). Geyer, W.R., Smith, J.D., 1987. Shear instability in a highly stratified estuary. Journal of Physical Oceanography 27, 1668e1679. Guan, W.B., Wolanski, E., Dong, L.X., 1998. Cohesive sediment transport in the Jiaojiang River estuary, China. Estuarine, Coastal and Shelf Science 46, 861e871. Hill, P.S., McCave, I.N., 2001. Suspended particle transport in benthic boundary layers. In: Boudreau, B., Jorgensen, B.B. (Eds.), The Benthic Boundary Layer: Transport Processes and Biogeochemistry. Oxford Press, pp. 78e103. Ippen, A.T., 1966. Estuary and Coastline Hydrodynamics. McGrawHill, New York, 744 pp.
762
W.B. Guan et al. / Estuarine, Coastal and Shelf Science 65 (2005) 747e762
James, I.D., 1996. Advection schemes for shelf sea model. Journal of Marine Systems 8, 237e254. Jiang, J., Mehta, A.J., 2000. Lutocline behavior in a high-concentration estuary. Journal of Waterway, Port, Coastal and Ocean Engineering 126 (6), 324e328. Jiang, J., Mehta, A.J., 2002. Interfacial instabilities at the lutocline in the Jiaojiang Estuary, China. In: Winterwerp, J.C., Kranenburg, C. (Eds.), Fine Sediment Dynamics in the Marine Environment. Elsevier, Amsterdam, pp. 125e137. Jianhua, J., Wolanski, E., 1998. Vertical mixing by internal wave breaking at the lutocline, Jiaojiang estuary, China. Journal of Coastal Research 14, 1426e1431. Kantha, L.H., Clayson, C.A., 2000. Numerical Models of Oceans and Oceanic Processes. Academic Press, San Diego, 940 pp. Krone, P.B., 1962. Flume Studies of the Transport in Estuarine Shoaling Processes. Final Report, Hydraulics Engineering Laboratory, University of Berkeley, California, USA, 110 pp. Le Hir, P., 1997. Fluid and sediment ‘‘integrated’’ modeling: application to fluid mud flows in estuaries. In: Burt, N., Parker, R., Watts, J. (Eds.), Cohesive Sediment. John Wiley, Chichester, UK, pp. 417e428. Le Hir, P., Bassoullet, P., Jestin, H., 2001. Application of the continuous modeling concept to simulate high-concentration suspended sediment in a macrotidal estuary. In: McAnally, W.H., Mehta, A.J. (Eds.), Coastal and Estuarine Fine Sediment Processes. Proceedings in Marine Science, vol. 3. Elsevier Science, Amsterdam, pp. 229e248. Le Hir, P., Cayocca, F., 2002. 3D application of the continuous modeling concept to mud slides in open sea. In: Winterwerp, J.C., Kranenburg, C. (Eds.), Fine Sediment Dynamics in the Marine Environment. Proceedings in Marine Science, vol. 5. Elsevier Science, Amsterdam, pp. 545e561. Le Normant, C., Peltier, E., Teisson, C., 1998. Three dimensional modeling of cohesive sediment transport in estuaries. In: Dronkers, J., Scheffer, M.B.A.M. (Eds.), Physics of Estuaries and Coastal Sea. Balkema, Rotterdam, pp. 65e71. Li, Y., Wolanski, E., Xie, Q., 1993. Coagulation and settling of suspended sediment in the Jiaojiang estuary, China. Journal of Coastal Research 9, 390e402. Luyten, P.J., Jones, J.E., Proctor, R., Tabor, A., Tett, P., Wild-Allen, K., 1999. COHERENS e A Coupled HydrodynamicaleEcological Model for Regional and Shelf Seas: User Documentation. MUMM Report, Management Unit of the Mathematical Models of the North Sea, 914 pp. Ly, L.N., O’Connor, W.P., Paduan, J.D., 2000. A barotropic simulation of semi-diurnal and diurnal tides for the Monterey Bay region. In: Blain, W.R., Brebbia, C.A. (Eds.), Hydraulic Engineering Software VIII. Proceedings of the Eighth International Conference on Hydraulic Engineering Software. Witpress, Southampton, Boston, pp. 155e169. Mehta, A.J., Srinivas, R., 1993. Observations on the entrainment of fluid mud in shear flow. In: Mehta, A.J. (Ed.), Nearshore Estuarine Cohesive Sediment Transport. American Geophysical Union, Washington, DC. Mellor, G.L., 2001. One-dimensional, ocean surface layer modeling, a problem and a solution. Journal of Physical Oceanography 31, 790e809. Mellor, G.L., Hakkinen, S., Ezer, T., Patchen, R., 2001. A generalization of a sigma coordinate ocean model and an intercomparison of model vertical grids. In: Pinardi, N., Woods, J.D. (Eds.), Ocean Forecasting: Conceptual Basis and Applications. Springer-Verlag, pp. 55e72.
Mellor, G.L., Yamada, T., 1982. Development of turbulence closure model for geophysical fluid problems. Reviews of Geophysics and Space Physics 20, 851. Odd, N.V.M., 1988. Mathematical modeling of mud transport in estuaries. In: Dronkers, J., van Leussen, W. (Eds.), Physical Processes in Estuaries. Springer Verlag, pp. 503e531. Odd, N.V.M., Bentley, M.A., Waters, C.B., 1992. Observations and analysis of the movement of fluid mud in estuaries. In: Mehta, A.J. (Ed.), Nearshore and Estuarine Cohesive Sediment Transport. Coastal and Estuarine Studies, vol. 42. AGU, Washington, pp. 430e446. Odd, N.V.M., Cooper, A.J., 1989. A two dimensional of the movement of fluid mud in a high energy turbid estuary. Journal of Coastal Research 5, 185e193. Owen, M.W., 1975. Erosion of Avonmouth Mud. Report No. INT 150, Hydraulics Research Station, Wallingford, United Kingdom. Partheniades, E., 1965. Erosion and deposition of cohesive soils. Journal of the Hydraulic Division, 105e139 (ASCE 91, No. HY1). Reid, R.O., Bodine, B.R., 1968. Numerical model for storm surges in Galveston Bay. ASCE, Journal of Waterways and Harbors Division 94, 33e57. Roberts, W., 1993. Development of a mathematical model of fluid mud in the coastal zone. Proceedings of the Institution of Civil Engineers e Water, Maritime and Energy 101, 173e181. Ross, M.A., 1988. Vertical structure of estuarine fine sediment suspensions. Ph.D. Thesis, University of Florida, Gainesville, 206 pp. Smagorinsky, J., Manabe, S., Holloway, J.L., 1965. Numerical results from a nine-level general circulation model of the atmosphere. Monthly Weather Review 93, 727e768. Teisson, C., Simonin, O., Galland, J.C., Laurence, D., October 1992. Turbulence and mud sedimentation: a Reynolds stress model and a two-phase flow model. In: Proceedings of the 23rd International Conference on Coastal Engineering. American Society of Civil Engineers, Venice, Italy, pp. 2853e2866. Thatcher, M.L., Harleman, R.F., 1972. A Mathematical Model for the Prediction of Unsteady Salinity Intrusion in the Estuaries, Massachusetts Institute of Technology, Report No. 144. Vinzon, S.B., Mehta, A.J., 2001. Boundary layer effects due to suspended sediment in the Amazon estuary. In: McAnally, W.H., Mehta, A.J. (Eds.), Coastal and Estuarine Fine Sediment Processes. Proceedings in Marine Science, vol. 3. Elsevier Science, Amsterdam, pp. 359e372. Winterwerp, J.C., 1999. On the dynamics of high-concentrated mud suspensions. Doctoral Thesis for the Technical University of Delft. Winterwerp, J.C., 2002. On the flocculation and settling velocity of estuarine mud. Continental Shelf Research 22, 1339e 1360. Winterwerp, J.C., Uittenbogaard, R.E., de Kok, J.M., 2001. Rapid siltation from saturated mud suspensions. In: McAnally, W.H., Mehta, A.J. (Eds.), Coastal and Estuarine Fine Sediment Processes. Proceedings in Marine Science, vol. 3. Elsevier Science, Amsterdam, pp. 1e18. Wolanski, E., Chappell, J.R., Ridd, P., Vertessy, R., 1988. Fluidization of mud in estuaries. Journal of Geophysical Research 93, 1535e 1545. Zhu, Y.K., 1986. Some characteristics of stream-like macro-tidal estuary (Jiaojiang). Geographical Study 5 (1), 21e29 (in Chinese).
3-D fluid-mud dynamics in the Jiaojiang Estuary, China W.B. Guan a,b,*, S.C. Kot b, E. Wolanski c a
Laboratory of Ocean Dynamic Processes and Satellite Oceanography, Second Institute of Oceanography, State Oceanic Administration, Hangzhou 310012, China b Department of Mechanical Engineering, Hong Kong University of Science and Technology, Clear Water Bay, Hong Kong, PR China c Australian Institute of Marine Science, P.M.B. No. 3, Townsville M.C., Queensland 4810, Australia Received 10 March 2004; accepted 27 May 2005 Available online 9 November 2005
Abstract A 3-D model has been developed for the muddy Jiaojiang Estuary and adjoining coastal waters, and verified against field observations. To simulate fluid-mud formation, the model uses a fine resolution grid near the bottom and involves coupling processes between hydrodynamics and fluid mud such as the sediment-induced buoyancy, increasing turbulent kinetic energy sink and kinematic viscosity, mixing by internal waves riding on the lutocline, and non-Newtonian properties of fluid mud. The effective hydrodynamic drag was reduced in the presence of fluid mud. It is shown that the estuary is infilled by tidal pumping and that longitudinal and transversal gradients of suspended sediment concentration, salinity, and currents control the formation of mud banks. Thus a 3-D model is necessary to estimate the fate of mud, although the model results are very sensitive to details of the parameterization of the hydrodynamics-mud feedback processes. Ó 2005 Elsevier Ltd. All rights reserved. Keywords: cohesive sediment; model; turbidity maximum; lutocline; fluid mud; China
1. Introduction The dynamics of fine-grained sediment is of major importance in many coastal areas and estuaries, because this sediment generates siltation and degrades the water quality. Fine-grained sediment (mud) is composed mainly of clay (particle size ! 2 mm) and silt (particle size ! 63 mm). Fluid mud is formed at suspended sediment concentrations greater than 5 g l1 and is the major mechanism of the transport of cohesive sediments in coastal zones and estuaries (Faas, 1984; Odd et al., 1992). Estuarine models that include fluid mud flows seldom deal with the internal structure of fluid mud layers; * Corresponding author. Laboratory of Ocean Dynamic Processes and Satellite Oceanography, Second Institute of Oceanography, State Oceanic Administration, Hangzhou 310012, China. E-mail address: [email protected] (W.B. Guan). 0272-7714/$ - see front matter Ó 2005 Elsevier Ltd. All rights reserved. doi:10.1016/j.ecss.2005.05.017
instead they assume a sharp interface between the mud layer and the water above, with semi-empirical exchange formulations between them for entrainment and settling. For example, Odd and Cooper (1989) developed a depth-averaged, 2-D fluid mud model to describe the formation, movement, re-entrainment and dewatering of fluid mud layers, and the fluid mud was assumed to be a Bingham plastic fluid with a uniform density. Roberts (1993) replaced the Bingham plastic assumption with a high viscosity Newtonian fluid model. Le Normant et al. (1998) coupled a horizontal 2-D fluid mud model to a 3-D model for simulating the turbidity maximum in the Loire Estuary. This model assumed a sharp discontinuity at water/fluid mud interface and no vertical velocity gradient within the fluid mud layer, and it failed to reproduce some observations (Le Hir et al., 2001). In fluid mud flows the suspended sediment concentration, C, is so high that interactions between particles
748
W.B. Guan et al. / Estuarine, Coastal and Shelf Science 65 (2005) 747e762
and water become strong, and a description of its mechanics needs a multiphase approach. This was achieved by Teisson et al. (1992) who modeled the fluid mud layer by using a two-phase flow model. This approach also enables a combination of hydrodynamics and soil mechanics, for which pore-water movement relative to settling particles is considered. However, the computations are very time-consuming; also, unsteady configurations cannot be readily simulated. In many cases of natural fine sediment transport, particles and water have similar horizontal velocity components, so that a unique set of horizontal momentum equations can be used for the water/ sediment mixture. For example, Le Hir (1997) presented an integrated modeling approach, simulating water- and mud dynamics by solving for the mass conservation and momentum equations over the entire column. In their vertical 2-D modeling study of the fluid-mud dynamics and the turbidity maximum in the Gironde Estuary, Le Hir et al. (2001) replaced the integrated formulations by a continuous one, which is probably more representative of a progressive transition between the water column and bottom sediment. Le Hir and Cayocca (2002) adopted the concept of continuous modeling to simulate 3-D turbidity currents on a slope. In their model, a zlevel coordinate system with a stepped bottom topography was chosen, which is perhaps suitable for their simplified problem but not for an actual case, especially a macro-tidal estuary. Based on simple analysis and historical measurements in the Amazon Estuary, Vinzon and Mehta (2001) found that it is essential to account for the strong interaction between flow and sediment dynamics within the lower boundary layer for modeling sediment discharge. In this paper these concepts were applied to a finescale, 3-D model of the topographically-complex, muddy, macro-tidal, Jiaojiang Estuary in China for which an extensive data set on the dynamics of water and suspended sediment was available. It is shown that the model can fairly reliably predict cohesive sediment dynamics, but that the results are very sensitive to details of the parameterization of these various processes. The model thus points to the need for processbased studies of these processes. The extremely turbid Jiaojiang Estuary (Fig. 1) is located about 200 km south of the Yangtze River (Changjiang River) mouth; it drains the Lingjiang and Yongning rivers. The Yongning River is dammed and has minimal runoff. The Lingjiang River drains a catchment area of 6520 km2. Its annual runoff is 6.66 ! 109 m3, corresponding to a mean discharge of 211 m3 s1 (Li et al., 1993). The Jiaojiang Estuary is about 35 km long. In the estuary, the mean width of the channel is about 1.2 km with a maximum of 1.8 km at the mouth. The estuary faces Taizhou Bay, which is shallow (the depth varies between 1 and 3 m at low
spring tide). Semi-diurnal macro-tides prevail with a mean tidal range of about 4 m and a maximum tidal range of 6.3 m. The vertically-averaged tidal current peaks at 2.0 m s1. At Haimen, the flood tide lasts 5.1 h while the ebb tide lasts 7.3 h (Dong et al., 1997). Off Haimen, the maximum flood and ebb currents are, respectively, 2.1 and 1.8 m s1 (Zhu, 1986). The ratio of the freshwater volume to the tidal volume is about 0.04, but more than 0.1 during periods of very high river discharge (Bi and Sun, 1984). A turbidity maximum zone is present where C peaks at 40 g l1 (Bi and Sun, 1984; Fu and Bi, 1989). These authors estimated that the fluvial sediment inflow is about 1.2 ! 106 tons year1, corresponding to a mean sediment concentration of 0.18 g l1. Most fine-grained sediment found in the estuary comes from coastal waters, and originates from the Yangtze River (Bi and Sun, 1984; Zhu, 1986). Siltation in the estuary is about 0.2 m year1, only 5% of which is due to fluvial material (Bi and Sun, 1984; Fu and Bi, 1989; Li et al., 1993; Dong et al., 1997; Guan et al., 1998).
2. Methods 2.1. Field observations An Inter-Ocean S4 electromagnetic current meter (Courtesy, Mark Geneau, Inter-Ocean) was deployed 1.75 m above the bed at mooring station M1 (see Fig. 1 for a location map) during April 12e18, 1991. At station C1 also three backscatterance nephelometers were deployed at elevations of 0.35, 0.55 and 1.67 m above the bed. These three sensors recorded data at an interval of 5 min. The nephelometers were calibrated in situ to yield C (Li et al., 1993; Dong et al., 1997; Guan et al., 1998). At stations C1 and C2, 25-h-long observational stations were maintained at spring tides and neap tides, during which hourly, vertical profiles of salinity, C and currents were obtained using the ship-born CTD cumnephelometer of Wolanski et al. (1988) and a rotor-type current meter (model SLC9-1, made by the Institute of Marine Instrument, Ocean University of China, Qingdao). Longitudinal and transversal transects of temperature, salinity and C were also obtained using the shipborn CTD. 2.2. 3-D model The generalized s-coordinate Princeton Ocean Model (POM) was used (Blumberg and Mellor, 1987; Mellor et al., 2001) to simulate the hydrodynamics in the estuary, and further a sediment transport module was developed to integrate with the hydrodynamics model. The incompressible, hydrostatic and Boussinesq approximations were used. The wateresediment mixture
W.B. Guan et al. / Estuarine, Coastal and Shelf Science 65 (2005) 747e762
749
(a) 28.85 N 28.8 Yongquan
0
2
4
6 km
28.75
Latitude (°N)
Lutou
Taizhou Bay
Shixianfu
M1, C2
28.7 Luocun
C1 LaoshuyuSongpuzha
Shanpu Haimen
Jiaojiang Estuary 28.65
T3
Korea
ina
28.6
Ch
Japan Study site
28.55
121.25
121.3
121.35
121.4
121.45
121.5
121.55
121.6
121.65
Longitude (°E)
(b)
Yongquan
Luocun
Songpuzha
Haimen
0
z (m)
-2
-4
-6
-8 0
5
10
15
20
25
30
35
40
45
Distance along transect (km)
(c)
C2
M1
Fig. 1. (a) Orthogonal curvilinear grid used in the simulations. M1 is the mooring stations; C1 and C2 are anchor stations; T3 is the tide stations. The transect connecting all thick dots will be used to show the model results below. The triangles denote locations of towns. (b) Vertical coordinates for the longitudinal transect indicated in Fig. 2. The triangles denote locations of towns. (c) Bathymetry (in m) for the estuarine part from Songpuzha to Yongquan. A dike built on the central shoal of the estuary has been indicated.
750
W.B. Guan et al. / Estuarine, Coastal and Shelf Science 65 (2005) 747e762
was treated as a single-phase fluid in which all particles follow water except for their settling velocity (Winterwerp et al., 2001). Some non-Newtonian effects of highconcentrated mud suspensions were included. The equations were written in a right-handed Cartesian coordinate system, (x, y, z, t), with the z-axis pointing vertically upwards, and t denoting time. u, v and w are the velocity components in the x-, y- and z-directions, respectively, r is the fluid density, r0 is the reference density, g is the gravitational acceleration, AM is the horizontal turbulent eddy viscosity, KM is the vertical turbulent eddy viscosity, KMB is the background vertical viscosity representing the mixing process due to the internal wave and the increment due to the suspended sediment, and f is the Coriolis parameter. The continuity and momentum equations were:
sediment concentration, C, and therefore a background turbulence value was imposed. The equations for salinity, S, and suspended sediment concentration, C, were:
vu vv vw C C Z0; vx vy vz
where ws is the settling velocity of suspended sediment particles (positive downward), AH and KH are the horizontal and vertical turbulent mass diffusivities, respectively, and KHB is the background value of the vertical diffusivity. The leapfrog scheme used in the original POM was replaced by a total-variation-diminishing (TVD) scheme with the superbee limiter as a weighting function between the upwind scheme and either the LaxeWendroff scheme in the horizontal or the central scheme in the vertical for the advective terms of scalar variables (James, 1996; Luyten et al., 1999). The drying and flooding algorithm of Flather and Heaps (1975) was used. Temperature differences are small in the Jiaojiang Estuary and were thus ignored (Guan et al., 1998). The density r (g l1) of the sedimentewater mixture was computed as (Odd, 1988):
vu vu2 vuv vuw C C C fv vt vx vy vz 1 vp v vu CFx C ðKM CKMB Þ Z ; r0 vx vz vz vv vvu vv2 vvw C C C Cfu vt vx vy vz 1 vp v vv CFy C ðKM CKMB Þ ; Z r0 vy vz vz
0Z
1 vp r g; r0 vz r0
ð1Þ
ð2Þ
ð3Þ
ð4Þ
with v vu v vu vv Fx Z C 2AM C AM ; vx vx vy vy vx v vv v vv vu C Fy Z 2AM C AM : vy vy vx vx vy In these equations, a differentiation is made between background turbulence and flow-generated turbulence. This was necessary because strong sediment-induced stratification results in zero or negligible vertical mixing, when the MelloreYamada level 2.5 scheme is used. However, field data reveal that background level mixing exists due to internal waves riding on the lutocline (Jianhua and Wolanski, 1998; Jiang and Mehta, 2000, 2002). This background turbulence affects also the conservation equations for salinity, S, and suspended
vS vuS vvS vwS C C C vt vx vy vz v vS v vS v vS Z AH C AH C ðKH CKHB Þ ; ð5Þ vx vx vy vy vz vz
vC vuC vvC vðw ws ÞC C C C vt vx vy vz v vC v vC v vC Z AH C AH C ðKH CKHB Þ ; ð6Þ vx vx vy vy vz vz
r rw rZrw CC s Z1000C0:78SC0:62C; rs
ð7Þ
where C is in g l1, rw Z 1000 C 0.78S is the water density (g l1), and rs Z 2650 (g l1) is the sediment density. Horizontal turbulent mixing, which describes the effect of sub-grid scale motion, was represented by a Smagorinsky-type formulation (Smagorinsky et al., 1965). Vertical turbulent mixing was calculated using the Mellor and Yamada (1982) turbulence closure model. This model calculates q2, the turbulence kinetic energy, and q2l, where q is the turbulence velocity scale and l is the length scale. The wall proximity function of Blumberg et al. (1992) was used, to account for the presence of free surface and to reproduce data that suggest that l is larger near the free surface than near a solid wall. Mellor (2001) introduced a correction function in the turbulence closure model based on experiments, and thus a critical constant GHc must be
751
W.B. Guan et al. / Estuarine, Coastal and Shelf Science 65 (2005) 747e762
tuned so that model calculations agree more closely with field data. At the bottom, ws CCðKH CKHB Þ
vC ZFD FE ; at zZ Hðx; yÞ; vz
ð12Þ
m
with ð8Þ
where FD is the deposition rate, and FE is the erosion rate of sediment. To calculate FE, the bed was discretized into a number (NLAY Z 10) of elements. The erosion rate equation developed by Ariathurai (1974) from the flume data of Partheniades (1965) was used: k E½t =t 1 for tb OtE FE Z e b E ; ð9Þ 0 for tb %tE where tb is the shear stress at the bottom, tE is the threshold shear stress, E is the erosion constant, and ke (ke Z 1 at the uppermost bed layer) is the remaining ratio of the capability of flow for erosion of the newly exposed layer. The threshold bed shear strength tE was: tE Zas Cbs s ;
for C%Chind max k1 Cm1 ; wsf m ; ws Z ð1C=Cgel Þ 2 ð1C=rs Þ wsr for Chind !C!Cgel 1Ca Cbm (
ð10Þ
in which Cs is the dry sediment density, as Z 6.85 ! 106 N m2 and bs Z 2.55 (Owen, 1975). The sediment budget calculations start from the top layer. If the dry sediment mass per unit bed area eroded during one time step (FEDt) did not exceed the total amount Fs contained in the top layer, the total erosion rate for this time step was set to be FE. Otherwise, the erosion rate from this layer was set to be Fs/Dt and the value of ke was decreased by Fs/FEDt. The same procedure was applied for the underlying bed layers until ke Z 0 or the bottom bed layer was reached. The total erosion rate for the present time step was the sum of the erosion rates for all 10 (NLAY) layers. The dry sediment densities for these 10 layers increased smoothly downward from 80 g l1 in the top layer to 152 g l1 in the bottom layer, while the maximum loads all were 5 kg m2. The deposition law followed Krone (1962): Cb ws b ½1 tb =tD for tb !tD FD Z ; ð11Þ 0 for tb RtD where FD is the mass deposited per unit area per time, Cb is the value C in the bottommost (‘b’) water layer, and ws b the corresponding settling velocity, and tD is the critical shear stress required for deposition. C was assumed to be smaller than a maximum value Cmax. At Cb O Cmax, FD was decreased so that the resulting Cb Z Cmax. The mud was assumed stationary when Cb O Cmax Z 90 g l1 (Odd et al., 1992). In water with salinity larger than a critical value Sfloc, the floc settling velocity is strongly dependent on C. Accordingly, for S O Sfloc, the settling velocity was
b
1 wsr Zk1 Cm hind
m 1Cam Chind m2 ; 1 Chind =Cgel ð1 Chind =rs Þ
where Chind is the value of C at which hindered settling starts, Cgel is the sediment concentration at which a space-filling network is formed, wsf is the free settling velocity, wsr is a reference settling velocity for hindered settling, k1, m1, am and bm are empirical constants, and m2 is an empirical exponent accounting for non-linear effects of return flow and wake due to the settling of sediment particle (Winterwerp, 2002). In the water column, C should be smaller than Cgel. The settling velocity ws is set to be the free settling velocity wsf for salinity !Sfloc and C ! Chind. Suspended sediment acts as a sink for turbulent kinetic energy. This decreases the values of the eddy viscosity and the sediment diffusivity. Thus, for a given shear stress the velocity gradient will be steeper in the stratified flow case and the drag will be reduced (Hill and McCave, 2001). The MonineObhukhov length was used to account for the effects of sediment stratification: Lz Z
U3b r0 ; kgr0 w0
ð13Þ
where U*b is the friction velocity at the bottom, k Z 0.4 is the von Karman constant and r0 w0 is the turbulent density flux: r0 w0 Z KH
vr r Z 0 PB ; vz g
ð14Þ
where PB is the buoyancy production (or sink), g vr PB Z KH : r0 vz Following Businger et al. (1971), the velocity profile in a stratified flow is given by: Ub HCz aðHCzÞ UZ C ; ð15Þ ln z0 Lz k where a is an empirical constant between 4.7 and 5.2. A quadratic friction law at the bed was assumed, tb =r0 ZU2b ZCD U2b ;
ð16Þ
where Ub is the velocity at bottommost grid point zb. Substituting Eq. (16) into Eq. (15) yields a relationship for CD,
752
W.B. Guan et al. / Estuarine, Coastal and Shelf Science 65 (2005) 747e762
1ZAðCD Þ1=2 B=CD ;
ð17Þ
with 1 HCzb AZ ln ; k z0
BZ
aðHCzb ÞPB : U3b
Thus A and B are known and independent of CD. The value of CD was calculated from Eq. (17). For the non-Newtonian fluid-mud mixture, the dependence of the apparent viscosity on the shear rate was calculated following Le Hir et al. (2001), m0 1Cam Cbm r0 1 aB CbB qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiCKIW ; C r0 2 2 Y0 C ðvu=vzÞ Cðvv=vzÞ
KMB Z
ð23Þ ð18Þ
where aB, bB, and Y0 are empirical coefficients, m0 is the dynamic viscosity of clear water, and in the fluid mud layer h 8 i3 < KMB0 1 Rig =Rigc 2 for 0!Rig !Rigc KIW Z 0 ; for Rig RRigc : KMB0 for Rig !0 ð19Þ where KMB0 is a empirical constant, Rig is the gradient Richardson number, 2 2 g vr vu vv Rig Z þ ; ð20Þ r0 vz vz vz and Rigc is the critical gradient Richardson number (Kantha and Clayson, 2000). Without any guidance from the literature, the background diffusivity coefficient KHB was assumed to be constant. The boundary conditions at the closed lateral boundaries are straightforward. The conditions of zero mass, momentum, salt, sediment, and turbulence fluxes must be satisfied at a solid lateral boundary. Along the shoreline where river inflow or outflow may occur, the depth-averaged velocities can be calculated. Assuming that the vertical distributions of velocities u and v are parabolic, they are specified as follows (Chen, 1994): uZ3u
ðHCzÞ2 ; D2
ð21Þ
2
vZ3v
ðHCzÞ ; D2
vertical mean velocity components used in the external (vertically integrated) mode of the dynamic equations. At the head of the river, zero salinity is applied, and C can approximately be specified as vertically uniform. Generally, zero normal gradient conditions are appropriate for q2 and l at the riverine boundary. At the seaward open boundary, when there is an inflow, S and C need to be prescribed. Immediately after low water slack, the salinity of the inflowing water will not be equal to the salinity of seawater. It will take a certain time before the maximum salinity Smax is found at the boundary. This lag effect in the salinity at the boundary was modeled following Thatcher and Harleman (1972) as: bðtÞSmax C½1 bðtÞSðtlws Þ for tlws %t%tlws CTTH SZ ; Smax for tlwsCTTH %t%thws
ð22Þ
where D Z H C h is the total depth of the water column, and u and v denote the depth-averaged velocities, i.e., the
where TTH is the ThatchereHarleman time lag (return time), tlws and thws are the times of the last low water slack and the next high water slack. For time t within tlws and tlws C TTH, the salinity S(t) is obtained by a weighted interpolation of S(tlws) and Smax. The salinity is equal to Smax during the remaining inflow period from t Z tlws C TTH to t Z thws. The time function b(t) is defined by: t tlws 2 bðtÞZsin 0:5p : ð24Þ TTH The following equation was used to calculate the seaward, open boundary conditions for C: ðCL CCH Þ ðCH CL Þ t tlws p C sin 2pmin ;1 CZ ; 2 2 2 TF ð25Þ where TF is the mean duration of flood current, CL is C at low water slack, and CH is C at peak flood current. At the seaward open boundary during an outflow, v v ðS; CÞCUN ðS; CÞZ0; vt vN
ð26Þ
where N denotes the coordinate normal to the boundary, and UN equals u or v according to the normal direction of the boundary. As far as q2 and l are concerned, any advection at the lateral boundaries can be assumed to be zero. For tidal calculations, either the free surface elevation or the tidal current velocity can be formulated. Observations of sea levels are more available than ocean currents. By modifying two classical open boundary conditions proposed by Ippen (1966) and Reid and Bodine (1968), Ly et al. (2000) gave a new formulation.
753
W.B. Guan et al. / Estuarine, Coastal and Shelf Science 65 (2005) 747e762 Table 1 Specification of some model parameters Parameter
Description
Value
Source
DtE
Time step for the external mode
2s
Dt
Time step for the internal mode
6s
r0 CM CH q
Reference density Dimensionless coefficient of horizontal eddy viscosity Dimensionless coefficient of horizontal eddy diffusivity Filter coefficient in the Asselin filter used in the solving the momentum equations Critical value of the Richardson number GH Critical elevation difference chosen so as to prevent the drying and flooding of a grid box at alternate time steps Minimum allowable value of the total water depth to avoid singularities in some relevant terms as a grid box dries Effective roughness length Empirical constant in the velocity profile equation for stratified flow Dynamic viscosity of clear water Empirical coefficient for the apparent viscosity of mudewater mixture Empirical exponent for the apparent viscosity of mudewater mixture Empirical coefficient for the Bingham stress of mudewater mixture Empirical exponent for the Bingham stress of mudewater mixture Empirical constant for the Moore rheological model of mudewater mixture Empirical constant reflecting the effect of internal wave on vertical mixing Critical gradient Richardson number above which the effect of internal wave should be taken into account Background diffusivity coefficient Critical salinity value above which flocculation is enhanced Critical C value above which the hindered settling occurs Free settling velocity Empirical coefficient for the settling velocity of floc Empirical exponent for the settling velocity of floc Empirical exponent accounting for non-linear effects of return flow and wake due to the settling of sediment Sediment concentration at which a space-filling network is formed Critical shear stress for deposition Erosion constant
1012 g l1 0.05 0 0.1
Determined by the CFL condition Determined by the CFL condition Dong et al. (1997) Mellor (2001) Mellor (2001) Mellor (2001)
120 0.02 m
Model calibration Flather and Heaps (1975)
0.2 m
Flather and Heaps (1975)
0.5 mm 4.9 1 ! 103 N s m2 0.2 1 8.7 ! 107 N m2 2.55 1 ! 103 s1 5 ! 104 m2 s1 0.33
Model calibration Businger et al. (1971) Molecular diffusion Ross (1988) Ross (1988) Ross (1988) Ross (1988) Le Hir et al. (2001) Model calibration Geyer and Smith (1987)
5 ! 105 m2 s1 0.01 10 g l1 0.1 mm s1 0.089 mm s1 1.3 1
Model Model Model Model Model Model Model
120 g l1 0.6 N m2 2.5 ! 103 kg m2 s1
Winterwerp (1999) Model calibration Model calibration
GHc e Dcrit z0 a m0 am bm aB bB Y0 KMB0 Rigc KHB Sfloc Chind wsf k1 m1 m2 Cgel tD E
According to Ly et al. (2000), the depth-averaged current normal to the open boundaries has the form: UN ZUNi Gðg=DÞ
1=2
ðhi hs Þ;
ð27Þ
where UNi and hi are the first interior depth-averaged velocity and sea level in the model domain, respectively. If the normal direction of the boundary is the same as the direction of the coordinate axis, Eq. (27) has a positive sign and otherwise a negative sign. The quantity hs is the specified time dependent tidal elevation using a harmonic expansion of the form: hs ðtÞZA0 C
NT X
Am cosðum tC40m 4m Þ;
ð28Þ
m
where A0 represents the residual sea level, Am, um, 40m, and 4m are, respectively, the amplitude, frequency, initial value of the phase umt C 40m and the epoch of the m-th tidal constituent, and NT is the number of tidal constituents.
calibration calibration calibration calibration calibration calibration calibration
The model domain (Fig. 1a) includes the entire Jiaojiang Estuary, the two tributaries, and coastal waters (Taizhou Bay). In the horizontal, an orthogonal curvilinear grid was used, with 102 ! 29 grid elements. In the vertical, 15 unequally spaced sigma-layers were used in the water, and 10 layers in the substrate, with a much higher resolution near the bottom. The bottommost grid was located 0.06 m above the bed. The simulation started at 0 h on April 10, 1991, and ended at 1200 hours on April 22, 1991. The model parameters are given in Table 1. The initial salinity distribution was assumed to be vertically uniform. From Songpuzha (Fig. 1) to the eastern boundary, the salinity was assumed to be constant along all cross-channels and to increase gradually from 0 to 24. For the other regions, the salinity was set to be 0. Initial values of C were 5 g l1. The initial surface elevations and velocities were set to zero everywhere. The stationary sediment in the bed was assumed to be uniformly 30 kg m2 everywhere.
754
W.B. Guan et al. / Estuarine, Coastal and Shelf Science 65 (2005) 747e762
Velocity (m s-1)
(a) 1 0 -1 -2
20
80
100
120
140
160
180
80
100
120
140
160
180
80
100
120
140
160
180
80
100
120
140
160
180
(b)
16 12 8
Suspended sediment concentration (g l-1)
4 0
24
(c)
20 16 12 8 4 0
50
(d)
40 30 20 10 0
Time (h) Fig. 2. Comparison of observed (solid lines) and predicted (dashed lines) (a) longitudinal velocity at 1.75 m above the bed, and the suspended sediment concentration (SSC, in g l1) at (b) 1.67 m, (c) 0.55 m and (d) 0.35 m above the bed, in April 1991 at station M1.
In addition, no wind, evaporation and precipitation are considered at the sea surface. Since the seaward open boundary is parallel to the depth contours, the surface elevation at this boundary is horizontal. The surface elevation was calculated from 11 tidal constituents (Q1, O1, P1, K1, N2, M2, S2, K2, M4, MS4 and M6) and zero residual sea level (A0). The outer salinity Smax was set to be 25 and the return time TTH was set to be 2.5 h all along the seaward open boundary. The time of the last low water slack tlws was recorded by the numerical model for each grid point of the open boundary. Eq. (25) was used to calculate the C boundary conditions during an inflow. At the sea surface, the C at the low water slack CL was set to 0.5 g l1, while the C at the peak flood current CH was set to 1.0 g l1. At the bottom, CL Z 2.0 g l1 and CH Z 4.0 g l1. In April 1991, the discharge of the Lingjiang River ranges from 71 to 954 m3 s1 with a mean value of
225 m3 s1. During the simulations, the instantaneous freshwater discharge was interpolated from these daily gauging data. The discharge in the Yongning River was zero. The salinity of the freshwater was 0, while its C was 0.18 g l1.
3. Results The model was verified against the field data including tidal levels at the T3 station, the longitudinal velocity and suspended sediment concentration (C ) data at station M1 (see Fig. 2 and see the standard experiment in Table 2) and longitudinal velocity, salinity and C data at stations C1 and C2. The currents at 1.75-m elevation, including the tidal asymmetry, were faithfully reproduced by the model. The flood current duration is shorter and the peak current is greater than the
755
W.B. Guan et al. / Estuarine, Coastal and Shelf Science 65 (2005) 747e762 Table 2 Root mean square (RMS) errors and other statistical values at mooring station M1 Elevation
1.75 m
1.67 m
Parameter
RMS error in longitudinal velocity (m s1)
RMS error in C (g l1)
Field data Standard experiment E1 E2 E3 E4 E5 E6
0.00 0.24 0.28 0.24 0.32 0.23 0.24 0.22
0.00 3.58 15.68 3.72 3.50 3.59 3.95 3.71
0.55 m
0.00 4.39 15.54 4.82 4.87 4.41 5.14 4.29
corresponding value at ebb tide. Ebb velocities peak at about 1 m s1, while flood velocity peak at 1.8 m s1. A similar asymmetry is also found, and faithfully reproduced by the model, at stations C1 and C2. Seaward, the tidal current asymmetry is smaller.
0.35 m
0.00 7.23 15.16 7.39 7.71 7.20 7.56 7.00
0.55 m Minimum of C (g l1)
Maximum of C (g l1)
Mean of C (g l1)
Standard deviation of C (g l1)
1.54 2.07 4.84 2.43 2.90 2.08 2.21 2.17
21.52 20.83 37.90 24.44 23.35 24.38 26.34 19.37
7.39 6.02 20.17 6.74 6.99 5.93 7.45 6.10
3.45 3.43 8.67 4.14 4.26 3.39 4.43 3.34
The data at the two lower elevations (0.55 and 0.35 m above the bed) of station M1 (Fig. 2c and d) show five flood tide peaks in C during the first 60 h of the records. These peaks were reproduced successfully by a vertical 2-D model under the assumption of the existence of
(a)
(b)
(c)
Fig. 3. Time series of the predicted (a) longitudinal velocity (m s1), (b) salinity and (c) suspended sediment concentration (SSC, in g l1) at station M1 during a spring tide.
756
W.B. Guan et al. / Estuarine, Coastal and Shelf Science 65 (2005) 747e762
Amount of Sediment (×109 kg)
4
Suspended particulate matter Deposited mud Total sediment mass
3
2
1
0 80
100
120
140
160
180
Volume-averaged SSC Volume-averaged kinetic energy
8
1.6
1.2 6 0.8 4 0.4
2
Kinetic Energy (m2 s-2)
Suspended Sediment Concentration (g l-1)
Time (h)
0 80
100
120
140
160
180
Time (h)
Fig. 4. Time series of sediment amount, volume-averaged SSC and volume-averaged kinetic energy over the estuarine part from Songpuzha to Yongquan.
a limited soft mud layer deposited at the bed (Guan et al., 1998). Since the 3-D model also includes a conceptually similar sediment bed module, it also successfully reproduces these peaks (Fig. 2bed). The comparison between measured and simulated surface elevations for the T3 station at the mouth near Haimen also shows a good agreement between field data and model results, both in amplitude and phase (not shown). The tidal wave at station T3 shows a tidal asymmetry, although a symmetric tidal wave is forced at the nearby seaward open boundary. The model predicts that during a tidal cycle, the isohaline of 1 migrates between Luocun and Songpuza; it also predicts that at station M1 (see Fig. 3) during peak currents the vertically-averaged suspended sediment concentration (C ) commonly reaches 5 g l1 in the turbidity maximum zone, that the near-bottom C peaks at 90 g l1 in the estuary, that at the beginning of the flood tide the values of C are higher than those at peak currents because of low tidal level, that the vertical profiles of velocity near the bed are not logarithmic when C O 20 g l1, and that there is a profound asymmetry in currents. The model also shows that the lutocline rises during accelerating tidal currents and falls during decelerating currents. A lutocline occurs during the whole tidal cycle except around high water slack. Defining fluid mud as the water with C O 20 g l1 (Mehta and Srinivas, 1993), the distribution and
thickness of fluid mud layer are predicted by the model. Accordingly, at the beginning of the flood tide, fluid mud occurs near Songpuzha and in Taizhou Bay, while there is no fluid mud yet in the estuary from Yongquan to Songpuzha. During swift tidal currents, fluid mud occurs also near Songpuzha and generally throughout the estuary. The maximum thickness of the fluid mud layer is about 0.60 m and this occurs near Haimen. All these predictions match well the observations of Guan et al. (1998). Fig. 4 shows that during the accelerating tidal currents the bottom sediment is eroded and entrained into the mobile sediment, to be deposited upstream at flood tide and downstream at ebb tide. There is a significant exchange of sediment between various areas of the estuary. Overall, the total deposited mud increases with time, as result of net advection from coastal waters. The upper limit of the total suspended sediment mass is about 1.2 ! 106 tons, which is comparable to the annual fluvial sediment discharge. During each tide cycle, there are two peaks in volume-averaged kinetic energy (Fig. 4), with the largest peak corresponding to the peak flood current and the smallest peak to the peak ebb current. The time-series plots of predicted, tidally-averaged salinity, currents and C (Figs. 5e7) suggest that there is little net gravitational circulation in the longitudinal direction and that it is the asymmetry in tidal currents
W.B. Guan et al. / Estuarine, Coastal and Shelf Science 65 (2005) 747e762
757
(a)
(b)
(c)
Fig. 5. Horizontal distributions of the predicted, vertically-averaged (a) velocity, (b) salinity and (c) suspended sediment concentration (SSC, in g l1) time-averaged over several tidal cycles.
that drives the outer sediment into the estuary from the coastal waters. The turbidity maximum is found to extend longitudinally about 16 km, and its center is located off Haimen (Figs. 5 and 6). There, cross-channel gradients in currents and C are apparent, with greater values occurring at the southern side than at the opposed side, indicating the existence of cross-channel residual circulation (Fig. 7). For the tidally-averaged mud fluxes the models predict that, off Haimen the northern side of the estuary is flood-dominant, while the southern side is slightly ebb-dominant. Since there is net flood dominance throughout the estuary mouth, mud is pumped into the estuary from coastal waters and is redistributed asymmetrically within the estuary, resulting in the formation of mud banks. This observation is in quantitative agreement with the location of the mud banks and with port authorities dredging logbooks. 4. Discussion Because of the asymmetry of currents, the eddy viscosity and eddy diffusivity are larger at flood than at
ebb tide. Large Richardson numbers occur in the nearbottom layers which have a strong stratification in C (Jianhua and Wolanski, 1998; Jiang and Mehta, 2000, 2002). Mixing occurs as a result of the breaking of internal waves riding on the lutocline. This process is not included in the turbulence closure model due to strong buoyancy damping at the lutocline. Therefore, background mixing coefficients must be added. Without this, the model would predict that the top of the water column would completely clear of sediment at slack water, a prediction at odds with observations (see also Run E6 discussed below). The numerical experiment described in the last section is called a standard experiment. Following these specifications of the standard experiment (Table 1), eight simulations using the 3-D model were completed to explore parameter sensitivity. Table 3 summarizes the objectives and specifically varied parameters in the eight simulations. Discussion on sediment dynamics associated with these simulations including the standard experiment will be expounded sequentially in this section.
758
W.B. Guan et al. / Estuarine, Coastal and Shelf Science 65 (2005) 747e762
(a)
(b)
(c)
Fig. 6. Longitudinal distributions of the time-averaged (a) along-transect velocity, (b) salinity and (c) suspended sediment concentration (SSC, in g l1) over several tidal cycles.
4.1. Buoyancy effects Time series results from Run E1 (no sedimentinduced buoyancy effects) are plotted in Fig. 8, while the corresponding results from Run E2 (no salinityinduced buoyancy effects) are shown in Fig. 9. Without sediment-induced buoyancy effects (Run E2), the predictions are unrealistic. Indeed, the predicted C values are unreasonably high and sometimes are uniformly equal to 30 g l1 in the vertical, which is unrealistic when compared to field observations (Guan et al., 1998). This finding results from the fact that, when there are no sediment-induced buoyancy effects, there is no damping of turbulence in the presence of fluid mud. The differences between Run E2 and the standard experiment are smaller. Nevertheless, without salinityinduced buoyancy effects, the vertical mixing coefficients are also unrealistically high especially in the upper layers, since all weak stratification near the sea surface has been ignored. 4.2. Hydraulic roughness and drag Run E3 was undertaken to demonstrate how sensitive the model is to the specific value of the
effective roughness length. When this length is increased by a factor of 10, the predicted surface elevation at low tide became unrealistically high, and the predicted ebb tidal currents become unrealistically small. The predicted tidal fluctuations of C also became unrealistic. Run E4 investigated the effect of the stratification on the effective hydraulic drag. At peak currents, stratification in C became stronger and thus the overall hydraulic roughness was reduced. However, this impact was moderate.
4.3. Sensitivity of a tuning constant in the turbulence closure model Guided by laboratory data, Mellor (2001) introduced a Richardson number cutoff and thus a tuning closure constant GHc to the turbulence dissipation in the MelloreYamada turbulence closure model. An empirical value of GHc Z 6.0 was adopted. In this case with GHc Z 6.0 (Run E5), the peaks in C during flood tide became unrealistically small. In the standard experiment, GHc Z 120.
759
W.B. Guan et al. / Estuarine, Coastal and Shelf Science 65 (2005) 747e762
(a)
Table 3 List of the additional experiments. The parameters are identical to those in the standard experiment unless otherwise noted Experiment Objective
Parameters different from the standard experiment
E1
rZ1000C0:78S (g l1)
E2
(b) E3
E4
E5
(c)
E6
E7
E8
Fig. 7. Tidally-averaged, cross-channel distributions in a transect past station M1 of the time-averaged (a) velocity (m s1); the longitudinal component is shown as a contour line, and the cross-channel component as a vector, (b) salinity, and (c) suspended sediment concentration (SSC, in g l1).
4.4. Background vertical mixing Generally, the relative, mud-induced increase of the eddy viscosity is much smaller than that for the eddy viscosity, except for flows with a high C and a low shear rate. Although increasing the mud-induced viscosity is necessary to keep the numerical stability of the model when high C values occur, it does not significantly change the overall result. On the other hand, the Richardson-number-dependent coefficient for the background vertical viscosity, which is considered to account for the effects of the internal wave, can influence the model results near the bed. Run E6 demonstrates that the time-averaged C value near the bed is reduced without the background vertical viscosity associated with internal wave. That is to say, internal waves contribute to the vertical transport of momentum and thereby indirectly affect the transport of mass (Uittenbogaard, pers. comm.).
To investigate sediment-induced buoyancy effects To investigate salinity-induced buoyancy effects To demonstrate how sensitive the model is to the effective roughness length To investigate the effect of stratification on the effective hydraulic drag To demonstrate the sensitivity of the tuning constant GHc in the MelloreYamada turbulence closure model To demonstrate the effect of background vertical mixing controlled by internal wave To understand the consequences of a small freshwater discharge of the Lingjiang River. To understand the consequences of a large freshwater discharge of the Lingjiang River.
rZ1000C0:62C (g l1)
z0 Z 5 mm
aZ0
GHc Z 6
KIW Z 0 m2 s1
Freshwater discharge is set to 71 m3 s1 Freshwater discharge is set to 954 m3 s1
4.5. Effects of the river discharge During a low river discharge (Run E7), the total mass of suspended mass decreases, and this results in increased sedimentation. Salty water moves further upstream, and the vertically-averaged and time-mean values of C are reduced. During high river discharge (Run E8), the sediment is gradually carried out of the estuary. The kinetic energy increases during ebb and decreases during flood. The masses of mobile sediment, stationary sediment, and total sediment are reduced, while turbidity maximum moves seaward to Taizhou Bay. 5. Conclusions A 3-D, POM-based model has been developed here for the highly turbid, macro-tidal, Jiaojiang Estuary and adjoining coastal waters. The model was verified against field observations of current and suspended sediment concentration, C. A fine resolution grid near the bottom was necessary to simulate fluid-mud formation. Sediment-induced buoyancy effects are found to be significant. Turbulence is significantly dampened at the high C. Because of
760
W.B. Guan et al. / Estuarine, Coastal and Shelf Science 65 (2005) 747e762
(a)
(b)
(c)
Fig. 8. Same as Fig. 3 for Run E1.
the strong stratification at the lutocline, the Mellore Yamada scheme tends to under-predict vertical mixing, thus requiring background mixing to be added. The MonineObukhov length scale appears successful at calculating the influence of sediment-induced stratification on the effective hydrodynamic drag. The hydraulic roughness is predicted to be measurably decreased in the very turbid estuary. The effect of mud in increasing the viscosity is small but found to be necessary for numerical stability. The TVD scheme avoids all negative values of salinity and C, which would occur frequently otherwise in the presence of large horizontal gradients of salinity and vertical gradients of C. A lutocline occurs almost over the entire tidal cycle, except near high water slack. During accelerating currents, the lutocline rises, and during slack water it falls. The total mass of mobile sediment trapped in the estuary is around 1.2 ! 106 tons, a value comparable to the annual river discharge. Asymmetry in tidal currents generates tidal pumping and estuarine infilling from sediment pumped from coastal waters. During high flow periods, the sediment is
carried to coastal waters; this sediment is advected back in the estuary during low flow periods. The 3-D model appears able to realistically predict the location of fluid mud banks as the areas where nearbottom currents advect mud either towards the deeper regions of the estuary or in coastal waters near the mouth. The simulations are very sensitive to details of the parameterization, suggesting that predictive 3-D modeling of mud dynamics remain an art rather than a hard science, and that improvements in the model requires process-based studies to improve the parameterization of these processes.
Acknowledgments This study was supported by the China Natural Science Foundation under contract No. 40206014, and by the National Key Basic Research Program from the Ministry of Science and Technology of China under contract No. G1999043702. The field data were collected under the Sino-Australia cooperation between the
W.B. Guan et al. / Estuarine, Coastal and Shelf Science 65 (2005) 747e762
761
(a)
(b)
(c)
Fig. 9. Same as Fig. 3 for Run E2.
Second Institute of Oceanography, People’s Republic of China and the Australian Institute of Marine Science. It is a pleasure to thank Dr. A. Mehta for editing this paper. Two anonymous reviewers helped improve this paper.
References Ariathurai, R., 1974. A finite element model for sediment transport in estuaries. Ph.D. thesis, University of California, Davis, 192 pp. Bi, A.H., Sun, Z.L., 1984. A preliminary study on the estuarine processes in the Jiaojiang River, China. Journal of Sediment Research 3, 12e26 (in Chinese). Blumberg, A.F., Galperin, B., O’Connor, D.J., 1992. Modeling vertical structure of open channel flows. Journal of Hydraulic Engineering 118, 1119e1134. Blumberg, A.F., Mellor, G.L., 1987. A description of a three-dimensional coastal ocean circulation model. In: Heaps, N. (Ed.), ThreeDimensional Coastal Ocean Models. Coastal and Estuarine Science, vol. 4. American Geophysical Union, Washington, DC, pp. 1e16. Businger, J.A., Wyngaard, J.C., Izumi, Y., Bradley, E.F., 1971. Fluxe profile relationship in atmospheric surface layer. Journal of Atmospheric Science 28, 190e201.
Chen, X., 1994. Effects of hydrodynamics and sediment transport process on nutrient dynamics in shallow lakes and estuaries. Ph.D. Dissertation, University of Florida. Dong, L.X., Wolanski, E., Li, Y., 1997. Field and modeling studies of fine sediment dynamics in the extremely turbid Jiaojiang River estuary, China. Journal of Coastal Research 13, 995e1003. Faas, R.W., 1984. Time and density-dependent properties of fluid mud suspensions, NE Brazilian Continental Shelf. Geo-Marine Letters 4, 147e152. Flather, R.A., Heaps, N.S., 1975. Tidal computations for Morecambe Bay. Geophysical Journal of the Royal Astronomical Society 42, 489. Fu, N.P., Bi, A.H., 1989. Discussion on the problems in the sediment transport of the Jiaojiang River. Journal of Sediment Research 8, 51e57 (in Chinese). Geyer, W.R., Smith, J.D., 1987. Shear instability in a highly stratified estuary. Journal of Physical Oceanography 27, 1668e1679. Guan, W.B., Wolanski, E., Dong, L.X., 1998. Cohesive sediment transport in the Jiaojiang River estuary, China. Estuarine, Coastal and Shelf Science 46, 861e871. Hill, P.S., McCave, I.N., 2001. Suspended particle transport in benthic boundary layers. In: Boudreau, B., Jorgensen, B.B. (Eds.), The Benthic Boundary Layer: Transport Processes and Biogeochemistry. Oxford Press, pp. 78e103. Ippen, A.T., 1966. Estuary and Coastline Hydrodynamics. McGrawHill, New York, 744 pp.
762
W.B. Guan et al. / Estuarine, Coastal and Shelf Science 65 (2005) 747e762
James, I.D., 1996. Advection schemes for shelf sea model. Journal of Marine Systems 8, 237e254. Jiang, J., Mehta, A.J., 2000. Lutocline behavior in a high-concentration estuary. Journal of Waterway, Port, Coastal and Ocean Engineering 126 (6), 324e328. Jiang, J., Mehta, A.J., 2002. Interfacial instabilities at the lutocline in the Jiaojiang Estuary, China. In: Winterwerp, J.C., Kranenburg, C. (Eds.), Fine Sediment Dynamics in the Marine Environment. Elsevier, Amsterdam, pp. 125e137. Jianhua, J., Wolanski, E., 1998. Vertical mixing by internal wave breaking at the lutocline, Jiaojiang estuary, China. Journal of Coastal Research 14, 1426e1431. Kantha, L.H., Clayson, C.A., 2000. Numerical Models of Oceans and Oceanic Processes. Academic Press, San Diego, 940 pp. Krone, P.B., 1962. Flume Studies of the Transport in Estuarine Shoaling Processes. Final Report, Hydraulics Engineering Laboratory, University of Berkeley, California, USA, 110 pp. Le Hir, P., 1997. Fluid and sediment ‘‘integrated’’ modeling: application to fluid mud flows in estuaries. In: Burt, N., Parker, R., Watts, J. (Eds.), Cohesive Sediment. John Wiley, Chichester, UK, pp. 417e428. Le Hir, P., Bassoullet, P., Jestin, H., 2001. Application of the continuous modeling concept to simulate high-concentration suspended sediment in a macrotidal estuary. In: McAnally, W.H., Mehta, A.J. (Eds.), Coastal and Estuarine Fine Sediment Processes. Proceedings in Marine Science, vol. 3. Elsevier Science, Amsterdam, pp. 229e248. Le Hir, P., Cayocca, F., 2002. 3D application of the continuous modeling concept to mud slides in open sea. In: Winterwerp, J.C., Kranenburg, C. (Eds.), Fine Sediment Dynamics in the Marine Environment. Proceedings in Marine Science, vol. 5. Elsevier Science, Amsterdam, pp. 545e561. Le Normant, C., Peltier, E., Teisson, C., 1998. Three dimensional modeling of cohesive sediment transport in estuaries. In: Dronkers, J., Scheffer, M.B.A.M. (Eds.), Physics of Estuaries and Coastal Sea. Balkema, Rotterdam, pp. 65e71. Li, Y., Wolanski, E., Xie, Q., 1993. Coagulation and settling of suspended sediment in the Jiaojiang estuary, China. Journal of Coastal Research 9, 390e402. Luyten, P.J., Jones, J.E., Proctor, R., Tabor, A., Tett, P., Wild-Allen, K., 1999. COHERENS e A Coupled HydrodynamicaleEcological Model for Regional and Shelf Seas: User Documentation. MUMM Report, Management Unit of the Mathematical Models of the North Sea, 914 pp. Ly, L.N., O’Connor, W.P., Paduan, J.D., 2000. A barotropic simulation of semi-diurnal and diurnal tides for the Monterey Bay region. In: Blain, W.R., Brebbia, C.A. (Eds.), Hydraulic Engineering Software VIII. Proceedings of the Eighth International Conference on Hydraulic Engineering Software. Witpress, Southampton, Boston, pp. 155e169. Mehta, A.J., Srinivas, R., 1993. Observations on the entrainment of fluid mud in shear flow. In: Mehta, A.J. (Ed.), Nearshore Estuarine Cohesive Sediment Transport. American Geophysical Union, Washington, DC. Mellor, G.L., 2001. One-dimensional, ocean surface layer modeling, a problem and a solution. Journal of Physical Oceanography 31, 790e809. Mellor, G.L., Hakkinen, S., Ezer, T., Patchen, R., 2001. A generalization of a sigma coordinate ocean model and an intercomparison of model vertical grids. In: Pinardi, N., Woods, J.D. (Eds.), Ocean Forecasting: Conceptual Basis and Applications. Springer-Verlag, pp. 55e72.
Mellor, G.L., Yamada, T., 1982. Development of turbulence closure model for geophysical fluid problems. Reviews of Geophysics and Space Physics 20, 851. Odd, N.V.M., 1988. Mathematical modeling of mud transport in estuaries. In: Dronkers, J., van Leussen, W. (Eds.), Physical Processes in Estuaries. Springer Verlag, pp. 503e531. Odd, N.V.M., Bentley, M.A., Waters, C.B., 1992. Observations and analysis of the movement of fluid mud in estuaries. In: Mehta, A.J. (Ed.), Nearshore and Estuarine Cohesive Sediment Transport. Coastal and Estuarine Studies, vol. 42. AGU, Washington, pp. 430e446. Odd, N.V.M., Cooper, A.J., 1989. A two dimensional of the movement of fluid mud in a high energy turbid estuary. Journal of Coastal Research 5, 185e193. Owen, M.W., 1975. Erosion of Avonmouth Mud. Report No. INT 150, Hydraulics Research Station, Wallingford, United Kingdom. Partheniades, E., 1965. Erosion and deposition of cohesive soils. Journal of the Hydraulic Division, 105e139 (ASCE 91, No. HY1). Reid, R.O., Bodine, B.R., 1968. Numerical model for storm surges in Galveston Bay. ASCE, Journal of Waterways and Harbors Division 94, 33e57. Roberts, W., 1993. Development of a mathematical model of fluid mud in the coastal zone. Proceedings of the Institution of Civil Engineers e Water, Maritime and Energy 101, 173e181. Ross, M.A., 1988. Vertical structure of estuarine fine sediment suspensions. Ph.D. Thesis, University of Florida, Gainesville, 206 pp. Smagorinsky, J., Manabe, S., Holloway, J.L., 1965. Numerical results from a nine-level general circulation model of the atmosphere. Monthly Weather Review 93, 727e768. Teisson, C., Simonin, O., Galland, J.C., Laurence, D., October 1992. Turbulence and mud sedimentation: a Reynolds stress model and a two-phase flow model. In: Proceedings of the 23rd International Conference on Coastal Engineering. American Society of Civil Engineers, Venice, Italy, pp. 2853e2866. Thatcher, M.L., Harleman, R.F., 1972. A Mathematical Model for the Prediction of Unsteady Salinity Intrusion in the Estuaries, Massachusetts Institute of Technology, Report No. 144. Vinzon, S.B., Mehta, A.J., 2001. Boundary layer effects due to suspended sediment in the Amazon estuary. In: McAnally, W.H., Mehta, A.J. (Eds.), Coastal and Estuarine Fine Sediment Processes. Proceedings in Marine Science, vol. 3. Elsevier Science, Amsterdam, pp. 359e372. Winterwerp, J.C., 1999. On the dynamics of high-concentrated mud suspensions. Doctoral Thesis for the Technical University of Delft. Winterwerp, J.C., 2002. On the flocculation and settling velocity of estuarine mud. Continental Shelf Research 22, 1339e 1360. Winterwerp, J.C., Uittenbogaard, R.E., de Kok, J.M., 2001. Rapid siltation from saturated mud suspensions. In: McAnally, W.H., Mehta, A.J. (Eds.), Coastal and Estuarine Fine Sediment Processes. Proceedings in Marine Science, vol. 3. Elsevier Science, Amsterdam, pp. 1e18. Wolanski, E., Chappell, J.R., Ridd, P., Vertessy, R., 1988. Fluidization of mud in estuaries. Journal of Geophysical Research 93, 1535e 1545. Zhu, Y.K., 1986. Some characteristics of stream-like macro-tidal estuary (Jiaojiang). Geographical Study 5 (1), 21e29 (in Chinese).