Modeling tidal circulation and stratification in Skagit River estuary using an unstructured grid ocean model

Ocean Modelling 28 (2009) 34–49

Contents lists available at ScienceDirect

Ocean Modelling journal homepage: www.elsevier.com/locate/ocemod

Modeling tidal circulation and stratification in Skagit River estuary using an unstructured grid ocean model Zhaoqing Yang *, Tarang Khangaonkar Battelle Marine Sciences Laboratory, 1100 Dexter Avenue North, Suite 400, Seattle, WA 98109, USA

a r t i c l e

i n f o

Article history: Received 19 April 2008 Received in revised form 12 July 2008 Accepted 21 July 2008 Available online 29 July 2008 Keywords: Unstructured Mesh Estuary Tidal flats Tidal Circulation Salinity Intrusion

a b s t r a c t Tidal circulation and river plume dynamics in shallow-water estuarine systems with large intertidal zones are complex. Strong asymmetries in tidal currents and stratification often occur in the intertidal zones and subtidal channels over a tidal cycle. The Skagit River is the largest estuary with respect to the discharge of a significant amount of freshwater and sediment into Puget Sound, Washington. It consists of a large intertidal zone with multiple tidal channels near the mouth of the estuary. To simulate the tidal circulation and salinity stratification accurately in the intertidal region, an unstructured grid numerical model with wetting–drying capability and the capability to accurately represent the bathymetry of tidal flats and the geometry of shallow distributary channels is necessary. In this paper, a modeling study for the Skagit River estuary using a three-dimensional unstructured grid, finite-volume coastal ocean model (FVCOM) supported by high-resolution LIDAR data is presented. The hydrodynamic model was validated with observed water surface elevation, velocity, and salinity data over spring and neap tidal cycles under low-river-flow and high-river-flow conditions. Wetting and drying processes in the intertidal zone and strong stratification in the estuary were simulated successfully by the model. Model results indicate that the Skagit River estuary is a highly stratified estuary, but destratification can occur during flood tide. Tides and baroclinic motion are the dominant forcing in the Skagit River estuary, but strong wind events can affect the currents in the intertidal zone significantly. Preliminary analysis also indicated that the salinity intrusion length scale is proportional to the river flow to the 1=4 power. Ó 2008 Elsevier Ltd. All rights reserved.

1. Introduction Intertidal zones of estuaries with large tidal flats play an important role on tidal circulation, salinity intrusion, and sediment transport. Intertidal zones often include different marine habitats that support a wide variety of marine species. The long-term stability and sustainability of these habitats is threatened by growing anthropogenic pressures and from potential changes brought about by climate change and sea level rise. Tidal circulation and salinity intrusion processes in intertidal regions generally are very complex and highly dynamic due to the effect of extreme shallow water depths in the tidal flats and tidal channels. This complexity is illustrated, for instance, by asymmetry of vertical mixing between ebb and flood currents, and dispersion of salinity gradients due to lateral advection over the intertidal zone (Geyer et al., 2000; Ralston and Stacey, 2005). Complexities can also arise because during low tides, some of the channels may be blocked, resulting in ponding of water, occurrence of transient shallow wave induced currents, and wind effects (Greenberg et al., 2005; Gorman and Neilson, 1999; Collins et al., 1998). However, because * Corresponding author. Tel.: +1 206 528 3056; fax: +1 206 528 3556. E-mail address: [email protected] (Z. Yang). 1463-5003/$ - see front matter Ó 2008 Elsevier Ltd. All rights reserved. doi:10.1016/j.ocemod.2008.07.004

of difficulties in conducting undisturbed field measurements, the lack of high-resolution bathymetric and topographic data, and computational limitations with simulation of wetting and drying process, the intertidal zone has not received much attention, and circulation processes are not understood as well. In recent years, some advanced numerical models have been developed for coastal oceans and estuaries with improved performance, including the capability of simulating wetting and drying process (Zhang et al., 2004; Ji et al., 2001; Chen et al., 2003; Oey, 2006; Casulli and Zanolli, 1998, 2002; Cheng and Casulli, 2004). Numerical models with an unstructured grid framework (Chen et al., 2003; Zhang et al., 2004; Cheng and Casulli, 2004) and the availability of high-resolution light detection and ranging (LIDAR) data in the tidal flats make it possible to simulate hydrodynamics in estuaries with large tidal flats with satisfactory accuracy. This capability is of considerable value to a large community of biologists, geomorphologists, ecologists, and restoration planners engaged in the activity of restoring tidal estuarine functions and nearshore habitat in support of fish stock recovery as in the case of in Puget Sound, Washington, and other estuarine and coastal regions around United States. In this paper, we present a modeling study of tidal circulation and salt intrusion in a highly stratified

Z. Yang, T. Khangaonkar / Ocean Modelling 28 (2009) 34–49

estuary with large tidal flats in Puget Sound, Washington, using an unstructured grid finite-volume coastal ocean model (FVCOM) and LIDAR bathymetric data. 2. The Skagit River estuary The Skagit River is the largest river in the Puget Sound estuarine system (Fig. 1). The drainage area of the Skagit River watershed is about 8000 km2, and the Skagit River alone is responsible for about 34–50% of the total freshwater flow into Puget Sound, depending on the season (Hood, 2006; Babson et al., 2006; Cannon, 1983). The river flow peaks both in winter (due to runoff), and again in late spring or early summer (due to snow melt), and is often at a minimum in September. The mean flow of the Skagit River at Mt. Vernon is 468 m3/s with records of maximum and minimum flows at 5100 and 78 m3/s, respectively (Wiggins et al., 1997). The Skagit River splits into the North Fork and the South Fork before it enters the Skagit Bay. The North Fork channel runs westerly as a dike-bounded channel between the marshlands and carries about 60% of the total Skagit River flow. The South Fork, instead, enters Skagit Bay through many small tidal channels. The area bounded by the North Fork and the South Fork is known as Fir Island. It was extensively diked in the late 19th and early 20th centuries to provide flood protection for the agricultural lands. The construction of levees has interrupted periodic sediment distribution from the river to Fir Island and the surrounding marshes (Hood, 2004). Skagit Bay is very shallow in most of the region. Water depths in most of Skagit Bay are less than 5 m. The deepest area is about 40 m near Deception Pass in the north of Skagit Bay. A large tide flat exists at

35

the mouth of the estuary, and most of the northeastern region of Skagit Bay is above the mean lower low water (MLLW) line. Unlike many estuaries and bays that often connect to open coastal waters, the Skagit River estuary is confined by the coast of Whidbey Island directly in front of the mouth of the Skagit River flow west of the intertidal zone. A relatively deep and narrow channel (25–30 m) exists between the east coast of Whidbey Island and the tidal flats of Skagit Bay. Skagit Bay connects to the Strait of Juan de Fuca through Deception Pass to the east and Padilla Bay through the Swinomish Channel to the north, which are two of the three connections of Puget Sound to the coastal waters enclosed by Whidbey Island, also known as Whidbey Basin. Deception Pass is a narrow and relatively shallow passage with strong tidal currents. Dredging of the Swinomish Channel began in the 1890s and continues today. The Swinomish Channel needs to be dredged every 3–4 years to prevent groundings. It is dredged by the U.S. Army Corps of Engineers to a depth of 3.6 m. The channel is used by boat builders, log-tow companies, fishermen, dryboat operators, and boaters on their way from Puget Sound to the San Juan Islands. A jetty was built in 1938 along the northern river bank of the North Fork to prevent sediment transport into the Swinomish Channel. Skagit Bay is subjected to tides from Puget Sound, primarily propagating from the south of the bay. Tides in Skagit Bay are mixed, semi-diurnal dominant tides and show large inequalities in tidal range and strong spring-neap tidal cycle. Density-induced baroclinic currents are important in Skagit Bay because of high salinity stratification and a strong freshwater front produced by Skagit River flow. The wind effect could be important in the tidal flats during strong wind events.

Fig. 1. Skagit River estuary in Puget Sound, Washington.

36

Z. Yang, T. Khangaonkar / Ocean Modelling 28 (2009) 34–49

3. Methodology 3.1. Data To simulate the tidal circulation and salinity stratification accurately in the Skagit River estuary and bay, high-resolution bathymetric data are required, particularly in the tidal flats area. The Skagit Bay hydrodynamic model was initially developed using the bathymetric data obtained from the University of Washington Puget Sound Digital Elevation Model (DEM) (Finlayson, 2005). The DEM bathymetric data had a horizontal grid resolution of 9 m by 9 m, which was sufficient to define bathymetry in most of the region in the Skagit River estuary and bay. However, the vertical positional errors of the DEM data were as high as 1.5 m at certain locations in the nearshore region where bathymetric and topographic data were combined from different data sources. The DEM data alone was not sufficient to describe in detail the tidal channels and tidal flats in the nearshore region. Elevations and geometry details of the intertidal zone with subtidal channels have been shown to play an important role in transport and exchange processes in estuaries (Simenstad, 1983; Hood, 2007; French and Spencer, 1993; Ralston and Stacey, 2005). After the initial model calibration against field data, it became apparent that modeling the hydrodynamics in the Skagit River estuary could be strongly affected by the accuracy of the bathymetry data over tidal flats. Therefore, the model bathymetry was updated with high-resolution LIDAR bathymetric data collected by Skagit River System Cooperative (SRSC) and U.S. Geological Survey (USGS) in the Skagit Bay tidal flats. The LIDAR data have a horizontal resolution of 1.8 m by 1.8 m and a vertical resolution of 0.15 m. Details of tidal channels and sea-bottom elevations in the intertidal zone were represented well in the LIDAR data collected by SRSC (Fig. 2). Using information provided from LIDAR data, the model grid was constructed such that unstructured grid cells were specified along tidal channels in the nearshore region. The accuracy of model bathymetry in the tidal flats was improved significantly and features of multiple tidal channels were better represented with LIDAR data. The river cross-section data in the Skagit River and the bathymetric data in the Swinomish Channel were obtained from the U.S. Amy Corps of Engineers, Seattle District. Skagit River inflow data were obtained from a USGS stream gage about 25 km upstream of the estuarine mouth, at Mt. Vernon, WA. The river flow is regulated by the upstream hydropower dams that result in diurnal peaking cycles due to power generation. Wind data were obtained from National Oceanic and Atmospheric Administration (NOAA), National Weather Service (NWS), Paine Field meteorological station. Tidal elevations at the model open boundaries at the mouth of Skagit Bay, Deception Pass, and the entrance of Swinomish Channel were obtained from the XTIDE prediction based on NOAA’s National Oceanic Service algorithms. To support calibration of the model, synoptic field-data collection was conducted from June 7–22, 2005. Two continuously monitoring stations equipped with Inter Ocean S4 current meters and Hydrolab conductivity-temperature depths (CTDs) near the bottom of the water column were deployed at the mouth of the North Fork branch of the Skagit River (NF) and near the deep channel of Skagit Bay (SB). Tidal elevations, velocities, and salinities were measured at these two stations. The average river flow for the model calibration period was about 325 m3/s, which corresponds to the low-river-flow condition. Oceanographic data were collected at a number of locations in Skagit Bay and the Swinomish Channel by USGS from May 1–31, 2006 (Grossman et al., 2007). A high-flow event occurred from May 17–22, 2006 with a daily peak flow of 1133 m3/s on May 19, 2006. One monitoring station equipped with bottom-mounted

Accoustic-Doppler Velocimeter (ADV) and CTD was deployed in Skagit Bay north of the Skagit River (S6) from May 16–31, 2006, which covered this high-river-flow event. The average flow for this field-data collection period was 714 m3/s. Tidal elevation, velocity, and salinity data from the S6 station were used for the model verification for the high-river-flow condition. 3.2. Hydrodynamic model FVCOM The three-dimensional (3D) FVCOM developed by Chen et al. (2003) was used in this study to simulate the tidal circulation and salt intrusion in the Skagit River Estuary. FVCOM simulates water surface elevation, velocity, temperature, salinity, sediment, and water-quality constituents in an integral form by computing fluxes between non-overlapping horizontal triangular control volumes. This finite-volume approach combines the advantages of finite-element methods for fitting complex boundaries and finite-difference methods for simple discrete structures and computation efficiency. A sigma-stretched coordinate is used in the vertical plane to better represent the irregular bottom topography. Unstructured triangular cells are used in the horizontal plane. The detailed theoretical aspects of the model can be found in Chen et al. (2006). The model has been used to simulate circulation and transport processes in many estuaries, coastal waters, and open oceans (Zheng et al., 2003; Chen and Rawson, 2005; Zhao et al., 2006; Weisberg and Zheng, 2006a,b; Isobe and Beardsley, 2006; Aoki and Isobe, 2007; Chen et al., 2008; Yang and Khangaonkar, 2008 ; Yang et al., submitted for publication). The three 3D governing equations for Reynolds-averaging turbulent flows with the Boussinesq approximation are presented in the flowing form in FVCOM:

ou ov ow þ þ ¼0 ox oy oz   ou ou ou ou 1 oP o ou þ u þ v þ w  fv ¼  þ Km þ Fu ot ox oy oz oz qo ox oz   ov ov ov ov 1 oP o ou þ Fv þ u þ v þ w þ fu ¼  þ K ot ox oy oz qo oy oz m oz oP ¼ qg oz

ð1Þ ð2Þ ð3Þ ð4Þ

where (x, y, z) are the east, north, and vertical axes in the Cartesian coordinates; (u, v, w) are the three velocity components in the x, y and z directions; (Fu, Fv) are the horizontal momentum diffusivity terms in the x and y directions; Km is the vertical eddy viscosity coefficient; q is density; P is pressure; and f is the Coriolis parameter. The 3D transport equations for temperature and salinity are solved simultaneously in FVCOM:

  oT oT oT oT o oT þ FT þu þv þw ¼ Kh ot ox oy oz oz oz   oS oS oS oS o oS þu þv þw ¼ Kh þ FS ot ox oy oz oz oz

ð5Þ ð6Þ

where T and S are temperature and salinity, Kh is the vertical eddy diffusivity coefficient, and (FT, FS) are the horizontal thermal and salt diffusivity terms. Temperature and salinity are related to density through the equation of state

q ¼ qðT; SÞ

ð7Þ

The bottom friction is described by the quadratic law with the drag coefficient determined by the logarithmic bottom layer as a function of bottom roughness. Wind stress is specified at the water surface. The model employs the Smagorinsky scheme for

Z. Yang, T. Khangaonkar / Ocean Modelling 28 (2009) 34–49

37

Fig. 2. Skagit Bay DEM and LIDAR bathymetric data.

horizontal mixing (Smagorinsky, 1963) and the Mellor Yamada level 2.5 turbulent closure scheme for vertical mixing (Mellor and Yamada, 1982). The effects of wind-driven surface waves breaking on the turbulent energy (Mellor and Blumberg, 2004) are also incorporated in FVCOM. The unstructured grid for Skagit Bay is shown in Fig. 3. The upstream end of the model domain was extended considerably beyond the split of the Skagit River into the North Fork and the South Fork. The river inflow data collected at the upstream USGS gage were used to specify the boundary condition of the model river inflow. To simulate the tidal wave propagation and salinity intrusion properly in the multi-channel and tide flat area, finer grid cells were specified in the North and South Forks and the river delta region. The model element size varied from around 15 m in the estuary to 400 m at the entrance of Skagit Bay. The model grid resolution was gradually reduced away from the estuarine delta to the open boundaries. The model consists of 9122 elements and 5496 nodes in the horizontal plane. To simulate salinity stratification accurately in the estuary, 25 uniform vertical layers were specified in the water column in a sigma-stretched coordinate system. The model was set up in Universal Transverse Mercator (UTM) North American Datum (NAD) 83 coordinates in the horizontal plane with reference to MSL in the vertical direction. The model consists of three open boundaries: the mouth of Skagit Bay, Deception Pass at the north of Skagit Bay, and the Swinomish Channel at the entrance of Padilla Bay. Model open

boundary conditions are specified by tidal elevations simulated using the XTIDE program based on National Oceanic Service algorithms. There were no salinity data available along the open boundaries. Salinity profiles along the open boundaries were initially estimated based on historical data in the Puget Sound area and further adjusted during the model calibration. Field data indicated that temperature variations in the study area are within 3 °C, and salinity variations are more than 20 psu. Therefore, the effect of temperature on density is not as significant as salinity. Therefore, temperature transport was not simulated in this study. Wind stress was specified at the water surface. Wind stress was applied uniformly to the entire model domain. 4. Results and discussion 4.1. Model calibration and verification Model calibration and verification were conducted from June 7–23, 2005, and May 16–31, 2006, corresponding to low and high river flows, respectively. Wetting and drying processes in the tidal flats of Skagit Bay were simulated in the model. A water depth of 5 cm was used as a dry cell criterion in the model. Model results were compared to observed data, and error statistics were calculated to measure the level of accuracy for model calibration and verification.

38

Z. Yang, T. Khangaonkar / Ocean Modelling 28 (2009) 34–49

Fig. 3. FVCOM model grid of Skagit Bay.

The model calibration and verification was conducted mainly through refining the model grid to accurately represent the bathymetry and geometry as well as adjusting the bottom roughness and background vertical eddy viscosity. The horizontal diffusion coefficient multiplier was specified as 0.2 for the Smagorinsky parameterization. The background vertical eddy viscosity was 105 (m2/s). The bottom roughness was set to 0.002 m, and the minimum bottom drag coefficient was specified as 0.0025. The model was run on a 36-processor Beowulf cluster in parallel mode with 5-second time steps. Comparisons of simulated and observed water surface elevation and velocity time-series for the calibration period at station NF in the North Fork Skagit River are shown in Fig. 4. Velocities were measured at about 1.5 m above the bottom. Velocity sensor malfunctioned after June 16, 2005, and therefore no data were available. Model results were generally in good agreement with observed data. Velocity results were projected along the riverchannel direction with positive in the upstream (flooding) and negative in the downstream (ebbing) directions. Salinity comparison is not shown in Fig. 4 because both simulated and observed salinities were zero for the entire period at station NF. Zero salinity at station NF indicated that the salinity intrusion point was located downstream of station NF. The water surface elevation at station NF was raised because of the effects of river backwater and high bed elevations in the intertidal zone. Water surface elevations at NF were about 1.0 m above mean sea level (Fig. 4a). The tidal ranges, therefore, were significantly reduced near the mouth of the North Fork Skagit River. The small phase differences between model prediction and observed data may be due to the lack of accurate river cross-section data. Both model results and field data indicated a consistent unidirectional downstream flow (negative) in the river

channel because of strong river flow. Sharp spikes of velocity corresponded to flood tides during each tidal cycle (Fig. 4b). Comparisons of simulated elevation, velocity, and salinity timeseries to field measurements at station SB near the deep channel of Skagit Bay are provided in Fig. 5. Velocities and salinities were measured at about 1.5 and 1.0 m above the bottom, respectively. Model results matched the observed data reasonably well. There were some uncertainties about the accuracy of the observed data during certain periods. For example, the tidal range after June 19, 2005, was expected to increase continuously, but a sudden drop was observed. Observed velocities during the period June 15, 2005, through June 19, 2005, were consistently close to zero, which were also questionable. Unlike station NF, strong tidal variations in water surface elevation, velocity, and salinity were observed at station SB. The spring-neap tidal cycle and the diurnal inequality were well reproduced in the model at station SB. The salinity in Skagit Bay increased up to 30 ppt during flood tide because of the intrusion of Puget Sound saline water and dropped below 10 ppt during ebb because the North Fork Skagit River freshwater plume was dispersed (Fig. 5c). Salinity profiles were collected near station SB during the deployment of instruments in the field-data collection program. Fig. 6 presents a comparison of simulated and observed salinity profiles on June 7, 2005, 12:38 PM. Strong salinity stratification in the bay shown in the measured data was also reproduced in the model. The pycnocline was at a very shallow depth of around 2 m. To illustrate the model capability of simulating the wetting and drying processes in the intertidal zone, simulated surface velocity distribution in Skagit Bay along with wetting and drying areas during low tide on June 22, 2005, is shown in Fig. 7. The gray areas in the figure indicate the cells that were dry at that time step. It is

Z. Yang, T. Khangaonkar / Ocean Modelling 28 (2009) 34–49

39

Fig. 4. Comparisons of simulated and observed water surface elevation (a) and along-channel velocity (b) at station NF.

seen that a significant portion of Skagit Bay became dry during low tide in the spring tide periods, which was consistent with the general observation of wetting–drying processes in the bay and the LIDAR data (Fig. 2). Model verification was conducted from May 17–31, 2006, corresponding to high-river-flow conditions. The model configuration and all parameters in the model during verification remained the same as in model calibration. Comparison of model results of water surface elevation, velocity, and salinity to field data at station S6 are presented in Fig. 8. Velocities and salinities were measured at about 1.4 and 0.7 m above the bottom, respectively. Model simulations of tidal elevation and velocity agreed with observed data reasonably well. The distribution of water surface elevation at station S6 was very similar to station SB. Tidal currents at S6 were stronger than those at SB because station S6 was at a deeper level closer to the center of the deep channel in Skagit Bay. Salinities at station S6 generally remained high and did not show large variations as observed in station SB. This is because S6 is located north of the North Fork Skagit River. The model over-predicted salinity after May 24, 2006, following the high-flow event. One of the reasons that accounts for this error is that the jetty that trains the river flow and sediment away from the mouth of Swinomish Channel was treated as a solid boundary in the model. In reality, however, the jetty is leaking in some locations and also has a weir for fish passage. Considerable fresh water, therefore, that leaked through the jetty to the northern side during high-flow conditions were not accurately represented in the model. During flood tide, the freshwater plume was being pushed back to the river channel while during low tide, the freshwater plume was blocked by the jetty near the mouth and was forced to move to the south. Small drops observed in salinity time histories in station S6 were caused in part by the freshwater being transported up to the north around the jetty during flood tides as well as leaking through the jetty (Fig. 8c). Error statistics of model calibration (Station NF and SB) and verification (Station S6) for water surface elevations, velocities, and salinities are provided in Table 1. AME and RMSE denote abso-

lute-mean error and root–mean–square error, respectively. In general, error statistics of model predictions are reasonable. Errors for water surface elevations and velocities at Station S6 are greater than those in Station NF and SB, probably due to the complexity of geometry and tidal flats near Station S6. Although relatively large salinity errors are present at Station SB, model results follow the measured distribution trend very well (Fig. 5c). The large errors are mainly contributed to the large variations and sharp drops of the salinities in a very short time. 4.2. Effect of river flow on circulation and salinity distribution Estuarine circulation and salinity distribution are generally affected by the magnitude of freshwater inflow. The effects of river flow on the circulation and salinity transport in Skagit Bay are investigated based on the model results for calibration (low-river-flow) and verification (high-river-flow) periods. The unique geometric setting of the Skagit River estuary and bay has significant influence on the distributions of the tidal currents and salinity in Skagit Bay. The deep channel along the coastline of Whidbey Island and the large tide flat play an important role in tidal mixing and the interactions between the freshwater plume and the saline water from Puget Sound. The fact that the orientation of the Skagit River plume is nearly normal to the orientation of the deep channel constrains the spreading of the river plume and results in strong horizontal gradients in tidal currents and salinities (fronts). Fig. 9 shows surface currents and salinity distributions during flood tide and ebb tide under low-river-flow conditions on June 22, 2005 (Fig. 9a and b) and high-river-flow conditions on May 17, 2006 (Fig. 9c and d). During flood tide, strong tidal currents are present along the deep channel, and currents over the tide flat area propagate up the bay front, pushing the river plumes from North and South Fork back to the mouth of the river channels (Fig. 9a). Salinities in the middle region of the bay front along the shoreline are generally higher than the north and south sides because the North and South Fork river plumes are separated by Fir Island. During ebb tide, freshwater primarily flows out through

40

Z. Yang, T. Khangaonkar / Ocean Modelling 28 (2009) 34–49

Fig. 5. Comparisons of simulated and observed water surface elevation (a), along-channel velocity (b), and salinity (c) at station SB.

Fig. 6. Comparisons of simulated and observed salinity profiles at station SB.

multiple tidal channels, and velocities tend to zero as the tidal flats become dry (Fig. 9b). Compared to low-flow conditions, the river

plumes were found to spread much further out to the tidal flats under high-river-flow conditions (Fig. 9c and d). In particular, the freshwater plume from the South Fork occupied most of Skagit Bay and extended further out of Skagit Bay during low tide. Horizontal gradients of salinity in the North Fork freshwater plume are much stronger than those in the South Fork plume because the North Fork branch carries about two-thirds of the total Skagit River flow (Pickett, 1997). The South Fork river plume occupies a much larger area in the tide flat than the North Fork river plume. This is because the North Fork Skagit River enters Skagit Bay through a single river channel while the South Fork enters the bay through multiple channels that spread the freshwater over a much larger area. The other factor limiting the size of the North Fork river plume is the coastline of Whidbey Island, which is closer to the North Fork than the South Fork. It is also seen that the river plume front terminates near the deep channel in the bay because of mixing and transport caused by the strong cross flow caused by tidal currents along the deep channel. To examine the salinity intrusion and stratification in the Skagit River estuary, salinity distribution along North Fork branch from the western shoreline of Skagit Bay to station NF is shown in Fig. 10 under low-flow (June 22, 2005) and high-flow (May 17, 2006) conditions. Salinity transects were generated for four different tidal phases (strong ebb, lower low water, strong flood, and higher

Z. Yang, T. Khangaonkar / Ocean Modelling 28 (2009) 34–49

41

Fig. 7. Simulated surface velocities and tidal flats during low tide (gray areas indicate dry cells).

high water) during spring tide. In low-river-flow conditions, strong salinity stratification was captured in the model during ebb tide as freshwater flowed downstream in the North Fork river channel and the tidal flats (Fig. 10a). Salinity vertical variation was as high as 16 ppt over a shallow water depth of 2 m. Salinity in the deep channel remained relatively high, indicating that the water mass in the deep channel was mainly dominated by saltwater even during ebb tide under low-flow conditions. During low tide, the freshwater plume was pushed further out of the tidal flats, and a strong horizontal gradient of salinity was developed near the bank of the deep channel in Skagit Bay (Fig. 10b). Salinity in the deep channel was partially mixed because some of the freshwater was transported to and mixed in the surface layer in the deep channel. As the tide reversed to flood tide, the freshwater front that was well-mixed vertically reversed back upstream (Fig. 10c). A signal corresponding to the upwelling of saltwater from the deep channel to the shallow tidal flats was observed near the sharp bank of the deep channel. At high tide, the salt intrusion point reached upstream as far as approximately 500 m downstream of station NF near River Mile 2 in the North Fork branch. A strong salt wedge developed as saltwater intruded upstream from the bottom of the water column (Fig. 10d). Under high-river-flow conditions, the river plume was seen to spread further out of the tide flat region and occupied the surface layer of the deep channel in Skagit Bay. Stratification in the deep channel became stronger during ebb and low tides (Fig. 10e and f). Salinity profiles were well-mixed in the water column during flood tide, similar to the low-river-flow conditions, except the river front was much further out in the bay (Fig. 10g). The salt intrusion point retreated downstream close to the tide flat region, about 500 m further downstream relative to the low-river-flow conditions (Fig. 10h). Both horizontal and vertical salinity gradients during flood and high tides were much stronger under high-river-flow conditions than those under low-river-flow conditions.

4.3. Effect of wind on circulation and transport processes Circulation in the Skagit Bay region is dominated by tides and river flow, and baroclinic motion can be important in the upstream river channels and in regions where strong stratification exists. However, over the large-tide-flat areas in Skagit Bay, currents are generally small because of the shallow water depths and strong wind events can significantly affect the tidal circulation. Jones and Davies (2008) applied an unstructured grid model to examine the effect of wind on tidal circulation in the nearshore region of the eastern Irish Sea. Study results showed that wind induced surge causes significant modification of tidal elevations, currents and wetting and drying process in the shallow coastal region. Collins et al. (1998) showed that flow structure in very shallow waters can be affected by wind action. In particular, current speeds measured within the upper part of the water column may cause deviations from the standard assumption of logarithmic profile. Dyer et al. (2000) showed that circulation and suspended sediment transport processes were very sensitive to the wind. A complete reversal of transport processes going from the conventional ebb dominance to flood dominance can occur. This is compounded because sediments can be eroded in the shallows by even small waves when onshore winds are present (Lumborg et al., 2006). Ridderinkhof et al. (2000) also pointed out that when the wind speed increases, the net transport over tidal flats increases drastically, and the direction is dependent on wind direction and the location of source material. During the simulation period from June 7–23, 2005, two episodic wind events occurred with wind speed greater than 10 m/s. Strong wind from the north occurred in the early morning of June 12, 2005 (Julian Day 163), and strong wind from the southeast was observed in the mid-day of June 13, 2005 (Julian Day 164) (Fig. 11). To investigate the wind effect on circulation in Skagit Bay,

42

Z. Yang, T. Khangaonkar / Ocean Modelling 28 (2009) 34–49

Fig. 8. Comparisons of simulated and observed water surface elevation (a), along-channel velocity (b), and salinity (c) at station S6.

Table 1 Model calibration and verification error statistics Station

NF SB S6

Water level (m)

Velocity (m/s)

Salinity (ppt)

AME

RMSE

AME

RMSE

AME

RMSE

0.093 0.097 0.148

0.118 0.115 0.189

0.042 0.094 0.150

0.057 0.116 0.181

0.00 3.03 1.62

0.00 4.76 2.03

model runs with and without wind effects for the same simulation period were conducted. Simulated surface velocities and salinities with and without wind effects at 4:00 AM of June 12, 2005, corresponding to ebb tide, are shown in Fig. 12. Without wind effects (Fig. 12a), the North Fork Skagit river plume moved to the southwest direction aligned with the river-channel orientation, and the South Fork Skagit River plume moved to the west across the tidal flats. Under strong north-wind conditions, ebb currents in the deep channel and tidal flats were enhanced in the south and southeast directions, which pushed the freshwater from the North Fork branch to the southeast along the bay front shoreline (Fig. 12b). The freshwater plume discharged from the South Fork branch was forced to turn to the southwest direction. Fig. 13 shows the surface velocity and salinity distributions in Skagit Bay with and without wind effects at 12:00 PM on June

13, 2005, which also corresponded to ebb tide. Current directions in the tide flat area are towards to the northwest into the north end of the bay. Freshwater from the South Fork branch was transported to the North Fork along the bay front coastline. At the same time, freshwater discharged from the North Fork was pushed back to the mouth of the South Fork. Without wind effects, the surface velocity and salinity distributions were somewhat similar to those without wind effects on June 12, 2005 (Fig. 12a). However, under strong southeast wind conditions, the freshwater plume of the North Skagit River was depressed back to the river mouth, and the South Fork plume was pushed to the northwest along the shoreline of Fir Island. Currents in the tide flat area mainly flowed to the northwest (Fig. 13b). The sensitivity analysis for the wind effect on current and salinity distribution in Skagit indicated that wind can alter the current and salinity distribution patterns significantly in the tide flat region. Under northwest and southeast wind conditions, salinity in the region along the coastline of Fir Island could be increased because of the transport of river plumes. 4.4. Residual circulation and mean salinity distribution To examine long-term mean velocities and salinity distribution in Skagit Bay, temporal averages of velocities and salinities over a spring-neap tidal cycle were calculated under low- and

Z. Yang, T. Khangaonkar / Ocean Modelling 28 (2009) 34–49

43

Fig. 9. Surface salinities and velocities at mid-flood and mid-ebb tides in low-flow conditions (a, b) and high-flow conditions (c, d).

high-river-flow conditions (Fig. 14). The depth-averaged velocities and salinities were used in this calculation. Apparently, a high salinity and low velocity regime exists near the bay front along the coastline of Fir Island between the North and South Forks. In general, residual currents flow towards the northwest and to the north end of the bay. Residual currents near the mouths of the North and South Forks of Skagit are mainly dominated by the river flows, and currents in the tidal flats are small. Mean current and salinity under high-river-flow conditions are very similar to lowriver-flow conditions, except salinities are reduced about 5 ppt in the tide flat region, and velocity magnitudes near the river mouths are increased proportional to the river flow. Salinity intrusion in estuaries is controlled by various mechanisms, including tidal dispersion, vertical mixing, geometry, and river inflow. Extensive studies have been conducted on this topic through analytical solutions or observed data analysis (Pritchard, 1952; Hansen and Rattray, 1965; MacCready, 1999; Monismith et al., 2002). At steady-state, or tidally averaged (spring-neap tidal

cycle) conditions, the salinity intrusion length scale Xs as a function to river inflow Q can be described by a power-law relation (Hansen and Rattray, 1965; Monismith et al., 2002):

Xs  Q a

ð8Þ

where a is the power dependence coefficient. Previous studies showed that the power dependence coefficient a varies in a wide range in different estuarine conditions. Abood (1974) found that a  ? for low-flows in the Hudson River, which matched the theory derived by Hansen and Rattray (1965), but a  1 for high-flow conditions. Oey (1984) showed a  1/5 based on observed data in the Hudson River for the highest flows. Monismith et al. (2002) calculated a  1/7, based on 21 years of observed salinity data in San Francisco Bay, and pointed out that the weaker dependence of salinity intrusion on flow is due to both the geometry of San Francesco Bay and the effects of stratification on vertical mixing. A preliminary estimate of salinity intrusion in the Skagit River was calculated using model results based on Eq. (8). The mean

44

Z. Yang, T. Khangaonkar / Ocean Modelling 28 (2009) 34–49

Fig. 10. Transect salinities along North Fork branch of Skagit River at ebb, low, flood and high tides in low-flow conditions (a, b, c, d) and high-flow conditions (e, f, g, h).

Fig. 11. Wind time-series from June 9, 2005 to June 15, 2005 at Paine Field, WA.

salinities averaged over a spring-neap tidal cycle along the North Fork transect branch under low- and high-river-flow conditions were calculated (Fig. 15). Model results show that the Skagit River estuary is stratified in both low- and high-river-flow conditions. The mean salt intrusion point moves up and down the river over a distance of around 1200 m between low- and high-river-flow conditions. Assuming relation (8) holds true in the low-flow and high-flow conditions in the Skagit River, the power-law dependence coefficient a can be estimated by

a  log









X sLow Q sLow  log X sHigh Q sHigh

1 ð9Þ

From Fig. 15, we can determine XsLow = 6600 m for low-flow QLow = 325 m3/s and XsHigh = 5400 m for high-flow QHigh = 714 m3/ s. Substituting these values into Eq. (9) gives a  1=4 , which is slightly smaller than the value of ? derived from Hansen and Rattray’s theory. 5. Summary In this paper, an unstructured grid coastal ocean model was applied to simulate the tidal circulation, salinity stratification, and freshwater plume dynamics in the intertidal region of the Skagit River estuary in Puget Sound. The unstructured grid model grid

Z. Yang, T. Khangaonkar / Ocean Modelling 28 (2009) 34–49

45

Fig. 12. Comparisons of surface velocities and salinities without wind (a) and with north wind (b) at 4:00 AM on June 12, 2005.

allowed variable grid resolution in the intertidal zone and in the multiple distributary tidal channels. The model reproduced the tidal currents, salinity stratification, and wetting–drying processes over tidal flats as observed in the field data and provided a good understanding of the overall physical process in the Skagit River estuary.

During ebb tide, the freshwater plume generally moved out a considerable distance downstream of the river because of the high river discharge of the Skagit River. Salinity intrusion in the main river branch (North Fork) reached an upstream distance of approximately River Mile 2 under low-river-flow conditions.

46

Z. Yang, T. Khangaonkar / Ocean Modelling 28 (2009) 34–49

Fig. 13. Comparisons of surface velocities and salinities without wind (a) and with southeast wind (b) at 12:00 PM on June 13, 2005.

During high-river-flow conditions, the salt intrusion point was further pushed out by the strong river plume, approximately 1000 m downstream from the river mouth (RM 0). Strong stratification and destratification processes occur in the estuary during a tidal cycle. While the river plume occupies a significant portion of the tidal

flats in Skagit Bay, the heavily diked Fir Island separates the river plumes from the North and South Forks and results in a relatively high salinity and low velocity zone in the bay front of Fir Island. Tidal residual currents in Skagit Bay are in the range of 0.1–0.2 m/s in most of the tide flat area under high-river-flow

Z. Yang, T. Khangaonkar / Ocean Modelling 28 (2009) 34–49

47

Fig. 14. Temporal averages of depth-averaged velocities and salinities over spring-neap tidal cycle in low-flow (a) and high-flow (b) conditions.

conditions and less than 0.1 m/s in low-river-flow conditions. Residual currents are generally in the northwest direction towards to the north end of the bay. Although circulation in Skagit Bay is dominated by barotropic tide and baroclinic motion, model simulations indicated that strong winds can affect the circulation in

the shallow tidal flats. Strong wind events can change the current directions completely and salinity distribution in the tidal flats of Skagit Bay. Preliminary modeling analysis indicated that the salinity intrusion length scale in the Skagit River estuary is proportional to the

48

Z. Yang, T. Khangaonkar / Ocean Modelling 28 (2009) 34–49

Fig. 15. Temporal averages of transect salinities along North Fork branch of Skagit River over spring-neap tidal cycle in low-flow (a) and high-flow (b) conditions.

river flow to the 1=4 power. To calculate the salinity intrusion more accurately as a function of river flow, tidal dispersion, and the fortnight timescale of spring-neap cycles, more detailed model simulations and data analysis are required for future studies. Acknowledgements This study was funded by the Skagit Watershed Council (SWC) and the Skagit River System Cooperative (SRSC) through grants from the Washington Salmon Recovery Funding Board. We thank Ms. Shirley Solomon of SWC and Dr. Greg Hood of SRSC for providing reviews and useful comments throughout the duration of this study. We also thank Dr. Eric Grossman of USGS for providing part of the LIDAR bathymetry data in Skagit Bay and oceanographic data for model verification. References Abood, K.A., 1974. Circulation in the Hudson estuary. In: Roels, O.A. (Ed.), Hudson River Colloquium. Annals of the New York Academy of Science, vol. 250. NY Academy of Science, pp. 38–111. Aoki, K., Isobe, A., 2007. Application of finite volume coastal ocean model to hindcasting the wind-induced sea-level variation in Fukuoka bay. Journal of Oceanography 63 (2), 333–339. Babson, A.L., Kawase, M., MacCready, P., 2006. Seasonal and interannual variability in the circulation of Puget Sound, Washington: a box model. Atmosphere-Ocean 44 (1), 29–45. Cannon, G.A., 1983. An Overview of Circulation in the Puget Sound Estuarine System. NOAA Technical Memorandum. ERL PMEL-48, 30 pp. Casulli, V., Zanolli, P., 1998. A three-dimensional semi-implicit algorithm for environmental flows on unstructured grids. In: Proc. of Conf. on Num. Methods for Fluid Dynamics. University of Oxford. Casulli, V., Zanolli, P., 2002. Semi-implicit numerical modelling of non-hydrostatic free-surface flows for environmental problems. Mathematical and Computer Modelling 36, 1131–1149. Cheng, R.T., Casulli, V., 2004. Modeling a three-dimensional river plume over continental shelf using a 3D unstructured grid model. In: The Proceedings of the 8th Inter. Conf. on Estuarine and Coastal Modeling, Monterey, CA, November 2003, pp. 1027–1043.

Chen, C., Rawson, M., 2005. An ecosystem management model system for the Satilla River Estuary, Georgia. Oceans, 2005. In: Proceedings of MTS/IEEE, vol. 1, pp. 622–632. Chen, C., Beardsley, R.C., Cowles, G., 2006. An unstructured grid, finite-volume coastal ocean model (FVCOM) system. Special issue entitled ‘‘Advance in Computational Oceanography”. Oceanography 19 (1), 78–89. Chen, C., Liu, H., Beardsley, R.C., 2003. An unstructured, finite-volume, threedimensional, primitive equation ocean model: application to coastal ocean and estuaries. Journal of Atmospheric and Oceanic Technology 20, 159–186. Chen, C., Xue, P., Ding, P., Beardsley, R.C., Xu, Q., Mao, X., Gao, G., Qi, J., Li, C., Lin, H., Cowles, G., Shi, M., 2008. Physical mechanisms for the offshore detachment of the Chanjiang diluted water in the East China Sea. Journal of Geophysical Research 113, C02002. doi:10.1029/2006JC003994. Collins, M.B., Ke, X., Gao, S., 1998. Tidally-induced flow structure over intertidal flats. Estuarine, Coastal and Shelf Science 46, 233–250. Dyer, K.R., Christie, M.C., Feates, N., Fennessy, M.J., Pejrup, M., van der Lee, W., 2000. An investigation into processes influencing the morphodynamics of an intertidal tide flat, the Dollard Estuary, The Netherlands. I: hydrodynamics and suspended sediment. Estuarine, Coastal and Shelf Science 50, 607–625. Finlayson, D.P., 2005. Combined Bathymetry and Topography of the Puget Lowland, Washington State. University of Washington. Available from: . French, J.R., Spencer, T., 1993. Dynamics of sedimentation in a tide dominated backbarrier salt marsh. Northfolk, UK. Marine Geology 110, 315–331. Geyer, W.R., Trowbridge, J.H., Bowen, M.M., 2000. The dynamics of a partially mixed estuary. Journal of Physical Oceanography 30, 2035–2048. Gorman, R.M., Neilson, C.G., 1999. Modelling shallow water wave generation and transformation in an intertidal estuary. Coastal Engineering 36, 197–217. Greenberg, D.A., Shore, J.A., Page, F.H., Dowd, M., 2005. A finite element circulation model for embayments with drying intertidal areas and its application to the Quoddy region of the Bay of Fundy. Ocean Modeling 10, 211–231. Grossman, E., Stevens, A., Gelfenbaum, G., Curran, C., 2007. Nearshore Circulation and Water Column Properties in the Skagit River Delta, Northern Puget Sound, Washington – Juvenile Chinook Salmon Habitat Availability in the Swinomish Channel: U.S. Geological Survey Scientific Investigations Report 2007-5120, 96 pp. Hansen, D.V., Rattray Jr., M., 1965. Gravitational circulation in straits and estuaries. Journal of Marine Research 23, 104–122. Hood, W.G., 2004. Indirect environmental effects of dikes on estuarine tidal channels: thinking outside of the dike for habitat restoration and monitoring. Estuaries 27, 273–282. Hood, W.G., 2006. A conceptual model of depositional, rather than erosional, tidal channel development in the rapidly prograding Skagit River Delta (Washington, USA). Earth Surface Processes and Landforms 31, 1824–1838.

Z. Yang, T. Khangaonkar / Ocean Modelling 28 (2009) 34–49 Hood, W.G., 2007. Scaling tidal channel geometry with marsh island area: a tool for habitat restoration, linked to channel formation. Water Resources Research 43, W03409. doi:10.1029/2006WR005083. Isobe, A., Beardsley, R.C., 2006. An estimate of the cross-frontal transport at the shelf break of the East China Sea with the Finite Volume Coastal Ocean Model. Journal of Geophysical Research 111, C03012. doi:10.1029/2005JC003290. Ji, Z.G., Morton, M.R., Hamrick, J.M., 2001. Wetting and drying simulation of estuarine processes. Estuarine, Coastal and Shelf Science 53 (5), 683–700. Jones, J.E., Davies, A.M., 2008. On the modification of tides in shallow water regions by wind effects. Journal of Geophysical Research 113, C05014. doi:10.1029/ 2007JC004310. Lumborg, U., Andersen, T.J., Pejrup, M., 2006. The effect of Hydrobia ulvae and microphytobenthos on cohesive sediment dynamics on an intertidal mudflat described by means of numerical modeling. Estuarine, Coastal and Shelf Science 68, 208–220. MacCready, P., 1999. Estuarine adjustment to changes in river flow and tidal mixing. Journal of Physical Oceanography 29, 708–726. Mellor, G.L., Yamada, T., 1982. Development of a turbulence closure model for geophysical fluid problems. Reviews of Geophysics and Space Physics 20, 851– 875. Mellor, G.L., Blumberg, A., 2004. Wave breaking and ocean surface layer thermal response. Journal of Physical Oceanography 34, 693–698. Monismith, S.G., Kimmerer, W., Burau, J.R., Stacey, M.T., 2002. Structure and flowinduced variability of salt intrusion in San Francisco Bay. Journal of Physical Oceanography 32, 3003–3019. Oey, L.Y., 1984. On steady salinity distribution and circulation in partially mixed and well mixed estuaries. Journal of Physical Oceanography 14, 629–645. Oey, L.Y., 2006. An OGCM with movable land–sea boundaries. Ocean Modelling 13, 176–195. Pickett, P.J., 1997. Lower Skagit River Total Maximum Daily Load Water Quality Study, Publication Number 97-326. Washington Department of Ecology, Olympia, WA. Pritchard, D.W., 1952. Salinity distribution and circulation in the Chesapeake Bay estuarine system. Journal of Marine Research 11, 106–123.

49

Ralston, D.K., Stacey, M.T., 2005. Stratification and turbulence in subtidal channels through intertidal tidal flats. Journal of Geophysical Research 110, C08009. doi:10.1029/2004JC002650. Ridderinkhof, H., van der Ham, R., van der Lee, W., 2000. Temporal variations in suspended sediment concentrations and current velocities above a tidal flat in the Ems-Dollard Estuary: variability on different timescales. Continental Shelf Research 20 (12–13), 1479–1493. Simenstad, C.A., 1983. The Ecology of Estuarine Channels of the Pacific Northwest: A Community Profile, Rep. FWS/OBS 83/05. U.S. Fish and Wildlife Serv., Washington, DC. Smagorinsky, J., 1963. General circulation experiments with the primitive equations. I: the basic experiment. Monthly Weather Review 91, 99–164. Weisberg, R.H., Zheng, L., 2006a. Circulation of Tampa Bay driven by buoyancy, tides, and winds, as simulated using a finite volume coastal ocean model. Journal of Geophysical Research 111, C01005. doi:10.1029/2005JC003067. Weisberg, R.H., Zheng, L., 2006b. Hurricane storm surge simulations for Tampa Bay. Estuaries and Coasts 29 (6A), 899–913. Wiggins, W.D., Ruppert, G.P., Smith, R.R., Reed, L.L., Courts, M.L., 1997. Water Resources Data, Washington, Water Year 1997. U.S. Geological Survey, Water Data Report WA-97-1, Tacoma, Washington. Yang, Z., Khangaonkar, T., 2008. Modeling of salt Intrusion, intertidal mixing, and circulation in a braided estuary. Journal of Coastal Research 52, 171–180. Yang, Z., Sobocinski, K.L., Khangaonkar, T., Thom, R., Fuller, R., Heatwole, D., submitted for publication. Hydrodynamic and ecological assessment of nearshore restoration: a modeling study. Ecological Modeling. Zhang, Y.L., Baptista, A.M., Myers, E.P., 2004. A cross-scale model for 3D baroclinic circulation in estuary-plume-shelf systems. I: formulation and skill assessment. Continental Shelf Research 24, 2187–2214. Zhao, L., Chen, C., Cowles, G., 2006. Tidal flushing and eddy formation in Mount Hope Bay and Narragansett Bay: an application of FVCOM. Journal of Geophysical Research 111, C10015. doi:10.1029/2005JC003135. Zheng, L., Chen, C., Liu, H., 2003. A modeling study of the Satilla River Estuary, Georgia. Part I: flooding/drying process and water exchange over the salt marsh-estuary-shelf complex. Estuaries 26 (3), 651–669.