Cross-Scale Baroclinic Simulation of the Effect of Channel Dredging in an Estuarine Setting

Abstract

Holistic simulation approaches are often required to assess human impacts on a river-estuary-coastal system, due to the intrinsically linked processes of contrasting spatial scales. In this paper, a Semi-implicit Cross-scale Hydroscience Integrated System Model (SCHISM) is applied in quantifying the impact of a proposed hydraulic engineering project on the estuarine hydrodynamics. The project involves channel dredging and land expansion that traverse several spatial scales on an ocean-estuary-river-tributary axis. SCHISM is suitable for this undertaking due to its flexible horizontal and vertical grid design and, more importantly, its efficient high-order implicit schemes applied in both the momentum and transport calculations. These techniques and their advantages are briefly described along with the model setup. The model features a mixed horizontal grid with quadrangles following the shipping channels and triangles resolving complex geometries elsewhere. The grid resolution ranges from ~6.3 km in the coastal ocean to 15 m in the project area. Even with this kind of extreme scale contrast, the baroclinic model still runs stably and accurately at a time step of 2 min, courtesy of the implicit schemes. We highlight that the implicit transport solver alone reduces the total computational cost by 82%, as compared to its explicit counterpart. The base model is shown to be well calibrated, then it is applied in simulating the proposed project scenario. The project-induced modifications on salinity intrusion, gravitational circulation, and transient events are quantified and analyzed.

1. Introduction

Many estuaries and coastal regions are influenced by anthropogenic activities [1,2,3,4]. Accordingly, water resources managers and engineers need to assess the impacts of any proposed activity prior to its inception. Numerical models are frequently required in understanding the effects of existing projects and predict the changes under future plans [5,6,7,8]. This demands good accuracy, efficiency and robustness from the modeling tool, which can be challenging in an estuarine setting because: (1) human activities often involve small-scale alterations on natural conditions, such as the construction of hydraulic structures or localized changes in shoreline and bathymetry features, which require a high-resolution model grid to resolve; (2) the spatial scales of such systems are intrinsically linked, which can range from kilometer-scale shelf dynamics to meter-scale local dynamics at the project site. Therefore, it is preferable to include all major processes in a single model domain. To reconcile the needs between local fine resolutions and a large model domain, a grid-nesting approach is often adopted, especially for models based on structured grids [9]. However, this may not be ideal because (1) gird-nesting is cumbersome; (2) the non-smooth transition at sub-model interfaces introduces errors; (3) a multi-model configuration may have difficulties where the connectivity from large to small spatial scales is the study object itself, e.g., fish larvae transport along the tributary-estuary-ocean continuum, or an engineering project that traverses different spatial scales as described in this paper.

Unstructured-grid (UG) models are promising for simulating cross-scale processes due to their superior ability in local refinement and boundary-fitting [10]; the latter is especially important in near-shore regimes as the underlying geometry and bathymetry can be very complex. Nevertheless, the large scale-contrasts in estuarine applications still pose challenges for the UG model in terms of efficiency and accuracy. The former is often constrained by the highest resolution in the model grid; the latter requires both high-order schemes and a faithful representation of the underlying bathymetry. As the developing team of the Semi-implicit Cross-scale Hydroscience Integrated System Model (SCHISM [11]), we have been making progress in dealing with these issues to improve the model’s cross-scale capability. In this paper, we demonstrate an efficient and accurate cross-scale baroclinic simulation that utilizes newly developed techniques, including a mixed quadrangle-triangle grid that resolves features in the horizontal dimension [11], a flexible vertical grid that faithfully follows bathymetry while reducing pressure gradient error and unphysical diapycnal mixing [11,12], and an implicit transport solver that maintains efficiency under contrasting scales [11,13]. The main purpose of this study is to use these new techniques to quantify the effect of channel dredging and land expansion on estuarine hydrodynamics and to assist decision makers.

The rest of the text is structured as follows: Section 2 briefly describes the study domain, observation assets, and the proposed engineering plans; Section 3 provides the details on model setup; Section 4 validates the baseline model (without project) against observation; Section 5 uses the validated model to simulate the proposed project scenario and analyzes its effects on hydrodynamics; Section 6 summarizes the paper.

2. Study Object

2.1. Spatial Domain

The domain of interest is a “subtributary-tributary-estuary-ocean” continuum of the Chesapeake Bay (abbreviated as “Bay” hereafter). As illustrated in Figure 1, the ocean part of the domain includes the continental shelf, the shelf break, and the shelf slope. The motivation for using this large domain is to adequately represent the shelf dynamics and incorporate some signals from the Gulf Stream through the ocean boundary. The estuarine part includes the entire Chesapeake Bay and its major tributaries. A complex configuration of deep channels and shallow shoals creates unique circulation patterns near the Bay mouth [14]. Salinity intrusion mainly occurs through the deep channels, some of which are artificially maintained for navigation and are of interest to the current project (labeled in Figure 2).

This domain poses a stringent test on the model’s cross-scale capability due to its complex geometry and bathymetry, as well as its highly variable stratification [15]. Moreover, the Bay processes are closely linked to both the small-scale processes occurring in the tributaries/sub-tributaries (such as the point/non-point sources of freshwater inputs) and the large-scale processes in the coastal ocean (such as the wind-induced Ekman transport) [14,16,17,18,19].

2.2. Project Specification

The proposed engineering project involves the deepening of the Atlantic Ocean Channel, the Thimble Shoal Channel, the Newport News Channel, and the Norfolk Harbor Channel (see locations in Figure 2). The changes in bathymetry are illustrated in Figure 3. The project also involves a land expansion near the Craney Island (Figure 3b,c), and a new channel will be dredged (referred to as “access channel” hereafter) between the eastern edge of the expanded land to similar depths of the adjacent Norfolk Harbor channel. Other minor changes include the removal of small shoreline structures (Figure 3c).

2.3. Observational Assets

The Chesapeake Bay is a well-instrumented system with extensive in-situ and remote sensing observational assets. EPA’s Chesapeake Bay Program Office (CBPO) maintains a network of stations throughout the Bay and regularly conducts surveys (usually bi-weekly in summer or monthly in winter) on both physical and biogeochemical variables. The modeled salinity and temperature profiles are compared to CBPO observations. The National Oceanic and Atmospheric Administration (NOAA) of United States maintains over 40 tide gauges in and around the Bay and the adjacent shelf. In addition, it also maintains Acoustic Doppler Current Profilers (ADCPs) at a few stations in the Bay during some recent years, which are used to validate our modeled current profiles. The locations of all monitoring stations used for model validation are marked in Figure 1b. The main bathymetry source is FEMA [20], with further modifications provided by the Norfolk District of the U. S. Army Corps of Engineers that reflect recent dredgings.

3. Model Setup

The unstructured-grid (UG) SCHISM model [11] is used to simulate the circulation near the project site, including the lower Chesapeake Bay, the James River and the Elizabeth River. The flexibility provided by the UG model allows us to simulate the processes in the entire Bay and continental shelf as a whole, thus greatly simplifying the boundary condition (B.C.) requirement.

The model grid has 24,520 nodes, 40,593 mixed triangular-quadrangular elements in the horizontal dimension, with resolutions varying from 15 m to 6.3 km; typical element sizes in different regions are illustrated in (Figure 4). Note that fine-sized quadrilateral elements are used to delineate the deep channels (Figure 4b–d). A flexible LSC${}^2$ (Localized Sigma Coordinates with Shaved Cells) vertical grid [12] is used, which efficiently covers depths from deep to shallow regions, with a maximum of 80 levels used at the maximum depth of ~3.6 km near the east ocean boundary, 29–41 levels for the shipping channel, and a minimum of 5 levels for the shallow area. On average, there are 28 levels used in the vertical dimension. Figure 5 illustrates the flexible configuration of the vertical grid on a representative transect from the deep ocean into the estuary. Note that the ocean part is refined near the surface to capture the freshwater plume and the estuary part is refined at mid-depths to resolve the thin pycnocline. Another nice feature of the vertical grid is the “degenerate elements” applied at large slopes (e.g., the near-bottom triangles in the zoomed-in view of Figure 5b). The benefits of this configuration include: (1) keeping vertical layers in upper levels more horizontal than those in a traditional terrain following grid, reducing pressure gradient errors; (2) following terrain, but effectively cutting off the unphysical diapycnal mixing between deep and shallow regions [12].

The use of flexible horizontal and vertical grids ensures an accurate representation of key bathymetric features as well as the project’s modifications, which are essential for simulating the gravitational circulation and assessing the project impact. Another important consideration is the computational cost. Although only one major scenario is presented in this paper, multiple scenarios were investigated that differ in dredging locations and depths, which resulted in a large number of model runs. Due to the large scale contrast in the model grid, the simulation speed was barely practical using traditional explicit transport schemes. This issue is fixed by the newly developed high-order implicit transport scheme in SCHISM [11,13], which in this case reduces the total computational cost to 1/6 of that with the explicit solver. Now the model runs 360 times faster than real time with 80 Xeon cores (on William and Mary’s Whirlwind cluster) at a time step of 120 s.

The model is forced by USGS-measured flows from the 7 major tributaries of the Bay (Susquehanna, Patuxent, Potomac, Rappahannock, York, James, and Choptank). Additionally, non-point sources of freshwater inputs are specified at all land boundaries, which is based on the operational watershed model [21] of the Chesapeake Bay program. At the air-water interface, the model is forced by the wind, atmospheric pressure, and heat fluxes predicted by North American Regional Reanalysis (NARR). At the outer ocean boundary, the elevation B.C. is obtained from inverse-distance interpolated values from two tide gauges at Duck, NC and Lewes, DE. The salinity and temperature B.C.’s are interpolated from HYbrid Coordinate Ocean Model (HYCOM). In addition, a 20-km nudging zone near the ocean boundary is used where the salinity and temperature are relaxed to the corresponding HYCOM values in order to prevent long-term drift, with a maximum relaxation period of 1 day. The model uses the turbulence closure of k-kl [22].

4. Validation of the Baseline Model

The simulation period is from 1 July 2009 to 1 January 2014. The first half-year period in 2009 is regarded as ramp-up; the remaining 4-year period is used for model-data comparison on elevation, velocity, salinity, and temperature. It should be noted that the model assessment in this paper only focuses on the project area; while a comprehensive model assessment for the entire Chesapeake Bay and coastal ocean will be presented in a separate paper [15].

The simulated surface elevation is compared to the observations on four NOAA stations near the project site (see locations in Figure 1). Both sub-tidal and tidal components of the elevation are well captured by the model. The correlation coefficients (CC) between the modeled and observed sub-tidal elevations are above 0.95, and the mean absolute error (MAE) are below 3.7 cm (Figure 6). The averaged CC and RMSE over the four stations are 0.97 and 3.5 cm respectively. The hurricane-induced set-ups are also captured by the model (i.e., Hurricane Irene on 27 August 2011). For brevity, the time series are only shown for one year (2011 with a storm event) since the skills are similar among years. The tidal component, i.e., the amplitudes and phases of the major constituents are also well modeled, which is shown by a harmonic analysis in Figure 7. Moreover, the skill of the simulated elevation has no bias in space from the Bay mouth to the sub-tributary.

The along channel velocity is validated against the ADCP measurements at four NOAA stations (see locations in Figure 1) at multiple (13–18) vertical layers.The correlation coefficients exceed 0.7 at all stations, and are mostly between 0.8–0.95. The averaged root mean square error (RMSE) at all stations is 8 cm/s. The time series of the depth-averaged along-channel velocity during a randomly chosen period are shown in Figure 8.

The simulated salinity and temperature are compared to the monthly/bi-weekly multi-layer samples at the CBP stations (see locations in Figure 1). The overall statistics are summarized in the form of target diagrams (Figure 9). The x-axis shows the un-biased RMSE (i.e., with the mean removed), with positive values (x > 0) indicating an over-estimation in the variability of the observations, and vice versa. The y-axis shows the model bias. The RMSE’s are all within 2 PSU and 2 °C, which are small compared to the standard deviation of the salinity and temperature data used to validate the model. The vertical salinity and temperature profiles are also compared to observation at each station. For brevity, only the salinity profiles are shown at one station (CB7.4) near the Bay mouth (Figure 10). This is the most challenging station to simulate due to its (1) large maximum stratification and (2) large temporal variations in stratification. The model nicely captures these features, e.g., the stratified-mixed-stratified cycle between 26 August and 9 December in 2013 (Figure 10). For a comprehensive model-data comparison on all stations, readers are referred to a parallel study on the same model domain [23].

5. Simulation of the Proposed Scenario

The analyses in this section focus on the project-induced changes in salinity intrusion and circulation patterns. These are essential for the water quality and biological processes around the project area. A parallel study on water quality is presented by our colleagues in [24]. It should be briefly mentioned that the tidal range throughout the whole domain is only marginally increased, with the largest increase of 0.7% found in the Elizabeth River. This is consistent with other studies using analytical solutions [5]. The change is subtle because the dredged channels are much narrower than the full river width.

5.1. Salinity

To assess the salinity change, time-averaged differences are presented both in plan view of the entire project area (Figure 11) and along the axis of the Elizabeth River where the largest change is found (Figure 12). In general, an increase in salinity is observed in the Elizabeth River, the James River, and the lower Chesapeake Bay. The largest change (~2.5 PSU) is near the land expansion site near the mouth of the Elizabeth River (red color in Figure 11a). This is mainly due to the deepening (from ~5 m to ~15 m) between the eastern edge of the expanded land and the dredging in the access channel (Figure 3b,c). In other areas, the salinity increase generally occurs near the deepened channels. In particular, a channelized pattern is visible at the bottom (Figure 11a). The time-averaged salinity difference for different regions is summarized in

Table 1. Salinity increase (PSU; temporally and spatially averaged) in different regions under the project scenario.

	Surface	Bottom
Lower Bay	0.15	0.14
James River	0.23	0.26
Elizabeth River	0.86	0.86

.

Figure 12 provides a more detailed view on the salinity increases along the axis of the Elizabeth River, where the largest increase is found. The increased salinity intrusion suggests that the gravitational circulation is generally enhanced with channel dredging, especially near the mesohaline regions. Since the strength of the gravitational circulation is a function of freshwater inflow, seasonal variabilities in the salinity change are seen in all transects. The largest increase (~2 PSU) is found in Elizabeth River during winter and spring months when large freshwater inputs are common. This has implications for the local biological communities that are sensitive to salinity changes.

5.2. Non-Tidal Circulation

The previously found salinity increase is closely linked to the changes in the gravitational circulation. According to the classic theory for an idealized estuary [25], the magnitude of the two-layer exchange flow ($u_E$) is related to the water depth (H), the longitudinal salinity gradient ($\overline{s_x}$), and the effective eddy viscosity ($K_M$) as:

$$u_E=\frac{\beta \overline{s_x}{H}^3}{48K_M},$$

(1)

where $\beta$ is the salinity contraction coefficient for the salinity-induced density change, which can be treated as a constant here.

Firstly, the changes in water depth due to dredging are obviously important for the exchange flow. In the Elizabeth River, the river bed within the project zone is typically deepened by 7–16%, translating to an amplification factor of 23–56% according to Equation (1). However, the expected increase in the exchange flow is much smaller than those values, since the dredged area only consists a small portion of the full river area.

Vertical mixing is another important factor for the gravitational circulation. Since the effective eddy viscosity can be viewed as a modification on the traditional eddy viscosity by incorporating the contributions from tidal averaging [25], both the mixing patterns and the tidal currents need to be examined. The overall mixing patterns are barely changed because the dredged area is relatively small and the artificially maintained channel bed is fairly smooth before and after the project (Figure 13). On the other hand, the amplitude of the tidal flow is not significantly changed in most areas, with the exception of a local reduction near the land expansion site at the Craney Island (which will be further discussed in Section 5.3). This reduction tends to locally strengthen the gravitational circulation.

The longitudinal salinity gradient (approximately the slope of the along-thalweg bottom salinity in Figure 12) is mostly unchanged in the Elizabeth River, except for a slight reduction near the mouth that tends to weaken the gravitational circulation locally. The likely cause is that the enhanced salinity intrusion through the dredged channels nudges the longitudinal salinity distribution in the Elizabeth River from the “central” regime (large-gradient zone in the mid-estuary) towards the “outer” regime (small-gradient zone near the mouth) as defined in Hansen and Rattray’s classic work [26]. Geyer and MacCready identified the opposite trend in which the retreat of salinity intrusion tends to increase the longitudinal salinity gradient [27].

Additionally, the land expansion near the Craney Island and the newly dredged access channel modify the local flow patterns near the mouth of the Elizabeth River. Under the baseline condition, the core of the bottom inflow tends to concentrate near the bed of the narrow central channel; while in the project scenario, the inflow core becomes more diffusive as it shifts westward into the newly vacated space. Similarly, the surface outflow core also shifts from the central channel towards the west. This trend is seen in Station #1 and #2 in Figure 14, the former is representative of the new access channel and the latter is representative of the pre-existing central channel.

The land expansion also modifies a bottom “double-gyre” pattern near the Elizabeth River mouth (highlighted in Figure 15). Under the baseline condition (Figure 15a), a part of the James River inflow jet detaches from the mainstream and forms a semi-closed cyclonic gyre to the north of the Craney Island, some of the cyclonic flow later merges into an anticyclonic gyre to its east. The similar pattern is hinted in the results of Shen et al.’s James River model [28], though its low spatial resolution at this location (which is not the focus of that study) masks many details. As indicated by the present model results, the two gyres become relatively independent in the project scenario, as less flow from the cyclonic gyre merges into the anticyclonic gyre and more re-integrates into the James River jet. The separation of the two gyres is likely due to the weakening of the anticyclonic gyre, which is in turn due to the restriction of the river mouth. In the project scenario, the bottom out-flow on the west side of the Elizabeth River is no longer present, as the anticyclonic gyre moves out of the river. The corresponding changes in the vertical exchange flow pattern are seen in Station #1 of Figure 14, where the out-flux is present in most of the water column under the baseline condition; whereas a normal 2-layer exchange pattern is seen in the project scenario.

To sum up, the net effect of the above changes is an increased gravitational circulation in most parts of the Elizabeth River. Despite the multiple project modifications near the Craney Island, the total sub-tidal volume flux through the entire 2D transect of the Elizabeth River mouth is barely altered. The magnitude of the influx is increased by 3.0% and the out-flux by 4.4%, i.e., a slight increase in the exchange flow. Although not shown here, the gravitational circulation generally increases in other dredged channels in the lower Bay and the James River.

5.3. Tidal Currents and an Extreme Event

In addition to the previous analyses that mostly focus on time-averaged properties, this section further examines some transient events, including tidal currents and a storm event in 2012 (Hurricane Sandy). Tidal currents during flood and ebb are expected to be modified by the changes in channel shape and bathymetry. However, since the dredgings are mostly confined in narrow channels, their effects on tidal currents are small (mostly less than 1 cm/s), with large changes confined in the access channel near the Craney Island expansion site (Figure 16a). This is not surprising considering the drastically changed channel shape. The relocation of bed material is expected to transform the gently-sloping shoal into a steep bank, which leads to a reduction in the surface channel width by ~36%. With the tidal range essentially unchanged, the reduced surface area translates into 13% and 29% reductions in the tidal current amplitude at 2 representative stations near the land expansion site (Figure 16b).

On the other hand, extreme events such as storm surges are also of interest, since they may disrupt navigation and local biota. Two hurricanes (Irene and Sandy) are present in the simulation period. Hurricane Sandy is analyzed here due to its larger local impacts. During the course of the storm surge, the peak elevation reaches 1.6 m above the mean sea level in the lower James River (Figure 17). The project-induced changes in total elevation (tidal and non-tidal signals) are negligible. Similarly, the transient currents remain mostly the same in the project scenario, except near the land expansion site. In this region, large differences in the horizontal flow patterns are found, which are illustrated in Figure 18. The model has some errors during the course of the surge but the peak value is well captured, so the snapshot of the horizontal currents is taken at a flood tide closest to the peak elevation’s timing. The most obvious change by the project scenario is that the cyclonic eddy near the Elizabeth River mouth in the baseline case is transformed into a much smaller-scale anticyclonic eddy. This is likely because the restricted river mouth in the project scenario tends to prevent eddy formation, similar to the case of the bottom residual currents discussed earlier. Additionally, the flow magnitude near the northern edge of the Craney Island increases, and the high flow zone along the east bank of the Elizabeth River shifts southward (Figure 18, right panel). These may have implications for shoreline management. It should be noted that the two hurricanes in the simulation period did not directly hit the Chesapeake Bay. Additional scenarios need to be conducted to fully assess the potential storm risk under the project scenario; a good candidate can be Hurricane Isabel in 2003.

In summary, the project changes lead to a net increase in salinity intrusion. This is supported by an overall increase in gravitational circulation. The region most sensitive to the project changes is the Elizabeth River, due to its small channel volume and the restricted exchange with upstream freshwater at channel-narrowing points (such as Great Bridge Lock in Figure 1). The land expansion at the Craney Island leads to localized changes on both non-tidal and transient currents.

6. Conclusions

In this paper, a baroclinic model was built for the Chesapeake Bay, covering different spatial scales from the small tributaries/sub-tributaries to the adjacent coastal ocean. The model utilized recently developed techniques (such as the flexible vertical grid design and the implicit transport scheme) to achieve good accuracy and efficiency. Near the project area, the simulated elevations, depth-averaged velocity, salinities, and temperatures of the lower Bay, James River and Elizabeth River regions have an averaged RMSE of 6 cm, 8 cm/s, <2 PSU, <2 °C respectively, and the model captured key processes (salt intrusion, periodic stratification-destratification) well. Then the model was successfully applied in quantifying the effect of channel dredging and land expansion on hydrodynamics. The increased water depths in the shipping channels enhance salinity intrusion and gravitational circulation. The maximum salinity increase is over 2 PSU in the mesohaline region of the Elizabeth River in spring months. However, the increases are generally small (~0.1 PSU) away from the project area (with a distance exceeding ~30 km from the project site). Moreover, the land expansion at the Craney Island also has local effects on gravitational circulation, tidal exchange, and eddy formation. The successful simulation of the mixed effects of the channel dredging (along the ocean-estuary-river axis) and the land expansion (localized in a sub-tributary) in a coastal system exemplifies the cross-scale capability of the model and the importance of adopting multi-scale approaches in coastal management.