Modeling bed-load transport of coarse sediments in the Great Bay Estuary, New Hampshire

Estuarine, Coastal and Shelf Science 58 (2003) 937e950

Modeling bed-load transport of coarse sediments in the Great Bay Estuary, New Hampshire A. Bilgilia,*, M.R. Swiftb, D.R. Lyncha, J.T.C. Ipa a Dartmouth College, Thayer School of Engineering, 8000 Cummings Hall, Hanover, NH 03755, USA University of New Hampshire, Ocean and Mechanical Engineering Departments, Durham, NH 03824, USA

b

Received 24 April 2003; accepted 7 July 2003

Abstract Current, sea level and bed-load transport are investigated in the Lower Piscataqua River section of the Great Bay Estuary, New Hampshire, USAda well-mixed and geometrically complex system with low freshwater input, having main channel tidal currents ranging between 0.5 and 2 m s
1. Current and sea level forced by the M2M4M6 tides at the estuarine mouth are simulated by a vertically averaged, non-linear, time-stepping finite element model. The hydrodynamic model uses a fixed boundary computational domain and accounts for floodingedrying of tidal flats by making use of a groundwater component. Inertia terms are neglected in comparison with pressure gradient and bottom friction terms, which is consistent with the observed principal dynamic balance for this section of the system. The accuracy of hydrodynamic predictions in the study area is demonstrated by comparison with four tidal elevation stations and two cross-section averaged current measurements. Simulated current is then used to model bed-load transport in the vicinity of a rapidly growing shoal located in the main channel of the lower system. Consisting of coarse sand and gravel, the shoal must be dredged every five to eight years. Two approaches are takendan Eulerian parametric method in which nodal bed-load flux vectors are averaged over the tidal cycle and a Lagrangian particle tracking approach in which a finite number of sediment particles are released and tracked. Both methods yield pathways and accumulations in agreement with the observed shoal formation and the long-term rate of sediment accumulation in the shoal area. Ó 2003 Elsevier Ltd. All rights reserved. Keywords: bedload; sediment transport; numerical model; vertically integrated; tidal currents; tidal estuaries; sand waves; New Hampshire coast

1. Introduction Bed-load sediment transport in estuaries is an important process that has the potential of causing negative long-term navigational and environmental effects. For cases where the tidal range is on the order of the mean depth, the net transport in estuaries is dominated by asymmetries in tidal currents that the strongly non-linear hydrodynamics and their complex interactions with basin geometry can create. These asymmetries can export (ebb dominant) or import (flood dominant) sediment from/to the domain (Harris and Collins, 1988; Larcombe and Ridd, 1996; Larcombe and Jago, 1996; Bryce et al., 1998), but may also cause * Corresponding author. E-mail address: [email protected] (A. Bilgili). 0272-7714/03/$ - see front matter Ó 2003 Elsevier Ltd. All rights reserved. doi:10.1016/j.ecss.2003.07.007

deposition and erosion in the form of converging and, respectively, diverging bed-load flux vectors (Ludwick, 1975; Bilgili et al., 1996; Signell et al., 2000). Advance knowledge of these bed-load processes is crucial for planners and scientists who need to make sound decisions concerning the long-term environmental management of both channels requiring dredging and dredged material. Also of concern are situations requiring immediate action like the dredging of navigational hazards (Bilgili, 1993) or environmentally damaging changes in the flow dynamics caused by depositional zones in narrow and shallow estuaries (Cxelikkol et al., 2002). The Great Bay Estuarine System, in New Hampshire, USA (Fig. 1), is a relatively shallow estuary with a tidal amplitude to mean depth ratio of 0.18. The typical tidal range is on the order of 2.5 m at the Gulf of Maine mouth and decreases to about 1.8 m in inner estuarine

938

A. Bilgili et al. / Estuarine, Coastal and Shelf Science 58 (2003) 937e950

Dover Salmon Falls River Cocheco River

NEW HAMPSHIRE

43d10min

74542.69 m Upper Piscataqua River

Bellamy River

MAINE

South Berwick

Durham D t in Po er ov

Oyster River

T-14

C-124

Adams Point

Little Bay

Frankfort Island T-14A SHOAL AREA

Newington

Furber Strait

43d05min

T-13 T-12

Eliot

DISPOSAL AREA

Spruce Creek

Boiling Rock

C-119 Lower Piscataqua River

Kittery

65136.65 m Newmarket

NEW HAMPSHIRE

Great Bay Lamprey River

Squamscott River

Greenland

Winnicut River

70d50min 367856.80 m

Portsmouth

Sagamore Creek

Portsmouth Harbor Gerrish Island New y Castle ar Island nd u o B n pe O el Gulf od of Maine M Odiornes 70d40min Point 381317.79 m

Fig. 1. Map of the Great Bay Estuarine System. The locations of the four tidal stations (T-) and two cross-sectionally averaged velocity transects (C-) from the Swift and Brown (1983) 1975 study are also shown. The coordinates are given in degrees and minutes, followed by the New Hampshire State Plane coordinate system which are in meters.

locations. The Gulf of Maine connects to the inner estuary by way of the Lower Piscataqua River, whose channel depth is on the order of 15 m with maximum currents ranging between 0.5 m s
1 and 2.0 m s
1. The inner estuary, consisting of Little Bay and Great Bay proper, has depths on the order of 10 m and currents on the order of 0.5 m s
1. The vertical variability of tidal currents is shown to be uniform by the field work of Swenson et al. (1977). Several freshwater tributaries exist in this section, but their overall input to the system is low, representing only 1% or less of the tidal prism under normal conditions. This makes the Great Bay Estuary a tidally dominated, well-mixed system. Extensive tidal flats in this part of the estuary cause more than 50% of the system’s surface to be exposed as mud flats at low tide. Swift and Brown (1983) showed that the M2 tidal constituent is dominant at least by an order of magnitude over the other semi-diurnal and diurnal constituents in the Great Bay Estuary. They also showed that the principal force balance is between the pressure gradient and total bottom stress, a characteristic also observed in other shallow and narrow tidal embayments (Friedrichs et al., 1992). They suggest that the narrow channel at Dover Point, which connects the Lower

Piscataqua River, Upper Piscataqua River and Little Bay acts as a hydraulic choke and separates the Great Bay into two different sections with contrasting dynamic regimes: a more dissipative lower section that extends from the mouth at the Gulf of Maine to Dover Point, which has a progressive wave characteristic; and a less dissipative upper section that extends from Dover Point into Great Bay, which has a standing wave characteristic. The Lower Piscataqua River is subjected to a high rate of bed-load transport and resulting shoaling in the navigational channel, under the combined effect of geometry, high current regime, strongly non-linear physics and coarse bottom sediment availability. This shoaling occurs in the form of three large sand waves (Fig. 2) with a height of about 3 m and a wavelength of about 49 m. As seen in Fig. 2, they are not perpendicular to the channel axis but rather have a bias towards the southwest (New Hampshire) side of the river. Historically, sedimentation rates have been about 0.3 m year
1 for the sand wave growth. The shoal has required frequent dredging operations on a cyclic basis to deepen the main channel to allow port access to terminals on the New Hampshire side of the river. During the dredging operation that took place in December of

A. Bilgili et al. / Estuarine, Coastal and Shelf Science 58 (2003) 937e950

939

5m 10m 13m 10m13m 5m

13m

13m 13m

18m

Fig. 2. Bed-load transport transects used in this study shown on top of the bottom features of the study area as derived from a side-scan sonar survey. A side-scan sonar picture of the actual sand waves is also shown in the inset (Bilgili et al., 1996).

1991, 15,300 m3 of dredged material was disposed of at a downriver disposal site shown in Fig. 2. The choice of an adjacent downriver riverine disposal site has lead to concerns about the fate of the disposed sediment and these concerns were addressed by many studies which showed that the sediment disposed of at the disposal site is most likely to be redeposited back at the dredge site (Bilgili et al., 1996; Clere, 1993; New Hampshire Fish and Game Department, Marine Fisheries Division, 1983; Maine Geological Survey Department, 1988). During the last dredging operation that took place in November 2000, 10,640 m3 of coarse bottom material was removed from an area of 18,000 m2 (Marconi, 2001). This corresponds to a yearly accumulation of 0.065 m year
1, assuming that the sediment is uniformly spread over the dredge area. In this study, the authors concentrate on the application of a depth-averaged, non-linear, time-stepping, kinematic finite element model with wetting and drying capabilities to the Great Bay Estuary with the purpose of investigating the bed-load sediment transport patterns in the Lower Piscataqua River. Hydrodynamic model performance in the study area is evaluated through comparison with tidal predictions based on a suite of surface elevation and cross-section averaged current observations from the 1975 Great Bay study described by Swift and Brown (1983). Bed-load sediment transport is modeled through the use of Eulerian (parametric bed-load sediment transport models) and Lagrangian

(particle tracking) approaches using currents predicted by the numerical model. While the Eulerian approach provides both quantitative and qualitative results, the Lagrangian approach concentrates on particle trajectories only in order to investigate overall sediment movement patterns in the study area. Assuming that the sediment is available to the flow, bed-load transport modeling performance is then evaluated through comparison with historic dredging data and with results of a field study described in Bilgili et al. (1996).

2. Numerical model description We adopt the fixed-domain, kinematic (that is, without local and advective accelerations), finite element model, ADAM, developed by Ip et al. (1998) to simulate the tides of the Great Bay Estuary. The absence of the accelerations is justified by the aforementioned finding that the force balance is dominated by the bottom stress and pressure gradient. This assumption is strengthened by the dynamical scaling work of McLaughlin (2001), which showed that this particular force balance holds true during high flow periods. The inertia was found to play a role near slack water only, affecting the timing of the shift and hence, the tidal phase. Our numerical model offers efficient and high-resolution numerics for simulating physically realistic flooding and dewatering of tidal flats in shallow embayments. In this approach

940

A. Bilgili et al. / Estuarine, Coastal and Shelf Science 58 (2003) 937e950

the water column is modeled as a composite system of open channel with a porous medium sub-layer between it and the bedrock, as shown in Fig. 3. We neglect Coriolis accelerations because of the relatively small length scales, and the use of a 2-D depth-averaged model is justified because of the observations that show the estuary to be vertically well-mixed. Thus the depthaveraged kinematic equations are employed to model the flow in the open channel. Flow in the sub-bottom porous medium is described by Darcy’s law. With the aforementioned assumptions, the governing equations are: Continuity:

for the porous medium sub-layer, the governing equations for the composite system are reduced to a scalar nonlinear diffusion equation of the form: S

vz
V  DVz ¼ 0; vt

in which the storage coefficient S and the non-linear diffusion coefficient D are specified as functions of H. For our system:
S¼

ð1Þ D¼

Horizontal momentum: ð2Þ

where, H is the total depth of the water column, v the depth-averaged velocity, t time, V the horizontal gradient differential operator, z the surface elevation relative to a horizontal datum, g the acceleration of gravity and Cd is the bottom drag coefficient. For better approximating a realistic bottom, we adopted the Manning approach in parameterizing the depth dependent bottom friction in the form: 2

Cd ¼

gn ; H1=3

1 if H > h0 ; 3 if H%h0

ð5Þ

and

vH þ V  Hv ¼ 0: vt

gHVz þ Cd jvjv ¼ 0:

ð4Þ

ð3Þ

where n is the Manning’s bottom roughness coefficient. Substituting the kinematic momentum balance into continuity and incorporating the Darcian description

8 < :

kh0 þ ðH
h0 Þ

3=2



g Cd jVzj

1=2

kH

if H > h0 ;

ð6Þ

if H%h0

where H ¼ h þ h0 þ z is the total depth of the water column, h the bathymetric depth, and h0 is the thickness of the porous medium with porosity 3 and hydraulic conductivity k. This non-linear diffusion equation (Eq. (4)) is discretized on linear triangles using the Galerkin weighted residual method in space. Conventional, implicit finite differences are applied in time: S

kþq



   zkþ1
zk
VDkþq V qzkþ1 þ ð1
qÞzk ¼ 0; Dt

ð7Þ

with time level k and temporal weighting factor, 0!q%1. The sea level solution is obtained iteratively in each time step, as described in Ip et al. (1998).

Fig. 3. Schematic geometry of the two-media system used by ADAM (Ip et al., 1998).

A. Bilgili et al. / Estuarine, Coastal and Shelf Science 58 (2003) 937e950

3. Computational setup The finite element domain boundary corresponding to the mean high water (MHW) datum is digitized from the United States Geological Survey (USGS) topographical maps. Tributary rivers in Great Bay are excluded because of their minimal effect to the overall estuarine circulation. A minimum element inside angle restriction of 33( and a maximum element area restriction of 4500 m2 are imposed during the grid generation. The resulting finite element grid consists of 22,140 nodes and 39,617 elements. The grid resolution, d, estipffiffiffiffiffiffiffiffi mated by d ¼ 2Ae , with Ae being the element area, ranges between 4.72 m and 204.66 m, with an average of 47.91 m. A sensitivity study applied to element size revealed that further refinement did not affect present results (Bilgili, 2000). The reader is referred to Ertu¨rk et al. (2002) for more information on the finite element grid, including detailed snapshots at four key areas. The nodal bathymetries of the grid are interpolated from a dataset digitized from the US National Oceanic and Atmospheric Administration (NOAA) navigational charts of the domain, supplemented by measurements in areas where NOAA data were not available (Bilgili, 2000). The depths range between
1.17 m and 24.69 m, with 0 defining the mean water level (MWL) datum. Ertu¨rk et al. (2002) give drawings of bathymetric cross-sections shown at six selective locations throughout the domain. To take into account the effects of changing depth, the Manning’s coefficient, n, is specified as a linear function of instantaneous total depth, H, according to n ¼ 0:040
0:000492H;

ð8Þ

in which the intercept and slope are tuned for both elevation and velocity match with field data. This Manning’s description yields smaller errors than the constant n ¼ 0:035 value used in the modeling study by Ertu¨rk et al. (2002), which concentrated on the tidal hydrodynamics of the Great Bay Estuarine System. The porous medium thickness is set to 1 m and the porous medium hydraulic conductivity to 3:162!10
4 m s
1 . The open channel hydrodynamic solutions are found to be not sensitive to the groundwater solution and sensitivity analyses done by doubling these two parameters yielded no noticeable difference (!1%) in the open channel results. A forcing consisting of the M2, M4 and M6 constituents of tides is specified as a Dirichlet elevation boundary condition across the open ocean boundary at the Gulf of Maine entrance, shown in Fig. 1. The elevation time series used to force the model is predicted using the following amplitude/Greenwich phase couples given, respectively, for M2, M4 and M6 constituents: (1.3012 m, 102(G), (0.014 m, 347(G) and (0.0046 m, 149(G). These are derived from two offshore stations

941

near Boon Island, Maine (70.410(W, 43.120(N) and Hampton, New Hampshire (70.785(, 42.900(N) using a weighted average based on distance. The simulation is started with fluid at rest at high tide with a time step of 27.945 s and four non-linear iterations per time step, and terminated after six M2 cycles, well after dynamical equilibrium is established. 4. Hydrodynamical model results The model’s ability to simulate the tidal conditions is already evaluated in great detail in Ertu¨rk et al. (2002) for the Great Bay Estuarine System. However, considering that a different Manning’s coefficient is used for the simulations presented in this paper and for completeness, we still evaluated the hydrodynamics results by comparisons with hindcast time series for tidal elevation stations (T-) and cross-section averaged current stations (C-), shown in Fig. 1, based on the M2, M4 and M6 harmonic constants given by Swift and Brown (1983). As shown in Fig. 4, the agreement in elevation amplitudes is very good but the model phases appear to lead the hindcast phase slightly, starting at station T-13. At station T-14, which is located at the bend in Dover Point, the model starts to under predict the tidal amplitudes, with the model phase still leading the hindcast phase. This mismatch, which also occurs at Little Bay and Great Bay proper, is possibly a cumulative effect that is due to non-linearities, introduced by bottom stress and geometry, and the lack of local accelerations. The latter are non-negligible terms of the standing wave solution for estuarine tides and have an effect around the slack water, as discussed in Ertu¨rk et al. (2002) and McLaughlin (2001). Table 1 summarizes the tidal elevation comparisons in the form of various error norms. The average L2 error is 0.0832 m for the study zone.

Fig. 4. Measured (dashed line) and simulated (solid line) tidal elevation time series at four Swift and Brown (1983) stations.

942

A. Bilgili et al. / Estuarine, Coastal and Shelf Science 58 (2003) 937e950

Table 1 Error statistics at elevation stations (T-) and current cross-sections (C-) in MKS units. C is the correlation coefficient, LN is the absolute maximum error, L2 is the RMS error between measured and model time series, nL2 is the RMS error normalized by the RMS of the observations and S (skill) is defined as S ¼ 1
ðnL2 Þ2 Station

C

LN

T-12 T-13 T-14A T-14

0.998 0.992 0.989 0.992

0.07 0.17 0.18 0.13

C-119 C-124

0.980 0.958

0.33 m s
1 0.80 m s
1

L2 m m m m

0.050 0.091 0.103 0.089

m m m m

0.111 m s
1 0.322 m s
1

nL2

S

0.070 0.133 0.154 0.132

0.995 0.982 0.976 0.982

0.218 0.299

0.952 0.910

The time series of simulated cross-section averaged currents, shown in Fig. 5, are calculated by spatially averaging the longitudinal component of current across transects C-119 and C-124 (Fig. 1). Overall, model currents capture the correct shape and amplitudes with a slightly leading phase. At station C-124, the current amplitudes are very slightly under predicted. Phase errors and small differences could be due to local conditions, lack of accelerations, as well as the model’s ability to resolve flow details. Nevertheless, the agreement between predicted and observed time series is good. Table 1 summarizes current time series comparison errors. The average L2 error is 0.2165 m s
1 for the study area. The large LN of 0.80 m s
1 at cross-section C-124 is mostly due to the phasing error between the model and measured time series. This is thought to be the result of lacking inertia terms, which play a more important role when the current levels are down. This, however, occurs around slack water and therefore affects the bed-load, which is the principal coverage of this paper, should be minimal.

5. Eulerian bed-load transport modeling Eulerian bed-load transport potential is calculated using the sediment-specific parametric models of Bag-

nold (BA) (Bagnold, 1966, 1973), BrowneEinstein (BE) (Bogardi, 1974; Einstein, 1942; Brown, 1950), Engelund and Fredsœ (EF) (Engelund and Fredsœ, 1976; Fredsœ and Deigaard, 1992) and Meyer-Peter and Mu¨ller (1948) (MPM). At the present time, there is no agreement as to which model is best, so several models are used for comparative purposes, as suggested by Nakato (1990) and Vanoni (1984). Assuming that the sediment is in continual supply, the models calculate instantaneous bed-load flux, defined as dry mass transport per time per distance normal to the flux vector, throughout the domain at each time step using nodal velocities predicted by the hydrodynamic simulation. The time step used in these calculations is the same as the hydrodynamic model, which is 27.945 s. Fig. 6 displays one tidal cycle long time series of tidal current and bedload discharges, calculated using various models at an element located in the main channel just downriver of transect 2. This shows that flood currents, given by negative values, have higher peaks and shorter durations than ebb currents, a trend shared by many other main channel elements. This figure also shows the BE model’s sensitivity to peak currents. The calculated bed-load fluxes are then averaged over an M2 cycle to yield Eulerian residual bed-load vectors. These residuals are qualitatively analyzed and related to geological features throughout the Lower Piscataqua River region. The reader is referred to Appendix A for brief descriptions of the parametric bed-load models. Sediment budget studies are performed over an M2 cycle using calculated bed-load at six transects. Three of these transects define the upriver and downriver boundaries of the disposal and shoal sites (Fig. 2), as described in Bilgili et al. (1996). The fourth and fifth transects are chosen in the immediate upriver and downriver vicinity of the sand waves and the sixth transect is chosen just south of the Boiling Rock constriction to investigate its effects on the disposal area (Fig. 2). As shown in Fig. 2, the transects are limited to the main channel where the bed-load

Fig. 5. Measured (dashed line) and simulated (solid line) cross-sectional velocity time series at two Swift and Brown (1983) cross-sections. Relationship between measured vertically averaged and cross-section averaged currents was established by making use of tidal prism and crosssection area considerations. The estimated uncertainty of the measured cross-section averaged currents is G10%.

943

A. Bilgili et al. / Estuarine, Coastal and Shelf Science 58 (2003) 937e950

Tidal Current Speed and Bed-Load Discharge Time Series 0.75

1.2

Tidal Current

0.8

Bagnold

0.6

Brown-Einstein Engelund & Fredsoe

0.4

Meyer-Peter & Muller

0.5

0.25

0.2 0

0

-0.2 -0.4

-0.25

-0.6 -0.8

Bed-Load Discharge (kg/m.sec)

Current Speed (m/sec)

1

-0.5

-1 -1.2

-0.75 62

64

66

68

70

72

74

Hours in 6th Tidal Cycle

Fig. 6. Time series of tidal current velocities and bed-load discharges as calculated by the parametric models. Positive and negative values correspond to ebb and flood currents and discharges, respectively.

transport is concentrated. The budgets are then used to predict yearly accumulations (or erosions) (YA) in kilograms for the shoal (between transects 2 and 3), disposal (between transects 1 and 2) and sand wave (between transects 4 and 5) sites, according to: Z T
Z LLW Z LSW YA ¼ 705  qLW dL
qSW dL dt; ð9Þ 0

0

0

where 705 is the approximate number of M2 cycles per year, T the M2 tidal period in seconds (z44712), LLW and LSW the lengths of the landward and seaward transects bounding the site in question (m), qLW and qSW the space and time varying bed-load fluxes along landward and seaward transects simulated by parametric bed-load models (kg/m s), L the length and t is time. Average sedimentation rates (ASR) in meters per year, are given by: ASR ¼ YA=rp A;

ð10Þ

where rp is the sediment density including porosity (kg/ m3) and A is the site area (m2), and are also calculated

for all sites and compared with field data. A constant total stress to skin friction stress ratio of 1/6 is used to be consistent with previous Great Bay bed-load transport work (Bilgili et al., 1996). This choice of a constant is based on the work by Armstrong (1974), which states that the main channel of the Piscataqua River consists mostly of coarse and unconsolidated sediments. These experience very long bi-directional excursions every tidal cycle under the influence of high velocities, so that roughness is expected to be fairly uniform along the main channel. Table 2 summarizes the parameters used in this study in coordination with the bed-load transport models mentioned above. Figs. 7e9 show residual bed-load flux vectors around the shoal site as calculated by the BA, BE and EF models, respectively. The vectors calculated by the MPM model are not shown because of their similarity to the ones calculated by the EF model, with only minor differences in amplitude. As shown in Fig. 7, the BA model predicts an overall upriver directed residual transport in the center channel. A stagnation and even

Table 2 Sediment input parameters for parametric bed-load models

Characteristic sediment diameter, d50 (mm)

Maximum sediment diameter, dg (mm)

Sediment density, rs (kg m
3)

Sediment density, including porosity, rp (kg m
3)

1

4.8

2800

1750

Bagnold’s efficiency factor, eb

BA model critical velocity after Levi (Bogardi, 1974) (m s
1)

Sand dynamic friction angle, c (degrees)

Total to skin stress ratio, Rs

Critical shear stress for EF and MPM models, tC (non-dimensional)

0.12

0.66

27

1/6

0.05

944

A. Bilgili et al. / Estuarine, Coastal and Shelf Science 58 (2003) 937e950 4

4

x 10

x 10

6.8

6.8

5 6.78

13

5 6.78

10

13

Transect 3 Transect 4

13

13

10

Transect 3 Transect 4

13

13

6.76

6.76

Transect 5

Transect 5

Transect 2

Transect 2 6.74

6.74

10

10

17 13

17 13

6.72

6.72

0.1 kg/m

1 kg/m

5

5

6.7

6.7 3.708

3.71

3.712

3.714

3.716

3.708

3.718

3.71

3.712

3.714

3.716

3.718 5

x 105

x 10

Fig. 7. Residual bed-load transport vectors in the Lower Piscataqua River region as calculated by the Bagnold model. X and Y axis coordinates are in NH state plane coordinate system in meters. Bedload transects 3e5 and 2 are also shown together with 5, 10, 13 and 17 m isobaths. Note different vector scales in Figs. 7e9.

Fig. 8. Residual bed-load transport vectors in the Lower Piscataqua River region as calculated by the BrowneEinstein model. X and Y axis coordinates are in NH state plane coordinate system in meters. Bedload transects 3e5 and 2 are also shown together with 5, 10, 13 and 17 m isobaths. Note different vector scales in Figs. 7e9.

a small reversal of bed-load transport is observed at the shoal site, contributing to the formation of the sand waves. The residual transport is from the disposal area towards the shoal area, a prediction which agrees with the field results of Bilgili et al. (1996). Similarly, the BE model shows an upriver directed overall residual transport. This upriver transport is reversed just upriver of the shoal, as a number of downriver directed vectors show in Fig. 8. This local reversal, combined with two small amplitude gyres on the shoal contribute to sediment deposition predicted by the BE model. Fig. 8 also shows the residual transport direction to be from the disposal area to the shoal area. The EF model prediction, given in Fig. 9, has residual transport vector directions very similar to the BE predictions, but the vector magnitudes are approximately three times smaller. The same is also true for MPM model predictions, in which the residual vector magnitudes are even smaller than the EF model. Just downriver of the disposal site, all models predict a bed-load divergence zone, which is consistent with patches of eroded bedrock, shown in Fig. 2. Quantitative Eulerian model results, given in Table 3, show scatter, as one would expect from parametric sediment models (White et al., 1973; Nakato, 1990; Vanoni, 1984; Soulsby, 1997; Bilgili, 1993). At the shoal site, bounded by second and third transects, all models show deposition in the main channel. This deposition occurs in the form of ebb-directed bed-load transport at transect 3 and flood-directed transport across transect 2.

Yearly sedimentation rates vary between a maximum of 0.170 m year
1 and a minimum of 0.040 m year
1, predicted by the BE and BA/MPM models, respectively. In the middle of these extreme values, the EF model nicely approximates the historical uniform yearly accumulation of 0.065 m year
1. The rate simulated by the BE model is about half of the 0.35 m year
1 value calculated using field measured velocities in Bilgili et al. (1996) during the pre-dredge and summer studies. This may be due to the fact that Bilgili et al. methodology used only three velocity point measurement stations across the channel while calculating the bed-load flux. In this study, the fluxes are calculated using velocities simulated at each of the grid elements across a transect (on the average, there are 10 elements across a transect). These are then weighted by the characteristic element size, d, and averaged over the width of the cross-section. Net transport results were compared instead of point measurements to reduce the inherent variability associated with these and because insufficient information on weather and local topographic conditions was available to precisely hindcast the Bilgili (1996) study. Between transects 4 and 5, which bound the actual sand waves, all models show increased deposition because of the reduced area of the site, although the bed-load transport into the area is reduced compared to the shoal site. This suggests that although there is likely to have deposition spread around the entire shoal site, the highest accumulation occurs in the exact location of

945

A. Bilgili et al. / Estuarine, Coastal and Shelf Science 58 (2003) 937e950 4

x 10

6.8

5 6.78

13

13

10

Transect 3 Transect 4

13

6.76

Transect 5 Transect 2 6.74

10

17 13

6.72

1 kg/m

5

6.7 3.708

3.71

3.712

3.714

3.716

3.718 5

x 10

Fig. 9. Residual bed-load transport vectors in the Lower Piscataqua River region as calculated by the Engelund and Fredsœ’s model. X and Y axis coordinates are in NH state plane coordinate system in meters. Bed-load transects 3e5 and 2 are also shown together with 5, 10, 13 and 17 m isobaths. Note different vector scales in Figs. 7e9.

the sand waves. The yearly sedimentation rates range between 0.31 m year
1 and 0.07 m year
1, calculated by the BE and MPM models, respectively. It is seen that the BE value nicely approximates the historical record of 0.3 m year
1, while the EF model predicts a yearly growth of 0.13 m, which is about half of the historical record. At the disposal site, bounded by the first and second transects, all models, except BE, predict erosion. This erosion occurs in the form of flood-directed transport from the disposal area into the shoal area through the upriver transect and ebb-directed transport across the downriver transect. Contrary to the other models,

the BE model predicts upriver transport across transect 1 also. This causes deposition as the amount of bed-load leaving the domain across transect 2 is approximately 4.5 times smaller than the portion entering it across the downriver boundary. This is due to the higher order flow velocity dependence (U6) of the BE model amplifying the velocity effects of the Boiling Rock constriction, which creates faster velocities by acting almost as a jet during flood tide and slows down the ebb velocities by acting as a choke during ebb tide. However, it is unlikely that this amount of sediment will be transported upriver through transect 1 due to the unavailability of the sediment in the vicinity of Boiling Rock, characterized by a relatively deeper channel with solid bedrock bottom (Allen et al., 1993). Across transect 6, all models predict ebb-directed transport with increased magnitudes. This is due to ebb current velocities being amplified by the narrow constriction. The noticeably large predicted values shown in Table 3 for this transect are a result of bedload models (especially BE) overpredicting the sediment discharges in high shear stress regimes, a behaviour observed earlier by other scientists (Nakato, 1990; Fredsœ and Deigaard, 1992; Dyer, 1986). It is also unlikely that this transport occurs across the transect because of the unavailability of sediment in the area. Although the transport rates are not realistic, the calculations still show that the overall sediment movement direction across this transect is downriver.

6. Lagrangian bed-load transport pathways Lagrangian particle trajectories are calculated using a time-stepping particle tracking algorithm (Blanton, 1995a,b), which uses a fourth order RungeeKutta integration scheme to compute particle movements in a given velocity field over a time interval. One hundred and ninety-eight particles are released at high tide and

Table 3 Quantitative parametric Eulerian bed-load modeling results. The upper six rows give the net sediment transport in dry mass through specific transects over an M2 cycle Bagnold

BrowneEinstein

Engelund and Fredsœ

Meyer-Peter and Mu¨ller

Transect 1 (+): ebb directed Transect 2 (+): ebb directed Transect 3 (+): ebb directed Transect 4 (+): ebb directed Transect 5 (+): ebb directed Transect 6 (+): ebb directed


1256 kg
1697 kg
427 kg
354 kg
1494 kg +68,006 kg


143,272 kg
31,285 kg +13,740 kg +24,038 kg
17,788 kg +7,310,151 kg

+8280 kg
7903 kg +13,657 kg +13,341 kg
4194 kg +277,985 kg

+2847 kg
3690 kg +8922 kg +8589 kg
1517 kg +273,890 kg

M2 accumulation (kg) (+): deposition

Disposal site:
442 kg Shoal site: +1270 kg Sand waves: +1140 kg

Disposal site: +111,987 kg Shoal site: +45,026 kg Sand waves: +41,826 kg

Disposal site:
16,183 kg Shoal site: +21,560 kg Sand waves: +17,535 kg

Disposal site:
6537 kg Shoal site: +12,613 kg Sand waves: +10,106 kg

Yearly sedimentation rate (m/year) (+): deposition

Disposal site:
0.001 kg Shoal site: +0.040 kg Sand waves: +0.080 kg

Disposal site: +0.270 kg Shoal site: +0.170 kg Sand waves: +0.310 kg

Disposal site:
0.040 kg Shoal site: +0.070 kg Sand waves: +0.130 kg

Disposal site:
0.020 kg Shoal site: +0.040 kg Sand waves: +0.070 kg

946

A. Bilgili et al. / Estuarine, Coastal and Shelf Science 58 (2003) 937e950

low tide throughout the estuary on two different simulations, with increasing density in the Lower Piscataqua River where the shoal is located. Levi’s (Bogardi, 1974) expression for the critical bed-load velocity of graded sediment is used to implement a sediment movement initiation threshold velocity. This expression is given by: sffiffiffiffiffiffiffiffiffi  1=7 pffiffiffiffiffiffiffiffiffi h dg Vc ¼ 1:4 gd50 ln ; ð11Þ 7d50 d50 where h is a characteristic depth, d50 the mean sediment diameter and dg is the maximum sediment diameter. With an h of 15 m, a d50 of 1 mm and a dg of 4.8 mm, Eq. (11) gives a threshold velocity of 66 cm s
1. A 33 cm s
1 threshold velocity is also implemented according to Eq. (11) to investigate the sediment paths of finer graded bed-load material with a d50 of 0.125 mm and a dg of 1.25 mm. The quadratic threshold stress calculated using the threshold velocities given by the equations from Bogardi (1974) are within the limits of threshold stresses given in many references, including that of Soulsby (1997), Miller et al. (1977), Buffington (1999), Paphitis (2001) and Lavelle and Mofjeld (1987) for coarse and fine sands. Considering that the Lagrangian simulations are only concentrating on predicting overall sediment movement patterns and excursions, we did not see the need to fine tune the threshold velocity or to use a threshold stress as the movement initiation criteria. Above the threshold velocities, all particles are assumed to travel at one-sixth of the depth-averaged current. This assumption is consistent with the following range, given by Bogardi (1974): 1 1 V R UB R V; 4 7

ð12Þ

where UB is the mean transport velocity of a bed-load particle and V is the depth-averaged flow velocity. The final location of each sediment particle after one M2 tidal cycle is saved and analyzed. As seen in Fig. 10, Lagrangian sediment tracks obtained using a 66 cm s
1 threshold velocity for high tide released particles show an upstream displacement from the disposal area into the shoal area over an M2 tidal cycle. This pattern can be explained in part by the depth-averaged tidal current being asymmetric over a tidal cycle, resulting in a net upriver transport in the case of coarse sediments, which are more sensitive to instantaneous peak velocities. This upriver transport is found to be independent of the particle release time, although a release at low tide, right before flood, increased the amount of upriver transport from the disposal area into the shoal area by 20%, compared to a high tide release. The displacements are more pronounced along the northeast (Maine) side of the channel and may be due to an abrupt flood current direction

change towards the north, located just southwest of Boiling Rock. This creates a tidal flow hugging the Maine side of the channel with increased velocities. The upstream transport peaks around the shoal and starts decreasing as one moves towards Frankfort Island. The majority of particles released upriver of the shoal and downriver of the Boiling Rock constriction, however, end their trajectories downriver of their initial positions after an M2 cycle. The tracks indicate that the current, and consequently the particles, follow the channel axis, staying almost on the same track during ebb and flood. Particle excursions are considerably larger in the center channel than along the sides because of higher velocities, with many off-channel particles not moving at all because the movement initiation threshold velocity is never reached. An animation of flood-released particles moving in the flow field, which can be accessed at http:// www-nml.dartmouth.edu/NML-02-1/dranim66.flc (Bilgili, 2002) shows particles accumulate at the shoal site both during flood and ebb. It is interesting to see that the disposal site is also a depositional site during ebb for particles initially positioned between the shoal site and the disposal site, acting as an immediate sediment source to the shoal site during flood. The sand wave bias to the southwest of the channel, shown in Fig. 2, is also replicated by the model at the end of an M2 cycle, in both low and high tide release simulations. The high tide release results obtained using a 33 cm s
1 threshold velocity show a net downriver transport for all particles (Fig. 11), suggesting this section of the estuary to be an exporter of finer sediments. Most particles end their trajectories downriver of their starting points after an M2 cycle with no spatial concentration anywhere throughout the domain, a result easily observable from the corresponding drogue animation, which can be seen at http://www-nml.dartmouth.edu/NML02-1/dranim33.flc (Bilgili, 2002). This trend is also found to be independent of particle release time.

7. Conclusion Hydrodynamic model result time series comparison with observations and error analyses indicates that the circulation model predicts the main features of the M2M4M6 tidesdsurface elevation and cross-sectionally averaged currentsdin the Lower Piscataqua River, where the kinematic assumption holds true. Using simulated currents for sediment transport calculations is an even more severe test of a numerical model than just a basic comparison of hydrodynamical variables since spatial resolution, residual of currents raised to various powers, and residuals of Lagrangian currents subject to thresholds must also be accurate. The general agreement of predicted flux with the formation of known geological features indicates that the

A. Bilgili et al. / Estuarine, Coastal and Shelf Science 58 (2003) 937e950

947

Fig. 10. Displacement vectors of sediment particles released at high tide and subjected to a threshold velocity of 66 cm s
1. The dots mark the initial positions and the tips of the lines mark the final positions of the particles after an M2 tidal cycle. Isobaths (5, 10, 13 and 18 m) and the three sand waves are also plotted for reference. Coordinates are in NAD83 NH state plane coordinates in meters.

hydrodynamic model does simulate these flow details well at the current grid resolution level. Although there is always going to be a better hydrodynamic model that covers a wider range of estuarine dynamics, we believe that the conclusions that we derive here are well-founded. This is due to the wellmixed nature of the domain and the scaling (both empirical and theoretical) of the dynamic terms for this particular section of the Great Bay Estuary, where the non-linear effects responsible for shoaling manifest themselves in the horizontal plane, rather than the vertical, and are strongly related to bottom topography and channel geometry. The model’s neglect of acceleration plays a role near slack water when current levels are down, but the effects of this on bed-load transport, which is more sensitive to peak currents, should be minimal. The vertical variability of residual currents is also assumed to be minimal and a vertically averaged velocity is used to drive the Eulerian and Lagrangian bed-load models. This is consistent with the parametric

methods used to derive these models, and consequently, with the input parameters that they require. The main advantage that we get from the model by making these assumptions is very good wettingedrying capability and, consequently for this system, very good predictions of horizontal distribution. This is due to the fact that current speeds in the Piscataqua River are directly affected by filling and emptying of Great Bay proper, which is a flooding and drying basin. All bed-load models tested in this study predict accumulation in the shoal area and sediment movement from the disposal site into the dredge site. The use of several Eulerian bed-load flux vector models, having different quantitative outcomes, reflects the state of the art in which there is no overall agreement as to which one is the best. Largest instantaneous bed-load fluxes are simulated by the BE model. This result is consistent with the calculations of Nakato (1990), who found that the bed-load predictions by the BE model can be more than 10 times greater than those predicted by the other models,

948

A. Bilgili et al. / Estuarine, Coastal and Shelf Science 58 (2003) 937e950

Fig. 11. Displacement vectors of sediment particles released at high tide and subjected to a threshold velocity of 33 cm s
1. The dots mark the initial positions and the tips of the lines mark the final positions of the particles after an M2 tidal cycle. Isobaths (5, 10, 13 and 18 m) and the three sand waves are also plotted for reference. Coordinates are in NAD83 NH state plane coordinates in meters.

including the MPM and EF models. Based on the results presented here, the BE model predictions are somewhat more consistent with the yearly growth of the sand waves while the EF predictions agree more with the uniform yearly accumulation over the entire dredge site. The BA and MPM models replicate the observed pathways but underestimate long-term accumulations and depositions. Although shoal formation can be explained by what we do here, future work should also investigate effects of episodic events and springeneap cycle. The Lagrangian study also yields net movement associated with the shoal buildup. This trend is found to be independent of the sediment particle release time in the tidal cycle. However, the apparent increase in the upriver transport after a low tide release suggests that a high tide or initial ebb stage disposal should be planned to keep the immediate movement of the disposed sediment into the shoal area at a minimum level when a downriver riverine disposal site is used. In addition, the Lagrangian model reveals that simulation results are

sensitive to particle size through the threshold velocity. This indicates a significant sediment sorting characteristic, which deserves further investigation throughout the estuary. Future work should also attempt to identify the original source of coarse sediments forming the sand waves. Although the downriver disposal area may act as an immediate source to the shoal site, the fact that the sand waves form in the absence of a riverine disposal suggests an original upriver source. This idea is backed up by the fact that coarse sediment availability is limited at downriver locations beyond Boiling Rock and the net bed-load transport is towards the estuarine mouth between Boiling Rock and Portsmouth Harbor (Ertu¨rk et al., 2002). Although this study was exploratory in its design, it shows that Lagrangian methods have great potential for sophisticated simulations of sedimentary processes (Davies et al., 2000). This is especially true considering the rapid growth in the computing power in recent years, which will allow greater statistical sophistication.

949

A. Bilgili et al. / Estuarine, Coastal and Shelf Science 58 (2003) 937e950

Acknowledgements We thank Prof. Wendell Brown and Prof. Barbaros Cxelikkol for many helpful suggestions, and Mr. Brian Blanton for his help in the particle tracking program setup. Field observations were obtained as part of the 1975 Great Bay Estuarine Field Program under the New Hampshire Sea Grant College Program and NOS. Financial support from the Cooperative Institute for Coastal and Estuarine Environmental Technology (CICEET), NOAA/UNH (grant numbers NA77OR0357, NA87OR0512, NA87OR338) is gratefully acknowledged. Appendix A. Eulerian parametric bed-load models The models described below were used to calculate the bed-load flux at each time step (every 27.945 s). Time varying bottom stress and vertically averaged velocities are provided by the two-dimensional hydrodynamic model. A.1. Bagnold’s model Bagnold’s model is given by:   gs
g  eb ; q#B tanðcÞ ¼ tB U g

ð13Þ

A.2. BrowneEinstein model

ð14Þ

where qB is the bed-load flux expressed as dry weight per time per unit distance across the channel, rs the density of sediment, K a time factor, and F is a parameter called the intensity of bed-load transport. The time factor K is given by: sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 36n2 36n2 þ 3 K¼
ð15Þ 3 3 gd 50 ðgs =g
1Þ gd 50 ðgs =g
1Þ; where n is the kinematic viscosity of water, taken as 1:9 ! 10
6 m2 s
1 . The intensity of bed-load transport, F, is a function of another parameter called the intensity of flow, j, and is given by:

ð16Þ

after Rouse (1950), where j is defined as: J¼

ðgs
gÞd50 : Rs tB

ð17Þ

In Eq. (17), Rs is the total bottom stress to skin stress ratio. A.3. Engelund and Fredsœ’s model Engelund and Fredsœ’s model is based on the assumption that the bed-load is the transport of certain fraction, p (probability), of the particles in a single layer of sediments: rffiffiffiffiffi  p tC qB ¼ 9:3 pd50 U 1
0:7 : ð18Þ 6 t# pffiffiffiffiffiffiffiffiffiffiffiffi Where, U* is the frictional velocity (¼ jtS j=r, with tS being the skin stress) and tC and t# are the Shields parameters (or critical non-dimensional shear stresses) given by: tC ¼

where gs and g are the specific weights of sediment and water, respectively, q#B the bed-load flux expressed as submerged weight per unit width per unit time, c the dy the vertically namic friction angle, tB the bottom stress, U averaged velocity, and eb is an efficiency factor. Bagnold (1966, 1973) presents two curves for the estimation of eb and c under various flow and sediment conditions. A sediment movement initiation threshold of 66 cm s
1 is  in this study. implemented on the magnitude of U

BrowneEinstein model is given by: sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi   ! gs qB ¼ rs K
1 d 350 F; g g

 3 1 ; F ¼ 40 J

U 2C ðrs =r
1Þgd50

and

t# ¼

U 2 ; ðrs =r
1Þgd50

ð19Þ

for critical and skin stresses, respectively, where U*C is the critical frictional velocity. The probability, p in Eq. (18) is given by: "  4 #
1=4 ðp=6ÞtanðcÞ p¼ 1þ : ð20Þ t#
tC

A.4. Meyer-Peter and Mu¨ller model Similar to the Engelund and Fredsœ’s, Meyer-Peter and Mu¨ller’s model relates bed-load transport to an effective shear stress and defines the beginning of motion, so that the transport becomes zero when t# is equal to tC: 
1=2  
1=2  r 1 3=2 qB ¼ 8ðt#
tC Þ gs : ð21Þ rs
r gd 350

References Allen, D.K., Celikkol, B., Swift, M.R., 1993. Bathymetric and sediment studies of the Piscataqua River using a Chirp Sonar. Technical report, New Hampshire Coastal Program, NH Office of State Planning, New Hampshire, USA, 160 pp. Armstrong, P., 1974. Copper, zinc, chromium, lead and cadmium in the unconsolidated sediments of Great Bay Estuary, New Hampshire. M.Sc. thesis, University of New Hampshire, Durham, New Hampshire, 85 pp. Bagnold, R.A., 1966. An approach to the sediment transport problem from general physics. U.S. Geological Survey Professional Paper, no. 422-I.

950

A. Bilgili et al. / Estuarine, Coastal and Shelf Science 58 (2003) 937e950

Bagnold, R.A., 1973. The nature of saltation and ‘bed-load’ transport in water. Proceedings of the Royal Society of London Series A 332, 473e504. Bilgili, A., 1993. A study of the bed-load sediment transport in the Piscataqua River navigation channel, New Hampshire. M.Sc. thesis, University of New Hampshire, Durham, NH, 131 pp. Bilgili, A., 2000. Finite element modeling of tidal flow and bed-load transport in the Great Bay Estuarine System, New Hampshire. Ph.D. dissertation. University of New Hampshire, Durham, NH, 255 pp. Bilgili, A., 2002. Animations of Lagrangian sediment transport in the Great Bay Estuary. Numerical methods laboratory internal report, no: NML-02-01. Dartmouth College, Hanover, New Hampshire. http://www-nml.dartmouth.edu/Publications/internal-reports/ NML-02-01/ (last visited April 18th, 2002). Bilgili, A., Swift, M.R., Cxelikkol, B., 1996. Shoal Formation in the Piscataqua River, New Hampshire. Estuaries 19(3), 518e525. Blanton, B.O., 1995a. Modifications to Quoddy2 time-stepping code to include Drog3ddt. Memorandum to GLOBEC team. Ocean Processes Numerical Methods Laboratory, University of North Carolina at Chapel Hill, NC, 4 pp. Blanton, B.O., 1995b. Drog3d: user’s manual for 3-dimensional drogue tracking on a finite element grid with linear finite elements. Ocean Processes Numerical Methods Laboratory, University of North Carolina at Chapel Hill, NC, 15 pp. Bogardi, J., 1974. Sediment Transport in Alluvial Streams. Akade´miai Kiado´, Budapest, 826 pp. Brown, C.B., 1950. Sediment transportation (Chapter 12). In: Rouse, H. (Ed.), Engineering Hydraulics. John Wiley and Sons, New York, pp. 796e799. Bryce, S., Larcombe, P., Ridd, P.V., 1998. The relative importance of landward-directed tidal sediment transport versus freshwater flood events in the Normanby River estuary, Cape York Peninsula, Australia. Marine Geology 149, 55e78. Buffington, J.M., 1999. The legend of shields. Journal of Hydraulic Engineering 125 (4), 376e387. Clere, J.C., 1993. Sediment transport modeling of the Piscataqua River navigation channel. M.S. thesis, University of New Hampshire, Durham, NH, 317 pp. Cxelikkol, B., Ward, L.G., Moynihan, R., Shevenell, T.C., Mahmutog˘lu, S., Scott, J.P., Adams, J., Bilgili, A., Aydınog˘lu, Z., 2002. Hydrodynamic investigation of Hampton/Seabrook Harbor, New Hampshire, report submitted to the Pease Development Authority, Division of Ports and Harbors, State of New Hampshire by the University of New Hampshire, Jere A. Chase Ocean Engineering Lab., Durham, NH. Davies, M., Serrer, M., Watson, D.A.W., 2000. Psed 2000: a Lagrangian sediment transport model. Technical report HYD-TR-051. Canadian Hydraulics Centre, Ottawa, Canada, 37 pp. Dyer, K.R., 1986. Coastal and Estuarine Sediment Dynamics. John Wiley and Sons, New York, 342 pp. Einstein, H.A., 1942. Formulas for the transportation of bed load. Transactions of the American Society of Civil Engineers, New York 107, 561e597. Engelund, F., Fredsœ, J., 1976. A sediment transport model for straight alluvial channels. Nordic Hydrology 7, 293e306. Ertu¨rk, S x.N., Bilgili, A., Swift, M.R., Brown, W.S., Cxelikkol, B., Ip, J.T.C., Lynch, D.R., 2002. Simulation of the Great Bay Estuarine System: tides with tidal flats wetting and drying. Journal of Geophysical ResearchdOceans 107(C5), 29 pp. (10.1029/ 2001JC000883). Fredsœ, J., Deigaard, R., 1992. Mechanics of Coastal Sediment Transport. World Scientific Publishing Co., Singapore, 369 pp. Friedrichs, C.T., Lynch, D.R., Aubrey, D.G., 1992. Velocity asymmetries in frictionally-dominated tidal embayments: longitudinal and lateral variability. In: Prandle, D. (Ed.), Dynamics and Exchanges in Estuaries and the Coastal Zone. Coastal and

Estuarine Studies 40. American Geophysical Union, Washington, DC, pp. 277e312. Harris, P.T., Collins, M., 1988. Estimation of annual bed-load flux in a macrotidal estuary: Bristol Channel, United Kingdom. Marine Geology 83, 237e252. Ip, J.T.C., Lynch, D.R., Friedrichs, C.T., 1998. Simulation of estuarine flooding and dewatering with application to Great Bay, New Hampshire. Estuarine, Coastal and Shelf Science 47, 119e141. Larcombe, P., Jago, C.F., 1996. The morphological dynamics if intertidal megaripples in the Mawddach Estuary, North Wales and the implications for palaeflow reconstruction. Sedimentology 43, 541e559. Larcombe, P., Ridd, P.V., 1996. Dry season hydrodynamics and sediment transport in mangrove creeks. In: Pattariatchi, C. (Ed.), Mixing Processes in Estuaries and Coastal Seas. Coastal and Estuarine Studies 46. American Geophysical Union, Washington, DC, pp. 388e404. Lavelle, J.W., Mofjeld, H.O., 1987. Bibliography on sediment threshold velocity. Journal of Hydraulic Engineering 113 (3), 389e393. Ludwick, J.C., 1975. Tidal currents, sediment transport and sand banks in Chesapeake Bay entrance, Virginia. In: Cronin, L.E. (Ed.), Estuarine Research, vol. 2. Academic Press, New York, pp. 365e380. Maine Geological Survey Department, 1988. Inter-departmental memorandum to Maine Environmental Protection Department about the maintenance dredging of the Piscataqua River. Marconi, G., 2001. Personal communication with Portsmouth, N.H., harbor master Geno Marconi about the amount of dredged sediment during the November 2000 Piscataqua River dredging. McLaughlin, J.W., 2001. Two dimensional finite element model for estuarine hydrodynamics with application to the Great Bay, NH. M.Sc. thesis, Thayer School of Engineering, Dartmouth College, Hanover, NH, USA, 158 pp. Meyer-Peter, E., Mu¨ller, R., 1948. Formulas for bed-load transport. Proceedings of Second Meeting of the International Association for Hydraulic Research, pp. 39e64. Miller, M.C., McCave, I.N., Komar, P.D., 1977. Threshold of sediment motion under unidirectional currents. Sedimentology 24, 507e527. Nakato, T., 1990. Tests of selected sediment-transport formulas. Journal of Hydraulic Engineering 116 (3), 362e379. New Hampshire Fish and Game Department, Marine Fisheries Division, 1983. Letter to the New Hampshire Office of State Planning, Concord, NH. Paphitis, D., 2001. Sediment movement under unidirectional flows: an assessment of empirical threshold curves. Coastal Engineering 43, 227e245. Rouse, H. (Ed.), 1950. Engineering Hydraulics. John Wiley and Sons, New York, pp. 769e857. Signell, R.P., List, J.H., Farris, A.S., 2000. Bottom currents and sediment transport in Long Island Sound: a modeling study. Journal of Coastal Research 16 (3), 551e566. Soulsby, R., 1997. Dynamics of Marine Sands. Thomas Telford Publications, London, United Kingdom, 249 pp. Swenson, E., Brown, W.S., Trask, R.P., 1977. Great Bay estuarine field program 1975 data report. Part 1: Currents and sea levels. University of New Hampshire Sea Grant technical report UNHSG-157, Durham, NH, 109 pp. Swift, M.R., Brown, W.S., 1983. Distribution of bottom stress and tidal energy dissipation in a well-mixed estuary. Estuarine, Coastal and Shelf Science 17, 297e317. Vanoni, V.A., 1984. Fifty years of sedimentation. Journal of Hydraulic Engineering, ASCE 110 (8), 1022e1057. White, W.R., Milli, H., Crabbe, A.D., 1973. Performance of theoretical method when applied to flume and field data. sediment transport: an appraisal of available methods. vol. 2, Report INT 119, Hydraulic Research Station, Wallingford, England.