Ocean Modeling and Prediction

Development of a coupled hydrodynamic/biogeochemical

for ecological forecasting in Chesapeake Bay

The Chesapeake Bay is America largest and most biologically diverse estuary, home to more than 3600 species of plants, fish and animals. The Chesapeake is a commercial and recreational resource for more than 15 million people who live in its basin. However, the water quality and living resources in the Bay have been in a dramatic decline over the last few decades. A major cause for this environmental degradation is an oversupply of nutrients entering the Bay, including those from point sources like wastewater treatment plants and nonpoint sources such as fertilizer running off farmland. Because nutrients are released into the Bay through river runoffs, the eutrophication problem is often compounded by the climatic/weather effects. A computer model which simulates the hydrodynamic and biogeochemical processes in the Chesapeake Bay can be of great help in simulating and analyzing nutrient-loading and climate variability scenarios and assist managers in designing nutrient reduction and restoration strategies.

Figure 1. A satellite image of Chesapeake Bay (lower) and Delaware Bay (upper) in the Middle Atlantic Bight.

Using the state-of-art Regional Ocean Modeling System (ROMS), we have developed a new hydrodynamic/biogeochemical model for the Chesapeake Bay. We have made hindcast simulations for the period between 1995 and 2000, covering years of highly variable climatic forcing conditions. The model results have been validated against a suite of observational data, including water quality data collected by the Chesapeake Bay Program, sea level measurements obtained at tidal gauges stations, temperature, salinity and current measurements obtained at the Chesapeake Bay Observing System (CBOS), surface chlorophyll distributions obtained from aircraft remote sensing. We are grateful for funding support from CICEET/NOAA.

Hydrodynamic Model

The Regional Ocean Modeling System (ROMS) is a state-of-the-art regional ocean model that has found wide-ranging applications, including basin-scale ocean circulation in the North Atlantic basin (Haidvogel et al., 2000), shelf circulation in the California Current System (Marchesiello et al., 2003) and the Hudson River estuary (Warner et al., 2005). We have configured ROMS for the Chesapeake Bay (Li et al, 2005). Figure 2a and 2b show the model’s bathymetry and grid. The model domain is made much wider than the main stem of the Bay in order to include all major tributaries and a part of the coastal ocean to facilitate free exchange across the Bay mouth. An orthogonal curvilinear coordinate system is designed to follow the general orientation of the deep channel and the coastlines of the main stem. Coastal boundaries are specified as a finite-discretized grid via land/sea masking. The grid spacing is less than 1 km in the cross-channel (latitudinal) direction and about 2- 3 km in the along-channel (longitudinal) direction. The total number of grid points is 120x80. High resolution is placed in the mainstem of the Bay to ensure adequate resolution over the deep channel which is the main conduit for landward salt transport. Bottom topographic irregularities with horizontal scales smaller than the grid sizes are truncated by the model. The model has 20 layers in the vertical direction.

Figure 2. An orthogonal curvilinear coordinate system is designed to follow the general orientation of the deep channel and the coastlines of the main stem. The model is forced by open ocean tides, freshwater inflows at river heads, wind and heat exchange across the water surface.

Biogeochemical model

The planktonic ecosystem in Chesapeake Bay is a complex system involving many species and multiple pathways of nutrient fluxes. Rather than attempting to reproduce the observed community structure, we shall instead explore a simplest possible biogeochemical model that can capture the most fundamental processes. Towards this end, we adopt a five-compartment (nitrate, ammonium, phytoplankton, zooplankton and detritus) biogeochemical model that has been developed for studying open-ocean biogeochemistry (Fasham et al., 1990). We use a description of nutrient kinetics that allows phytoplankton to uptakes ammonium preferentially over nitrate. In addition to nutrient limitation, phytoplankton growth is subject to light and temperature limitation. We introduce detritus so that plankton mortality and zooplankton sloppy feeding are converted to detritus which sink to the bottom. This nitrogen-based biogeochemical model has been widely used in open-ocean and coastal biogeochemical investigations (Fasham, 1995; Mosian and Hofmann, 1996; Fennel et al. 2005). To apply it to shallow coastal and estuarine waters, however, we need to consider additional biogeochemical processes. Sediment loading is a major factor in determining water clarity in Chesapeake Bay. In the present model, we prescribe surface PAR values using Fisher et al. (2002) measurements and a light-attenuation coefficient using Harding et al. (2002) optical measurements. Dead plankton bodies sink into the bottom and are remineralized into ammonium by bacteria function. This remineralization process is greatly accelerated during the summer when the water temperature is high and bacteria activity is intense (e.g. Kemp et al., 1990, Kemp and Boynton, 1992). To parameterize this process, we specify a temperature-dependent remineralization rate. To account for the effects of sediment biogeochemistry on the overlying water column, we specify a denitrification rate as a linear function of the ambient nitrate concentration in bottom water, according to the laboratory experiments by Kana et al. (1998). Figure 3 shows a schematic diagram of the biogeochemical model.

Figure 3. A schematic diagram of the biogeochemical model.

Model validation against observations

In order to develop a forecasting model for Chesapeake Bay, we have validated the coupled hydrodynamic-biogeochemical model against a suite of measurements. For the hydrodynamic model, the data used for comparison include: (1) tidal elevations at tidal gauge stations; (2) tidal current measurements; (3) long-term salinity time series at the monitoring stations maintained by Chesapeake Bay Program (EPA CBP); (4) real-time current velocity measurements obtained at the buoys of Chesapeake Bay Observing System (CBOS). In addition we have used three-dimensional synoptic salinity maps (along-channel and cross-channel sections) acquired during NSF-funded hydrographic surveys. For the biogeochemical model, the data used for comparison include: (1) nitrate, ammonium, chlorophyll measurements at CBP monitoring stations; (2) surface chlorophyll maps acquired by remote-sensing aircrafts (Chesapeake Bay Remote Sensing Program CBRSP); (3) depth-integrated primary productivity and euphotic-layer chlorophyll concentrations inferred from observations and statistical productivity models (Harding et al., 2002).

(1) Tidal prediction

Because tidal currents provide a major source of mechanical energy to Chesapeake Bay, we first examine model’s prediction for tidal elevations and tidal currents. In Figure 4 we compare the observed and predicted sea-surface elevations at six tidal stations in Chesapeake Bay. Tochester Beach and Annapolis stations are in the upper Bay, Solomons Island and Lewisetta stations in the middle Bay, Windmill Point and CBBT in the lower Bay. The predicted tidal elevations are in excellent agreement with the observed tidal records. The root-mean-square (rms) error averaged over 12 stations in the Bay is 5 cm, the relative average error is 5%, the correlation coefficient is 95% and the model predictive skill is 0.95 (Li et al., 2005).

Figure 4. Comparison between the predicted (blue) and observed (red) tidal elevations at 6 tidal gauge stations in Chesapeake Bay.

(2) Salinity prediction

Salinity largely determines density distribution in Chesapeake Bay. Salinity distribution along the center-axis of Chesapeake Bay is a key measure of estuarine dynamics. Figure 5a and 5b show a comparison of the along-channel salinity distribution between the model prediction and hydrographic data obtained from Scanfish- a towed undulating vehicle equipped with multiple sensors. The model reproduces a well-mixed surface layer as shown in the data. The predicted salinity distribution beneath the surface layer also shows a good agreement with the hydrographic data. The 13 isohaline intercepts the bottom boundary at about 38.8 latitude; the 15 isohaline reaches a location near 38 latitude. As a partially-mixed estuary, the Chesapeake Bay features a two-layer circulation. The hydrodynamic model reproduces this circulation pattern, as shown in the tidally-averaged along-channel velocity field (Figure 5c). Water moves seaward in the surface layers and landward in the bottom layers. Maximum residual velocity is about 0.2 ms^-1. The level of no motion separating the two counter-flowing layers lies somewhere between 5 and 8 meters (Li et al., 2005).

Figure 5. Comparison of along-channel salinity distribution on November 1, 1996 between high-resolution Scanfish measurements (a) and ROMS model (b). Tidally-averaged residual velocity in the along-channel section (c).

To evaluate how the model captures temporal salinity variability over seasonal and interannual time scales, we located four stations in the main stem: Station CB3.3C, CB4.4, CB5.4 and CB8.1. They occupy different salinity regimes along the main axis of the Bay, ranging from nearly-fresh water in the upper Bay to shelf salinity at the Bay mouth. In addition, we located one station in each of the two large tributaries: LE2.3 in the Potomac River and LE5.5 in the James River. At each monitoring station, we obtained time series of the surface (about 1 m depth) and bottom (1 to 2 m above the bottom boundary) salinity from the CBP database. For the normal runoff year of 1997, the model appears to have captured the seasonal salinity variations well (see Figure 6). Most of the skill scores exceed 0.8 and five salinity predictions produce scores above 0.9. Only one score falls below 0.7. The model skill scores for the high runoff year of 1996 are lower, possibly because the turbulent mixing parameterization is inaccurate in simulating mixing under strong stratification (Li et al., 2005a).

Figure 6. A comparison of modeled (solid lines) and observed (dots) surface/bottom salinity at CBP monitoring stations for the normal runoff year of 1997: CB3.3c (a)/(b), CB4.4 (c)/(d), CB5.4 (e)/(f), CB8.1 (g)/(h), LE2.3 (i)/(j) and LE5.5 (k)/(l). Surface salinity is in the left column and bottom salinity in the right column.

(3) Current prediction

Another important metric for evaluating the model results is the subtidal residual velocity. We carried out a comparison between the observed and predicted velocity over a forty-day period during the spring of 1997 (Figure 7). Both observations and model records were passed through a 34-hour Lanczos filter to remove tidal fluctuations. There were a series of strong wind events, each lasting for 2 to 5 days. The Chesapeake Bay responded to this local longitudinal wind forcing, with amplitude and duration of currents matching those in the wind record reasonably well. Figure 7 shows that the predicted currents track the observed current, but there were small but visible phase differences. The model skill is 0.82 and 0.84 for the surface and bottom subtidal velocities, respectively. Although wind-driven currents dominate at the daily time scales, the monthly-averaged velocity shows a net seaward flow at 2.4 m depth but a net landward flow at 19.6 m depth (Li et al, 2005).

Figure 7. Comparison of the subtidal velocity between the ROMS model and CBOS observations at the mid-Bay station during the spring of 1997. (a) Time series of the wind speed used to drive the ROMS model. (b) Low-pass filtered wind speed. Low-pass filtered along-channel current speed (positive for the seaward direction) at 2.4 (c) and 19 (d) meter depth. The solid line represents CBOS observations and the dashed line the ROMS predictions.

(4) Biophysical model predictions

We used the coupled hydrodynamic-biogeochemical model to simulate the plankton dynamics during the normal runoff year of 1997. The fresh-water plume moves seaward, carrying nitrate down the Bay. By the end of May, however, much of the nitrate is depleted due to the uptake by phytoplankton production. A spring phytoplankton bloom develops around mid-May, covering a region between the upper and mid Bay. Ammonium concentration is much lower than nitrate concentration, confirming that the spring bloom is fuelled by the external nutrient drained from the watershed. In contrast, summer plankton production appears to be largely supported by regenerated nutrients. Warm water temperature stimulates bacteria production and accelerates the remineralization of organic particles into ammonium.

The Chesapeake Bay Remote Sensing Program (CBRSP) has used measurements of ocean color from light aircraft to estimate chlorophyll concentrations in Chesapeake Bay for the past thirteen years (Harding, 1994; Harding and Perry, 1997). Flights are conducted at a frequency of up to twice per week on the main stem of the Bay and monthly on two selected tributaries using ocean-color sensors similar to the Sea-viewing Wide Field-of-view Sensor (SeaWiFS) instrument that is currently in earth orbit. Figure 8b shows a snapshot of surface chlorophyll distribution obtained from the aircraft remote sensing at the end of July. Figure 8a shows the predicted chlorophyll distribution at the same time. The model appears to be doing a good job in reproducing the observed spatial variability in phytoplankton populations.

Figure 8. Comparison between the predicted (a) and observed (b) surface chlorophyll map at the end of July of 1997.

To examine whether the model captures the ecosystem productivity as a whole, we plot the annual time series of depth-integrated primary productivity and euphotic-layer chorlophyll in Figure 9. We examine not only the averages over the whole Bay but also the regional averages in the lower, mid and upper Bay. Although phytoplankton biomass reaches a peak during the spring bloom, the primary productivity reaches a broad peak during the summer. This phase lag between the spring peak of euphotic chl a and the summer peak of net primary productivity is a distinct characteristic of Chesapeake Bay planktonic ecosystem and has been observed consistently in previous observational studies (Malone, 1992). The predicted values are broadly consistent with Harding et al. (2002)’s observational estimates but model over-predicts productivity in upper bay, possibly as a result of neglecting phosphorus limitation there.

Figure 9. Annual time series of depth-integrated productivity (red) and euphotic-layer Chla (black) averaged for the whole Bay (a), lower Bay (b), mid-Bay (c) and upper bay (d).

Publications:

Li, M., L. Zhong, and W. C. Boicourt. 2005. Simulations of Chesapeake Bay estuary: Sensitivity to turbulence mixing parameterizations and comparison with observations, J. Geophys. Res., 110, C12004, doi:10.1029/2004JC002585.

Li, M., L. Zhong, W. C. Boicourt, S. Zhang and D.-L. Zhang. 2005. Hurricane-induced storm surges, currents and destratification in a semi-enclosed bay. Geophys. Res. Lett., In press.

Zhong, L. and M. Li. 2005. Tidal energy fluxes and dissipation in the Chesapeake Bay. Cont. Shelf Res., In review.