ABSTRACT Title of Dissertation / Thesis: MODELING THE PHYSICAL, OPTICAL AND BIOLOGICAL PROPERTIES OF CHESAPEAKE BAY Jiangtao Xu Dissertation / Thesis Directed By: Shenn-Yu Chao, professor, UMCES Raleigh R. Hood, associate professor, UMCES This thesis describes a relatively simple biogeochemical model that I developed and coupled with a three-dimensional circulation model of the Chesapeake Bay. To improve the performance of the physical model I attempted to assimilate high-resolution salinity data using a Newtonian relaxation scheme. In general, the simple assimilation scheme leads to visibly improved density structures in the Bay. However, the injection of high-resolution salinity data produces transient gravitational readjustment, which can have a significant impact on the biogeochemical properties and processes in the estuary. Therefore, this approach cannot be directly applied in biogeochemical modeling studies. Instead, I show that adjusting the salinity at open-ocean boundaries is also able to improve the density structure of the inner estuary. To obtain a relatively simple but effective way to model light attenuation variability in the coupled physical-biological model, I adopted a simple, non-spectral empirical approach. Surface water quality data and light measurements from the Chesapeake Bay Program were used to determine the absorption coefficients in a linear regression relationship. The resulting model between light attenuation coefficient (K d ) and water quality concentrations (chlorophyll, TSS and salinity as a proxy for CDOM) gives generally good estimates of K d in most parts of Chesapeake Bay. I also discuss the feasibility and caveats of using K d converted from Secchi depth in the empirical method. To develop the relatively simple biogeochemical model for Chesapeake Bay, I adopted a simple NPZD-type biological model and added in necessary additional components and simple parameterizations of the important processes for estuarine applications. The coupled model is then run under very different conditions: a dry year (1995) and a very wet year (1996). Observations of DIN, chlorophyll, total suspended solids (TSS), dissolved oxygen (DO), and light attenuation coefficient (K d ) obtained from Chesapeake Bay Program are used to validate the model. I demonstrate that this simple biological model is capable of reproducing the major features in nutrient, phytoplankton, DO, TSS and K d distributions in a complex ecosystem like Chesapeake Bay, and the model is robust enough to generate reasonable results under both wet and dry conditions. Sensitivity studies on selected parameters are also reported. MODELING THE PHYSICAL, OPTICAL AND BIOLOGICAL PROPERTIES OF CHESAPEAKE BAY By Jiangtao Xu Dissertation submitted to the Faculty of the Graduate School of the University of Maryland, College Park, in partial fulfillment of the requirements for the degree of Doctor of Philosophy 2005 Advisory Committee: Associate Professor Raleigh R. Hood, Chair/Advisor Professor Shenn-Yu Chao, Co-Advisor Professor William C. Boicourt Professor James A. Carton Professor Thomas R. Fisher Associate Professor Harry V. Wang ? Copyright by Jiangtao Xu 2005 Dedication To my parents. ii Acknowledgements I am deeply grateful to my advisors Drs. Shenn-Yu Chao and Raleigh Hood for their guidance, support, and patience. I could not have finished this project without their encouragement and help. I am also grateful to Drs. Bill Boicourt, Tom Fisher and Harry Wang for their valuable advice and comments to this research. Special thanks go to Dr. Harry Wang for his generous help in the initial tuning of the physical model and Dr. James Carton for agreeing to sit in my committee as Dean?s representative on a short notice. I would like to thank Cathy Lascara and Russ Burgett of Old Dominion University who kindly helped me with the IDL code for retrieving and visualizing the Scanfish data. This work has also benefited from many hours of discussion on statistics with my friend Yi Huang at Johns Hopkins University. In addition, I am thankful for the two years of fellowship and five months bridge fund I received from Horn Point Laboratory. Finally, I wish to express my sincere thanks to Horn Point community for the support, help and friendship. iii Table of Contents Dedication..................................................................................................................... ii Acknowledgements......................................................................................................iii Table of Contents......................................................................................................... iv List of Tables ............................................................................................................... vi List of Figures............................................................................................................. vii Chapter 1: General Introduction and Motivation of Current Research......................... 1 1. Introduction........................................................................................................... 1 2. Background........................................................................................................... 4 2.1 Physical perspective of the Bay ...................................................................... 4 2.2 Nutrient loading to the Bay............................................................................. 6 2.3 Factors affecting phytoplankton growth ......................................................... 8 2.4 Seasonal and inter-annual changes in distribution of phytoplankton biomass and production .................................................................................................... 10 2.5 The role of marshes and wetlands on nutrient loading ................................. 12 2.6 Sediment biogeochemistry............................................................................ 14 3. Motivation and approach of current research ..................................................... 15 3.1 Importance of correctly modeling the physical processes ............................ 15 3.2 Underwater light field................................................................................... 17 3.3 Complexity of the biological model ............................................................. 18 References............................................................................................................... 19 Chapter 2: Assimilating High-Resolution Salinity Data Into a Model of a Partially Mixed Estuary............................................................................................................. 32 Abstract................................................................................................................... 32 1. Introduction......................................................................................................... 33 2. Oceanographic Setting........................................................................................ 36 3. The hydrodynamic model ................................................................................... 38 4. Data assimilation scheme.................................................................................... 45 5. Gravitational readjustment.................................................................................. 47 6. Model improvement through data assimilation .................................................. 50 7. Discussion and conclusions ................................................................................ 56 Appendix A: A shorter assimilation time scale ...................................................... 58 Appendix B: Modification of salinity at open-ocean boundaries ........................... 59 Acknowledgements................................................................................................. 60 References............................................................................................................... 61 Figures..................................................................................................................... 64 Chapter 3: A Simple Empirical Optical Model for Simulating Light Attenuation Variability in a Partially Mixed Estuary ..................................................................... 81 Abstract................................................................................................................... 81 1. Introduction......................................................................................................... 82 2. Derivation of empirical light model.................................................................... 86 2.1 Empirical linear light model derived using direct light measurements ........ 86 2.2 Empirical light model derived using Secchi Depth ...................................... 90 iv 2.3 Dual empirical linear light models with direct light measurement............... 92 3. The Role of each component in light attenuation ............................................... 96 4. Summary and conclusions .................................................................................. 97 Acknowledgements............................................................................................... 100 References............................................................................................................. 100 Figures................................................................................................................... 106 Chapter 4: Modeling Biogeochemical Cycles in Chesapeake Bay with a Coupled Physical-Biological Model........................................................................................ 112 Abstract................................................................................................................. 112 1. Introduction....................................................................................................... 113 2. Model Description ............................................................................................ 119 2.1 The physical model..................................................................................... 119 2.2 The biogeochemical model......................................................................... 120 2.3 Boundary conditions for the biogeochemical model .................................. 123 2.4 Parameterization in the model .................................................................... 128 3. Results and Discussion ..................................................................................... 134 3. 1 Main run results ......................................................................................... 135 3.2 Lateral variations of chlorophyll distribution ............................................. 145 3.3 Sensitivity studies ....................................................................................... 149 4. Summary and Conclusions ............................................................................... 153 Acknowledgements............................................................................................... 157 References............................................................................................................. 157 Figures................................................................................................................... 167 Chapter 5: Summary and Conclusions...................................................................... 191 References............................................................................................................. 203 References................................................................................................................. 204 v List of Tables Table 2.1 The time period and corresponding values of salinity change made at the open-ocean boundaries. Table 3.1 The coefficient, P value and partial R 2 for each term in equation 6. Table 3.2 The coefficient and P value for each term in equation 12. Table 4.1 Model parameters. Table 4.2 Initial values for variables in the biogeochemical model vi List of Figures Figure 2.1. Bathymetry of the Chesapeake Bay and adjacent coastal area. Major tributaries are marked. Depth scales are in the unit of meters. Figure 2.2. Daily freshwater discharge in units of 10 3 cubic meters per second from Susquehanna, Potomac, James and Rappahannock Rivers in 1995. Figure 2.3. Scanfish transects and selected salinity stations in 1995. Salinity stations were maintained by Chesapeake Bay Program (CBP); selected stations identities mentioned in the text are given. Zonal dashed lines delineate upper, middle and lower reaches of the Bay. Figure 2.4. Boundary fitted grid of the Chesapeake Bay hydrodynamic model. Figure 2.5. Model-derived time series of surface salinity and observed salinity data at (a) station CB3.3C and (b) station CB5.1 in 1995. Figure 2.6. Model-derived and observed time series of water surface temperature at (a) Tolchester Beach and (b) Solomons Island in 1995. Figure 2.7. Vertical sections of salinity (psu) from the model (upper panels) and Scanfish measurements (lower panels) along (a) transect 13 and (b) transect 20 in July. Horizontal axes are time (in hours) in 1995. Figure 2.8. As in Fig. 2.7 except along (a) transect 2 and (b) transect 16 in October. Figure 2.9. Longitudinal-vertical sections of monthly averaged currents and salinity in (a) February and (b) October of 1995, before data are assimilated. The longitudinal section follows the center axis of deep channel southward to the vii mouth of Rappahannock River and thereafter extrapolates farther southward to the southern land boundary. Contour intervals for salinity are 2 psu. Figure 2.10. Longitudinal-vertical sections of instantaneous currents and salinity (psu) at: (a) 11:30 am on October 23, (b) 11:30 am on October 24 and (c) 11:30 pm on October 27 of 1995. Panels (a), (b) and (c) also correspond respectively to the times right before the beginning of data assimilation, at 8 hours after the beginning of active data assimilation, and at 10 hours after the completion of data assimilation. Contour intervals for salinity are 2 psu. Figure 2.11. The root-mean-square errors (in psu) between modeled and observed salinity values with and without data assimilation in the (a) upper, (b) middle and (c) lower reaches of the Chesapeake Bay. With reference to Fig. 2.3, the ensemble contains all data points collected in the upper reaches in (a), only data at station CB4.3C in the middle reaches in (b), and all data in the lower reaches in (c). Temporally, data are grouped by the month except in (b). Arrows indicate the approximate time of data insertion. Figure 2.12. Zonal sections of model-derived salinity in October midway between Scanfish transects 3 and 4 in Fig. 2.3. Top panel is shortly before data injection in the vicinity of this section (at day 295.72). The middle panel is shortly after data injection at day 298.72. The bottom panel is at day 301.72. Fast changes of salinity at three-day intervals indicate a quick loss of memory after data injection in the lower Bay. The view is northward. Figure 2.13. Comparisons among observed salinity profiles (solid), modeled salinity profiles without data assimilation (dashed), and modeled salinity profiles with viii data assimilation (dash-dotted) at selected stations. Panels (a) and (b) are at station CB3.3C in the upper Bay, while panels (c) and (d) are at stations CB5.1 in the mid-Bay, Times of comparison as indicated in (a) and (c) are right after data assimilation in July and times in (b) and (d) are about one and a half months after active data insertions in October. Figure 2.14. Biweekly averaged flow fields without data assimilation [left panels (a) and (c)] and corresponding changes caused by data assimilation [right panels (b) and (d)]. Time averaging is over the second and third weeks of November 1995. Panels (a) and (b) are surface features, while (c) and (d) are taken at 10m below mean water level. Zonal length scales and zonal velocities are stretched for clarity. Figure 2.15. Biweekly averaged section of density (? t ) anomalies induced by the data assimilation in October along the main axis of the deep channel. Time averaging is from day 295 to day 310.55 in panel (a), and from day 310.55 to day 326.1 in panel (b). The longitudinal section is the same as in Figs. 2.9 and 2.10. The density anomaly is obtained by subtracting the model result without data assimilation from that with data assimilation. Solid and dashed contours correspond to density deficit and density surplus, respectively. Contours intervals are ?? t =0.2. Figure 2.16. As in Fig. 2.11 except T is reduced threefold to 2 hours and K is increased threefold to (5 hr) -1 . ix Figure 2.17. Model-derived time series of bottom salinity and observed data at CB3.3C (a and c) and CB5.1 (b and d) before (a and b) and after (c and d) adjustment of open-ocean boundary conditions in 1995. Figure 3.1. Locations of Chesapeake Bay Monitoring Program stations. Figure 3.2. Vertical light attenuation coefficient (K d ) (m -1 ) obtained from the uniform linear regression relation (shown in the figure) versus K d derived from direct light measurements. Figure 3.3. K d (m -1 ) derived from direct light measurements versus Secchi depth (m) and a comparison of three different conversion relationships between K d and Secchi depth. Figure 3.4. Vertical light attenuation coefficient (K d ) (m -1 ) obtained from the uniform linear regression relation (shown in the figure) versus K d converted from measured Secchi depth. Figure 3.5. A down Bay comparison between K d obtained from linear regression relations 11.1 and 11.2 and from direct light measurements. Figure 3.6. The contributions of water plus dissolved organic matter (K w ?+ K s [SAL]), phytoplankton (K c [CHL]+K t [Phy]), and seston (K t [TSS-Phy]) to total light attenuation coefficient for all data points used in the regression model. Figure 4.1. A schematic diagram of the ecosystem model. Figure 4.2. Daily freshwater discharge in units of 10 3 cfs from major tributaries of Chesapeake Bay in 1995 (dashed line) and 1996 (solid line). Figure 4.3. Selected Chesapeake Bay Program monitoring stations. Zonal dashed lines delineate upper, middle and lower reaches of the bay. x Figure 4.4. The modeled and observed seasonally averaged surface DIN concentration (?M) at the mainstem stations in 1995 and 1996. Data points represent the mean and error bars the 95% confidence interval. Figure 4.5. As in Fig. 4.4 except for chlorophyll (?g/l). Figure 4.6. As in Fig. 4.4 except for light attenuation coefficient (K d ) (m -1 ). Figure 4.7. As in Fig. 4.4 except for total suspended solids (mg/l). Figure 4.8. The modeled and observed spatially averaged surface chlorophyll concentration (?g/l) for upper, mid, and lower bay in 1995 and 1996. Figure 4.9. The modeled and observed seasonally averaged vertically profiles of chlorophyll concentration (?g/l) at (a) CB3.3C, (b) CB5.1 and (c) CB6.3 in 1995. Figure 4.10 As in Fig. 4.9 except for DIN concentration (?M). Figure 4.11. As in Fig. 4.9 except for total suspended solids concentration (mg/l). Figure 4.12. As in Fig. 4.9 except for dissolved oxygen concentration ((?M). Figure 4.13. Planview of modeled averaged surface DIN (?M), chlorophyll (?g/l) and light attenuation coefficient (K d ) (m -1 ) in 2 nd and 3 rd quarters in 1995. Figure 4.14. As in Fig. 4.13 except for 1996. Figure 4.15. Planview of modeled averaged surface DIN (?M) and chlorophyll (?g/l) in the 2 nd quarter, 1995 with implementation of bivalve filtration at shoal regions. Figure 4.16. Seasonally averaged surface chlorophyll (?g/l) in 2nd and 3rd quarters 1995 from observations, modeled with and without implementation of bivalve xi filtration at shoal regions at two pairs of stations: (a) CB4.2C and EE2.1, (b) CB5.3C and EE3.4. xii Chapter 1: General Introduction and Motivation of Current Research 1. Introduction In recent decades, coupled physical-biological models have been widely applied to the marine environment to simulate both the physical and biogeochemical processes and study the interactions between them, especially the effect of physical factors on biological communities. The complexity of the physical models ranges from box (Li et al., 2000) and 1-D models (cf. Doney et al., 1996; Hood et al., 2001; Marra and Ho, 1993) to fully 3-D hydrodynamic models (cf. Lima and Doney, 2004; Skogen et al., 1995). The biological models range from simple NPZ (nutrient, phytoplankton, zooplankton) (cf. McClain et al., 1996) or NPZD (nutrient, phytoplankton, zooplankton, detritus) model (cf. Doney et al., 1996; Oschlies and Garcon, 1999; Hood et al., 2003) to multi-nutrient, multi-species and size-structured ecosystem models (Moore et al., 2002; Lima and Doney, 2004). When such models are applied to estuarine and coastal waters, they can provide a means of assessing the potential impacts of local management strategies and hence provide useful information to decision-makers. Compared to the open-ocean, the estuarine environment is much more variable and complex due to the confluence of both fresh and oceanic waters. Elevated terrestrial and anthropogenic nutrient inputs and the hydrographic characteristics of estuaries promote the retention of both materials and organisms (cf. Taft et al., 1978; Malone, 1992), which support high productivity in estuaries. The 1 high primary production in these systems leads to disproportionately large yields from higher trophic levels including fisheries (Houde and Rutherford, 1993). Estuaries also provide irreplaceable habitats for living marine resources and wildlife and support recreation, tourism and other industries. These properties make estuaries very valuable natural resources. However, similar to other estuarine systems (e.g., Lapointe and Clark, 1992; Pitkanen et al., 1993; Nagy et al., 2002), Chesapeake Bay, the largest estuary in the United States, has been suffering from degradation of water quality due to increased environmental stresses (cf. Carpenter et al., 1969; Malone 1992). Eutrophication in Chesapeake Bay has caused serious economic, aesthetic and ecological problems: harmful algae blooms (Bowers et al., 2000), loss of submerged aquatic vegetation (SAV) (Orth and Moore, 1983), and hypoxia and anoxia at deep waters in summer (Cooper and Brush, 1991; Kemp et al., 1992). Increased loads of suspended solids from the surrounding land directly reduces water clarity and deposits of this material on the bottom can have detrimental impacts on benthic organisms and production (Airoldi, 2003; Miller et al., 2002). Efforts have been made to reduce the N and P inputs from point and non-point sources and land-based sediment runoff with the goal of restoring the Bay to conditions observed in the early 1950s (Chesapeake Bay Agreement 1983, 1987, 2000). Numerical models have been used as a key analytic tool to provide guidelines in setting goals of nutrient and sediment reduction to achieve water quality standards. In an effort to model the complexity of the ecosystem the Chesapeake Bay Program (CBP) has developed a whole package of models, which includes an airshed model (Regional Acid Deposition Model (RADM)) (Chang et al., 1990; Dennis, 2 1996), a watershed model (Hydrological Simulation Program-Fortran (HSPF)) (Bicknell et al., 1996; Greene and Linker, 1998), a hydrodynamic model (WES- CH3D) (Johnson et al., 1991; Hood et al., 1999; Sheng, 1986) and a water quality model (CE-QUAL-ICM) (Cerco and Cole, 1994; Cerco and Noel, 2004) coupled with a sediment (DiToro and Fitzpatrick, 1993) and living resources (including SAV and benthos) model (Madden and Kemp, 1996; Wetzel and Neckles, 1986). Even though this modeling system is extremely complicated, many processes, which could be important in the nutrient budgets, are not modeled or fully accounted for. These include, for example, the processes in marshes and wetlands, nutrient inputs through ground water and atmospheric dry deposition, etc. Among these, marshes and wetlands could be a key missing component. Marshes have high tolerance for increased nutrient loading and function as a buffer zone for intercepting nutrients before they reach the adjacent estuarine and coastal waters (Nixon, 1980). In addition, the water quality model in this package itself has 24 state variables including two physical variables (temperature and salinity), multiple algal groups, two zooplankton groups, and multiple forms of nitrogen, phosphorus, carbon and silica. The numerous processes and huge number of parameters in this model make it difficult to identify the essential components and diagnose the potential problems. In contrast to the extremely complicated, all-inclusive CBP modeling package, in this study we use the same hydrodynamic model but incorporate a relatively simple ecosystem model with the goal of capturing the lower order variability in the system. Using the simple ecosystem model we attempt to address the following questions: 1) What are the essential components in modeling the biological 3 processes in Chesapeake Bay? 2) Can a relatively simple model capture the 1 st order variability in nutrient cycles, oxygen concentration and phytoplankton biomass? 2. Background 2.1 Physical perspective of the Bay Chesapeake Bay, the largest estuary in the United States, stretches for about 320km from the mouth of the Susquehanna River to its seaward end at Cape Charles and Cape Henry. As a partially mixed estuary, the circulation is driven by the freshwater inflow, tides and wind forcing. There are 50 major rivers which discharge into the Bay. The total freshwater input to the system average about 2280 m 3 /s. Nearly half of the total (48.2%) is supplied by the Susquehanna River. Freshwater sources along the western shore account for 43.6% of the total average freshwater inflow. The remaining 8.2% of the total average freshwater inflow comes from the rivers at the eastern shore. All of the major freshwater sources entering the Bay exhibit considerable month-to-month and year-to-year variability. Typical of mid-latitude rivers, they have high discharge in spring followed by low to moderate flow throughout the rest of the year, with lowest flow in late summer and early fall. However, the year-to-year variation in the monthly averaged flow can be relatively large. Over long time scales the range of monthly averaged flows considerably exceeds the median flow for each month (Schubel and Pritchard, 1987). Tides in Chesapeake Bay have a predominant lunar semidiurnal component with a significant spring and neap tide cycle. Tidal waves progress from the Bay?s 4 entrance to the mouth of Susquehanna River in approximately 14 hours (Hicks, 1964), so that the Bay is able to contain a complete semidiurnal tidal wave at all times. The bay is also wide enough to show the earth?s rotational effects: the mean tidal range is significantly larger on the eastern shore. Close to the head of the Bay, the friction and reflection effects on the propagation of tidal waves become significant (Schubel and Pritchard, 1987). Averaging over a long period, the residual circulation in an estuary is predominantly density driven. In the Bay, the spatial variation in density results primarily from spatial variation in salinity. However, the vertical variation in temperature can account for as much as 20% of the vertical variation in density in the middle and upper reaches of the mainstem during late spring and early summer. The density distribution in the Bay depends primarily on the amount and timing of the freshwater inflow. When averaged over a sufficient time interval, the residual current field of Chesapeake Bay exhibits the classical two-layer estuarine circulation pattern of a partially mixed estuary. However, certain parts of the mainstem and some tributaries do, at times, approach conditions of a vertically homogeneous estuary. Wind-induced circulation in the Bay can be locally driven or remotely driven. The upper layer of the estuary responds directly to the north-south component of the local wind forcing, and the bottom layer responds in the opposite direction with a time lag, presumably in response to the downwind setup of the sea surface (Wang, 1979). The circulation also reflects the Ekman transport effect produced by the cross- bay component of the winds in the middle and lower reaches of the Bay (Elliott et al., 5 1978) or the north-south wind in the adjacent coastal ocean which results in sea level fluctuations (Elliott and Wang, 1978; Wang and Elliott, 1978). 2.2 Nutrient loading to the Bay The nutrient enrichment of lakes, estuaries and coastal systems from anthropogenic activities has been a concern for the last several decades (e.g.: Nixon and Pilson 1983, Turner and Rabalais 1994). Chesapeake Bay has shown many symptoms of eutrophication (Heinle et al. 1980) and high nutrient concentrations in the Bay cause rapid growth of phytoplankton and algae. In shallow areas, the excess algae block the sunlight which is important for the growth of submerged aquatic grasses. This degrades the habitat and causes the eventual loss of these grass beds. In deeper areas, the decomposition of dead algae uses up available oxygen in the water, causing prolonged anoxic or hypoxic condition in warmer summer months. Excessive nutrient loading is now a serious problem in managing the water quality of the Bay. Nutrients enter Chesapeake Bay from rivers, atmospheric deposition, point sources, ground water discharge and diffusive sources along the shoreline. Of the nutrients from the nine major tributaries of Chesapeake Bay, three rivers have the highest flow: the Susquehanna, the Potomac, and the James Rivers, and they contribute the largest nutrient loads to the tidal part of the Bay. The Susquehanna River is the largest river entering the Bay and contributes about 60% of the total streamflow, 62% of the total nitrogen load, and 34% of the total phosphorus load to the nontidal part of the Chesapeake Bay. The Potomac River is the second largest river draining to the Bay and the second largest source of nitrogen and phosphorus. It 6 contributes about 20% of the total streamflow, 28% of the total nitrogen load and 33% of total phosphorus load. Of the nine major tributaries, the James River contributes about 12% of the stream flow, 5% of the total nitrogen load and 20% of total phosphorus load. Collectively, the Susquehanna, Potomac and James Rivers contribute about 95% of the nitrogen load and 87% of the phosphorus load from the nine major tributaries draining to the Bay (Belval and Sprague, 2000). Atmospheric deposition is also a significant source of nutrient input for the Bay. Phosphorus has no significant atmospheric sources (ref) and anthropogenic deposition of atmospheric nitrogen is largely produced from the burning of fossil fuels. Studies have shown that atmospheric deposition of nitrogen is an important element of total nitrogen load to the Chesapeake Bay (Meyers et al., 2000, Fisher and Oppenheimer, 1991). Meyers et al. (2000) estimated that the total nitrogen deposition to the watersheds of Chesapeake Bay is about 13 kg N ha -1 yr -1 , of which 35% originated from dry deposition and the total direct nitrogen deposition to the water surface is about 11x10 6 kg yr -1 . Including riverine nitrogen loadings attributable to atmospheric sources, atmospheric nitrate deposition makes up as much as 25% of the anthropogenic nitrogen loading to the Bay, and atmospheric ammonium deposition contributes another 14% of the total (Fisher and Oppenheimer, 1991). However, only about 10 percent of the Bay's nitrogen load is the result of airborne nitrogen, deposited directly on the water surface of the mainstem and the tidal portions of its tributaries. Point sources also make important local contribution to nutrient inputs. Discharges from municipal wastewater treatment plant contribute the majority of 7 point source loadings. Nitrogen and phosphorus loads from point sources were reduced 33% and 56% between 1985 and 2001, respectively (The Chesapeake Bay Program). Chesapeake Bay Program data from 1995 and 1996 show that the total ?end-of-pipe? discharge, including industrial sources of nitrogen and phosphorus to all receiving water bodies are about 38 x 10 6 kg yr -1 and 2.7 x 10 6 kg yr -1 , respectively. Point sources are the second largest contributor to the total nutrient loadings. The diffusive nutrient sources entering the Bay from the coastal plain portion of the basin are not accounted for in the river inputs. Due to the proximity to the water body, with potentially less retention, nutrient inputs from land runoff from the adjacent watershed may be significant (e.g. Lee et al., 2001). The below fall line diffuse sources contribute about 5% and 10%, respectively, of total nitrogen and phosphorus input to the Maryland mainstem bay (Boynton et al., 1995). The nutrient loading through ground water discharge could be an additional large input (Staver and Brinsfield, 1996). However, only a little is known about the amount of nutrient from ground water entering Chesapeake Bay. 2.3 Factors affecting phytoplankton growth Nutrient concentrations in Chesapeake Bay have been undesirably high in recent years. It is widely accepted that N is the essential limiting nutrient for phytoplankton biomass on a baywide scale (Malone et al. 1996, Harding 1994). However, the relative nutrient availability exhibits large seasonal and spatial variations, leading to possible seasonal shifts in the limiting nutrient. This is 8 profoundly important when considering the control of the algae growth by reducing the nutrients entering the Bay. Rivers supply large amounts of N relative to P due to the greater percentage of N (particularly NO 3 - ) transported via nonpoint runoff. Atmospheric deposition adds solely nitrogen to the Bay. During late winter and spring, Chesapeake Bay receives maximum fresh water runoff, and nutrient loadings from rivers dominate all sources. However, in summer, the freshwater discharges decrease, and regenerated nutrients become more important. Fisher et al. (1992, 1999) found that along the main axis of the Bay DIN/PO 4 in the surface water was typically greater than 100 in spring, indicating potential P- limitation of the phytoplankton growth rate and the accumulation of algae biomass. In summer, DIN/PO 4 was generally <10:1, indicating potential N-limitation. Their study, which includes assessments of nutrient turn over time, alkaline phosphatase activity and nutrient enrichment bioassays, strongly suggests that the limiting nutrient shifts from P (possibly Si) in spring to N in summer. This study also shows that the likelihood of P-limitation is greatest in the upper Bay where the impact of high N:P ratio non-point source loading is greatest. Similar seasonal changes in the limiting nutrients were also found in the Patuxent and York subestuaries of Chesapeake Bay (D?Elia et al., 1986, Webb 1988). The relative availability of P is also closely linked to sedimentation processes. During winter-spring, the high-flow season, the total suspended solid input is also high. When suspended solids settle, they remove phosphorus to the sediments by absorption (Stumm and Morgan 1981, Fox et al. 1985). In warm seasons the DIP 9 fluxes across the sediment-water interface (to the water column) are highest (Boynton et al., 1980). In summer time, the anoxic conditions in deep water appear to result in substantial release of PO 4 from sediments (Boynton and Kemp, 1985). Besides nutrient limitation, the phytoplankton distribution in Chesapeake Bay also shows strong light limitation in nutrient-rich, high turbidity zones. This will be emphasized in the next section. 2.4 Seasonal and inter-annual changes in distribution of phytoplankton biomass and production Fresh water inflow, winds and tides determine the circulation features, nutrient and sediment loading of an estuary. These factors in turn determine the distribution of phytoplankton in the system. The interplay of these properties produces a spatially and temporally heterogeneous distribution of phytoplankton growth and the accumulation of biomass. For example, a persistent patch of high phytoplankton concentration between the mouths of the Potomac and Rappahannock Rivers is probably induced by the sill off the Rappahannock River. Annual phytoplankton production in Chesapeake Bay appears to be more sensitive to nitrogen than to phosphorus loading. The Susquehanna River is the major external source of nitrogen to Chesapeake Bay, and the fresh water discharge has a spring maximum and a fall minimum. As a result, about 50-60% of the annual nitrogen input to the upper Bay occurs during the spring freshet (Schubel and Pritchard 1986). Most of this nitrogen supply is assimilated downstream of the turbidity maximum in the mesohaline region of the Bay where phytoplankton productivity and chlorophyll a concentration are highest (Harding et al. 1986, Fisher 10 et al. 1988). In the turbid zone light limits nutrient uptake and phytoplankton productivity even though nutrients are abundant (Harding et al. 1986). Consequently, the timing and magnitude of the spring bloom in Chesapeake Bay largely depend on the timing and magnitude of the spring freshet, which controls the delivery of the nutrients and the estuarine circulation. In a high flow year, the sediment loading is high and the water at the upper reaches of the Bay is more turbid and has lower salinity The turbidity maximum locates further downstream due to the higher flow. Consequently, the high production zone in a high flow year exhibits a greater downstream displacement. The annual cycle of riverine nutrient input is in phase with phytoplankton biomass in the mesohaline reach of the Bay, but out of phase with phytoplankton productivity in this region (Malone et al. 1988). The accumulation of phytoplankton biomass is correlated with the external nutrient input while the phytoplankton productivity is correlated with light and temperature. Therefore, phytoplankton biomass peaks in spring, but phytoplankton productivity peaks in summer. These two peaks are coupled via the sedimentation of phytoplankton biomass during spring and subsequent recycling of regenerated nitrogen into the euphotic zone during summer. The fluxes of regenerated nutrients from bottom water to the euphotic zone in summer are mainly controlled by the vertical stratification. The water column is generally stratified in summer but the vertical density gradients can be weakened by tidal mixing in shallower regions. Also, occasional large and often rapid lateral oscillations of the pycnocline can change the vertical density structure and lateral transport of nutrients which may be responsible for the higher phytoplankton 11 production over the flanks of the main channel relative to production in the channel itself (Malone et al. 1986). Furthermore, sporadic mixing events in summer can break down the stratification in the mainstem and bring the regenerated nutrients from deep waters to euphotic zone where the nutrients support high phytoplankton production. The species composition of phytoplankton in the Bay also shifts seasonally. The classic view is that in spring a diatom bloom accounts for the annual biomass peak, but in summer flagellates and dinoflagellates make up most of the phytoplankton population (Malone et al., 1988; 1991). Accordingly, the food web of the system is dominated by diatom-mesozooplankton in spring, but microbially dominated in summer. 2.5 The role of marshes and wetlands on nutrient loading The amount of nutrients entering an estuarine or coastal system determines, to a large extent, the biogeochemical activities in the system. Consequently, great care must be taken concerning nutrient loading when modeling such systems. However, processes such as ground water flow, stream-edge land runoff, nutrient retention/release in marshes and wetlands and atmospheric dry deposition are still very difficult to fully represent in numerical models. Among these uncertain factors, marshes and wetlands may be the most important factor for nutrient sinks/sources due to their special role in nutrient budgets. Intertidal marshes are generally described as net exporters of organic material. Export occurs in the form of detritus and through the activities/feeding and migration of animals (Teal 1962, Deegan and Garritt 1997) and the extent of export or outwelling from tidal marshes is related to the level of productivity, marsh coverage, 12 tidal amplitude and the geomorphology of the estuarine landscape. The occurrence is often intermittent and largest during rainstorms and high spring tides (Odum, 2000). However, marshes have been described both as sinks and sources of inorganic nutrients (Nixon 1980). The function of marshes to intercept land-derived nutrient and hence act as a buffer zone between land and the adjacent estuarine and coastal waters has been of great interest in light of increasing nutrient loading to these systems (Valiela et al, 1976). Nitrogen entering the marshes is removed from the biologically active systems primarily either by being trapped in refractory organic matter or through loss to the atmosphere as N 2 by denitrification. In oxic water, phosphorus is generally found as insoluble salts and is transported to marshes attached to particles. Therefore, phosphorus is trapped and buried in marsh sediments in both organic and inorganic forms. Marshes at the upper Rhode River trap 700 moles d -1 N and 34 moles d -1 P in sediments, respectively (Jordan et al., 1983, 1991). A recent study of Patuxent River tidal freshwater marshes shows that the marshes retain 35% of the nitrogen and 81% of phosphorus inputs to the tidal fresh portion of the river and remove approximately 10% of nitrogen by denitrification (Merrill, 1999). Teal and Howes (2000) noted that the nitrogen removal through denitrification is more effective on inputs through ground water than through surface/tidal waters. Valiela et al. (2000), summarizing studies on 19 salt marshes, found that mature marshes export more NH 4 and NO 3 than young marshes. Despite these studies, there is still not enough information on nutrient exchange among linked watershed, estuary, and marsh systems to clearly understand the different roles of marshes. 13 The nutrient uptake or release from marshes and wetlands in temperate climates usually exhibits a seasonal pattern or at least differences in retention rate. Jordan et al.(1991) showed that dissolved PO 4 3- production peaked in summer when watershed discharge was lowest, but NO 3 - consumption peaked in spring when watershed discharge was highest in the upper Rhode River. Ammonium release in summer due to the mineralization of organic matter and ammonification of nitrate taken up by the sediments under anoxic conditions was observed in Patuxent River tidal marshes (Merrill, 1999). Chesapeake Bay tidal marshes cover approximately 1.7 million acres. Further studies and new techniques are needed to assess the role of these marshes and other sinks/sources and processes in nutrient balances. 2.6 Sediment biogeochemistry Sediment processes are important for nutrient cycling in estuaries. Organic nitrogen and phosphorus are remineralized in sediments. The inorganic nitrogen and phosphate in interstitial waters can then be released to the water column and become available to the plankton community. The regeneration of both ammonium and phosphate are temperature dependent with maximum fluxes from sediment in summer (Boynton et al. 1980). Sediment regeneration of phosphate and ammonium can provide an average of 28% of phytoplankton P requirement (Fisher et al. 1982) and 13 to 40% of phytoplankton N requirement with a higher percentage in the summer period (Boynton and Kemp, 1985). In the sediments some of the ammonium released from the organic matter is nitrified to nitrate, and denitrification converts a substantial proportion of the nitrogen to N 2 in hypoxic and anoxic strata (Kemp et al., 1990). 14 Hence the coupled process of nitrification-denitrification represents a pathway for nitrogen loss. Sediment nitrification rates are generally regulated by availabilities of O 2 and NH 4 + . The regeneration of NH 4 + has the highest rate in summer; however, the O 2 penetration into sediments declines in summer due to increased temperature and organic inputs. In the mesohaline region of Chesapeake Bay, a relatively high rate of nitrification and denitrification in spring and fall and virtual elimination of both processes in summer have been observed (Kemp et al. 1990). The sediment redox potential also affects phosphate precipitation and dissolution. In anaerobic sediments bacteria reduce ferric iron (Fe 3+ ) to ferrous iron (Fe 2+ ) in the presence of H2S. Ferrous iron is much less effective in adsorbing phosphate than ferric iron (Krom and Berner, 1980), which makes dissolved phosphate available under anaerobic conditions. Dissolved phosphate may leave anoxic sediments, but some of the phosphate may reprecipitate as FePO 4 and/or sorb to oxyhydroxides at the oxic-anoxic interface (Krom and Berner 1980). Boynton and Kemp (1985) suggested that the relatively low observed sediment flux of DIP in some Chesapeake Bay sediments is due to the presence of O 2 in the overlying water. 3. Motivation and approach of current research 3.1 Importance of correctly modeling the physical processes In aquatic systems, biological and chemical processes are, to a large degree, controlled by physical processes. For example, rates of phytoplankton growth, zooplankton grazing, and organic matter remineralization are temperature dependent, and species composition of the biological community are affected by salinity and 15 temperature. Recent open ocean modeling studies have emphasized the importance of improving the representation of physical processes and variability in order to improve the performance of biogeochemical models (Oschlies and Garcon, 1999; Hood et al., 2003; Friedrichs et al., 2004). In Chesapeake Bay, river flow during the spring freshet, which determines nutrient and sediment input, largely dictates the timing, magnitude and location of the spring bloom (Malone et al., 1988; Fisher et al., 1988), which consequently determines the available nitrogen for recycling in summer months (Malone et al., 1988). In summer, stratification and mixing events control the amount of nutrient delivered to the euphotic zone and affect the extent of hypoxia and anoxia conditions in the Bay. Due to the paramount importance of reproducing the correct physical conditions in a coupled physical-biological model for modeling biogeochemical variability, my first task in this thesis (Chapter 2) was to validate the physical model for the studied period: 1995 and 1996. These two years were chosen because of their very contrasting river flow conditions, i.e., 1995 was a low flow year and 1996 was a very high flow year. Because computer resources are limited, one inevitable problem faced in all numerical modeling studies is resolution. Even with today?s powerful computers, it is still impossible to resolve all processes at all relevant scales in a model. Moreover, there are still some processes, such as turbulence, whose mechanisms are not fully understood. Data assimilation has emerged as a powerful tool to improve model performance (e.g., Ezer amd Mellor, 1994; Forbes and Brown, 1996; Wu et al., 1999) 16 and/or infer uncertain parameters (Bogden et al., 1996; Smedstad and Obrien, 1991; Ullman and Wilson, 1998). In Chapter 2 I also made an attempt to assimilate the high resolution salinity data obtained by Scanfish to improve the density structure, and hence the velocity structure, of the 3-D hydrodynamic model. It is noted, however, that instabilities generated by the data assimilation can have a significant impact on the circulation and mixing in the estuary, and therefore may interfere with biogeochemical modeling studies. 3.2 Underwater light field Light is essential for photosynthetic plants and algae. However, due to the rapid attenuation in water, light is often a limiting factor in primary production in the aquatic environment (e.g. Fisher et al., 1999). The degree of light attenuation also varies tremendously in aquatic systems due to the variable presence of ?chromophoric? organic matter, such as phytoplankton, DOM and detritus. Therefore, reproducing the correct underwater light field is a key problem in modeling the biogeochemical processes in these systems. Because the incoming photosynthetically active radiation (PAR) at the air- water interface can be measured or calculated quite accurately (e.g. Fisher et al., 2003), the main issue in calculating the underwater light field is to have correct estimate of the vertical light attenuation coefficient (K d ). For monochromatic light the vertical light attenuation can be decomposed as a set of partial attenuation coefficients, each characterizing absorption and scattering by a different waterborne material. Strictly speaking, a complete spectrum of K d (K d (?)) is needed to obtain the average K d for the whole photosynthetic waveband and for each narrow band it is 17 necessary to know the wavelength-specific absorption and scattering coefficients for each waterborne material. Spectral bio-optical models have been developed and applied to different kinds of water bodies (e.g.: Arrigo and Sullivan 1994; Gallegos et al. 1990; Platt and Sathyendranath1988). However, due to the optical complexity of estuarine waters, and because our goal is to keep our biological model as simple as possible, we developed a simple, empirical, non-spectral bio-optical model for the Chesapeake Bay for estimating K d variability. This model is described and validated in Chapter 3 of this thesis. 3.3 Complexity of the biological model As we discussed above, the water quality model developed by the Chesapeake Bay Program for management purpose is extremely complex. However, a recent biogeochemical model intercomparison study has shown that increasing model complexity may not lead to increased skill or predictive ability (Friedrichs et al., 2004). It is not clear, however, whether or not these results, which were derived from a study of open-ocean models, is applicable in a complex system like Chesapeake Bay. We therefore set out to develop a simple biogeochemical model for Chesapeake Bay that includes only the essential components that are necessary for modeling nutrient cycling, oxygen and phytoplankton biomass variability. We also developed and incorporated into this model a suite of simple paramaterizations that account for some of the key sources of higher order variability, such as phosphorus limitation, temperature effects and seasonal changes in ecosystem structure. To bypass the effect of marshes and wetlands we use a simple nudging scheme to ?push? the model towards observations at the head of tributaries. We then couple this model with our 18 3-D hydrodynamic model and run them for both 1995 and 1996 to test whether or not this simple configuration is robust enough to reproduce the tremendous seasonal and interannual variability in Chesapeake Bay. We conclude that it is, although with some significant caveats. An overview of the dissertation is as follows. The effort of improving the 3-D hydrodynamic model through data assimilation is described and discussed in chapter 2. Chapter 3 deals with the development and validation of a simple empirical, non- spectral light model for Chesapeake Bay, and examines the issue of calculating K d using direct light measurement versus using K d derived from Secchi depth (SD). Chapter 4 provides a detailed description of the biogeochemical model we developed and we validate it by comparing the model results with Bay Program monitoring data. Sensitivity studies using selected parameters are also reported in chapter 4. Finally, a general summary and conclusions are given in chapter 5. References Arrigo, K. R. and C. W. Sullivan, 1994. A high resolution bio-optical model of microalgal growth: tests using sea-ice algal community time-series data. Limnology and Oceanography 39 (3): 609-631. Belval, D. L. and L. A. Sprague, 2000. Monitoring nutrients in the major rivers draining to Chesapeake Bay, USGS water resources investigations report 99- 4238. 19 Bicknell, B., J. Imhoff, J. Kittle, A. Donigian Jr., R. Johanson and T. Barnwell. 1996. Hydrologic Simulation Program-Fortran user?s manual for release 11. U.S. Environmental Protection Agency Environmental Research Laboratory, Athens, GA. Bogden, P. S., P. Malanotte-Rizzoli, and R. Signell, 1996. Open-ocean boundary conditions from interior data: Local and remote forcing of Massachusetts Bay, Journal of Geophysical Research 101: 6487-6500. Boynton, W. R., J. H. Garber, R. Summers and W. M. Kemp. 1995. Inputs, transformations, and transport of nitrogen and phosphorus in Chesapeake Bay and selected tributaries. Estuaries 18 (1B): 285-314. Boynton, W. R., W. M. Kemp, 1985. Nutrient regeneration and oxygen consumption by sediments along an estuarine salinity gradient. Marine Ecology Progress Series 23: 45-55. Boynton, W. R., W. M. Kemp and C. G. Osborne, 1980. Nutrient fluxes across the sediment-water interface in the turbid zone of a coastal plain estuary. Estuarine Perspectives, Academic Press. Bowers H. A., E. Tengs, H. B. Glasgow, J. M. Burkholder, P. A. Rublee and D. W. Oldach, 2000. Development of real-time PCR assays for rapid detection of Pfiesteria piscicida and related dinoflagellates. Applied and Environmental Microbiology 66 (11): 4641-4648. Carpenter J. H., D. W. Pritchard and R. C. Whaley, 1969. Observations of eutrophication and nutrient cycles in some coastal plain estuaries. In: Eutrophication: causes, consequences, correctives-proceedings of a 20 symposium. US national academy science publ 1700, Washington DC, pp210- 221. Chang, J., P. Middleton, W. Stockwell, C. Walcek, J. Pleim, H. Lansford, S. Madronich, F. Binkowski, N. Seaman and D. Stauffer. 1990. The regional acid deposition model and engineering model, NAPAP SOS/T report 4. In National Acid Precipitation Assessment Program: State of science of technology, Vol. 1, National Acid Precipitation Assessment Program, Washington, D.C. Cooper S. R. and G. S. Brush, 1991. Long-term history of Chesapeake Bay anoxia. Science 254: 992-996. Deegan, L.A. and R.H. Garritt, 1997. Evidence for spatial varability in estuarine food webs, Marine Ecology Progress Series 147: 31-47. D?Elia, C. F., D. M. Nelson, W. R. Boynton, 1986. Nutrient enrichment studies in a coastal plain estuary: phytoplankton growth in large-scale, continuous cultures. Canadian Journal of Fishery and Aquatic Science 43: 397-406. Dennis, R. 1996. Using the Regional Acid Deposition Model to determine the nitrogen deposition airshed of the Chesapeake Bay watershed. In Atmospheric Deposition to the Great Lakes and Coastal Waters. Ed.: Joel Baker. Society of Environmental Toxicology and Chemistry. Di Toro, D. and J. Fitzpatrick. 1993. Chesapeake Bay sediment flux model. U.S. Army Corps of Engineers Waterways Experiment Station, Vicksburg, MS. Contract report EL-93-2. 21 Doney, S. C., D. M. Glover and R. B. Najjar, 1996. A new coupled, one-dimensional biological-physical model for the upper ocean: Applications to the JGOFS Bermuda Atlantic time-series study (BATS) site. Deep-Sea Research II 43 (2- 3): 591-624. Elliott, A.J. and D.P. Wang, 1978. The effect of meteorological forcing on the Chesapeake Bay: the coupling between an estuarine system and its adjacent coastal waters, In Hydrodynamics of Estuaries and Fjords, edited by NCN. Nihoul, Elsevier, Amsterdam, The Netherlands. Ezer, T., and G, K, Mellor, 1994. Continuous assimilation of Geosat altimeter data into a three-dimensional primitive equation Gulf Stream model. Journal of Physical Oceanography 24: 832-847. Forbes, C., and O. Brown, 1996. Assimilation of sea surface height data into an isopycnic ocean model. Journal of Physical Oceanography 26: 1189-1213. Fisher, D.C. and M. Oppenheimer, 1991. Atmospheric nitrogen deposition and the Chesapeake Bay estuary, AMBIO, vol. 20, No.3-4. Fisher, T.R., P.R. Carlson and R.T. Barker, 1982. Sediment nutrient regeneration in three North Carolina estuaries. Estuarine, Coastal and Shelf Science 14: 101- 116. Fisher, T. R., A. B. Gustafson, G. M. Radcliffe, K. L. Sundberg and J. C. Stevenson, 2003. A long-term record of photosynthetically available radiation (PAR) and total solar energy at 38.60 degrees N, 78.2 degrees W. Estuaries 26 (6): 1450- 1460. 22 Fisher, T. R., A. B. Gustafson, K. Sellner, R. Lacouture, L. W. Haas, R. L. Wetzel, R. Magnien, D. Everitt, B. Michaels and R. Karrh, 1999. Spatial and temproral variation of resource limitation in Chesapeake Bay. Marine Biology 133 (4): 763-778. Fisher, T. R., L. W. Harding, D. W. Stanely, L.G. Ward, 1988. Phytoplankton, nutrients, and turbidity in the Chesapeake, Delaware, and Hudson estuaries. Estuarine, Coastal and Shelf Science 27 (1): 61-93. Fisher, T. R., E. R. Peele, J. W. Ammernan, L. W. Harding, 1992. Nutrient limitation of phytoplankton in Chesapeake Bay. Marine Ecology Progress Series 82: 51- 63. Fox, L.E., S.L. Sager and S.C. Wofsy, 1985. Factors controlling the concentrations of soluble phosphorus in the Mississippi estuary. Limnology and Oceanography 30: 826-832. Friedrichs, M. A. M., R. R. Hood and J. D. Wiggert, submitted. Ecosystem model complexity versus physical forcing: quantification of their relative impact with assimilated Arabian Sea data. Deep Sea Research II. Gallegos, C. L., D. L. Correll, and J. W. Pierce. 1990. Modeling spectral diffuse attenuation, absorption, and scattering coefficients in a turbid estuary. Limnology and Oceanography 35: 1486-1502. Greene, K. and L. Linker. 1998. Chesapeake Bay Watershed Model application and calculation of nutrient and sediment loadings-phase IV Chesapeake Bay Watershed Model-appendix a: model hydrology calibration results. EPA 903- 23 R-98-004, CBP/TRS 196/98, Chesapeake Bay Program Office, Annapolis, MD. Harding, L. W., 1994. Long-term trends in the distribution of phytoplankton in Chesapeake Bay: role of light, nutrients and streamflow. Marine Ecology Progress Series 104: 267-291. Harding, L. W., B. W. Meeson, T. R. Fisher, 1986. Phytoplankton production in two east coast estuaries: photosynthesis-light functions and patterns of carbon assimilation in Chesapeake and Delaware Bays. Estuarine, Coastal and Shelf Science 23: 773-806. Heinle, D.R., C.F. D?Elia, J.L. Taft, J.S. Wilson, M. Cole-Jones, A.B. Caplins and L.E. Cronin, 1980. Historical review of water quality and climate data from Chesapeake Bay with emphasis on effects of enrichment, USEPA Chesapeake Bay Program final report, Grant No. R806189010, Chesapeake Research Consortium, Inc. Publ. No. 84, Annapolis, MD. Rep. No. TR-002E. UMCEES 80-15CBL. Hicks, S.D., 1964. Tidal wave characteristics of Chesapeake Bay. Chesapeake science 5, 103-113. Hood, R.R., N.R. Bates, D.G. Capone and D.B. Olson, 2001. Modeling the effect of nitrogen fixation on carbon and nitrogen fluxes at BATS. Deep-sea research II 48: 1609-1648. Hood, R. R., K. E. Kohler, J. P. McCreary, S. L. Smith. 2003. A four-dimensional validation of a coupled physical-biological model of the Arabian Sea. Deep- Sea Research II 50 (22-26): 2917-2945. 24 Hood, R. R., H. V. Wang, J. E. Purcell, E. D. Houde and L. W. Harding Jr.. 1999. Modeling particles and pelagic organisms in Chesapeake Bay: Convergent features control plankton distributions. Journal of Geophysical Research 104: 1223-1243. Houde, E. D. and E. S. Rutherford. 1993. Recent trends in estuarine fisheries ? predictions of fish production and yield. Estuaries 16 (2): 161-176. Johnson, B. H., K. W. Kim, R. H. Heath, H. L. Butler and B. B. Hsieh 1991 Development and verification of a three-dimensional numerical hydrodynamic, salinity, and temperature model of the Chesapeake Bay, Technical Report HL-91-7, US Army Engineer Waterways Experiment Station, Vicksburg, MS. Jordan, T. E, D.L. Correll, J. Miklas and D.E. Weller, 1991. Nutrients and chlorophyll at the interface of a watershed and an estuary. Limnology and Oceanography 36(2): 251-267. Jordan, T. E, D.L. Correll and D.F. Whigham, 1983. Nutrient flux in the Rhode River: tidal exchange of nutrients by brackish marshes. Estuarine, Coastal and Shelf Science 17: 651-667. Kemp, W.M., P. A. Sampou, J. Caffrey, M. Mayer, K. Henriksen and W.R. Boynton, 1990. Ammonium recycling versus denitrification in Chesapeake Bay sediments. Limnology and Oceanography 35(7): 1545-1563. Kemp, W.M., P. A. Sampou, J. Garber, J. Tuttle and W. R. Boynton. 1992. Seasonal depletion of oxygen from bottom waters of Chesapeake Bay: roles of benthic 25 and planktonic respiration and physical exchange processes. Marine Ecology Progress Series 85: 137-152. Krom, M.D and R.A Burner, 1980. The diffusion coefficients of sulfate, ammonium and phosphate ions in anoxic marine sediments. Limnology and Oceanography 25: 327-337. Lapointe, B. E. and M. W. Clark, 1992. Nutrient input from the watershed and coastal eutrophication in the Florida Keys. Estuaries 15 (4): 465-476. Lee, K. Y., T. R. Fisher and E. Rochelle-Newall, 2001. Modeling the hydrochemistry of the choptank river basin using GWLF and Arc/Info: 2. Model validation and application. Biogeochemistry 56 (3): 311-348. Li, M., A. Gargett and K. Denman. 2000. What determines seasonal and interannual variability of phytoplankton and zooplankton in strongly estuarine systems? Application to the semi-enclosed estuary of Strait of Georgia and Juan De Fuca Strait. Estuarine Coastal and Shelf Science 50 (4): 467-488. Lima, I. D. and S. C. Doney, 2004. A three-dimensional, multinutrient, and size- structured ecosystem model for the North Atlantic. Global Biogeochemical Cycles 18 (3), Art. No. GB3019. Madden, C. and W. M. Kemp. 1996. Ecosystem model of an estuarine submersed plant community: calibration and simulation of eutrophication responses. Estuaries 19 (2B), 457-474. Malone, T.C., 1992. Effects of water column processes on dissolved oxygen, nutrients, phytoplankton and zooplankton, in oxygen dynamics in the Chesapeake Bay, a synthesis of recent research, edited by D.E. Smith, M. 26 Leffler and G.E. Machiernan, Maryland Sea Grant College, Univ. of Maryland, College Park, Maryland, pp61-112. Malone, T. C., D. J. Conley, T. R. Fisher, P. M. Glibert, L. W. Harding, K. G. Sellner, 1996. Scales of nutrient-limited phytoplankton productivity in Chesapeake Bay. Estuaries 19: 371-385. Malone, T. C., L.H. Crocker, S.E. Pike and B.W. Wendler, 1988. Influences of river flow on the dynamics of phytoplankton production in a partically stratified estuary. Marine Ecology Progress Series 48: 235-249. Malone T. C., H. W. Ducklow, E. R. Peele and S. E. Pike. 1991. Picoplankton carbon flux in Chesapeake Bay. Maine Ecology Progress Series 78 (1): 11-22. Malone, T. C., W.M. Kemp, H.W. Ducklow, W.R. Boynton, J. H. Tuttle and R. B. Jonas, 1986. Lateral variation in the production and fate of phytoplankton in a partically stratified estuary. Marine Ecology Progress Series 32: 149-160. Marra, J. and C. Ho, 1993. Initiation of the spring bloom in the northeast Atlantic (47- degrees-N, 20 degrees-W)- A numerical-simulation. Deep-Sea Research II 40 (1-2): 55-73. McClain, C. R., K. Arrigo, K. S. Tai and D. Turk, 1996. Observations and simulations of physical and biological processes at ocean weather station P, 1951-1980. Journal of Geophysical Research-Oceans 101 (C2): 3697-3713. Merrill, J.Z., 1999. Tidal freshwater marshes as nutrient sinks: particulate nutrient burial and denitrification, Ph. D dissertation, Univ. of Maryland. Meyers T., J. Sickles, R. Dennis, K. Russell, J. Galloway and T. Church, 2000. Atmospheric nitrogen deposition to coastal estuaries and their watersheds, in 27 nitrogen loading in coastal water bodies: an atmospheric perspective, edited by Richard A. Valigura, American Geophysical Union. Moore, J. K., S. C. Doney, J. A. Kleypas, D. M. Glover and I. Y. Fung. 2002. An intermediate complexity marine ecosystem model for the global domain. Deep Sea Research II 49 (1-3): 403-462. Nagy, G. J., M. Gomez-Erache, C. H. Lopez and A. C. Perdomo, 2002. Distribution patters of nutrients and symptoms of eutrophication in the Rio de la Plata river estuary. Hydrobiologia 475 (1): 125-139. Nixon, S.W., 1980. Between coastal marshes and coastal waters-a review of twenty years of speculation and research on the role of salt marshes in estuarine productivity and water chemistry, in estuarine and watershed processes, edited by P. Hamilton and K.B. Macdonald, pp. 437-525, plenum press, New York. Nixon, S.W. and M.G. Pilson, 1983. Nitrogen in estuarine and coastal marine ecosystems, in nitrogen in the marine environment, edited by E.J. Carpenter and D.G. Capone, academic press, New York, pp565-648. Odum, E. P., 2000. Tidal marshes as outwelling/pulsing systems, in Concepts and controversies in tidal marsh ecology, edited by M.P. Weinstein and D.A. Kreeger, Kluwer Academic Publishers, The Netherlands. Orth R. J. and K. A. Moore, 1983. Chesapeake Bay-an unprecedented decline in submerged aquatic vegetation. Science 222 (4619): 51-53. Oschlies, A. and V. Garcon, 1999. An eddy-permitting coupled physical-biological model of the North Atlantic- 1. Sensitivity to advection numerics and mixed layer physics. Global Biogeochemical Cycles 13 (1): 135-160. 28 Pitkanen H., T. Tamminen, P. Kangas, T. Huttula, K. Kivi, H. Kuosa, J. Sarkkula, K. Eloheimo, P. Kauppila and B. Skakalsky, 1993. Late summer trophic conditions in the northeast gulf of Finland and the river Neva estuary, Baltic Sea. Estuarine Coastal and Shelf Science 37 (5): 453-474. Platt, T., and S. Sathyendranath, 1988. Ocean primary production: estimation by remote sensing at local and regional scales. Science 241: 1613-1620. Price J. F., R. A. Weller and R. Pinkel, 1986. Diurnal cycling: Observations and models of the upper ocean response to diurnal heating, cooling, and wind mixing. Journal of Geophysical Research 91 (C7): 8411-8427. Schubel, J.R. and D.W. Pritchard, 1986. Responses of upper Chesapeake Bay to variations in discharge of the Susquehanna River. Estuaries 9: 236-249. Schubel, J.R. and D.W. Pritchard, 1987. A brief physical description of the Chesapeake Bay, in contaminant problems and management of living Chesapeake Bay resources, edited by S.K. Majumdar, I.W. Hall and H.M. Austin, the Pennsylvania Academy of Science. Skogen, M. D., E. Svendsen, J. Berntsen, D. Aksnes and K. B. Ulvestad, 1995. Modeling the primary production in the North-Sea using a coupled 3- dimensional physical-chemical-biological ocean model. Estuarine Coastal and Shelf Science 41 (5): 545-565. Sheng, Y. P. 1986. A three-dimensional mathematical model of coastal,estuarine and lake currents using a boundary fitted grid, Rep. 585. ARAR Group, Titan System, Princeton, NJ. 29 Smedstad, O. M. and J. J. Obrien, 1991. Variational data assimilation and parameter- estimation in an equatorial Pacific-Ocean model. Progress in Oceanography 26 (2): 179-241. Staver, K.W, and R.B Brinsfield, 1996. Seepage of groundwater nitrate from a riparian agroecosystem into the Wye River Estuary. Estuaries 19: 359-370. Stumm W. and J.J. Morgan, 1981. Aquatic chemistry, 2nd edition, Wiley. Taft, J. L., A. J. Elliot and W. R. Taylor, 1978. Box model analysis of Chesapeake Bay ammonium and nitrate fluxes. In: Wiley, M. L. (ed.) Estuarine interactions. Academic Press, New York, pp115-130. Teal, J.M., 1962. Energy flow in the salt marsh ecosystem of Georgia. Ecology 43: 614-624. Teal, J.M. and B.L. Howes, 2000. Salt marsh values: retrospection from the end of the century, in Concepts and controversies in tidal marsh ecology, edited by M.P. Weinstein and D.A. Kreeger, Kluwer Academic Publishers, The Netherlands. Turner, R.E. and N.N. Rabalais, 1994. Coastal eutrophication near the Mississippi river delta. Nature 368: 619-621. Ullman, D. S. and R. E. Wilson, 1998. Model parameter estimation from data assimilation modeling: temporal and spatial variability in bottom drag coefficient. Journal of Geophysical Research-Oceans 103 (C3): 5531-5549. Valiela, I., M.L. Cole, J. Mcclelland, J. Hauxwell, J. Cebrian and S.B. Joye, 2000. Role of salt marshes as part of coastal landscapes, in Concepts and 30 controversies in tidal marsh ecology, edited by M.P. Weinstein and D.A. Kreeger, Kluwer Academic Publishers, The Netherlands. Valiela, I., J. M. Teal and N. Y. Persson, 1976. Production and dynamics of experimentally enriched salt-marsh vegetation-belowground biomass. Limnology and Oceanography 21 (2): 245-252. Wang, D. P., 1979. Wind-driven circulation in the Chesapeake Bay, winter 1975, Journal of physical oceanography 9: 564-572. Wang D.P. and A. J. Elliott, 1978. Non-tidal variability in the Chesapeake Bay and Potomac River: evidence for non-local forcing, Journal of physical oceanography 2: 225-232. Webb, K. L., 1988. Comment on ?Nutrient limitation of phytoplankton growth in brackish coastal ponds? by Caraco, Tamse, Boutros and Valiela (1987). Canadian Journal of Fishery and Aquatic Science 45: 380. Wetzel, R. and H. Neckles. 1986. A model of Zostera marina L. photosynthesis and growth: simulated effects of selected physical-chemical variables and biological interactions. Aquatic Botany 26: 307-323. Wu, C. R., P. T. Shaw and S. Y. Chao, 1999. Assimilating altimetric data into a South China Sea model. Journal of Geophysical Research-Oceans 104, C12, 29987- 30005. 31 Chapter 2: Assimilating High-Resolution Salinity Data Into a Model of a Partially Mixed Estuary Abstract A three-dimensional circulation model of the Chesapeake Bay is used to validate a simple data assimilation scheme, using high-resolution salinity data acquired from a ship towed undulating vehicle (a Scanfish). The simulation period spans the entire year of 1995 during which the high-resolution Scanfish data were available in July and October, lasting a few days each. Since Scanfish data were irregularly distributed in time and space, only salinity fields are nudged in the model for simplicity. Model improvements through data assimilation are evaluated from a pair of experiments: one with data assimilation and one without. Data from scattered Chesapeake Bay Program monitoring stations and a few stations maintained by the National Ocean Service inside the Bay are used independently to check the model performance. In general, the simple assimilation scheme leads to visibly improved density structures in the upper and middle reaches of the Bay. The improvement in the lower Bay is equally pronounced after data assimilation but diminishes in a shorter time scale because of faster flushing from the adjacent coastal ocean. Moderate to weak nudging normally enhances the gravitational circulation. Strong nudging may produce transient overshooting, during which the gravitational circulation is renewed vigorously. 32 1. Introduction In modeling partially mixed estuaries, a major difficulty to overcome is the numerical damping. The numerical representation of three-dimensional flow and density fields by a finite number of computation cells invariably increases friction. Part of the enhanced friction arises from grid-scale mixing, because friction coefficients must be made proportional to a power of grid spacing to achieve computational stability. Numerical form drag also enhances friction when irregular coastlines and bottom topographies are approximated by groups of computation cells. While these problems are common to all ocean models, they become particularly acute in models of long and narrow estuaries with excessive coastline and topography irregularities. The bottom inflow must follow a long and often sinuous path to enter, upwell and return seaward. To overcome numerical damping, it is often necessary to enhance bottom inflow of seawater from the mouth region in order to produce a realistic two-layered circulation well inside an estuary. There are two ways to further enhance model realism. One way is to improve grid resolution at the expense of computation speeds. The other way is through assimilation of high-resolution data. With the availability of satellite altimeter and climatological data sets, data assimilation is now widely used in large-scale ocean models and dynamic principles have been developed for the purpose of nudging several variables simultaneously (e.g., Ezer and Mellor, 1994; Forbes and Brown, 1996; and Wu et al., 1999). Similar efforts in shallow reaches of the ocean are deliberately simplified for lack of climatological data sets and reliable altimeter data. For example, Spitz and Klinck (1998) used an adjoint variational method along with 33 tide gauge data from a few stations to improve sea level predictions in a two- dimensional tidal model of the Chesapeake Bay. Similar methods were also used to assimilate current velocity data into shallow water equation models of Massachusetts Bay (Bogden et al., 1996) and Long Island Sound (Bogden and O?Donnell, 1998). Because of the resolution problem, highly nonlinear phenomena such as internal bore intrusion and sill-induced hydraulic jumps in estuaries are often smeared out by friction in numerical models. In this light, successful assimilation of high- resolution data would be highly desirable for it allows modelers to reproduce these physical processes in more realistic settings. This remains as a lofty goal at the present time. Recent advances in Undulating Oceanographic Recorders (such as Scanfish manufactured by Danish company Geological & Marine Instrumentation) offer an alternative. Through rapid vertical undulations, a ship-towed vehicle can provide a reasonably synoptic, three-dimensional view of the density structure in a large body of shallow waters. The high-resolution data, though irregularly distributed in space and time, may be assimilated into numerical models. Ideally, one would like to derive climatological data sets for shallow bodies of waters from repeated sampling over many years, and assimilate climatological data into models in a dynamically consistent fashion. This option is presently not feasible for obvious reasons. An attempt is made below to assimilate the Scanfish data into a Chesapeake Bay three-dimensional circulation model. Year 1995 was chosen because it was the first year the high-resolution salinity data became available. Further, the hydrodynamic model simulation has not been carried out and verified beyond 1995 at 34 the time this research was initiated. The simulation period spans the entire year of 1995, during which two rapid sampling cruises covering the main stem of the Bay were made in July and October. In the same year, 49 mainstem monitoring stations were maintained by Chesapeake Bay Program of the Environmental Protection Agency. Each station was sampled 16 to 20 times in the year at irregular time intervals. Only Scanfish data from the July and October cruises are assimilated. Selected salinity data from the 49 fixed stations are used to evaluate the performance of the data assimilation scheme. Intuitively, direct assimilation of hydrographic data into a model seems like an effective way to improve model realism. The ideal scenario is that assimilation improves the density structure, and the improved density field supports a more realistic circulation field. While this is generally true, the improvement does not come without penalties. In the subject at hand, it is found that quick injections of data may trigger brief moments of readjustment in gravitational circulation. Circulation during brief periods of gravitational readjustment may be unrealistic. In this light, the speed of data injection must be optimized, so that the gain will outweigh the loss. The oceanographic setting and data availability are described in section 2. In section 3, a hydrodynamic model of the Chesapeake Bay is described. Section 4 discusses the data assimilation scheme. The undesirable consequence of data assimilation, i.e., the gravitational readjustment, is elaborated in section 5. Section 6 summarizes benefits brought about by data assimilation. Section 7 concludes this work. 35 2. Oceanographic Setting A deep channel running north-south more or less along the western side of the main stem dominates the bathymetry in the middle reaches of the Chesapeake Bay (Fig. 2.1). The main channel, being bounded to the south by a sill at about 37.6?N, is completely closed below the sill depth of about 14 m. South of 37.6?N, the deep channel becomes somewhat shallower and ill-defined, often branching in multiple directions. Between 37.6?N and 39?N, the main deep channel harbors a rather persistent, river-induced two-layered circulation although the gravitational circulation is often influenced by winds and stratification (Goodrich et al., 1987). Along the main stem of the Bay, drainage from eight major tributaries (Susquehanna, Patapsco, Patuxent, Potomac, Rappahannock, York, James and Choptank) contributes to most of the river input. The Susquehanna River in the northern extreme of the Bay provides the largest freshwater influx among the eight, approximately one half of the total freshwater input. Tidal forcing is modest inside the Bay with tidal range rarely exceeding 1 m. Winds are generally episodic with dominant periods of 2-7 days. In the upper and middle reaches of the Bay, northwesterly winds dominate in winter months (November-February), but are more frequently disrupted by southerly winds of several days each in summer. Fig. 2.2 shows daily discharge rates from the four largest tributaries (Susquehanna, Potomac, James and Rappahannock) for the entire year of 1995, which was perceived as an abnormally dry year. Discharge from other rivers was considerably lower. Discharge was generally high from mid-January to April, further enhanced by several peaks of 5-10 days duration. It decreased markedly in summer 36 and was relatively high again in late fall. When averaged over sufficiently long periods to filter out wind and tidal effects, the annual variation of subtidal circulation in the Chesapeake Bay is expected to be dominated by the strength of freshwater discharge. See Goodrich et al. (1987) for some observations. Fig. 2.3 shows 24 Scanfish transects across the main stem of the Chesapeake Bay. Temperature and salinity were sampled along these transects from July 19 to July 22 and from October 24 to October 26 in 1995. Starting from transect 1 across the mouth of Chesapeake Bay, transects were visited sequentially as the ship moves up the estuary. Moving at an average speed of 6 knots, the ship covers each transect in about half an hour to four hours, and it takes about 3 days to complete a basinwide survey. The Scanfish follows slanted paths up and down the water column with inclination angles around 6?, sampling at time intervals of 0.5 seconds. During the July cruise, data along transect 1, transect 4 and transect 7 to 12 in the lower Bay were unfortunately lost. In consequence, effects of data assimilation in July come mostly from inner reaches of the Bay. Selected salinity stations maintained by Chesapeake Bay Monitoring Program are marked by dots in Fig. 2.3. These stations were visited 16 to 20 times in 1995 at somewhat irregular intervals. At each station, salinity profiles were measured with a vertical resolution of 1~2 m. These salinity profiles provide an independent data set that can be used to assess how well the data assimilation scheme works. Two temperature stations at Tolchester Beach and Solomons Island, marked by ?X? in Fig.2.3, were maintained by the National Ocean Service of NOAA in 1995 (http://www.co-ops.nos.noaa.gov/data-res.html). These time series can also be used 37 for model verification, although our major emphasis is on salinity as the major indicator of water density in the Bay. 3. The hydrodynamic model The model, originally developed by Sheng (1986), was subsequently modified extensively by the US Waterways Experiment Station (WES) for application to Chesapeake Bay (Johnson et al., 1991; Wang and Chapman, 1995). Under Boussinesq and hydrostatic approximations, the hydrodynamic model solves for salinity, temperature, water-level elevation and velocities in three dimensions. There are up to 19 layers in the vertical with a uniform layer thickness of 1.52 m, except that the top layer thickness fluctuates with sea level. Horizontally, the governing equations in the Cartesian coordinate system are recast in a boundary-fitted curvilinear coordinate system (Fig. 2.4) to cope with the irregular shoreline configuration and deep channel orientation. The model domain extends offshore to include a piece of coastal ocean with coarse resolution. The coastal ocean is included mainly as a buffer zone to facilitate free exchanges across the Bay mouth. Inside the Bay, typical grid size ranges from 1 to 5 km in the main stem of the Bay; bottom topographic irregularities with horizontal scales in and below this range are truncated by the model. Further, the prominent sill bounding the main deep channel in the south, located between the mouths of Rappahannock and Potomac Rivers, is marginally resolved. In spite of the coarse resolution, essential circulation features, such as the two-layered circulation in the main channel and major tributaries, can be reproduced by the model (Johnson et al., 1991; Hood et al., 1999). Similar to a host of primitive-equation models such as Blumberg and Mellor (1987), the staggered 38 Arakawa-C grid system is used in both horizontal and vertical directions of the computation domain. The vertical eddy viscosity and diffusivity are computed from mean flow and stratification characteristics using the second-order k-? turbulent closure scheme (see, for example, Kundu, 1980; Launder and Spalding, 1974). A quadratic stress is exerted at the bottom, assuming the bottom boundary layer is logarithmic over a bottom roughness height of 0.05 cm. Coefficients of horizontal eddy viscosity and diffusivity are set to 10 4 cm 2 /s. Originally, the initial three-dimensional salinity and temperature fields were constructed using the historical field data in January averaged over many years. Since initialization, this model simulation has been extended from 1985 to 1994 by WES. We used the model output at the end of 1994 as the initial condition for salinity and temperature fields. The initial velocity field was taken to be zero and the water surface was initially set at the mean sea level. The model is subsequently forced by open ocean tides, winds, freshwater inflows and the heat exchange at the water surface through 1995. Further, salinity and temperature fields were also prescribed on offshore open boundaries using monthly Levitus climatology data (Levitus, 1982) combined with field data at Duck, North Carolina (36.1833?N, 75.7467?W) acquired daily (with occasional lapses) by the Field Research Facility of the US Army Corps of Engineers. Daily freshwater inflow with zero salinity and time-varying temperature was prescribed for the eight major tributaries; arrows in Fig. 2.1 mark inflow locations. On each inflow cross-section, the incoming current is uniform with time-varying speeds regulated by the daily freshwater discharge rate. 39 Hourly wind stress in the lower and middle reaches of the Bay was linearly interpolated from data at the Norfolk International Airport (NIA), Patuxent River Naval Station (PRNS) and the Baltimore-Washington International (BWI) Airport. Their locations are marked by solid squares in Fig. 2.3. North of BWI, wind stress is assumed to be identical to that at BWI. Empirical factors for different regions were used to extrapolate winds over land to winds over water. Daily air-water heat exchange was computed from data taken at the Patuxent meteorological station using the formulation of Edinger, Bradly and Geyer (1974). Ideally, meteorological stations over the water are desirable, but few offshore stations were available in 1995. In constructing the wind field for the Bay model, it should be noted that longitudinal winds are much more effective than lateral winds in driving circulation along the main stem of the Bay (Wang, 1979a and b). The linear interpolation among the three meteorological stations (NIA, PRNS and BWI) is intended to improve spatial resolution of longitudinal winds along the main axis of the Bay. One reviewer of this manuscript pointed out two additional records at Solomons Island and Tolchester Beach, maintained by the National Ocean Service of NOAA. From a basin-wide perspective, the Solomons Island station and PRNS are practically at the same location. The inclusion of Solomons Island, therefore, will not improve the spatial resolution of winds along the main axis of the Bay. Tolchester Beach station and BWI are also at about the same latitude. If the two wind records differ substantially, the inclusion of Tolchester Beach will likely improve lateral resolution of wind forcing in the upper reaches of the Bay, but does little to enhance longitudinal resolution of the wind field. If winds at Tolchester Beach and BWI do not differ substantially, then 40 there is no reason to include it. To investigate further, we have performed cross- spectral analysis to document the relationship between winds at Solomons Island and PRNS, and at Tolchester Beach and BWI. Of major concern are low-frequency longitudinal winds with periods longer than a few days along the main axis of the Bay. High-frequency winds and cross-estuary winds are basically noise generators, ineffective to drive basin-scale subtidal circulation. In periods longer than about 2.5 days, the coherence squared between Solomons Island and PRNS is about 0.8 for longitudinal winds and approaches 0.7 for lateral winds. The corresponding coherence squared between Tolchester Beach and BWI approaches 0.76 for longitudinal winds and 0.68 for lateral winds. Further, the phase lag between each pair of stations is generally less than a few hours for longitudinal winds. Therefore, the subtidal circulation in the main stem of the Bay will not be significantly impacted whether the additional wind records are included or not. As a consistency check, we have also compared the modeled and observed surface water temperature at Tolchester Beach and Solomons Island in 1995. The result, to be shown later in Fig. 2.6, shows reasonable agreement even without the inclusion of the two additional wind records, lending support to the foregoing argument. Open-ocean boundary sea level was updated using data from stations at Wachapregue, VA (37.6067?N, 75.6867?W) and Duck, NC (36.1833?N, 75.7467?W) obtained from National Ocean Service (NOS), NOAA. These coastal sea-level data were first extrapolated offshore based on Green?s Law (Ippen, 1966). Namely, tidal amplitudes are assumed to be inversely proportional to the ? power of local water depth. In the along-shore direction, tidal amplitudes are linearly interpolated on the 41 offshore open boundary. While water level fluctuations are prescribed on open-ocean boundaries, the incoming and outgoing currents are induced by water level gradients normal to these boundaries. The model solves external and internal mode equations separately. The external mode consists of equations for the water surface elevation and vertically averaged flows in two horizontal directions. The internal mode computes the vertical shear of horizontal velocities, vertical velocity, temperature and salinity. Time steps for the external and internal modes are both set at 300s. The larger-than-normal time step for the external mode is made possible by an implicit solver which relaxes the stringent requirement for small time steps set by the Courant-Friedrichs-Levy computational stability criterion. Before data assimilation, the hydrodynamic model was tuned to reproduce observed surface salinity in the upper and middle reaches of the deep Channel. The initial tuning includes minor adjustment in the bottom topography, vertical mixing parameters and salinity on open-ocean boundaries. Fig. 2.5 shows the model- produced time series of surface salinity at stations CB3.3C and CB5.1 in the upper and middle reaches of the deep channel, respectively. Superimposed are corresponding data points that agree with the model reasonably well. The quality control at the two points over the deep channel ensures comparable model performance in the vicinity, at least near the water surface. Figure 2.6 compares the modeled and observed water surface temperature at Tolchester Beach and Solomons Island in 1995. The model outputs were retrieved at half-hour intervals while the observed time series were at hourly intervals. In general, 42 the model reproduces the seasonal trend of water surface temperature reasonably well, although the model tends to overestimate surface temperature in winter months. Note that the hydrodynamic model does not have an ice layer component and therefore cannot simulate occasional winter freezing events in shallow reaches of the basin. This deficiency is likely to cause some discrepancies in winter. The problem is not serious because the discrepancy diminishes quickly in warmer months. Leaving the model-data agreement aside, sizeable discrepancies still exist at depths and laterally. Figure 2.7 illustrates the general pattern of discrepancies by comparing model results with Scanfish data along two selected transects (sections 13 and 20 in Fig. 2.3) in July. The lower panels show sections of salinity fields derived from Scanfish data, while upper panels show corresponding sections retrieved from the model, following the same tracks and sampling intervals of the Scanfish. Section 13 (left panels) and section 20 (right panels) are in the middle and upper reaches of the Bay, respectively. The comparison points out a dominant trend. Namely, the model tends to overestimate salinity at depths and the deviation increases toward lower reaches of the Bay. A similar comparison in October (Fig. 2.8) leads to the same conclusion. Because of the availability of transects in lower reaches of the Bay, section 2 (left panels) and section 16 (right panels) are chosen to facilitate the comparison in lower and middle reaches of the Bay, respectively. In the lower Bay (left panels), the model-derived salinity is considerably higher than observations, and the discrepancy increases with depth. In the middle reaches (right panels), the modeled salinity structures compare favorably with those observed. 43 Figure 2.9 shows longitudinal sections of monthly averaged circulation and salinity field derived from the model in February (top panel) and October (bottom panel). The vertical slice follows the center axis of the deep channel southward to the mouth of Rappahannock River and thereafter extrapolates farther southward to the southern land boundary. The monthly averaging removes most of the wind and tidal influences and the residual circulation is mostly gravitational. The peak discharge from tributaries in January results in a pronounced two-layered circulation in February. The bottom inflow is visibly much stronger than the surface outflow because of the lateral confinement at depths. October is in the end of a long dry period and the two-layered circulation is much weaker. Waters in the deep channel also becomes saltier in the dry period. Conceivably, further tuning of the model will further reduce the discrepancies as illustrated in Figs. 2.7 and 2.8, but the point of diminishing return will be reached soon if the model resolution remains the same. As we stated earlier, the model resolution is ultimately responsible for this type of discrepancies. With coarse resolution, the bottom inflow is partially choked by numerical damping near the estuary mouth, and therefore must be enhanced in order to reproduce observed salinity structures in the middle and upper reaches of the Bay. In consequence, the model-produced bottom inflow becomes saltier, especially in the lower reaches of the Bay. The intention of data assimilation is therefore to reduce modeled salinity at depths and in lower reaches of the Bay. 44 4. Data assimilation scheme The hydrodynamic model receives irregularly spaced time series of Scanfish data through the salinity equation, using a Newtonian relaxation scheme (Anthes, 1974). Since salinity is the major indicator of water density in this region, temperature data are excluded from assimilation for simplicity. Let x be the longitudinal axis, y be the lateral axis and z be the vertical axis. At a given time (t = t 0 ), a model grid point at (x 0 ,y 0 ,z 0 ) may receive Scanfish data from distributed points (x i , y i , z i , t i ) in a four dimensional neighborhood. The governing equation for salinity (S) is )1(]}/)(/)( /)(/)({exp[)(][/ 22 0 22 0 22 0 22 0 TttZzz YyyXxxMaxSSKdiffusionDtDS ii ii obs ??? ??????+= where D/Dt is the substantial differential operator and [diffusion] accounts for turbulent mixing in three dimensions. In addition to advection and diffusion, the nudging term in (1) restores observed salinity at a fixed rate K. The Gaussian dependence in space ensures that the influence of a data point on a model grid point decays with distance away from the data point. The appropriate length scale of spatial decay is X in the longitudinal direction, Y in the lateral direction and Z in the vertical. A Gaussian time dependence ensures active data injection in a timescale (T) before and after the data arrival. At any given time step of data insertion, a model grid point must choose a point among distributed Scanfish data in a four-dimensional neighborhood to receive data. The winning data point at (x i , y i , z i , t i ) is the point making maximum contribution to a model grid point at (x 0 , y 0 , z 0 , t 0 ). In this light, the 45 observed salinity (S obs ) in eq. (1) represents the salinity value of the winning data point. The assimilation procedure as outlined in the proceeding paragraph is computationally demanding. At any given time step of integration, a winning data point must be chosen among millions for every grid point inside the model domain. The searching procedure is cumbersome and arises solely because the irregular distribution of Scanfish data is highly incompatible with modeled salinity fields. A few measures can be taken to speed up the search. For example, one can switch off the search routine if a model grid point is sufficiently away from Scanfish data in time or space. Further, the resolution of Scanfish data is unnecessarily high in terms of model needs. To enhance the efficiency of searching, the Scanfish data were sub- sampled at intervals of 5 seconds before they were used for data assimilation. For the assimilated results shown below, the salinity restoration rate (K) is chosen to be (15 hr)-1. The value of K needs to be large enough to make an impact while being small enough to avoid excitation of gravity waves. Haltiner and Williams (1980) suggested that the timescale for K should be smaller than the dominant timescale contained in observations. In anticipation of a fast changing estuarine environment, our timescale for K is considerably shorter than typical values used in open-ocean settings (Sarmiento and Bryan, 1982). While the nudging improves the modeled salinity fields, it also triggers brief moments of readjustment in gravitational circulation. The readjustment may occur during and shortly after the data insertion period, and brings unrealistic features into the model for a short period of time. In this light, K is optimized to maximize the gain and minimize the loss caused by data 46 injection. The choice of e-folding timescale, T = 6 hr, is comparable to the timescale of semidiurnal tides. Nudging length scales (X, Y and Z) have also been optimized. For the solution shown below, longitudinal and lateral (X and Y) scales are set to 40 km and 3 km, respectively. Vertical scale (Z) is considerably shorter; chosen to be 2 m. Model sensitivity to K, T, X, Y, Z will be discussed later, after the discussion of main results. 5. Gravitational readjustment The data assimilation may trigger renewed gravitational circulation because the density structure is significantly altered. Since each model or assimilation period varies in oceanographic setting and data availability, it is difficult to predict the timing and longitudinal extent of the readjustment. Nevertheless, the readjustment process documented below is likely to be encountered in a variety of models of long estuaries during periods of strong or moderate data injection. Prognostic models of long estuaries, if properly tuned to produce realistic features in the inner reaches, are likely to overestimate the salinity of bottom inflow especially near the estuary mouth. When data are inserted, the salinity is reduced in the mouth region and/or deeper reaches of the basin. The buoyant outflow may intrude farther seaward in response to the altered density structure. The strength of the renewed seaward expansion depends on the data injection speed. Further data injection in time will eliminate the undesirable transient and move the solution back to reality. In this model, high-resolution data were available only briefly in July and October. The data injection must be sufficiently strong to make a lasting impact. As a 47 result of the strong nudging, the gravitational readjustment occurs preferably in the early stage of data assimilation in October, soon after the salinity data in lower reaches of the Bay are inserted. The readjustment will not occur if either July or October is excluded from data assimilation. The combination of the two assimilation periods is necessary to trigger it. Figure 2.10 shows flow and salinity fields on the longitudinal-vertical section as in Fig. 2.9 before, during and after the gravitational readjustment in October. These snapshots are instantaneous, so that wind- and tide-induced currents are included. Locations of 16 psu isohalines are marked by arrows on top of each panel to highlight the gravitational readjustment. Shortly before data arrival (top panel), brackish water is confined in upper reaches and waters in lower reaches are quite saline. The middle panel shows the same vertical section 24 hours later. The time corresponds to 8 hours after the beginning of active data assimilation or 2 hours after the Scanfish was deployed in October. Recall that active data assimilation begins in an e-folding time scale (T = 6 hr) before data arrival. At this time, data are inserted only in regions around and seaward of transect 5 in Fig. 2.3. Nevertheless, the limited amount of data insertion is able to trigger a gravitational readjustment. As indicated by the 16-psu isohaline, the buoyant surface layer expands seaward by about 70 km in one spurt. Further, waters in lower reaches are freshened by about 2 psu or so. Thereafter, continuous data insertion would eliminate the artificial seaward expansion. The bottom panel shows the same section 10 hours after the October data assimilation ends. The snapshot is taken 84 hours after the middle panel. The artificial seaward expansion of buoyant layer is eliminated and the basin-wide salinity structure is 48 moved closer to observations. Our analysis of root-mean-square errors in section 6 will confirm this point. It is worth pointing out that the salinity restoration rate (K) is a crucial parameter controlling the strength and longitudinal extent of gravitational readjustment. With a larger K, the seaward expansion of the buoyant layer is greater, but subsequent data injection in middle and upper reaches of the Bay also eliminates the readjustment at a faster rate. In the other extreme, the readjustment process can be eliminated if K becomes exceedingly small. Leaving K aside, the longitudinal extent of data assimilation also influences the strength of gravitational readjustment. For example, one could limit the data insertion to upper and middle reaches of the Bay only. The consequent gravitational readjustment would be weaker. It is conceptually useful to interpret the renewed gravitational circulation in terms of pressure changes. In the model-derived two-layered circulation, the bottom inflow upwells and returns as a surface outflow. The gyre is maintained by proper pressure gradients. When data are assimilated in lower reaches of the Bay, pressure is reduced near the mouth. The consequent increase in the seaward pressure gradient triggers the seaward expansion of the buoyant layer. If data are inserted only in middle and upper reaches, the effect is essentially to reduce salinity at depths. The consequent pressure deficit at depths enhances the bottom inflow of saline water from the mouth region. Thereafter the pressure field is temporarily reduced in lower reaches due to the sudden loss of salinity. This may also cause the layer of brackish water to expand seaward. Following this line of reasoning, the data insertion in lower Bay has the immediate effect to encourage seaward expansion of buoyant layer. The 49 data injection in middle and upper reaches may also encourage the seaward expansion, but only after the pressure field in lower Bay is reduced. Thus, data assimilation in the lower Bay has a more profound effect than insertions elsewhere, insofar as the renewed gravitational readjustment is concerned. Our preliminary experiments by varying the region of assimilation lend support to the foregoing conclusion. 6. Model improvement through data assimilation Despite the undesirable consequence of gravitational readjustment, the agreement between model and the observations after data assimilation is generally improved. Salinity measurements at scattered stations in 1995 (Fig. 2.3) provide an independent data set to evaluate the model performance. Each station was visited 16~20 times in 1995 at somewhat irregular intervals with a vertical resolution of 1~2 m. Discrepancies between the model and data are first evaluated in terms of root- mean-square (RMS) errors. Figure 2.11 illustrates the RMS error as a function of time in the upper, middle and lower reaches of the Bay. In the top panel, each ensemble contains all data points collected north of the Choptank River in each month. The RMS error, ranging up to several psu, may have under-represented the model?s prognostic skill because salinity stations are fixed in space. A slight longitudinal shift of salinity patterns produced by the model may be seen as a large error at a fixed station. Leaving the magnitude of RMS errors aside, model improvements through data assimilation are apparent. Results from the pair of experiments, one with data assimilation and one without, contrast the difference brought about by data injections. 50 The reduction in RMS errors is maximum in the month immediately following the July and October assimilation, and decreases slowly thereafter. The model performance is similar at all salinity stations in the middle reaches of the Bay; therefore, only one station (CB4.3C) is chosen to illustrate the RMS error (middle panel). At this station, data are not grouped for each month and each ensemble consists of only a vertical profile of salinity with 1~2 m resolution. The improvement brought about by data assimilation is generally more profound in middle reaches than in upper reaches, except during the brief period of gravitational readjustment. At the end of July, the RMS error decreases twofold as a result of data insertions. The improvement diminishes slowly in time thereafter. After October data assimilation, the RMS error actually increases over a brief period in November as a result of the gravitational readjustment. Thereafter the RMS error decreases again after the adjustment is over. Data insertions generally also enhance the model performance in the lower Bay, but the improvement does not persist for a long time because the adjacent coastal sea is excluded from data assimilation. As in the top panel, each ensemble in the bottom panel of Fig. 2.11 contains all data points collected from stations south of Potomac River (Fig. 2.3) in each month. The RMS errors decrease little after data injections. On a longer time scale, marginal improvements resulted in the lower Bay despite massive injections of data with high spatial resolution. Shorter memories of flushing time scales in the lower Bay are responsible for the deficiency. Fig. 2.12 shows the model-produced variations of salinity structures in a zonal section midway 51 between Scanfish transects 3 and 4 as indicated in Fig. 2.3. Before the data arrival in October, the salinity field (top panel) shows considerable stratification. In reality, waters in the lower Bay are typically less stratified in winter. Destratification occurs shortly after the data injection (middle panel), bringing the model closer to reality. Thereafter the stratification returns in time (bottom panel). Apparently, the bottom inflow from the coastal ocean tends to reestablish the stratification. More illustrations of short flushing time scale in the lower Bay will be given later in Fig. 2.15. There are two ways to improve the model performance in the lower Bay. One way is to continuously inject data with high temporal resolution in the lower Bay. This measure, however, would defeat the purpose of the data assimilation scheme, which is meant to find ways to make lasting model improvements through occasional data injection. A more reasonable way would be to extend the assimilation areas to the adjacent coastal ocean. Without data assimilation, the salinity of the coastal ocean is highly constrained by climatology; subsequent intrusion into the Bay tends to offset the data injection effort especially in the lower Bay. If high-resolution salinity data are available in the coastal ocean during the data assimilation period, assimilation of these data in the coastal ocean would sustain the effect of data assimilation in the lower Bay for a longer period of time. This recommendation is not heeded herein for lack of qualified data off the Bay mouth. In general, model improvements through the data assimilation are not depth- sensitive. Figure 2.13 shows observed and modeled salinity profiles at selected stations. Top panels are derived from station CB3.3C in the upper Bay, while bottom panels are from station CB5.1 in the mid-Bay. Surface salinity data from these two 52 stations have been used in the initial tuning of the model before data assimilation (Fig. 2.5). Figure 2.13 indicates that the model improvements after July assimilation (left panels) and October assimilation (right panels) do not favor a particular depth. With data assimilation, the model-derived salinity profiles generally shift toward observed profiles at all depths with few exceptions. In the absence of concurrent flow measurements, it is a bit uncertain whether the data assimilation actually improves the model-produced flow fields. However, the collective wisdom from previous modeling experiences suggests that a more realistic density structure often supports a more realistic flow field. It is highly likely that the flow fields after the data assimilation are more realistic. Figure 2.14 illustrates the changes in the flow field induced by the data assimilation. Left panels are biweekly averaged flow fields at the surface (top panel) and 10 m below mean water level (bottom panel) without data insertions. The time average is over the middle two weeks in November (from day 310.5 to day 324.5). Since wind and tidal effects are filtered out through time averaging, the patterns are dominated by surface outflow and bottom inflow. The right panels show the corresponding difference caused by data assimilation. Leaving minor variations aside, it is clear that the data assimilation essentially enhances both the surface outflow and bottom inflow. The speed enhancement ranges up to about 4 cm/s. Similar changes in the circulation pattern were also found in August and September (not shown here) after data assimilation in July. The result is not surprising in light of the fact that most ocean models tend to underestimate the strength of density-driven currents because of grid-scale mixing. 53 Thus, assimilation of hydrographic data appears to offer a remedy to offset numerical damping. Dynamically, the enhanced bottom inflow and surface outflow can be regarded as a renewed adjustment under gravity. To illustrate this, Fig. 2.15 shows the time-averaged longitudinal section of density anomalies induced by data assimilation in October. The longitudinal section is the same as in Figs. 2.9 and 2.10. Further, the time averaging is from day 295 to day 310.55 in the top panel, and from day 310.55 to day 326.1 in the bottom panel. The averaging period in the top panel covers the time span of active data assimilation (roughly from day 296 to day 298). The density anomaly is obtained by subtracting model results without data assimilation from that with data assimilation. Since the time span extends to well after the period of active data assimilation, the density anomaly in the top panel of Fig. 2.15 does not correspond to a static change brought about solely by data assimilation. Dynamic adjustments also set in to change the density structure. Despite the complication, simple analyses below suggest that the density anomaly is becoming an integral part of the two-layered circulation. The data assimilation essentially induces density deficits that intensify toward the bottom of the deep channel. Data injections also produce a patch of density surplus in upper depths, confined mostly in the upper and middle reaches of the Bay. By comparison, the bottom-trapped density deficit is the most dominant signal. Roughly speaking, the density deficit in the top panel of Fig. 2.15 is characterized by ?? t = 0.8 over the bottom 10 m of the deep channel. The associated speed of internal gravity waves (c 0 ) is about 28 cm s -1 . The density deficit is mostly confined inside the 54 estuary because the adjacent coastal ocean is excluded from data assimilation. In consequence, this density deficit would induce bottom inflow from the coastal ocean. Let u 0 be the characteristic inflow speed. In the absence of topography drag, mixing and friction, inviscid theories such as Benjamin (1968) would suggest an internal Froude number (u 0 /c 0 ) of order one and the two-layered circulation would be enhanced by 30 cm s -1 or so. The actual enhancement is only about 5 cm s -1 in Fig 2.14, suggesting a characteristic Froude number much below one in this partially mixed estuary. The low Froude number governs not only the perturbation field induced by data assimilation, but also the mean circulation in the Chesapeake Bay as well. Relative to the seawater density near the Bay mouth, the density deficit (?? t ) in the upper reaches of the Bay is often in excess of 10 [see Goodrich et al. (1987) or Fig. 2.9]. The characteristic c 0 associated with this density deficit is about 100 cm s -1 . Mean speeds of bottom inflow are generally below 20 cm s -1 [see Goodrich et al. (1987) or Fig. 2.14]. Thus, the basin-scale mean circulation is also governed by a similarly low Froude number. The bottom panel of Fig. 2.15 provides an alternative to illustrate the faster loss of memory in the lower Bay. Generally speaking, the density anomalies decrease slowly in time after the data injection. The enhancement of two-layered circulation by data injections also decreases in time accordingly. However, the density deficit diminishes much faster in the lower reaches of the Bay, disappearing almost completely in the bottom panel of Fig. 2.15. As commented on earlier, the fast 55 disappearance arises because the adjacent coastal ocean is excluded from data assimilation. 7. Discussion and conclusions Using a Chesapeake Bay model as a test case, assimilation of high-resolution Scanfish data proved to be a useful tool to enhance model performance if certain precautions are taken to minimize volatile transients induced by fast data injections. If nudging is strong, the consequent transient may manifest as renewed gravitational circulation. Subsequent data assimilation will eliminate the transient overshooting. Brief lapses of model accuracy may be inconvenient if one desires to obtain a continuous quality output. If this is the major concern, one can blend in the model result without data assimilation using a time varying weighting function to eliminate the undesirable transient. Given a few narrow windows of high-resolution data in a year, the nudging must be strong enough to make a difference but also weak enough to minimize possible gravitational readjustment. The precaution is necessary because of the limited availability of high-resolution data. Ideally, the restoration rate of hydrographic data can be reduced to a bare minimum if the data string is more or less continuous in time. In this idealized setting, continuous nudging in time will minimize model discrepancies and the nudging rate need not be large because the discrepancy is kept small at all times. It is highly likely that the undesirable overshooting can be eliminated in this limit. Leaving the idealized scenario aside, the choice of nudging rate (K) must be optimized to maximize the gain and minimize the 56 loss brought about by the data assimilation. When this is done, the gravitational circulation is normally enhanced after the data insertion, and the enhancement will last for months. In theory, the amount of salinity anomaly received by the model from a single data point is linearly related to a four-dimensional volume (TXYZ) by the constant K. Since data continue to enter the model from different locations, the real relationship is quite complex. Nevertheless, K and T should be chosen to be inversely proportional to each other in order to approximately maintain the same intensity of nudging. Leaving K and T aside, the assimilation scheme still involves choices of proper length scales (X, Y and Z) in the longitudinal, lateral and vertical directions. The model is generally not sensitive to these choices as long as we maintain proper aspect ratios pertaining to the estuary basin. As a rule of thumb, the choice of X must be commensurate with the tidal excursion length in the longitudinal direction. Once X is fixed, Y and Z can be chosen proportionally to maintain the aspect ratios of the basin. After these choices are made, moderate variations in parameter space do not profoundly impact the model response. In the long run, repeated acquisitions of high-resolution hydrographic data would lead to climatological data sets for the basin. Climatological data are regularly distributed in time and space, and therefore can be assimilated more efficiently into models. Additional gains from regularly spaced climatological data can also be anticipated in the future if we borrow similar experiences from the open ocean modeling community. Through assimilation of regularly spaced climatological data, we may be able to adjust other variables such as sea level, vorticity and currents in a 57 dynamically consistent manner to maximize the gain. Similar methodology has been advanced considerably in the open ocean setting; see, for example, Ezer and Mellor (1994), Forbes and Brown (1996) and Wu et al. (1999) for several interesting applications. At the present time, it is not clear what dynamic constraints should be enforced when nudging several variables simultaneously in a tidally dominated estuarine environment. While the methodology still awaits future development, the simple assimilation scheme presented herein draws attention to this issue and makes a modest start. Appendix A: A shorter assimilation time scale One reviewer suggested that a shorter e-folding time scale (T = 1~2 hours) should be used for data assimilation in estuaries. Obviously, the choice of T should be constrained by the dominant tidal period (12.42 hours). Our numerical results indicate that this is a loose constraint. As long as T is not completely decoupled from the tidal period, the assimilation scheme achieves similar results if the restoration rate (K) and T are inversely proportional to each other. Taking the reviewer?s comment as an example, we can reduce T threefold (from 6 to 2 hours) and increase K threefold (from (15 hr) -1 to (5 hr) -1 ) to achieve similar effects. Figure 2.16 illustrates resulting RMS errors as functions of time in the upper, middle and lower reaches of the Bay. This figure is produced following the same recipe as that of Fig. 2.11; a comparison between the two figures highlights the insensitivity of the assimilation scheme to the e-folding assimilation time scale. In the upper Bay, the new e-folding time scale leads to similar reduction in RMS errors. In 58 the middle reaches, improvements after data assimilation and discrepancies induced by the brief gravitational readjustment after October are comparable to that in Fig. 2.11. Improvements in the lower Bay are as marginal as before. Since the difference between Figs. 2.11 and 2.16 is small, other details will not be presented. Appendix B: Modification of salinity at open-ocean boundaries Even though the gravitational readjustment that is induced as a result of the data injection is transient, it is problematic for carrying out biogeochemical calculations. The seaward expansion of the buoyant layer may wash out the biological organisms and produce an unrealistic distribution of biological variables, which will not be restored to reality by further physical data assimilation. Therefore, we made another attempt to improve the modeled density structure by adjusting the salinity at the open-ocean boundaries to set the stage for the biogeochemical modeling work described in the subsequent chapters. Without assimilation the model tends to overestimate salinity at depth in both years. Beside the numerical damping problem we discussed above, another possible reason for this could reside in the specification of the open-ocean boundary conditions from monthly Levitus climatological data, i.e., we use long term averaged seasonal salinity at the open-ocean boundaries that may not represent the specific year very well. We therefore set out to determine how the prescribed salinity values at the open-ocean boundary might be modified to give the best fit to the observed salinity values in the inner estuary. Because of the complexity of the three-dimensional numerical model and irregularity of available observations (the time series data are 59 available at scattered stations and the high resolution (scanfish) data are available only over short periods of time), the process of adjusting the boundary conditions was done manually and the optimization was achieved by try-and-error. We first estimated that the time lag for boundary salinity to affect the upper bay bottom salinity is about 40 days. Then corrections were made gradually at open-ocean boundaries to minimize the discrepancies between observed and modeled bottom salinity at the upper and mid bay (CB3.3C and CB5.1). The final adjustments we made to the boundary conditions are listed in table 2.1. These adjustments significantly improved the salinity structure at depth, especially in the upper and mid reaches of the bay shown here in the bottom salinity comparisons in both years (Fig. 2.17). This, in turn, resulted in substantial improvements in the agreement between the modeled and observed salinity profiles in both 1995 and 1996 without the need for data assimilation. Table B1: The time period and corresponding values of salinity change made at the open-ocean boundaries. Time Period Day 1-70, 1995 Day 130- 180, 1995 Day 10-60, 1996 Day 140- 200, 1996 Day 220- 290, 1996 Salinity Change -4.0 -2.0 +1.0 -1.2 -4.0 Acknowledgements Cathy Lascara and Russ Burgett of Old Dominion University kindly helped us with the IDL code for retrieving and visualizing the Scanfish data. Authors JX and RRH were supported by NASA grant NAG5-6286 and National Science Foundation grant OCE96-28888. Author SYC was supported by the National Science Foundation under grant DEB-94-12113. 60 References Anthes, D. L. T., Data assimilation and initialization of hurricane-predicting models, J. Atmos. Sci., 31, 702-719, 1974. Benjamin, T. B., Gravity currents and related phenomena, J. Fluid Mech., 31, 209- 248, 1968. Blumberg, A. F., and G. L. Mellor, A description of a three-dimensional coastal ocean circulation model, in Three-dimensional coastal ocean models, Coastal Estuary Sci., Vol. 4, edited by N.S. Heaps, 1-16, AGU, Washington, D. C., 1987. Bogden, P. S., P. Malanotte-Rizzoli, and R. Signell, Open-ocean boundary conditions from interior data: Local and remote forcing of Massachusetts Bay, J. Geophys. Res., 101, 6487-6500, 1996. Bogden, P. S., and J. O?Donnell, Generalized inverse with shipboard current measurements: Tides and nontidal flows in Long Island Sound, J. Mar. Res., 56, 995-1027, 1998. Edinger, J. E., D. K. Brady, and J. C. Geyer, Heat exchange and transport in the environment, Report 14, EPRI Publication No. 74-049-00-3, prepared for electric power research institute, Palo Alto, CA, 1974. Ezer, T., and G, K, Mellor, Continuous assimilation of Geosat altimeter data into a three-dimensional primitive equation Gulf Stream model, J. Phys. Oceanogr., 24, 832-847, 1994. Forbes, C., and O. Brown, Assimilation of sea surface height data into an isopycnic ocean model, J. Phys. Oceanogr., 26, 1189-1213, 1996. 61 Goodrich, D. M., W. C. Boicourt, P. Hamilton and D. W. Pritchard, Wind-induced destratificaiton in Chesapeake Bay, J. Phys. Oceanogr., 17, 2232-2240, 1987. Haltiner, G. J., and R. T. Williams, Numerical Prediction and Dynamic Meteorology, 477pp, John Wiley, New York, 1980. Hood, R. R., H. V. Wang, J. E. Purcell, E. D. Houde, and L. W. Harding Jr., Modeling particles and pelagic organisms in Chesapeake Bay: Convergent features control plankton distributions, J. Geophys. Res., 104, 1223-1243, 1999. Ippen, A. T., Tidal dynamics in Estuaries, part I: Estuaries of rectangular section, in Estuary and coastline hydrodynamics, edited by A. T. Ippen, 744pp, 1966. Johnson, B. H., K. W. Kim, R. H. Heath, H. L. Butler and B. B. Hsieh, Development and verification of a three-dimensional numerical hydrodynamic, salinity, and temperature model of the Chesapeake Bay, Technical Report HL-91-7, US Army Engineer Waterways Experiment Station, Vicksburg, MS, 1991. Launder, B. E., and D. B. Spalding, The numerical calculation of turbulence flows, Computer methods in applied mechanics and engineering, vol. 3, 269-289, 1974. Levitus, S., Climatological atlas of the World ocean, NOAA Professional Paper, No. 13, US Government Printing Office, Washington, DC, 173pp, 1982. Kundu, P. K., A numerical investigation of mixed layer dynamics, J. Phys. Oceanogr.,10, 220-236,1980. Sarmiento, J. L., and K. Bryan, An ocean transport model for the North Atlantic, J. Geophys. Res., 87, 394-408, 1982. 62 Sheng, Y. P., A three-dimensional mathematical model of coastal, estuarine and lake currents using a boundary fitted grid, Rep. No. 585, ARAR Group of Titan Systems, Princeton, N. J., 1986. Spitz, Y. H., and J. M. Klinck, Estimate of bottom and surface stress during a spring- neap tide cycle by dynamical assimilation of tide gauge observations in the Chesapeake Bay, J. Geophys. Res-Oceans, 103, C6, 12761-12782, 1998. Wang, D.-P., Subtidal sea level variations in the Chesapeake Bay and relation to atmospheric forcing, J. Phys. Oceanogr., 9, 413-421, 1979a. Wang, D.-P., Wind-driven circulation in the Chesapeake Bay, winter 1975, J. Phys. Oceanogr., 9, 564-572, 1979b. Wang, H. V., and R. S. Chapman, Application of vertical turbulence closure schemes in the Chesapeake Bay circulation model ? A comparative study, ASCE Estuarine and coastal modeling, 283-297, 1995. Wu, C. R., P. T. Shaw and S. Y. Chao, Assimilating altimetric data into a South China Sea model, J. Geophys. Res.-Oceans 104, C12, 29987-30005, 1999. 63 Figures Fig. 2.1 Bathymetry of the Chesapeake Bay and adjacent coastal area. Major tributaries are marked. Depth scales are in the unit of meters. 64 65 66 67 68 69 70 71 72 73 74 75 76 Fig. 2.14 Biweekly averaged flow fields without data assimilation [left panels (a) and (c)] and corresponding changes caused by data assimilation [right panels (b) and (d)]. Time averaging is over the second and third weeks of November 1995. Panels (a) and (b) are surface features, while (c) and (d) are taken at 10m below mean water level. Zonal length scales and zonal velocities are stretched for clarity. 77 Figure 2.15. Biweekly averaged section of density (? t ) anomalies induced by the data assimilation in October along the main axis of the deep channel. Time averaging is from day 295 to day 310.55 in panel (a), and from day 310.55 to day 326.1 in panel (b). The longitudinal section is the same as in Figs. 2.9 and 2.10. The density anomaly is obtained by subtracting the model result without data assimilation from that with data assimilation. Solid and dashed contours correspond to density deficit and density surplus, respectively. Contours intervals are ?? t =0.2. 78 Fig. 2.16 As in Fig. 2.11 except T is reduced threefold to 2 hours and K is increased threefold to (5 hr) -1 . 79 80 Chapter 3: A Simple Empirical Optical Model for Simulating Light Attenuation Variability in a Partially Mixed Estuary Abstract The representation of the submarine light field is a crucial component of pelagic ecosystem and water quality models. Modeling the light field in estuaries is a particularly complicated problem due to the significant influence of high concentrations of dissolved and particulate matter that are derived from both terrestrial and estuarine sources. The goal of this study was to develop a relatively simple but effective way to model light attenuation variability in a turbid estuary (Chesapeake Bay, USA) in a coupled physical-biological model. In this effort we adopted a simple, non-spectral empirical approach. Surface water quality data (salinity was used as a proxy of CDOM) and light measurements from the Chesapeake Bay Program were used to determine the absorption coefficients in a linear attenuation model using regression methods. This model predicts K c (specific attenuation due to phytoplankton/chlorophyll), K t (specific attenuation due to total suspended solids) and K s (a function of specific attenuation coefficients of CDOM in relation to salinity). The bay-wide fitted relation between light attenuation coefficient and water quality concentrations gives generally good estimates of total light attenuation, K d . However, the direct inclusion of salinity in the relationship has one disadvantage: it can yield negative values for K d at high salinities when applied in a numerical model. We therefore developed two separate models for two different 81 salinity regimes. This approach, in addition to solving the negative K d problem, also accounts for some changes in specific light absorption by chlorophyll, seston (non- phytoplankton particulate matter) and CDOM that apparently occur in different salinity regimes in Chesapeake Bay. The resulting model predicts the statistical characteristics (i.e., the mean and variance) of K d quite accurately in most part of Chesapeake Bay. We also discuss in this paper the feasibility and caveats of using K d converted from Secchi depth in the empirical method. 1. Introduction The intensity and spectral composition of light in aquatic systems change greatly with depth. These changes arise from the absorption and scattering by water and substances that are either suspended or dissolved. As a result, except in very shallow systems, light tends to limit primary production in deep water, but the depth at which this limitation occurs varies tremendously depending upon the concentrations of ?chromophoric? (optically active) dissolved and particulate matter in the water. Thus, the degree of light limitation, and therefore rates of primary production, in aquatic systems are a strong function of these constituents as well. Light availability also influences many other biological and chemical processes, including, among other things, species composition (Rijstenbil 1987; Jones and Gowen 1990), behavior of organisms (e.g., Gal et al. 1999; Graham et al. 2001; Dieguez and Gilbert 2003), phytoplankton physiological response (Cullen and Lewis 1987) and photochemical degradation (Anning et al. 2000). Therefore, reproducing the light field variability is a key problem in modeling biogeochemical processes in aquatic ecosystems. 82 Light intensity diminishes approximately exponentially with depth in water, so that I(Z) = I(0) EXP(-K d Z) (1) where I(Z) is the downward irradiance at depth Z; I(0) is the downward irradiance just beneath the air-water interface and K d is the vertical light attenuation coefficient. However, to determine precisely the amount of light available to phytoplankton, the spectral distribution of underwater light is needed because of the differential light absorption by water and chromophoric matter. If the quality of light is taken into consideration, then Eq. (1) can be expressed as: I(Z,?) = I(0,?) EXP(-K d (?)Z) (2a) and I(Z) = ? I(Z,?) d? (2b) Given that I(0) or I(0,?) can usually be directly measured or estimated at the sea surface, the major challenge is how to model K d or K d (?) variability in the water. For monochromatic light the vertical light attenuation can be decomposed as a set of partial attenuation coefficients, each characterizing absorption and scattering by a different waterborne material. Taking all wavelengths into consideration, spectral bio-optical models have been developed to calculate each partial attenuation coefficient for each narrow band of the spectrum. These models have, among other things, been applied to the interpretation of remote sensing data on ocean color (Platt and Sathyendranath1988; Sathyendranath and Platt 1989a) and development of numerical models for primary production (Smith et al. 1989; Sathyendranath and Platt 1989b; Arrigo and Sullivan 1994). 83 In a spectral light model, it is often necessary to know the wavelength-specific absorption and scattering coefficients for each waterborne material; this is not a trivial task in a complex and variable estuarine environment like Chesapeake Bay. Therefore, a simple approach expressed in eq. (1) is appealing, especially when the goal is to predict diffuse attenuation from numerically simulated concentrations of chromophoric substances in the water. Strictly speaking, a complete spectrum is still needed to determine the average K d for the whole photosynthetic waveband. However, as a simplification it is commonly assumed that the average K d for photosynthetically active radiation (PAR) can be decomposed as a set of partial attenuation coefficients in the same manner as for the monochromatic light. To further simplify the problem, K d (PAR) is often modeled as a linear function of water quality concentrations (Smith 1982; Stefan et al. 1983). Phytoplankton (and, if any, macrophytes), seston and chromophoric dissolved organic matter (CDOM) contribute to light attenuation (Kirk, 1994). Consequently, K d (PAR) can be approximated as: K d (PAR) = K w + K p [PHY] + K s [SES] +K o [CDOM] (3) where K w is the attenuation due to water, and K p , K s , and K o are the specific attenuation coefficients of phytoplankton, seston and CDOM, respectively. Estimation of K p , K s and K o is quite challenging in estuaries. The estuarine environment is more optically complex and variable than either open ocean or coastal waters due to the confluence of river water and sea water, leading to a broad suite of optically active constituents from both terrestrial and aquatic sources. The temporal and spatial distributions and compositions of phytoplankton, seston and CDOM vary 84 considerably. Theoretically, the corresponding specific coefficients, K p , K s and K o, will vary in time and space as well. Chesapeake Bay is a large partially mixed estuary in the United States. There are 50 major rivers discharging into the Bay. Among all tributaries, the Susquehanna River at the head is the primary source of freshwater as well as nutrients, dissolved organic matter and sediments. Generally speaking, concentrations of these chromophoric constituents are high in the upper Bay and low in the lower Bay. These gradients are mainly affected by the freshwater input. Typical of mid-latitude rivers, the discharge is high in spring followed by low to moderate flow throughout the rest of the year. The seasonal and spatial change in nutrient concentrations and turbidity greatly affect the light attenuation in the Bay (Harding 1994). A complex spectral light model has been developed and applied to Chesapeake Bay in studies related to restoring the Submerged Aquatic Vegetation (SAV) habitat (Gallegos et al. 1990; Gallegos 2001). In this paper we describe our efforts to develop a simple and suitable light model for calculating light penetration in a coupled physical-biological model of Chesapeake Bay. We use an empirical approach to estimate the specific attenuation coefficients for chlorophyll, total suspended solids (TSS), and CDOM (using salinity as a proxy). Specifically, water quality data obtained from Chesapeake Bay Program are used in a linear regression model to obtain a relation between K d and water quality concentrations, chlorophyll, TSS and salinity. We show that this method yields a simple bay-wide optical model which reproduces the observed K d variability remarkably well, and can be easily adopted in a numerical biogeochemical model. 85 2. Derivation of empirical light model 2.1 Empirical linear light model derived using direct light measurements The empirical, non-spectral approach is not a new method, having been used to study and model the relation between light attenuation and water quality concentrations in a number of different marine systems (e.g., Smith 1982; McMahon et al. 1992; Wang et al. 1996, Gallegos and Moore 2000). The basic premise behind this approach is that if one can simultaneously measure both K d and the concentrations of the optical constituents that determine K d , then multiple linear regression methods can be used to ?back out? the values of the specific attenuation coefficients. The regression method also has the advantage that it is simple and will generate coefficients that are specific to the particular water body from where the data were obtained. The 1995 and 1996 data from 70 stations maintained by Chesapeake Bay Program in the main stem of Chesapeake Bay and its tributaries (Fig. 3.1) were used to develop a set of empirical models for Chesapeake Bay (www.chesapeakebay.net/data/). At each station, underwater light intensity was measured every 0.25 m or 0.5 m, with the surface measurement at 0.1m. In addition, salinity was measured every 1 or 2 m; chlorophyll, dissolved organic carbon and total suspended solids (TSS) were usually measured at the surface (0.5 m or 1m) , above the pycnocline, under the pycnocline and at bottom. Each station was visited 12-20 times a year at somewhat irregular intervals. For the development of the empirical models presented in this paper, TSS data were used as a proxy for seston in equation 3. There was no direct measurement of CDOM and DOC cannot be used as a proxy 86 for CDOM because the proportion of CDOM in DOC varies considerably (Rochelle- Newall and Fisher 2002; Siegel et al. 2002). However, it has been shown that CDOM behaves conservatively like salinity in Chesapeake Bay (Rochelle-Newall and Fisher 2002) as in some other estuarine and coastal systems (Monahan and Pybus 1978; Bowers et al. 2004), with the primary source deriving from terrestrial freshwater inputs. As a result, CDOM varies inversely with salinity and salinity data can be used as a proxy for CDOM in equation 3. The linear representation of K d can then be written as K d (PAR) = K w ? + K c [CHL] + K t [TSS] + K s [Sal] (4) where K c and K t are the specific attenuation coefficients of chlorophyll and total suspended solid, respectively. K s is a function of specific attenuation coefficients of CDOM in relation to salinity, and K w ? is attenuation due to pure water and CDOM in fresh water. The combination of these latter two terms characterizes everything that impacts K d (PAR) except the effect of chlorophyll and total suspended solids. Including salinity in the equation makes the physical interpretation of the first and fourth term a little awkward. However, there are at least two advantages: firstly, the slope and intercept of the inverse relationship between CDOM and salinity are basically determined by the CDOM concentrations in fresh water. By using salinity directly in the equation we avoid the uncertainty and variability of the slope and intercept in the relation. Secondly, CDOM is usually not an explicit compartment in coupled physical-biological models while salinity is carried universally. Therefore, equation 4 is more readily applied in biogeochemical models for estimating K d . The light attenuation coefficient was calculated from Eq. (1) as: 87 )( )/ln( 12 21 ZZ II K d ? = (5) Where I 1 and I 2 are measured underwater light intensity at depth Z 1 and Z 2 , respectively. Z 1 was taken to be depth closest to the surface, usually at 0.1m and Z 2 the depth of 1.5m. If the water is shallower than 1.5m or the measurement at 1.5m was missing then a measurement closest to 1.5m was used. Using the entire database the regression of this calculated K d against measured surface chlorophyll, total suspended solids and salinity yields: K d = 1.932 - 0.004765[CHL] + 0.059[TSS] - 0.08667[Sal] (6) where [CHL] and [TSS] have the units of mg/m 3 (?g/l) and g/m 3 (mg/l), respectively. The value and probability for each coefficient are listed in Table 3.1. Also listed in Table 3.1 is the variability explained by each variable from a stepwise statistical model. TSS is by far the most important factor in controlling light attenuation in Chesapeake Bay. Alone it explains about 58% of the total variability in K d . CDOM variation (expressed by salinity here) is the second most important, which explains an additional 14% of K d variability. Chlorophyll enters the model only at 5% level and does not improve the R 2 . Thus, unlike oceanic waters, phytoplankton absorption plays a minor role in controlling the light field in Chesapeake Bay, though it is shown Table 3.1. The coefficient, P value and partial R 2 for each term in equation 6. Variables K w ? K t K s K c Value 1.932 0.059 -0.0867 -0.004765 Pr <0.001 <0.001 <0.001 0.0275 Partial R 2 NA 0.585 0.138 0.001 88 below that the importance of phytoplankton absorption increases down estuary. The negative sign in front of salinity is expected because CDOM is inversely related to salinity. However, the small negative but significant (at 0.05 level) coefficient for chlorophyll is somewhat puzzling given that increases in chlorophyll concentration should lead to increases in light attenuation. Gallegos and Moore (2000) thought that a negative coefficient for chlorophyll was a dubious result and so and dismissed these instances in their discussion. However, we suspect that this statistically significant coefficient reflects that fact that in the tributaries and in the upper bay (see next section) light often limits phytoplankton growth. That is, even though chlorophyll contributes to light attenuation, the light control of phytoplankton growth is so strong in some areas due to the influence of TSS and CDOM that it results in an inverse relation between K d and chlorophyll concentration. When K d is big, i.e. very turbid water, chlorophyll concentration is relatively low because phytoplankton cannot grow and vice versa. However, since phytoplankton is also component of TSS its net contribution to light attenuation can still be positive. The fit between K d obtained from eq. (6) (K d _predicated) against the calculated K d (K d _observed) appears to be linear (Fig. 3.2) and the R 2 = 0.72 is remarkably high, i.e., the empirical model explains about 72% of the observed K d variability. The variability that is not accounted for can, at least to some degree, be attributed to the fact that the composition of optically active constituents in the Bay changes considerably in time and space. Light absorption and scattering by phytoplankton changes with the species composition, pigment composition (Stuart et al. 1998), cell shape and size (Ciotti et al. 2002, Lorenz et al. 2003), and light 89 scattering also depends on the composition, geometric shape and size of TSS (Baker et al. 2001; Richardson 1987; Risovic 2002). Given the demonstrated importance of TSS in controlling K d , we conclude that the latter is responsible for most of the unresolved K d variability. 2.2 Empirical light model derived using Secchi Depth Secchi depth (SD) is routinely measured as a simple means of assessing water clarity in estuarine, coastal and open ocean waters. In the Chesapeake Bay Program, it was measured much more often (along with water quality measurements) than the direct light measurements that we used in our analysis above. For 1995 and 1996, 3428 SD measurements were made from 129 stations. Therefore, it is tempting to use attenuation coefficients derived from these measurements for our optical model development. However, the relationship between SD and attenuation coefficient is not fixed, i.e., it can vary by as much as sevenfold in waters with large variations in CDOM and turbidity (Koenings and Edmundson 1991). One must, therefore, be very careful when converting SD to K d in estuarine applications or highly inaccurate values may be obtained. Indeed, our attempts to derive an optical model using SD- derived K d values were unsuccessful. The following analysis illustrates the problem: There are a total of 1345 simultaneous measures of SD and K d in 1995 and 1996. A hyperbolic fit of the dataset gives us 01856.06.0 1 + = SD d Z K (R 2 = 0.78) (7) In contrast, a conversion obtained from Choptank River 90 1.0 97.0 ? = SD d Z K (8) gives a similar fit to the dataset (R 2 = 0.72) if data points of Z SD =0.1 are discarded, but Fig. 3.3 shows that Eq. (8) tends to give higher K d values when SD is less than 0.3 and lower K d values when SD is greater than 0.3. Also shown in Fig. 3.3 is the widely used conversion: SD d Z K 7.1 = (9) Using these three different relations to convert SD to K d we repeated the regression of this derived K d to chlorophyll, TSS and salinity: K d = 2.69 + 0.005[CHL] + 0.024[TSS] ? 0.108[SAL] (10.1) K d = 1.37 + 0.004[CHL] + 0.061[TSS] ? 0.082[SAL] (10.2) K d = 2.90 + 0.005[CHL] + 0.031[TSS] ? 0.122[SAL] (10.3) All three of these optical models describe much less of the observed Kd variability (R 2 = 0.56, R 2 = 0.55, R 2 = 0.55, respectively) compared to the model derived using the direct light measurements (Eq. 6). Moreover, they all tend to underestimate high K d and overestimate low K d , thus creating significant biases at both extremes (Fig. 3.4). Because all three models give similar fit of the dataset, only K d obtained from Eq. 10.1 vs. K d calculated from Eq. 7 is illustrated. The severity of these biases suggests that K d values derived from SD data cannot be used for optical model development. 91 2.3 Dual empirical linear light models with direct light measurement The uniform empirical linear light model described by Eq. 6 produces a reasonably good fit of the data. However, there is one significant problem: K d can go negative in high salinity regions. To surmount this problem and account for some geographic variability in specific attenuation coefficients as well, we divided the data into two groups based on salinity: one group of data with salinity less and equal to 15 and the other group with salinity greater than 15. If only the first and fourth terms in Eq. 6 are considered, only salinity value greater than about 22 could result in negative K d . The dividing criterion of salinity was chosen to obtain a statistically significant relationship for high salinity regime. Because the data are fewer and more scattered in the high salinity regime, a lower salinity value than the one could give negative K d is used here. The regression model (Eq. 4) leads to the following two empirical relations: K d = 1.80 - 0.0044[CHL] + 0.0673[TSS] - 0.096[Sal] (S ? 15, n=785) (11.1) K d = 1.17 + 0.024[CHL] + 0.006[TSS] - 0.0225[Sal] (S >15, n=563) (11.2) where S is salinity and n the number of available data points. The overall R 2 is 0.78 and thus is not dramatically improved compared to Eq. 6. To test the significance of the change in regression coefficients between 11.1 and 11.2, we used an indicator variable (A), or dummy variable, based on the grouping criterion-salinity (A=0 for S ? 15 and A=1 for S > 15) (for more detailed description of the method, please refer to Weisberg (1985) or Rosner (1995)). A new model is built as K d = K w ? + K c [CHL] + K t [TSS] + K s [SAL]+ 92 K w ? 1 A + K c1 [CHL?A] + K t1 [TSS?A] +K s1 [SAL*A] (12) Therefore, when A = 0, Eq. 11 collapses to 11.1 and when A = 1, it is equivalent to 11.2. Furthermore, the changes in K w ? , K c , K t and K s from 11.1 to 11.2 are given by K w ? 1 , K c1 , K t1 and K s1. Therefore the significances of changes in K w ? , K c , K t and K s in 11.1 and 11.2 are given by the significances of K w ? 1 , K c1 , K t1 and K s1 , respectively. The value and probability for each term are listed in Table 3.2. All the coefficients in eqs. 11.1 and 11.2 are significantly different (P<0.01 for intercept Kw? and P < 0.001 for the others, Table 3.2). The smaller intercept and the smaller coefficient for salinity in 11.2 are expected. They show, respectively, that in the lower Bay, the higher salinity water has less CDOM and CDOM has less influence on light attenuation. The decrease of Kt from low to high salinity in Eqs. (11.1-11.2) is not intuitive. Theoretically, the scattering process of light by suspended solids is size-dependent with smaller particles more optically active. Because bigger particles sink faster to the bottom the mineral particles in the water column get finer down the Bay. Following this reasoning, one would expect lower TSS but higher Kt as salinity increases. On the other hand, the lower Kt in higher salinity water agrees with our general knowledge that TSS has much less influence on attenuation in coastal and open ocean waters. Moreover, the weaker relationship between TSS and Table 3.2. The coefficient and P value for each term in equation 12. Variable K w ? K c K t K s K w1 ? K c1 K t1 K s1 Value 1.80 - 0.0044 0.0673 -0.096 -0.63 0.0284 - 0.0609 0.0735 Pr <0.001 0.0368 <0.001 <0.001 0.007 <0.001 <0.001 <0.001 93 K d is expected from scrutinizing the data, i.e., the values for TSS in the lower bay are highly variable, but this variability does not appear to be reflected in the K d variability. One other interesting change in the coefficients between Eqs. 11.1 and 11.2 is that the specific light attenuation coefficient for chlorophyll changes sign from negative to positive going from the low to high salinity regime. We interpret this change as reflecting a change from light controlling phytoplankton biomass in the turbid waters of the upper Bay (negative sign), to phytoplankton biomass controlling light in the clearer waters of the lower Bay (positive sign). We also speculate that the magnitude of K c in Eq. 11.2 is reduced compared to literature values because of these opposing effects. We see this same change when treating 1995 and 1996 data separately (analysis not shown). 1996 was a very high flow year while 1995 was a slightly below average flow year. As a result, the water was much more turbid in 1996 and phytoplankton growth was light limited over a considerably larger portion of the Bay. As we can infer from argument above we indeed have an inverse relationship between K d and chlorophyll in 1996, which is manifest in the lower Bay as well as the upper Bay. This leads to a general problem in the empirical method. In environment, such as estuaries and turbid inland water, where light could be a limiting factor in phytoplankton growth, field data will contain information from the two competing effects: Firstly, light limits phytoplankton growth. Hence phytoplankton grows more in clearer (lower K d ) conditions. Secondly, more phytoplankton biomass absorbs more light and attenuates light more quickly (higher K d ). Therefore, the empirical 94 light model obtained from field data will tend to underestimate the net contribution of phytoplankton to light attenuation (the second effect alone). Another problem revealed by the division of the data into two different salinity regimes is that the model explains a much smaller fraction of the optical variability in the higher salinity regime (R 2 = 0.15). It explains even less when a higher salinity criterion is chosen. Thus, it appears that in the higher salinity water, something other than chlorophyll, CDOM and TSS is substantially controlling variability in light attenuation. Alternatively, it is possible that the light controlling dependences change considerably in time and space so that the bulk relationship, fitted over many months and a large expanse of the lower Bay, fails. Fitting each month?s data in the high salinity group shows that there is a quite big range for each regression coefficient. Despite the problems discussed above, the empirical light model expressed in Eqs. 11.1 and 11.2 performs remarkably well over the entire Bay. Figure 3.5 compares K d derived from direct light measurements with K d estimated and from Eqs. 11.1 and 11.2 at all available main stem stations in 1995. The model-estimated K d s compare favorably with measurements and the model also captures the high K d variability at the upper bay and low variability in the mid bay. However, as we discussed above, the light model cannot explain a large part of the variability in K d in the lower bay (i.e., the high salinity region). 95 3. The Role of each component in light attenuation To estimate the role of each component in total K d we calculated the percentage contribution of each component using relations 11.1-11.2 for each data point. To simplify the interpretation we combined the first (K w ?) and fourth (K s [Sal]) terms, which yields the contributions from water itself and CDOM. We also separated phytoplankton from other suspended solids in the third term (K t [TSS ]) and added it to the second term (K c [CHL]) to give the total contribution of phytoplankton. Figure 3.6 shows the contribution of water + CDOM (K w ?+ K s [Sal]) , phytoplankton (K c [CHL]+K t [Phy]) and seston (K t [TSS-Phy]) to total K d . The phytoplankton in weight in the TSS term is derived from chlorophyll data assuming the ratio of chlorophyll to carbon is 1:50 and total organic weight is about twice the weight of carbon. We multiplied the carbon weight by 2.5 to get the total phytoplankton weight to account for the inorganic material. In the low salinity region (S ? 15), light attenuations by seston and water + CDOM are equally important. The contribution from phytoplankton is almost negligible. However, it does not necessarily mean that there is little phytoplankton accumulation in the water column or they don?t contribute to light attenuation. As we discussed above, the empirical model will tend to underestimate the net contribution of phytoplankton due to the competing effects of light attenuation by phytoplankton with light control of phytoplankton growth in turbid water. In the high salinity region (S > 15), phytoplankton plays a bigger role in attenuating light. The averaged contribution is below 20% but can be up to 50% at times. Light attenuation due to water and CDOM dominates. Because CDOM concentration is, presumably, much lower than in the low salinity region, a large part 96 of this term has to be attributed to water itself, with relatively small fluctuations in CDOM accounting for a large fraction of the K d variability. The contribution of seston becomes much smaller. Recall, however, that in the high salinity region the model explains only 15% of the observed K d variability. Thus, it is possible that all three of these sources of optical variability are small compared to some other unknown chromophoric constituents or some nonlinear effects. Overall, these results lead us to conclude that the variability of light attenuation due to phytoplankton in this partially mixed estuary is small compared to that due to seston and CDOM. In estuaries where CDOM behaves conservatively CDOM distribution can be adequately represented by salinity. Therefore, in modeling the light field in such environments, the first order importance is to reproduce both the mean distribution and variability of TSS. In contrast, in the open ocean and coastal waters, sediment loading and resuspension are negligible. Hence the role of TSS in light attenuation is insignificant. 4. Summary and conclusions This paper describes an effort to develop a relatively simple optical model for estimating the diffuse light attenuation coefficient (K d ) from variations in the concentrations of optically active constituents (i.e., chlorophyll-a, TSS and CDOM) in a partially mixed estuary. The estuarine environment is optically complex and variable, which poses a potential challenge in reproducing the underwater light field with a simple light model. Nonetheless, we have demonstrated that an empirical linear model derived 97 using in situ observations of K d , Chl-a ,TSS and salinity (as a proxy of CDOM) works remarkably well, i.e, the simple light models derived here generate reasonable K d values across the estuarine gradient. One caveat that arises in an optical model derived from using K d and Chl-a, TSS and salinity data for the entire Chesapeake Bay is that it can give rise to negative K d values in high salinity regions. Therefore, we divided the data into two groups by salinity (one for S ? 15 and one for S > 15). This approach resolved the negative K d problem in our model and it further demonstrated how the role of different optically active constituents can change over a wide range of salinity in an estuarine environment. In particular, models developed for the two different salinity regimes show that the specific absorption coefficient for chlorophyll changes sign (becomes positive) in high salinity waters. We believe this indicates that there are two competing factors controlling the relationship between K d and Chl-a. In turbid waters where constituents other than phytoplankton strongly influence K d , light controls phytoplankton growth and biomass, which will tend to give rise to a negative correlation between K d and Chl-a. In contrast, in clearer waters where phytoplankton growth and biomass are controlled by factors other than light (i.e., nutrients), chlorophyll strongly influences K d which will tend to give rise to a positive correlation between K d and Chl-a. Our empirical optical model derived using data from the entire Chesapeake Bay reveals that the former effect tends to dominate, i.e., the specific absorption coefficient for Chl-a is negative because light attenuation is strongly controlled by CDOM and seston. However, because the field data always 98 contain information from these two competing factors, we have to take precautions when analyzing the contribution of each optically active constituent to total light attenuation: An empirically obtained model will tend to underestimate the effect of Chl-a variability on K d and vice versa. We also show that Secchi depth-derived K d values cannot be used in the derivation of an empirical light model for Chesapeake Bay, i.e., when we try to use SD-derived K d values in the regression analysis the resulting optical model does a poor job of predicting K d and has strong biases at both high and low K d values. These problems arise because the equations used to convert SD to K d cannot describe all of the K d variability and they introduce biases at extreme K d values. Furthermore, SD cannot be used in the derivation of an optical model directly because the linear regression assumes that K d is linearly related to the water quality concentrations whereas SD is not linearly related to K d . Finally, the analysis of each optical constituent?s contribution to total light attenuation shows that in modeling the light field in estuaries, the first order importance is to reproduce both the mean distribution and variability of TSS and CDOM because they are often the dominant determinants of K d variability. In estuaries where CDOM behaves conservatively, CDOM concentration can be adequately represented by salinity. Therefore, in such systems, to correctly model the underwater light field the main problem is reproducing TSS variability. Because the light attenuation properties of phytoplankton, TSS and CDOM vary in different environments, the specific attenuation coefficients obtained in this study will tend to be locality-specific. None-the-less, we have demonstrated that the simple optical 99 model derived using linear regression yields reasonably good results in a variable and complex estuarine system. Acknowledgements We would like to thank Dr. Tom Fisher for helpful discussions and Ms. Yi Huang for consultation on statistics. This research was supported by funds from NASA, the National Science Foundation and Horn Point Laboratory?s Education Committee. References Anning,T., H. L. MacIntyre, S. M. Pratt, P. J. Sammes, S. Gibb and R. J. Geider. 2000. Photoacclimation in the marine diatom Skeletonema costatum. Limnology and Oceanography 45(8): 1807-1817. Arrigo, K. R. and C. W. Sullivan. 1994. A high resolution bio-optical model of microalgal growth: tests using sea-ice algal community time-series data. Limnology and Oceanography, 39 (3): 609-631. Baker, E. T., D. A. Tennant, R. A. Feely, G. T. Lebon and S. L. Walker. 2001. Field and laboratory studies on the effect of particle size and composition on optical backscattering measurements in hydrothermal plumes. Deep-Sea Research Part I 48 (2): 593-604. Bowers, D. G., and D. Evans, D. N. Thomas, K. Ellis, P. J. le B. Williams. 2004. Interpreting the colour of an estuary. Estuarine, Coastal and Shelf Science 59: 13-20. 100 Ciotti, A. M., M. R. Lewis and J. J. Cullen. 2002. Assessment of the relationships between dominant cell size in natural phytoplankton communities and the spectral shape of the absorption coefficient. Limnology and Oceanography 47 (2): 404-417. Cullen, J. J. and M. R. Lewis. 1988. The kinetics of algal photoadaptation in the context of vertical mixing. Journal Of Plankton Research 10 (5): 1039-1063. Dieguez, M.C. and J. J. Gilbert. 2003. Predation by Buenoa Macrotibialis (Insecta, Hemiptera) on zooplankton: effect of light on selection and consumption of prey. Journal Of Plankton Research 25 (7): 759-769. Ferrari, G. M, and M.D. Dowell. 1998. CDOM absorption characteristics with relation to fluorescence and salinity in coastal areas of the southern Baltic Sea. Estuarine, Coastal and Shelf Science 47: 91-105. Gal, G., E. R. Loew, L. G. Rudstam and A. M. Mohammadian. 1999. Light and diel vertical migration: spectral sensitivity and light avoidance by Mysis relicta. Canadian Journal Of Fisheries And Aquatic Sciences 56(2): 311-322. Gallegos, C. L. 2001. Calculating optical water quality targets to restore and protect submerged aquatic vegetation: overcoming problems in partitioning the diffuse attenuation coefficient for photosynthetically active radiation. Estuaries 24: 381-397. Gallegos, C. L., D. L. Correll, and J. W. Pierce. 1990. Modeling spectral diffuse attenuation, absorption, and scattering coefficients in a turbid estuary. Limnology and Oceanography 35: 1486-1502. 101 Gallegos, C. L., and K. A. Moore. 2002. Factors contributing to water-column light attenuation, p. 16-27. In R. A. Vatiuk, P. Bergstrom, W. M. Kemp, E. Koch, L. Murray, J. C. Stevenson, R. Bartleson, V. Carter, N. B. Rybicki, J. M. Landwehr, C. Ballegos, L. Karrh, M. Naylor, D. Wilcox, K. A. Moore, S. Ailstock, and M. Teichberg (eds.), Chesapeake Bay Submerged Aquatic Vegetation Water Quality and Habitat-based Requirements and Restoration Targets: A Second Technical Synthsis. U.S. Environmental Protection Agency, Chesapeake Bay Program, Annapolis, Maryland. Graham, W. M., F. Pages and W. M. Hammer. 2001. A physical context for gelatinous zooplankton aggregations: a review. Hydrobiologia 451: 199-212. Harding, L. W. 1994. Long-term trends in the distribution of phytoplankton in Chesapeake Bay: roles of light, nutrients and streamflow. Marine Ecology Progress Series 104: 267-291. Jones, K. J. and R. J. Gowen. 1990. Influence of stratification and irradiance regime on summer phytoplankton composition in coastal and shelf seas of the British- isles. Estuarine, Coastal and Shelf Science 30: 557-567. Keith, D. J., J. A. Yoder, and S. A. Freeman. 2002. Spatial and temporal distribution of coloured dissolved organic matter (CDOM) in Narragansett Bay, Rhode Island: implications for phytoplankton in coastal waters. Estuarine, Coastal and Shelf Science 55: 705-717. Kirk, J.T.O. 1994. Light and Photosynthesis in Aquatic ecosystems. Cambridge University Press, Cambridge. 102 Koenings J. P., and J. A. Edmundson. 1991. Secchi disk and photometer estimates of light regimes in Alaskan Lakes ? effects of yellow color and turbidity. Limnology and Oceanography 36: 91-105. Lohrenz, S. E., A. D. Weidemann and M. Tuel. 2003. Phytoplankton spectral absorption as influenced by community size structure and pigment composition. Journal of Plankton Research 25 (1): 35-61. McMahon, T. G., R. C. T. Raine, T. Fast, and J. W. Patching. 1992. Phytoplankton biomass, light attenuation and mixing in the Shannon Estuary, Ireland. Journal of Marine Biological Association of the United Kingdom 72 (3): 709-720. Monahan, E.C., Pybus, M. J. 1978. Colour, ultraviolet absorbance and salinity of the surface waters off the west coast of Ireland. Nature 274: 782-784. Platt, T., and S. Sathyendranath. 1988. Ocean primary production: estimation by remote sensing at local and regional scales. Science 241: 1613-1620. Richardson, M. J. 1987. Particle-size, light-scattering and composition of suspended particulate matter in the north-atlantic. Deep-Sea Research Part A 34 (8): 1301-1329. Rijstenbil, J. W. 1987. Phytoplankton composition of stagnant and tidal ecosystems in relation to salinity, nutrients, light and turbulence. Netherlands Journal Of Sea Research 21: 113-123. Risovic, D. 2002. Effect of suspended particulate-size distribution on the backscattering ratio in the remote sensing of seawater. Applied Optics 41 (33): 7092-7101. 103 Rochelle-Newall, E. J. and T. R. Fisher. 2002. Chromophoric dissolved organic matter and dissolved organic carbon in Chesapeake Bay. Marine Chemistry 77: 23-41. Rosner, B. 1995. Fundamentals of Biostatistics, 4th ed. Duxbury Press. Sathyendranath, S., and T. Platt. 1989a. Remote sensing of oceanic primary production: computations using a spectral model. Deep-Sea Research 36: 431- 453. Sathyendranath, S., and T. Platt, 1989b. Computation of aquatic primary production: extended formalism to include effects of angular and spectral distribution of light. Limnology and Oceanography 34: 188-198. Siegel, D. A., S. Maritorena, N. B. Nelson, D. A. Hansell and M. Lorenzi-Kayser. 2002. Global distribution and dynamics of colored dissolved and detrital organic materials. Journal of Geophysical Research-Oceans 107 (C12), Art. No. 3228. Smith, Jr., W.O. 1982. The relative importance of chlorophyll, dissolved and particulate material, and seawater to the vertical extinction of light. Estuarine, Coastal and Shelf Science 15: 459-465. Smith, R.C., B.B. Prezelin, R. R. Bidigare, and K. S. Baker. 1989. Bio-optical modeling of photosynthetic production in coastal waters. Limnology and Oceanography 34: 1524-1544. Stefan, H.G., J.J. Cardoni, F. R. Schiebe, and C.M. Cooper. 1983. Model of light penetration in a turbid lake. Water Resources Research 19: 109-120. 104 Stuart, V., S. Sathyendranath, T. Platt, H. Maass and B. D. Irwin. 1998. Pigments and species composition of natural phytoplankton populations: effects on the absorption spectra. Journal of Plankton Research 20 (2): 187-217. Wang, M., D. R. Lyzenga and V. V. Klemas. 1996. Measurement of optical properties in the Delaware Estuary. Journal of Coastal Research 12 (1): 211- 228. Weisberg, S. 1985. Applied Linear Regression, 2 nd ed. John Wiley & Sons, Inc. 105 Figures 106 107 108 109 110 111 Chapter 4: Modeling Biogeochemical Cycles in Chesapeake Bay with a Coupled Physical-Biological Model Abstract In this paper we describe the development and validation of a relatively simple biogeochemical model in Chesapeake Bay. This model (which is coupled to a three- dimensional hydrodynamic model of the Bay) was adapted from a NPZD-type open- ocean ecosystem model, which was then modified by adding additional compartments and parameterizations of biogeochemical processes that are important in estuarine systems. These modifications include an empirical optical model for predicting K d , compartments for representing oxygen and suspended sediment concentrations, and parameterizations of phosphorus limitation, denitrification, and seasonal changes in ecosystem structure and temperature effects. To show the overall performance of the coupled physical-biological model, the modeled dissolved inorganic nitrogen, phytoplankton, dissolved oxygen, total suspended solids and light attenuation coefficient in 1995 (a dry year) and 1996 (a very wet year) are examined and compared to observations obtained from the Chesapeake Bay Program. We demonstrate that this relatively simple model is capable of producing the general distribution of each field (both the mean and variability) in the mainstem of the Bay. And the model is robust enough to generate reasonable results under both wet and dry conditions. Some significant discrepancies are also observed, such as overestimation of phytoplankton concentrations in shoal regions and overestimation oxygen 112 concentrations in deep channels, which reveal some deficiencies in the model formulation. Some potential improvements and remedies are suggested. Sensitivity studies on selected parameters are also reported. 1. Introduction In recent decades, coupled physical-biological models have been widely applied to the marine environment to simulate both the physical and biogeochemical processes and study the interactions between them, especially the effect of physical factors on biological communities. The complexity of the physical models ranges from box (Li et al., 2000) and 1-D models (e.g., Doney et al., 1996; Hood et al., 2001; Marra and Ho, 1993) to fully 3-D hydrodynamic models (e.g., Lima and Doney, 2004; Skogen et al., 1995). The biological models range from simple NPZ (nutrient, phytoplankton, zooplankton) (e.g., McClain et al., 1996) or NPZD (nutrient, phytoplankton, zooplankton, detritus) models (e.g., Doney et al., 1996; Oschlies and Garcon, 1999; Hood et al., 2003) to multi-nutrient, multi-species and size-structured ecosystem models (e.g., Moore et al., 2002; Lima and Doney, 2004). When such models are applied to estuarine and coastal waters they can provide a means of assessing the potential impacts of local management strategies and hence provide useful information to decision-makers. Chesapeake Bay is the largest and most productive estuary in the United States. Similar to other estuarine systems (e.g., Lapointe and Clark, 1992; Pitkanen et al., 1993; Nagy et al., 2002), Chesapeake Bay has been suffering from degradation of water quality due to increased environmental stresses (Carpenter et al., 1969; Malone 113 1992). Eutrophication in Chesapeake Bay has caused serious economic, aesthetic and ecological problems: harmful algae blooms (Bowers et al., 2000), loss of submerged aquatic vegetation (SAV) (Orth and Moore, 1983), hypoxia and anoxia at deep waters in summer (Cooper and Brush, 1991; Kemp et al., 1992), among other things. And increased load of suspended solids directly reduces water clarity and when they deposit at the bottom they can have a detrimental impact on benthic organisms and production (Airoldi, 2003; Miller et al., 2002). Efforts have been made to reduce the N and P inputs from point and non-point sources and land-based sediment runoff with the goal of restoring the Bay to conditions observed in the early 1950s (Chesapeake Bay Agreement 1983, 1987, 2000). Numerical models have been used as a key analytic tool to provide guidelines in setting goals of nutrient and sediment reduction to achieve water quality standards. The Chesapeake Bay Program has developed a modeling system that is a state-of-the-art package of models that has been expanded and refined over more than a decade through the combined efforts of both scientists and managers. In an effort to model the complexity of the real world this package includes an airshed model (Regional Acid Deposition Model (RADM)) (Chang et al., 1990; Dennis, 1996), a watershed model (Hydrological Simulation Program-Fortran (HSPF)) (Bicknell et al., 1996; Greene and Linker, 1998), a hydrodynamic model (WES-CH3D) (Johnson et al., 1991; Hood et al., 1999; Sheng, 1986; Xu et al., 2002) and a water quality model (CE-QUAL-ICM) (Cerco and Cole, 1994; Cerco and Noel, 2004) coupled with a sediment (DiToro and Fitzpatrick, 1993) and living resources (including SAV and benthos) model (Madden and Kemp, 1996; Wetzel and Neckles, 1986). This 114 modeling system, which has been developed explicitly for management applications, is extremely complicated. The water quality model in this package alone has 24 state variables including two physical variables (temperature and salinity), multiple algal groups, two zooplankton groups, and multiple forms of nitrogen, phosphorus, carbon and silica. And there has been tendency toward including more and more complexity in an effort to simulate all of the potentially relevant biogeochemical processes. However, recent studies have demonstrated that more complexity in an ecosystem model does not necessarily improve model performance (Denman, 2003; Hood et al., 2003; Friedrichs et al., 2004). Friedrichs et al. (2004) has shown that a simple NPZD model can reproduce as much of the observed variability as more complicated models in an open ocean system, and that more complex model formulations can lead to reduced predictive ability if they are not adequately constrained with data. In addition, simple models have many advantages in terms of identifying the most important processes and parameters that drive the observed variability. It is not clear, however, whether or not these results, which were derived from an open ocean model intercomparison, are applicable in a complex estuarine system like Chesapeake Bay. Nitrogen, silica and iron are the major limiting nutrients in the open ocean. In estuaries, iron is not likely to be an important limiting nutrient due to the close proximity of terrestrial Fe sources. Rather, in estuaries nitrogen, silica and phosphorus limit phytoplankton growth, with the latter becoming particularly important during periods of high freshwater runoff (Fisher et al., 1992). Variations in freshwater flow can therefore lead to seasonal and regional shifts in these limiting factors (e.g., D?Elia et al. 1986; Fisher et al., 1992). 115 Benthic processes play a far more important role in biogeochemical cycling in estuarine systems due to the closer proximity of the bottom (Boynton and Kemp, 1985; Boynton et al., 1995; Seitzinger, 1988). Under different conditions, sediments can be either a sink or source of nutrients. Sediment regeneration of phosphate and ammonium can provide a significant portion of phytoplankton N and P requirement (Fisher et al. 1982; Boynton and Kemp, 1985; Malone et al., 1988). The coupled nitrification and denitrification process represents an important pathway for removing nitrogen from the system (Boynton et al., 1995). These benthic influences are particularly important in coastal plain estuaries like Chesapeake Bay which are very shallow. Another important difference between open-ocean and estuarine systems is the influence of suspended sediments on light transmission in the water column. High sediment loads in estuaries, which are also associated with periods of high freshwater flow, can lead to very rapid light attenuation in estuarine waters which limits primary production (see Xu et al., 2004 and references therein). Seasonal and interannual variability in river flow into the Chesapeake Bay is extremely large, with annual flow varying between about 20 ? 60 x 10 9 m 3 yr -1 (Harding, 1994). Because nutrient (and sediment) loads vary in direct proportion to flow, so does stimulation of phytoplankton growth, resulting in large seasonal and interannual variations in primary production (ranging from ~ 200 ? 600 gC m -2 yr -1 , Harding et al., 2002) and oxygen demand. Variations in river flow also impact sediment load / light and stratification which, in turn, controls the resupply of oxygen to bottom waters and regenerated nutrient in deep water to euphotic zone. 116 The species composition of phytoplankton in the Bay also shifts seasonally. The classic view is that in spring the diatom bloom accounts for the annual biomass peak but in summer flagellates and dinoflagellates make up most of the phytoplankton population. Accordingly, the food web of the system is dominated by diatom-mesozooplankton in spring but microbially in summer (Malone et al., 1988, Malone et al., 1991). Another unique phenomenon observed in Chesapeake Bay is the accumulation of phytoplankton biomass at depth in the deep channels of the mainstem Bay and tributaries (Chesapeake Bay Program database: www.chesapeakebay.net/data/ and our Fig. 4.9 as an example). In winter and spring, maximum chlorophyll concentrations are often observed at the bottom. In the open ocean, deep chlorophyll maximum are routinely observed near the bottom of the euphotic zone in association with the nutricline due to a combination of photoadaption and phytoplankton growth (e.g., Gundersen et al., 1998; Varela et al., 1992; Venrick, 1991), but this is a very different mechanism from the accumulation of phytoplankton biomass at the bottom in Chesapeake Bay. In the Bay this accumulation of chlorophyll happens well below the euphotic zone and appears to be associated with enhanced phytoplankton sinking and deposition. It has been hypothesized that high concentrations of phytoplankton and TSS combined with strong vertical and horizontal T-S gradients promote the flocculation and sinking of phytoplankton in Chesapeake Bay (and in estuarine systems in general, W. M. Kemp, personal communication). The shallow depth and two-layered estuarine circulation help to retain the biomass in the system so that the subsequent recycling of regenerated nitrogen into 117 the euphotic zone supports high phytoplankton production in summer (Malone et al. 1988). Therefore, in addition to correctly modeling the physical processes that deliver nutrient from the bottom to euphotic zone, it is also important to reproduce the bottom accumulation of phytoplankton biomass in spring in order to have a reasonable phytoplankton production and deep oxygen demand in summer. Clearly, we face numerous additional challenges when we attempt to model biogeochemical cycles in an estuarine system like Chesapeake Bay. But does this necessarily mean that we must employ a vastly more complicated model, like the CBP modeling system (which may have reduced predictive skill) in order to reproduce the observed variability? Or is it feasible to use a more simplified model and parameterize the impacts of these additional complexities? In this study, we adapted a simple NPZD-type biological model, that was originally developed for the open ocean, and coupled it with a numerical hydrodynamic model of Chesapeake Bay, with the goal of capturing the first order biogeochemical variability in the system. Using this coupled model we explore the possibility of using a simple ecosystem model to simulate the complex estuarine environment and simple ways to parameterize the important processes, such as phosphorus limitation, suspended sediment effects upon light penetration, benthic biogeochemical impacts and enhanced sinking and recycling of phytoplankton biomass. Observations are used to evaluated the model performance and analyze existing problems in the model. 118 2. Model Description The combined model consists of an eight-compartment biogeochemical model coupled to a numerical hydrodynamic model of the Chesapeake Bay. In the following subsections we describe these two models and how they are linked, as well as the forcing, and boundary conditions. Because the physical model has been discussed in detail previously (Xu et al., 2002), here we will focus on the biogeochemical part of the model. 2.1 The physical model The physical model, originally developed by Sheng (1986), was subsequently modified extensively by the US Waterways Experiment Station (WES) for application to Chesapeake Bay (Johnson et al., 1991; Wang and Chapman, 1995). Under Boussinesq and hydrostatic approximations, the hydrodynamic model solves for salinity, temperature, water-level elevation and velocities in three dimensions. There are up to 19 layers in the vertical with a uniform layer thickness of 1.52 m, except that the top layer thickness fluctuates with sea level. Horizontally, the governing equations in the Cartesian coordinate system are recast in a boundary-fitted curvilinear coordinate system (see Xu et al., 2002, their figure 4) to cope with the irregular shoreline configuration and deep channel orientation. The model domain extends offshore to include a piece of coastal ocean with coarse resolution. The coastal ocean is included mainly as a buffer zone to facilitate free exchanges across the Bay mouth. Inside the Bay, typical grid size ranges from 1 to 5 km in the main stem of the Bay. Similar to a host of primitive-equation models such as Blumberg and Mellor (1987), the staggered Arakawa-C grid system is used in both horizontal and 119 vertical directions of the computation domain. The vertical eddy viscosity and diffusivity are computed from mean flow and stratification characteristics using the second-order k-? turbulent closure scheme (see, for example, Kundu, 1980; Launder and Spalding, 1974). A quadratic stress is exerted at the bottom, assuming the bottom boundary layer is logarithmic over a bottom roughness height of 0.05 cm. Coefficients of horizontal eddy viscosity and diffusivity are set to 10 4 cm 2 /s. The model is forced by open ocean tides, winds, freshwater inflows and the heat exchange at the water surface. Further, salinity and temperature fields were also prescribed on offshore open boundaries using monthly Levitus climatology data (Levitus, 1982) combined with field data at Duck, North Carolina (36.1833?N, 75.7467?W) acquired daily (with occasional lapses) by the Field Research Facility of the US Army Corps of Engineers. Daily freshwater inflow with zero salinity and time-varying temperature was prescribed for the eight major Chesapeake Bay tributaries. Additional details about the physical model implementation and forcing can be found in Xu et al. (2002). In spite of the coarse resolution of the model, essential circulation features, such as the two-layered circulation in the main channel and major tributaries, and reasonable temperature and salinity structures can be reproduced by the model (Johnson et al., 1991; Hood et al., 1999; Xu et al., 2002). 2.2 The biogeochemical model The biological model (Fig. 4.1) was adapted for application in Chesapeake Bay from an open ocean model described in Hood et al. (2001). It has 6 compartments and a total of 8 state variables. The three nutrient pools include 120 dissolved inorganic nitrogen (DIN), dissolved inorganic phosphorus (DIP) and dissolved organic nitrogen (DON). The model also contains one compartment for phytoplankton (P), one compartment for heterotrophs (H) and one detrital compartment (D). The H compartment is considered to represent all heterotrophic groups including bacteria, microzooplankton, mesozooplankton etc. And in this model the heterotrophic processes are microbially dominated. Two additional state variables are carried in this model: dissolved oxygen (DO) and inorganic suspended solid (ISS). Dissolved oxygen is included to simulate anoxia and hypoxia in the Bay and it also serves as a natural trigger to slow down the respiratory processes of the heterotrophs under hypoxic and anoxic conditions. The inorganic suspended solid does not participate in the biological cycles. However, its existence has a great effect on the light attenuation, which in turn modifies the phytoplankton growth. The symbols and values of all model parameters are list in Table 4.1. Changes in H due to biological processes are determined as follows: DONHCgenhHHChgehDHCgedhPHCgeph t H mdonmhmdmp +?++= ? ? )1( (1) Where ??=??=??=??= /,/,/,/ hhdondonddpp hhhh And shdondp HKHDONDP +?+?+?+?=? , The HK s , the half-saturation constant, is assumed here to be the same for all substrates. The DIN compartment represents the sum of all forms of dissolved inorganic nitrogen: NO 3 , NO 2 and NH 4 + . Changes in DIN, due to biological processes, are determined by the heterotrophic remineralization of particulate and dissolved organic nitrogen to dissolved inorganic nitrogen and the uptake of DIN by phytoplankton: 121 PUDONHChgenaen HHChgehaehDHChgedaedPHChgepaep t DIN pmdon mhmdmp ??+ ?+?+?= ? ? )( )()()( ??? (2) Where U p = ),min()1( / / PN II II p KDIP DIP KDIN DIN ee k ++ ? ? ? ? ? . The phytoplankton DIN uptake is mediated by light and nutrients (N, P) availability. The equation for phytoplankton growth and mortality is: PHChPSPU t P mppP ??= ? ? ? (3) Similar to DIN, changes in DON, due to biological processes, have corresponding components of the heterotrophic remineralization in addition to contribution to the DON pool from phytoplankton due to direct exudation and natural mortality and consumption of DON by heterotrophs (i.e., bacteria): DONHChPSPUHHChgehaeh DHChgedaedPHChgepaep t DON mdonppmh mdmp ??+????+ ??+??= ? ? )1()1())(1( ))(1())(1( ??? ?? (4) The modeled biological processes for the detritus pool include the contributions due to egestion by heterotrophs and natural mortality of phytoplankton and removal due to consumption by heterotrophs: DHChPSDONHChaen HHChaehDHChaedPHChaep t D mdpmdon mhmdmp ?+?+ ?+?+?= ? ? ?)1( )1()1()1( (5) DO is not represented in the modeled biological processes directly. The changes in DO are based on the oxygen production in DIN uptake by phytoplankton due to photosynthesis and oxygen demand in respiratory processes using the photosynthetic quotient for oxygen and nitrogen (pqn). 122 pqnPUDONHChgenaen HHChgehaehDHChgedaedPHChgepaep t DO pmdon mhmdmp ])( )()()([ ??+ ?+?+?= ? ? ??? (6) DIP is not dynamically modeled in this study. Its parameterization will be discussed in the next section. In addition, there is a sinking term for phytoplankton, detritus and ISS in our model. The change of these variables due to sinking is represented by z CW t C s ? ? = ? ? (7) where C is the concentration of phytoplankton, detritus or ISS, W s is the corresponding sinking rate and z is the vertical coordinate. This model runs ?on line? in the hydrodynamic model. To save computing time, the biogeochemical model currently has a time step of 1 hour while the physical model has a time resolution of 5 minutes. 2.3 Boundary conditions for the biogeochemical model 2.3.1 Data availability Biological data from stations at the major tributaries were downloaded from Chesapeake Bay Program (CBP). At each station, DIN, DON, Chlorophyll, and TSS were measured once or twice every month and they were linearly interpreted in the model. Heterotrophic biomass was assumed to be ? of the phytoplankton in N unit. TSS was assumed to include everything particulate in the model: phytoplankton, heterotrophs, detritus and inorganic suspended solids (ISS). To estimate the detritus Table 4.1: Model parameters. Description Symbol Value Units 123 Growth efficiency for H on P Gep 0.2 Growth efficiency for H on DON Gen 0.2 Growth efficiency for H on D Ged 0.2 Growth efficiency for H on H Geh 0.2 Assimilation efficiency for H on P aep 0.7 Assimilation efficiency for H on DON aen 1.0 Assimilation efficiency for H on D aed 0.7 Assimilation efficiency for H on H aeh 0.7 Partition of P senesence ? 0.25 Partition of P production ? 0.7 Partition of excretion to DIN ? 0.75 Maximum phytoplankton growth rate ? p 0.96 (Sal ? 3) 3.22 (Sal > 3) d -1 Phytoplankton light saturation paramenter I k 40.0 W m - 2 Phytoplankton photoinhibition paramenter I ? 400.0 W m - 2 Half-sat. const. for DIN uptake by P K N 0.5 ?M Half-sat. const. for DIP uptake by P K P 0.015 ?M Phytoplankton natural mortality rate S p 0.01 d -1 Heterotrophic maximum consumption rate C m 0.8 (T ? 10 ?C) 6.4 (T > 10 ?C) d -1 Half-sat. const. for heterotrophic consumption HK s 0.8 ?M Heterotrophic preference for P ? p 0.3 (T ? 10 ?C) 0.1 (T > 10 ?C) Heterotrophic preference for D ? d 0.2 (T ? 10 ?C) 0.3 (T > 10 ?C) Heterotrophic preference for H ? h 0.3 Heterotrophic preference for DON ? don 0.2 (T ? 10 ?C) 0.3 (T > 10 ?C) Stoichiometric O:N ratio pqn 8.625 Stoichiometric P:N ratio pnr 0.0625 Detritus sinking rate w d 2 m d -1 Phytoplankton sinking rate w p 2.5 (Jan. - May) 1 (Jun. ? Dec.) m d -1 DIP concentration DIP 0.1 (Jan. ? May) 1.0 (Jun. ? Dec.) ?M Dinitrification loss rate at bottom layer R dnf 0.01 124 and ISS from TSS and chlorophyll data it was assumed that ISS contributes 50% of TSS. 2.3.2 Initial conditions and boundary conditions at tributaries and Open-Ocean All biological variables except DIN are assigned uniform values for the whole Bay at values that roughly approximate wintertime concentrations (Table 4.2). The along bay gradient of DIN concentration was estimated by fitting data from mainstem stations in winter to a power function and the concentration was hold constant vertically and laterally. Because the loadings and biological processes rapidly change and dominate the distributions of biogeochemical properties, the initial values have little impact on the solution. The amount of nutrients entering an estuarine or coastal system determines, to a large extent, the biogeochemical cycling in the system. Consequently, great care must be taken in specifying nutrient loading when modeling such systems. However, processes such as ground water flow, stream-edge land runoff, nutrient retention/release in marshes and wetland and atmospheric dry deposition are still very difficult to fully represent in numerical models. To obtain a correct estimate of the nutrient loading to the Bay and account for the uncertain nutrient sinks/sources we Table 4.2: Initial values for variables in the biogeochemical model DIN DON DIP P H D DO ISS Max{2, 2.71*1.63**(0.072*I a )- 4.88} ?M 0.15 ?M 0.1 ?M 6.0 ?M 1.0 ?M 10 ?M 300 ?M 10 mg/l a: I is the along-bay grid index. 125 nudged the biological variables? values towards CBP observations in the upper reaches of the tributaries. At the upper reaches of the tributaries the equation for each biological variable (C) becomes: )(log CCprocessesicalbiodiffusion Dt DC obs ?++= ? (8) In addition to the advection, diffusion and biological production/consumption, the nudging term in (8) restores observed value of C at a fixed rate ?. To obtain Cobs in each grid cell, available data from stations in the upper reaches of the tributaries was linearly interpolated along the tributaries. The nudging is not strong (? = 24 hours -1 ) and it is applied only at the tidal fresh and oligohaline regions of the tributaries so that the vast majority of the estuarine system still functions dynamically. At the open-ocean boundary, seasonally averaged data for DIN, phytoplankton and DO were obtained from NOAA (nodc.noaa.gov). The data were linearly interpolated to each boundary grid cell and the boundary conditions change seasonally. When data were not available zero gradient boundary conditions were used. 2.3.3 Atmospheric wet and dry deposition of DIN & DON For simplicity, atmospheric wet and dry deposition of DIN and DON are considered to be uniform across the Bay. Monthly wet DIN deposition data were obtained from National Atmospheric Deposition Program (NADP) at station MD13 (- 76.1525, 38.9131). And wet DON deposition is estimated to be 20% of total wet 126 deposition (Meyers et al., 2000). Dry deposition of DIN has been estimated to attribute 30 to 63% of the total atmospheric nitrogen deposition in the eastern US (Levy and Moxim, 1987, Logan, 1983, Meyers et al. 2000). In this study we assume that the dry deposition is 50% of total deposition. 2.3.4 Air-sea exchange of O 2 In addition to advection and diffusion and biological production and consumption, oxygen also exchanges between air and water at the estuary surface. For the surface layer, the change of oxygen concentration induced by this process is modeled by: Hs D[O 2 w]/ Dt = K(PO 2 w ? PO 2 a) (9) where Hs is the depth of the surface layer and K is the exchange rate. And K = kL Where k is the piston velocity and L is the solubility. Then equation (9) can be expressed as: Hs D[O 2 w]/ Dt = k L (PO 2 w ? PO 2 a) = k ([O 2 w] ? LPO 2 a) (10) The piston velocity and solubility are computed by using equations from Wanninkhof (1992): k = 0.31 U 2 (Sc / 660) 1/2 (11) Sc(T) = 1638 ?81.83T + 1.483T 2 ?0.008004T 3 (12) Where U is the wind speed in m/s and Sc is the Schimit number and T is temperature. Solubility L is a function of temperature and salinity and calculated as in the table A2 in Wanninkhof (1992). 127 For simplicity, the O 2 concentration in the air is taken to be constant: PO 2 a = 0.23 atm. 2.3.5 Sediment deposition and resuspension There is no explicit sediment layer in our model. Instead, the bottom layer of the physical model is treated as a reservoir of sediment as well as part of the water column. Inorganic and organic matter is treated differently at the bottom. For inorganic suspended sediments (ISS), we assume a 100% deposition rate, i.e., all particles that hit the bottom will deposit to the bottom. Many studies have been done on sediment deposition rate and there are a wide range of values. When bottom stress exceeds some critical shear stress sediments are resuspended, which is modeled by: Resuspension flux = M(?-? c ) where M is the resuspension rate, ? is the shear stress and ? c is the critical shear stress. In contrast, the organic detritus (D) accumulates in the bottom layer and is regenerated by heterotrophs in the same way as in the water column with a feedback from oxygen concentration, i.e., when DO is lower than 60 mg/l, the respiration rate decreases by 25% (T. R. Fisher, personal communication). 2.4 Parameterization in the model Except for those discussed below and in the previous sections, all of the biogeochemical model parameters are set as described in Hood et al. [2001, see their Table 1]. 128 2.4.1 P-Limitation Phosphorus limits phytoplankton growth over a significant extent of the upper Bay during high flow conditions (Fisher et al., 1992; 1999). As discussed in the introduction, proper representation of the phosphorus cycle in the Chesapeake Bay requires explicit modeling of the phosphorus cycle, i.e., at the very least, each compartment of the model would also have to be expressed in phosphorus units and cycled using different rules and parameterizations. In order to avoid this level of complexity, and considering that P-limitation effects are significant only during springtime and that there are only small differences in P concentration in the surface waters (where it matters) over the the whole bay (Fisher et al., 1992) we took a very simple approach: We assign a uniform value for the P concentration over the whole Bay, and then specify how this value varies seasonally relative to a fixed half- saturation constant for phosphorus uptake. In so doing, we can invoke P-limitation manually by simply specifying a low P concentration during the time when P- limitation is expected and a high value for other time when it is not. The specific seasonal values that we set for the Bay wide P concentration are specified in Table 4.1. The value for half saturation coefficient for phosphorus uptake (K p ) was set a posteriori to invoke a moderate level of P-limitation in the spring with P concentration set to be 0.1?M (Fisher et al., 1992). Note that what matters in the model is not the actual values but the relative magnitude of these two terms. 129 2.4.2 Light attenuation To calculate the under water light field we use a simple non-spectral light attenuation model: I(Z) = I(0) EXP(-K d Z) (14) where I(Z) is the photosynthetically available radiation (PAR) below the surface of the water, I(0) is PAR at the surface, K d is the diffuse attenuation coefficient for PAR, and Z is depth. Spatial and temporal variability in K d is specified using a simple empirical optical model that was derived specifically for the Chesapeake Bay (Xu et al., submitted). In developing the empirical optical model, surface water quality data and direct light measurements from the CBP were divided into two groups by salinity. Each group of data was then fitted against a linear attenuation model to determine the specific attenuation coefficients for chlorophyll, TSS and salinity (used as a proxy for CDOM). The two empirical relations are: K d = 1.80 - 0.0044[CHL] + 0.0673[TSS] - 0.096[Sal] (S ? 15, n=785) (15a) Kd = 1.17 + 0.024[CHL] + 0.006[TSS] - 0.0225[Sal] (S > 15, n=563) (15b) It is shown in Xu et al. (submitted) that this model can explain more than 70% of the observed K d variability in Chesapeake Bay. The solar radiation at the water surface is calculated from daily meterological data: cloud coverage, dew-point temperature, air temperature and wind speed, from the Patuxent Naval station in the same way that it is used to force the physical model (see Xu et al., 2002 for details) and the variability in the solar radiation flux is considered to be uniform across the Bay. This computed incoming radiation includes both long and short wave radiations. To obtain PAR at the surface (I(0)) a conversion 130 factor of 0.47 (estimated from the spectral energy distribution of solar radiation at sea level) is used. Below the water surface, the averaged light intensity for each physical model layer is applied in the biological model. 2.4.3 Denitrification Denitrification is a very important process in coastal marine and estuarine systems, where N losses via denitrification may account for a significant portion of total N input from terrestrial sources and regeneration (Boynton et al., 1995; Seitzinger, 1988; Seitzinger and Giblin, 1996). In Chesapeake Bay, denitrification loss of N has been estimated to be 24 % of total N budget in Maryland mainstem and ranges from 13-79% at other study sites (Boynton et al., 1995). Indeed, without any representation of denitrification processes in the model, it tends to greatly overestimate DIN concentrations at depth. Denitrification rates are largely determined by the availability of NO 3 - . However, our model does not differentiate different types of DIN. Therefore, a simple loss term is added in each grid cell at the bottom, which is proportional to the DIN concentration in the bottom layer. The loss was then set a posteriori at a rate that gives reasonable bottom water DIN concentrations (Table 4.1). This denitrification loss term yields a annual DIN sink of about 60 ?mol/m 2 /h in 1995 and 100 ?mol/m 2 /h in 1996, which is comparable to the denitrification loss estimated from field data (Boynton et al., 1995). 131 2.4.4 Seasonal adjustment of parameters Zooplankton grazing and bacterial remineralization are temperature dependent, i.e., lower in winter and higher in summer. In order to account for this effect, we lower the maximum heterotrophic consumption rate (C m ) by a factor of 6 when temperature drops below 10 ? C (Table 4.1). The choice of the temperature criterion and the magnitude of the rate decrease is consistent with the temperature dependence of mesozooplankton grazing rate described by Huntley and Lopez (1992). Along with grazing rate we also changed the grazing preferences in different seasons (Table 4.1). Specifically, the grazing preferences were adjusted so that the heterotrophs prefer to graze more on phytoplankton and less on detritus and DON in the wintertime. This change in grazing preference is intended to crudely parameterize broadscale seasonal changes in foodweb structure that are known to occur in Chesapeake Bay, i.e., in summer and fall, the Bay is dominated by microbial consumers, while in winter and spring, mesozooplankon grazing is considered to be more important (Malone and Ducklow, 1990; Malone et al., 1991). The dominant primary producers also change seasonally in Chesapeake Bay. In spring, diatoms dominate phytoplankton production whereas in summer flagellates and dinoflagellates dominate (Malone and Ducklow, 1990; Malone et al., 1991). In order to parameterize the effect that these floristic shifts have upon export flux, we use a higher sinking rate for phytoplankton in later winter to spring (Table 4.1). 132 2.4.5 Regional adjustment of parameters Finally, in order to reproduce the relatively low CBP-measured phytoplankton biomass observed at the head of the Bay it was necessary to lower the net growth rate of phytoplankton in our model. To achieve this change we simply lowered the growth rate of phytoplankton in the fresh region of the Bay (Table 4.1). This adjustment was necessary in spite of the fact that the model reproduces the observed K d values quite accurately in these waters. Although this adjustment may seem somewhat arbitrary, it can be justified for a number of reasons. In the estuarine environment the confluence of river water and sea water produces a salinity gradient which ranges from 0 to more than 30 across the system. Thus, biogeochemical cycling is driven by a more of a fresh water ecosystem at the head of the estuary compared to the mouth. Obviously, the species composition of both phytoplankton and zooplankton communities will be very different in marine and freshwater systems, and may therefore require different phytoplankton growth and grazing rate parameters. We can also speculate that the transition from fresh to more saline water at the head of the Bay will result in a decline in the phytoplankton growth rate due to senescence in this transitional zone. Finally, it also likely that in the Susquehanna flats at the head of the Bay, benthic bivalves graze heavily on phytoplankton. In the model, different approaches can be taken to lower the net primary production, i.e., lower phytoplankton growth rate, increase the grazing rate, increase the sinking rate of phytoplankton or any combination which will decrease the production term and/or increase the loss terms in equation 3. Here we chose the simple solution of lowering phytoplankton growth rate (Table 4.1). 133 3. Results and Discussion In this section we discuss our main run solution and carry out some selected parameter sensitivity studies. In so doing we point out and discuss some of the major successes and deficiencies of our model. Wherever possible, we suggest potential reasons for deficiencies and possible means of correcting them. In order to provide a robust test of the overall performance of this coupled model system, we have forced the model over two contrasting years in terms of flow: 1995 and 1996. The year1995 was a below normal flow year, while 1996 was a very high flow year. Daily discharge rates from the major tributaries in both years are shown in Fig. 4.2. The averaged fresh water discharge in 1996 is 112,000 cfs, comparing to an average of 49,000 cfs in 1995. During the first major freshet event in 1996, the discharge rate reached up to 1,123,000 cfs. The high nutrients and TSS input coming with the high fresh water discharge in 1996 resulted in a very different biogeochemical response bay-widely. As a starting point, we first examine the seasonally (quarterly) averaged characteristics of the nutrient cycles and phytoplankton distributions generated by the model and compare them directly to observations collected by the CBP. Given the fact that we do not expect our model to be able to reproduce short time and space- scale variability, and the general paucity of observations that are available for comparison with the model, we consider this level of comparison to be an appropriate starting point. 134 3. 1 Main run results 3.1.1 Seasonal comparison with CBP data at main stem stations along a longitudinal transect In order to compare our model results quantitatively with Chesapeake Bay Program data, we sampled the model in the same manner as the observations and calculated the seasonal mean and variability at each station for both modeled and observed fields. We selected the stations along the main axis of the Bay to present the comparisons (Fig. 4.3). (For clarification we used the same scale for all seasons and all stations. Some error bars are out of the plotting area. We didn?t rescale the plot because we think they do not interfere with the presentation of the idea and different scales may be confusing and misleading.) 3.1.1.1 DIN and Chlorophyll Figure 4.4 shows the seasonally averaged surface DIN in both 1995 and 1996 for each quarter. The data points represent the mean and the error bars are the 95% confidence interval. Data points with no error bars mean there is only one observation available for that period. The modeled DIN mean field follows the observed pattern quite closely in space and time in both years. For the most part, the error bars indicate the modeled and observed fields are statistically indistinguishable. However, these error bars are very wide due to the long averaging time periods. In the analysis that follows we point out differences between the model and the observations even though in most cases these differences are not statistically significant. 135 It appears that the model has a tendency to underestimate DIN in the first quarter of both years. Also in 1996, the model has a tendency to overestimate DIN in the upper Bay in the second, third and fourth quarters, which corresponds to the underestimation of chlorophyll (Fig. 4.5, see below). Thus, in 1996, when freshwater and nutrient loads were particularly high, the phytoplankton in the model often did not take up enough DIN in the surface layer. DIN concentrations are highest in the spring and lowest in the summer in both years in the model and the observations. The available DIN in the surface water decrease abruptly going down the bay due to the phytoplankton consumption except in winter and spring of 1996 due to the extremely high fresh water discharge rate. During these two seasons DIN at the surface layer in the lower Bay can still exceed 20 ?M, enough to support substantial phytoplankton growth. The variability is generally highest in the upper Bay because the available DIN mainly varies with the fresh water discharge in that region. In the lower part of the Bay DIN is scarce in the surface water and it varies little except in winter and spring of 1996. Our modeled DIN variability follows the observed patterns of variability quite well. However, in some cases the modeled variability appears to substantially exceed that in the observations (e.g., second quarter of 1996 in the upper Bay). Figure 4.5 shows the seasonally averaged surface chlorophyll concentration in 1995 and 1996. Again, our modeled chlorophyll concentrations compare reasonably well with the observations. In 1995, the chlorophyll-a maxima in all seasons occur to the north of Bay Bridge and there is a well defined decreasing trend going down the Bay (Fig. 4.5a) in both the model and the observations. In contrast, in 1996 the 136 surface chlorophyll is higher throughout the Bay expect in winter as observed (Fig. 4.5b). In spring of 1996, the high surface chlorophyll extends from Choptank River to the Bay mouth, and the chlorophyll maximum actually occurs in the mid to lower Bay. At the very head of the Bay phytoplankton growth in the model is light limited due to high light attenuation from TSS (see section 3.1.1.2 below). In the summer of 1996, even though DIN in the surface layer is scarce the chlorophyll concentration is still high in both the observations and the model, presumably because phytoplankton growth is supported by the regenerated nutrients during this period (Malone et al., 1988). From these seasonal comparisons we conclude that our simple model not only reproduces the bay-wide patterns, but also successfully captures the seasonal and wet versus dry year variability in DIN and Chlorophyll concentrations in Chesapeake Bay. Our model?s tendency to sometimes underestimate chlorophyll at the head of the bay in spring, summer and fall, is partly the result of specifying a lower phytoplankton growth rate in the fresh regions of the Bay. Because salinity changes dynamically with river discharge and other physical conditions, the parameterization has different effects under different conditions and perhaps too much effect during some seasons. As we will discuss later, the phytoplankton biomass in the upper Bay may be more grazing controlled. If this grazing control comes from bivalves in the Susquehanna River flats, then it will tend to occur at fixed locations, which will improve our model?s simulation of chlorophyll. 3.1.1.2 K d and TSS 137 Proper representation of the light field and, in particular, the diffuse attenuation coefficient (K d ) variability, is crucial for modeling primary production in estuarine systems where light attenuation varies tremendously in both space and time. Reproducing K d variability is particularly difficult in turbid estuaries like Chesapeake Bay where suspended sediments play an important role in controlling K d (Xu et al., submitted). The seasonally averaged, observed and modeled surface light attenuation coefficient (K d ) is shown in Fig. 4.6, where the observed K d is calculated from direct light measurements obtained from the CBP. (Although the number of K d values that can be derived from direct light measurements in the Bay is relatively sparse compared to those that can be derived form secchi-derived measurements, we use the former because the conversion equations that must be applied to calculate K d from secchi depth produce systematic biases in the derived K d values (Xu et al., submitted)). Figure 4.6 shows that our simple empirical light model is able to reproduce the observed K d patterns reasonably well, reproducing the spatial, seasonal and year- to-year variations where the data allow comparison. Both the model and the observations show that the variability in K d tends to be high in the upper Bay and low in the mid-Bay. However, the variability in the modeled K d in the lower bay is much less than that of the observed. This happens in part because modeled TSS variability is generally lower than observed in the lower lower Bay and partly because the light model itself does not explain as much of the observed K d variability in the lower bay (see Xu et al., submitted). 138 Figure 4.7 shows the observed and modeled surface TSS for 1995 and 1996. These plots reveal that the modeled TSS means and variability compare fairly well with observations in the mid and upper Bay regions. However, our model always tends to underestimate both the mean and variability in the lower Bay. The high TSS values and large variability is a persistent feature in the observations in the lower bay. The failure of our model to simulate TSS in the lower Bay is almost certainly related in some way to the simplicity of our TSS representation: Modeled TSS is the sum of the four biological model compartments: phytoplankton, zooplankton, detritus plus inorganic suspended solids (ISS). In Chesapeake Bay, ISS can constitutes a large part of TSS and our parameterization of processes controlling ISS distribution is oversimplified, i.e., we use a uniform and constant sinking speed, resuspension rate and critical shear stress. Alternatively, in the lower bay, there may be wave-current interactions, which are more effective in resuspending the sediments. There may also be some lateral TSS loading to this region that is not represented in our model, which loads ISS only in the tributaries and does not include shoreline erosion. 3.1. 2 Synoptic comparison with CBP data at different bay areas By seasonally averaging the data at each station we sacrifice the temporal resolution in our model, which makes it impossible to judge how well the model reproduces finer-scale (sub-seasonal) temporal variability. In this section, we divide the Bay into three different regions: upper bay, mid bay and lower bay (Fig. 4.3). For each region, spatially averaged data are presented for each cruise to get a synoptic view, i.e, we now calculate means and standard deviations that represent spatial means and spatial variability (Fig. 4.8). However, because each station was not 139 sampled simultaneously, the sampling time for each region can span a time period covering as much as 2 days. As in the previous plots, the model results were sampled exactly the same way as the observations. In 1995, both the mean and variability of chlorophyll are high in upper Bay and decrease down the Bay (Fig. 4.8a). Our model reproduces the temporal evolution of both the mean and the variability quite well in mid to lower bay. Note that there is not much seasonality in the mid to lower Bay in 1995 in the observations and the model, which we attribute to persistent nutrient limitation in these regions in a low flow year. It appears that the relatively weak spring bloom in 1995 did not deposit enough organic nitrogen on the bottom to fuel a substantial summer increase in phytoplankton growth in the mid to lower Bay (Malone et al., 1988). In the upper bay the model does not reproduce the high observed chlorophyll concentrations in May and tends to overestimate chlorophyll concentration in September and October. In 1996, chlorophyll concentrations in the mid to lower bay are higher than in 1995 and our model generally captures this phenomenon, which we attribute to generally higher nutrient loading and more recycling from organic nitrogen deposited on the bottom. In 1996 in the upper bay, our model simulates chlorophyll concentration and seasonality quite well. Note that there is much more seasonality in 1996 compared to 1995, especially in the mid the lower bay which is also generally captured by the model, although point-to-point discrepancies are apparent. In 1996, there was a strong summer/fall chlorophyll increase in the upper bay and spring/summer increase in mid to lower bay. It appears that the upper Bay bloom in 140 1996 was delayed by high turbidity and low water clarity associated with high freshwater flow and ISS loading (Fig. 4.6 and 4.7). 3.1.3 Seasonal comparison of vertical profiles at selected stations So far we have only examined the model performance at the surface. In this section we show some seasonally averaged vertical profiles at selected stations: CB3.3C in the upper bay, CB5.3 in the mid bay and CB6.3 in the lower bay. The positions of the stations are shown in Fig. 4.3. The observations were linearly interpolated in the vertical direction to the modeled layers. At each depth, both the modeled and observed values are presented as means with 95% confidence intervals. Figure 4.9 shows the seasonally averaged chlorophyll distributions at the three stations. Generally, our modeled vertical profiles compare favorably with the observed profiles. Note that in the upper Bay (station CB3.3C, Fig. 4.9a) both the model and the observations reveal substantially higher chlorophyll concentrations near the bottom in the first quarter, but the reverse tends to be true in the third and fourth quarters. The second quarter appears to be transitional. At this station (Fig. 4.9a), the model reproduces the observed profiles quite well in winter and fall, but tends to underestimate the chlorophyll concentration in the upper layer in spring and summer. In contrast, in mid to lower bay (Fig. 4.9b and 4.9c), the model reproduces the observed chlorophyll profiles remarkably well in all but the winter season (first quarter) where, there is considerable underestimation of chlorophyll. We speculate that this may have something to do with our increased sinking speed in the first and second quarters and the fact that our sinking rate is uniform across the bay. It is 141 possible that, in reality, phytoplankton sink faster in the upper bay due to flocculation because both phytoplankton and TSS concentrations are higher. In general, Figure 4.9 shows that the model reproduces the broad scale seasonal and spatial changes in chlorophyll concentrations at depth as well as at the surface, and also the pronounced bottom enhancement of chlorophyll concentration that is sometimes observed. We emphasize here that this bottom enhancement was achieved in the model through the direct application of a sinking term to the phytoplankton compartment, which is not necessary in open ocean models (e.g., Hood et al., 2001; 2004). Indeed, as we discussed in the introduction, it appears that the combination of high chlorophyll and ISS concentrations, perhaps combined with strong spatial gradients in temperature and salinity, accelerate phytoplankton sinking and export to the bottom in productive and turbid estuaries like Chesapeake Bay. It is possible that even better results could be obtained through the application of some kind of aggregation model or parameterization which substantially enhances the rate of chlorophyll sinking under high chlorophyll and ISS concentrations. The seasonally averaged vertical profiles of DIN at the same three stations are shown in Fig. 4.10. Reproducing the vertical distribution of DIN is more problematic. In general, the model reproduces the observed profiles better in the second half of the year at all three locations. In the first quarter, at all three stations the observed DIN is high at the surface and decreases at depth. However, our model produces a reversed curve with lower DIN at surface and higher DIN at the bottom. We speculate that this problem may have something to do with physical rather than biological processes, i.e. perhaps DIN is not transported down the bay fast enough or there is not enough 142 vertical mixing in the physical model to replenish DIN in the upper layer from the lower layers. In the second quarter, the model generally overestimates DIN at depth. This may also be explained by insufficient vertical mixing. Alternatively, there could be too much reminerization at depth during this time period. As with chlorophyll, the model reproduces the observed DIN profiles better at the two lower Bay stations. In spite of these discrepancies, Fig. 4.10 shows that, as with chlorophyll, the model reproduces the broad scale seasonal and spatial changes in DIN concentrations at depth as well as at the surface. Vertical profiles of TSS are illustrated in Fig. 4.11. In our model TSS matters most where it can affect the light attenuation hence the phytoplankton growth. Therefore, our first goal is to correctly estimate the TSS concentration in the upper layer. Given the simplicity of our parameterization of TSS production and burial processes, the model does a remarkable job of reproducing the observed TSS profiles in the upper and mid Bay regions (Fig. 4.11a, b). The model does not, however, reproduce TSS profiles very well in the lower Bay (Fig. 4.11c), where it generally underestimates the observed concentrations. As discussed above, we suspect that this has something to do with lateral loading and/or a resuspension problem. Finally, in Fig. 4.12 we show DO concentration profiles at these same three stations (Fig. 4.12). DO in our model is cycled in proportion to DIN using a fixed ratio and is therefore coupled tightly with primary production and remineralization except at the surface layer where there is air-sea exchange. DO concentration profiles are reproduced very well in spring and fall and in all seasons in the lower bay. In winter, the model tends to underestimate DO at depth in the upper to mid bay. This 143 could result from the underestimation of primary production in winter or, perhaps more likely, insufficient wintertime ventilation/vertical diffusion in the physical model. The model consistently overestimates DO at depth in the summer (third quarter), especially in the deep channel (Figs. 4.12a,b). Thus, it appears that there is not enough organic matter to fuel remineralization at depth during this time. This idea is consistent with our speculation that we have too much remineralization in the spring, which gives rise to the high DIN concentration at depth in spring. Perhaps a more sophisticated thermal regulation of remineralization could improve these aspects of the solution, as opposed to the step function that is currently employed to capture temperature control of heterotrophic processes. Even though 1996 is very different year in terms of freshwater forcing and nutrient loading, we see similar levels of agreement in the vertical profiles between the modeled and observed fields (results not shown). 3.1.4 Seasonal spatial plots for chlorophyll and DIN in 1995 and 1996 The modeled surface DIN and chlorophyll concentration in the second (spring) and third (summer) quarters in 1995 and 1996 are shown in Fig. 4.13 and Fig. 4.14, respectively. DIN concentration is usually highest in spring in the upper Bay. In 1995 high DIN concentration was confined to north of Bay Bridge (Fig. 4.13a), whereas in 1996, due to the large amount of fresh water discharge and DIN loading, the high DIN concentrations extended well into the mid Bay (Fig. 4.14a). Phytoplankton biomass, shown here as chlorophyll a concentration, was higher in 1996 throughout the Bay (Figs. 4.13b and 4.14b). Along with high nutrient input also came high ISS loading. As a result, there was much more intense light attenuation in 144 the upper Bay in 1996 compared to1995 (Figs. 4.13c and 4.14c). Consequently, the chlorophyll a concentration maximum was further south of the Bay in 1996 (Figs. 4.13b and 4.14b). In summer, surface DIN concentrations are generally low (Figs. 4.13d and 4.14d). However, in 1996 there was still a lot of DIN available in the water column in the upper Bay (Fig. 4.14d), which appears to have contributed to an extensive phytoplankton bloom in the summer of 1996 (Fig. 4.14e) while in 1995 the bloom was less intense and more restricted to the upper reaches of the Bay. The overall picture in our model matches favorably with our understanding of Chesapeake Bay plankton dynamics, i.e., the idea that seasons and years with high flow and nutrient loading tend to produce more intense phytoplankton blooms that extend further down Bay compared to low flow /load years (Harding, 1994). One persisting feature in our model that is not entirely consistent with observations is the very low chlorophyll in regions of very high turbidity in the upper Bay (see section 3.1 above). Another striking feature in these figures is the enhanced phytoplankton growth near the shores, especially in spring (Figs. 4.13b and Fig. 4.14b). This enhancement is not observed in either satellite maps of near surface chlorophyll (L. Harding, personal communication). In the following section we take a closer look at this problem and suggest some possible solutions. 3.2 Lateral variations of chlorophyll distribution Theoretically, we would expect higher chlorophyll concentrations at shoal areas than in the deep main stem of the Bay. At the shallow flanks, the mixing depth will be shallower than the euphotic zone more often due to the shallow depth which will tend to relieve light limitation. Moreover, organic matter trapped by the bottom 145 and regenerated there by benthic organisms will tend to resupply nutrients directly to the euphotic zone. As a result of these kinds of processes, our model generates unrealistically high chlorophyll concentrations in the shoal areas along the shores, especially in spring (Figs. 4.13b and Fig. 4.14b). The fact that we see (and expect) this enhancement in our model, it suggests that we are either mis-calculating something like light attenuation rate over the shoals, or that some processes is missing from the model that prevents phytoplankton growth and/or chlorophyll accumulation in the shoal regions. Here we put forward a set of hypotheses that could explain this discrepancy and test them with our model. One possible mechanism that could prevent the development of enhanced chlorophyll concentrations in the shoals is lateral transport of chlorophyll and organic matter from the shoals to deeper water where production is light limited. i.e, there may not be sufficient lateral transport and export from the shoals in the model. We tested this idea by increasing the horizontal diffusion in the model by tenfold, which has the effect of transporting organic matter from regions of high concentration (shoals) to low concentration (deep channel). Increasing the horizontal diffusion in this manner results in a general smoothing out of chlorophyll gradients in the Bay, but it does not correct the problem (results not shown). Of course, this does not rule out the possibility that we are missing some form of horizontal transport that moves chlorophyll and/or organic matter into the mainstem. For example, sediment transport from the shoal areas towards the main stem may play an important role. Because we do not have a sediment layer explicitly in the model there is no simple way to further explore this possibility. 146 The high biomass of phytoplankton concentrations that develop near the shoal regions in the model is related to the application of direct sinking losses to the phytoplankton compartment, which was done to reproduce the elevated chlorophyll concentrations at the bottom in the deep channel. Phytoplankton sinking has different effect between shallow versus deep waters in the model. In deep water, it transports the biomass out of the euphotic zone to lower layers where the phytoplankton stops growing. In contrast, in shallow water even the phytoplankton at the bottom may still continue to grow. This difference leads directly to enhanced primary production in shoal regions in the model when phytoplankton sinking is invoked. When phytoplankton sinking is turned off, the enhancement disappears. One possible solution to this problem would be to assume that the mortality rate of phytoplankton at the bottom increases. When we invoked this assumption in the model we found that it has a significant effect on production in both the deep channels and the shoals and so does not substantially alleviate the problem (results not shown). Another possible solution is that there is, in fact, enhanced phytoplankton mortality at the bottom, but that it occurs only over the shoals, perhaps due to grazing or some other source of mortality that is only occurs in shallow water. Indeed, Malone et al. (1986) found that even though phytoplankton biomass was low in shoal waters along the eastern shore the growth rate was high, which indicates grazing control. We therefore implemented a simple bivalve filtration loss term at the bottom for waters no deeper than 5m. As a first attempt, we adopted a high bivalve concentration of 100g/m 2 with a filtration rate of 0.25 l/hr/g (Newell and Koch, in press) and the filtration was turned on for spring and summer to crudely reproduce 147 increased bivalve filtration that occurs in summer (Apr. to Sep.). This filtration does, indeed, remove most of the high chlorophyll concentration near the shores (Compare Fig. 4.15b with Fig. 4.14e). To quantitatively show the effect, we selected two pairs of stations: one in the main stem and one in the shoal (CB4.2C and EE2.1, CB5.3 and EE3.4, positions shown in Fig. 4.3), and at each station we seasonally averaged the observed surface chlorophyll and compared it to average chlorophyll concentrations from model runs with and without bivalve filtration (Fig. 4.16). The biggest effect can be seen in the mid to lower bay in spring (left panel in Fig. 4.16b) where we see a substantial lowering of chlorophyll in the model run with filtration, which compares much more favorably with observed concentrations. Due to the uniformly high bivalve concentration we adopted in the model the chlorophyll concentration at EE3.4 in the summer even went too low. Another interesting fact is that even though bivalve filtration only happens in shallow water it has a significant effect on chlorophyll concentrations in the main stem as well. From this suite of tests we conclude that some form of enhanced mortality of phytoplankton may be occurring in shoal waters in Chesapeake Bay that prevents accumulation of phytoplankton biomass where it otherwise would have a tendency to increase relative to the mainstem Bay. One possibility is that this mortality is due to some form of benthic grazing that occurs in shallow waters, such as bivalve filtration. However, we cannot rule out the possibility that lateral transport from the shoals into the deep channel, perhaps associated with sediment transport, is occurring that is not properly represented in the model. 148 3.3 Sensitivity studies In this section we briefly describe how the model response to adjustments in selected parameters that are involved specifying P-limitation, light response, grazing and export/detrital sinking. Most of the sensitivity results discussed below were derived from parameter adjustments applied to our 1995 model run because 1995 represent the more typical year in terms of flow and loading. Plots are not shown. In its current configuration the model is surprisingly insensitive to small perturbations to its parameter values, which seems to indicate that the seasonal and spatial patterns that we see in the model (and by analogy the Chesapeake Bay) arise as a result of the strong influence of the freshwater forcing and nutrient and ISS loading. 3.3.1 Sensitivity to P-limitation Many publications have pointed out that phosphorus is important limiting nutrient in estuarine systems and particularly in the spring, especially in a wet year, in Chesapeake Bay (e.g., D?Elia et al. 1986; Fisher et al., 1992; Glibert et al., 1995). Adjusting the degree of P limitation in a model like ours provides an excellent means of illustrating biogeochemical effects. We explored the impacts of P-limitation by adjusting the value of the P concentration (Table 4.1) in summer and fall to a value that is much greater than Kp (Table 4.1) to turn the P-limitation off. As expected, without P-limitation the model produces a bigger spring bloom in both 1995 and 1996. The seasonally averaged chlorophyll concentration increases about 50% in upper to mid bay in spring without P-limitation. P-limitation needs to 149 be invoked in both years in order to get the correct bloom intensity and down Bay extent. It seems, therefore, that P-limitation was occurring during the springtime during both years even though they were very different in terms of flow and nutrient loading. However, a larger extend of the bay appears to be P-limited in 1996. In 1995 P-limitation extends from Baltimore Harbor to the mouth of Potomac River while in 1996 it extends from the Bay bridge to south of Rappahannock River. In addition, without P-limitation the seasonally averaged chlorophyll concentration in winter can be twice as high in the mid Bay in 1996. In contrast, P-limitation has little effect in 1995. The P-limitation effect is most pronounced in upper to mid bay. In the far upper reaches at the head of Bay, invoking P-limitation has little or no effect because in these waters light is the most important controlling factor. In the lower bay, DIN becomes depleted and is the primary limiting factor and so the impacts of P-limitation tend to be reduced there as well. 3.3.2 Sensitivity to I k Unlike many other model parameters, the 1995 solution is fairly sensitive to the half saturation parameter for light (I k ) because this parameter directly dictates how deep phytoplankton growth can occur in the water column. This, in turn, has a direct impact on the DIN distribution in the model. In winter and spring, the effects of changing I k can extend from the head of the Bay to south of Potomac River, which indicates that this large part of the Bay is light, rather than nutrient, limited due to the low surface irradiance, high turbidity, deeper winter mixing and high nutrient loading at these times of year. In summer and fall, the effects of changing I k are confined to 150 the upper Bay because most of the middle and lower Bay switch to nutrient limitation. Adjustments to I k show very clearly where the region of transition from light to nutrient limitation occurs in the Bay, i.e, it occurs where the impacts of adjusting I k goes to zero. For example, the transition occurs close to the mouth of Potomac River in winter while near bay bridge in summer. 3.3.3 Sensitivity to grazing parameters As in all NPZD-type models, the grazing terms are the closure of the model. Therefore, the specification of grazing rate and grazing preferences is very important. The model results show that changes in the grazing parameters (G m and the grazing preferences) affect the model to the largest extent in the upper and mid Bay in spring. While in other seasons, the effect is only apparent north of Bay bridge. This indicates that phytoplankton biomass is more grazing controlled in the upper bay and that this control is exerted further down Bay in spring when nutrient loads are higher. Conversely, this shows that the lower Bay is more subject to ?bottom up? control, i.e., nutrient, rather than grazing, limitation. 3.3.4 Sensitivity to sinking parameters The sinking rate of phytoplankton is an important parameter in our model. It transports phytoplankton from the surface to the bottom layers. Without sinking of phytoplankton, there is no accumulation of phytoplankton near the bottom in deep water anywhere in the model domain. Therefore, specification of the sinking rate of phytoplankton largely determines the redistribution of biomass in the water column. 151 However, as we mentioned in the previous section, phytoplankton sinking also results in the enhanced chlorophyll concentration on the shoals, especially along the eastern shore from upper to mid Bay. With phytoplankton sinking invoked, the model is not sensitive to the sinking rate of the detritus in the range we tested (1-4 m/d). However, when phytoplankton does not sink it can be very important especially in summer. The model is most sensitive to the sinking rate of detritus in summer and least sensitive in winter. As we mentioned before, the phytoplankton growth is mostly fueled by regenerated nutrients in summer and fueled by input from rivers at other times. Therefore, the amount of detritus that is deposited on the bottom directly determines the availability of DIN in the water column during summer as it is recycled through the heteotrophic compartment. Thus, these sensitivity results show that the model dynamics are consistent with our conceptual model of how detritus deposition associated with the spring phytoplankton bloom fuels summertime primary production through recycling (Malone et al., 1988). The sinking rate of ISS is very important in determining the ISS distribution, which largely determines the TSS distribution because ISS constitutes a large portion of TSS in the Bay. And TSS concentration in the upper layers greatly affects K d and the light available for phytoplankton growth. When the ISS sinking rate is high, it drives most of the mineral particles out of the surface layers. As a result, TSS concentrations in the upper layers drop and become much less variable. This, in turn, causes water clarity to increase (K d goes down) and become less variable. These kinds of effects are most pronounced in the upper Bay where the ISS loading is 152 highest and where TSS is the primary factor that controls light penetration. The current sinking rate used in the model was obtained by trial-and-error to achieve reasonable TSS concentrations and variability in the upper layers and hence correct K d values and variability. 4. Summary and Conclusions In this study we incorporated a simple NPZD biological model with a representation of mineral particles and dissolved oxygen into a 3-D numerical hydrodynamic model of Chesapeake Bay. To keep the ecosystem model as simple as possible we parameterized some important biochemical processes in the Bay: (1) instead of implementing a fully uncoupled P-currency model we enforced P- limitation using a simple parameterization where we specify uniform seasonal P concentrations with low values in winter and spring relative to a fixed half-saturation constant for P uptake; (2) underwater light attenuation is calculated using a simple empirical optical model that was developed for Chesapeake Bay which can account for more than 70% of the observed K d variability; (3) denitrification loss of DIN is accounted for by adding a loss term at the bottom which is proportional to DIN concentration in the bottom layer; (4) the temperature-dependent grazing rate of heterotrophs is simplified by a step function and the seasonal food-web structure change is represented as switch in grazing preferences; and (6) a sinking rate is added to the phytoplankton compartment that is modified seasonally to account bottom accumulation of chlorophyll and seasonal variations in phytoplankton species composition and export. 153 With its simple configuration, our model succeeds in producing the observed patterns in DIN, phytoplankton, DO, TSS and K d in the main stem in a dry year (1995) and a very wet year (1996) (Figs. 4.4-14). Modeled surface DIN distributions along the main stem compare quite favorably with the observed concentrations (Fig. 4.4). And the vertical DIN profiles in summer and fall also compare reasonably well with the observed. However, significant discrepancies remain in winter and spring, perhaps related to a combination of insufficient vertical mixing and too much remineralization at depth (Fig. 4.10). For phytoplankton, both the surface and subsurface distributions in the model generally compare well with the observations (Figs. 4.5 and 4.9). Discrepancies are mainly in the upper bay in 1996. In general, the model reproduces the observed concentrations better at locations with less variability (generally mid to lower bay). Even though DO simply cycles with DIN in the model (except for air-sea exchange), the vertical profiles of DO compare favorably with observations at most time and locations. The biggest problem is the overestimation of DO near bottom over the deep channel in summer, which may related to insufficient organic matter accumulation or remineralization in the spring that is too rapid. Correctly modeling the TSS distribution is as complex as modeling the biogeochemical system, but must be included at some level because TSS concentrations have a strong influence on light penetration. Despite the simplicity of our TSS parameterization, we obtain reasonable TSS distributions across the Bay except in lower bay, where the model tends to underestimate both the mean and variability. As a result, the model also reproduces the observed K d variability quite well except in the lower Bay. 154 The spatial picture of chlorophyll distribution in both years reveals that the model produces unrealistically high chlorophyll concentrations on the shoal regions, especially along the eastern shore, from the upper to mid bay. Using different test runs, we suggest that this problem may be due to grazing control at the shoals and/or physical processes that transport phytoplankton laterally from the shoals to the deep channels. We test the former hypothesis by implementing bivalve filtration which succeeds in reducing the high chlorophyll concentrations in these shoal regions. However, the lateral transport hypothesis is difficult to test with our current model setting. In addition, we also carried out some sensitivity studies with the simple biological model. Varying P-limitation in the model shows that P-limitation is important in both years but to a larger extend in 1996 when freshwater flow was much higher. The model is quite sensitive to the light saturation parameter I k for phytoplankton growth in winter and spring. While in summer and fall, the effect is confined to the upper bay only, revealing the seasonal and spatial extent of light limitation. As the closure of the model, we expect the grazing parameters are very important. Interestingly, the effects of adjusting the grazing parameters is mainly manifested in the upper bay except in spring when they extend much further south, which indicates that at most of time from mid to lower bay phytoplankton biomass is not grazing limited, but rather nutrient limited in the main stem. We have a sinking flux specified for three variables in our model: phytoplankton, detritus and ISS. Phytoplankton and ISS sinking is necessary to transport mass to deeper water and the 155 rate greatly affects their vertical distributions. With phytoplankton sinking, the model becomes less sensitive to the detritus sinking rate. In this study, we have demonstrated that a relatively simple biological model is capable of reproducing the major features in nutrient concentrations, phytoplankton biomass, oxygen concentration and underwater light attenuation in a complex biogeochemical system like Chesapeake Bay. And the model is robust enough to generate reasonable results under extreme conditions, i.e., in both a dry year (1995) and a very wet year (1996). Nonetheless, we also uncovered a number of significant discrepancies which suggest possible future improvement to the model. For example, inclusion of an explicit sediment layer could help to improve our solution in several ways. That is, by allowing a different zone for remineralization other than in the water column, and a more realistic interface for nutrient flux. An explicit sediment layer would also allow representation of sediment transport across the bed, which may help reduce unrealistically high chlorophyll concentrations in shoal regions, and a better representation of remineralization and denitrification, which might produce better DIN and oxygen profiles. Other potential improvement includes implementation of phytoplankton and TSS biomass dependent sinking rate for phytoplankton that can capture the enhanced sinking rates and export flux that seem to occur in turbid, productive estuaries. It may also be necessary to include grazing control in shoal regions in Chesapeake Bay to help reduce the high chlorophyll concentration in shoal regions at both the head of bay and along the shores. 156 Acknowledgements We would like to thank Drs. Tom Fisher, Larry Harding and Shenn-Yu Chao for helpful discussions and valuable input. This research was supported by funds from NASA, the National Science Foundation and Horn Point Laboratory?s Education Committee. References Airoldi, L. 2003. The effects of sedimentation on rocky coast assemblages. Oceanography and Marine Biology 41, 161-236. Bicknell, B., J. Imhoff, J. Kittle, A. Donigian Jr., R. Johanson and T. Barnwell. 1996. Hydrologic Simulation Program-Fortran user?s manual for release 11. U.S. Environmental Protection Agency Environmental Research Laboratory, Athens, GA. Boynton, W. R., W. M. Kemp, 1985. Nutrient regeneration and oxygen consumption by sediments along an estuarine salinity gradient. Marine Ecology Progress Series 23, 45-55. Boynton, W. R., Garber, J. H., Summers, R. and Kemp, W. M. 1995. Inputs, transformations, and transport of nitrogen and phosphorus in Chesapeake Bay and selected tributaries. Estuaries 18 (1B): 285-314. Bowers H. A., E. Tengs, H. B. Glasgow, J. M. Burkholder, P. A. Rublee and D. W. Oldach. 2000. Development of real-time PCR assays for rapid detection of Pfiesteria piscicida and related dinoflagellates. Applied and Environmental Microbiology 66 (11), 4641-4648. 157 Carpenter J. H., D. W. Pritchard and R. C. Whaley. 1969. Observations of eutrophication and nutrient cycles in some coastal plain estuaries. In: Eutrophication: causes, consequences, correctives-proceedings of a symposium. US national academy science publ 1700, Washington DC, pp210- 221. Cerco, C. F. and T.M. Cole 1994 Three-dimensional eutrophication model of Chesapeake Bay, Volume I: Main Report. U.S. Army Corps of Engineers Waterways Experiment Station,Vicksburg, MS. Technical Report EL-94-4. Cerco, C. F. and Noel M. R. 2004. The 2002 Chesapeake Bay eutrophication model. EPA 903-R-04-004. Prepared for Chesapeake Bay Program Office, U.S. Environmental Protection Agency. U.S. Army Corps of Engineers Waterways Experiment Station,Vicksburg, MS. Chang, J., P. Middleton, W. Stockwell, C. Walcek, J. Pleim, H. Lansford, S. Madronich, F. Binkowski, N. Seaman and D. Stauffer. 1990. The regional acid deposition model and engineering model, NAPAP SOS/T report 4. In National Acid Precipitation Assessment Program: State of science of technology, Vol. 1, National Acid Precipitation Assessment Program, Washington, D.C. Cooper S. R. and G. S. Brush. 1991. Long-term history of Chesapeake Bay anoxia. Science 254, 992-996. Denman K. L. 2003 Modelling planktonic ecosystems: parameterizing complexity. Progress in Oceanography 57 (3-4), 429-452. 158 Dennis, R. 1996. Using the Regional Acid Deposition Model to determine the nitrogen deposition airshed of the Chesapeake Bay watershed. In Atmospheric Deposition to the Great Lakes and Coastal Waters. Ed.: Joel Baker. Society of Environmental Toxicology and Chemistry. D?Elia, C. F., D. M. Nelson, W. R. Boynton, 1986. Nutrient enrichment studies in a coastal plain estuary: phytoplankton growth in large-scale, continuous cultures. Canadian Journal of Fishery and Aquatic Science 43, 397-406. Di Toro, D. and J. Fitzpatrick. 1993. Chesapeake Bay sediment flux model. U.S. Army Corps of Engineers Waterways Experiment Station, Vicksburg, MS. Contract report EL-93-2. Doney, S. C., D. M. Glover and R. B. Najjar, 1996. A new coupled, one-dimensional biological-physical model for the upper ocean: Applications to the JGOFS Bermuda Atlantic time-series study (BATS) site. Deep-Sea Research II 43 (2- 3), 591-624. Fisher, T.R., P.R. Carlson and R.T. Barker, 1982. Sediment nutrient regeneration in three North Carolina estuaries. Estuarine, Coastal and Shelf Science 14: 101- 116. Fisher, T. R., E. R. Peele, J. W. Ammernan, L. W. Harding, 1992. Nutrient limitation of phytoplankton in Chesapeake Bay. Marine Ecology Progress Series 82, 51- 63. Friedrichs, M. A. M., R. R. Hood and J. D. Wiggert, 2004. Ecosystem model complexity versus physical forcing, quantification of their relative impact with assimilated Arabian Sea data. Deep Sea Research II. 159 Glibert, P. M., D. J. Conley, T. R. Fisher, L. W. Harding, T. C. Malone. 1995. Dynamics of the 1990 winter spring bloom in Chesapeake Bay. Marine Ecology Progress Series 122 (1-3), 27-43. Greene, K. and L. Linker. 1998. Chesapeake Bay Watershed Model application and calculation of nutrient and sediment loadings-phase IV Chesapeake Bay Watershed Model-appendix a: model hydrology calibration results. EPA 903- R-98-004, CBP/TRS 196/98, Chesapeake Bay Program Office, Annapolis, MD. Gundersen, J. S., W. D. Gardner, M. J. Richardson and I. D. Walsh 1998 Effects of monsoons on the seasonal and spatical distributions of POC and chlorophyll in the Arabian Sea. Deep-Sea Research 45, 2103-2132. Harding, L. W. Jr. 1994 Long-term trends in the distribution of phytoplankton in Chesapeake Bay: role of light, nutrients and streamflow. Marine Ecology Progress Series 104, 267-291. Harding, L. W. Jr., A. Magnuson, and M. E. Mallonee. 2004. SeaWiFS retrievals of chlorophyll in Chesapeake Bay and the mid-Atlantic bight. Estuarine Coastal and Shelf Science, in press. Harding, L. W. Jr., M. E. Mallonee and E. S. Perry 2002 Toward a predictive understanding of primary productivity in a temperate, partially stratified estuary. Estuarine Coastal and Shelf Science 55 (3), 437-463. Harding, L. W., B. W. Meeson, T. R. Fisher. 1986. Phytoplankton production in two east coast estuaries: photosynthesis-light functions and patterns of carbon 160 assimilation in Chesapeake and Delaware Bays. Estuarine, Coastal and Shelf Science 23, 773-806. Hood, R. R., H. V. Wang, J. E. Purcell, E. D. Houde and L. W. Harding Jr.. 1999. Modeling particles and pelagic organisms in Chesapeake Bay: Convergent features control plankton distributions. Journal of Geophysical Research 104, 1223-1243. Hood, R.R., N.R. Bates, D.G. Capone and D.B. Olson 2001 Modeling the effect of nitrogen fixation on carbon and nitrogen fluxes at BATS. Deep-Sea Research II, 48, 1609-1648. Hood, R. R., K. E. Kohler, J. P. McCreary, S. L. Smith. 2003. A four-dimensional validation of a coupled physical-biological model of the Arabian Sea. Deep- Sea Research II 50 (22-26), 2917-2945. Huntley, M. E. and M. D. G. Lopez. 1992 Temperature-dependent production of marine copepods: a global synthesis. The American Naturalist 140 (2), 201- 242. Johnson, B. H., K. W. Kim, R. H. Heath, H. L. Butler and B. B. Hsieh 1991 Development and verification of a three-dimensional numerical hydrodynamic, salinity, and temperature model of the Chesapeake Bay, Technical Report HL-91-7, US Army Engineer Waterways Experiment Station, Vicksburg, MS. Kemp, W.M., P. A. Sampou, J. Garber, J. Tuttle and W. R. Boynton. 1992. Seasonal depletion of oxygen from bottom waters of Chesapeake Bay: roles of benthic 161 and planktonic respiration and physical exchange processes. Marine Ecology Progress Series 85, 137-152. Lapointe, B. E. and M. W. Clark, 1992. Nutrient input from the watershed and coastal eutrophication in the Florida Keys. Estuaries 15 (4), 465-476. Levy, H. and W. J. Moxim 1987. Fate of U.S. and Canadian combustion nitrogen emissions. Nature 328, 414-416. Li, M., A. Gargett and K. Denman. 2000. What determines seasonal and interannual variability of phytoplankton and zooplankton in strongly estuarine systems? Application to the semi-enclosed estuary of Strait of Georgia and Juan De Fuca Strait. Estuarine Coastal and Shelf Science 50 (4), 467-488. Lima, I. D. and S. C. Doney, 2004. A three-dimensional, multinutrient, and size- structured ecosystem model for the North Atlantic. Global Biogeochemical Cycles 18 (3), Art. No. GB3019. Linker L. C., G. W. Shenk, P. Wang, K. J. Hopkins and S. Pokharel. 2001. A short history of Chesapeake Bay modeling and the next generation of watershed and estuarine models. U.S. EPA Chesapeake Bay Program, Annapolis, MD. Logan, J. A. 1983. Nitrogen oxides in the troposphere: global and regional budgets. Journal of Geophysical Research 88, 10785-10807. Madden, C. and W. M. Kemp. 1996. Ecosystem model of an estuarine submersed plant community: calibration and simulation of eutrophication responses. Estuaries 19 (2B), 457-474. 162 Magnuson, A., L. W. Harding Jr., M. E. Mallonee, and J. E. Adolf. 2004. Bio-optical model for Chesapeake Bay and the Middle Atlantic Bight. Estuarine Coastal and Shelf Science 61, 403-424. Malone, T.C. 1992. Effects of water column processes on dissolved oxygen, nutrients, phytoplankton and zooplankton, in oxygen dynamics in the Chesapeake Bay, a synthesis of recent research, edited by D.E. Smith, M. Leffler and G.E. Machiernan, Maryland Sea Grant College, Univ. of Maryland, College Park, Maryland, pp61-112. Malone, T. C., L.H. Crocker, S.E. Pike and B.W. Wendler, 1988. Influences of river flow on the dynamics of phytoplankton production in a partically stratified estuary. Marine Ecology Progress Series 48, 235-249. Malone T. C. and H. W. Ducklow. 1990. Microbial biomass in the coastal plume of Chesapeake Bay: Phytoplankton-bacterioplankton relationships. Limnology and Oceangraphy 35 (2), 296-312. Malone T. C., H. W. Ducklow, E. R. Peele and S. E. Pike. 1991. Picoplankton carbon flux in Chesapeake Bay. Maine Ecology Progress Series 78 (1): 11-22. Malone T. C., W. M. Kemp, H. W. Ducklow, W. R. Boynton, J. H. Tuttle and R. B. Jonas. 1986. Lateral variation in the production and fate of phytoplankton in a partially stratified estuary. Marine Ecology Progress Series 32, 149-160. Marra, J. and C. Ho, 1993. Initiation of the spring bloom in the northeast Atlantic (47- degrees-N, 20 degrees-W)- A numerical-simulation. Deep-Sea Research II 40 (1-2), 55-73. 163 McClain, C. R., K. Arrigo, K. S. Tai and D. Turk, 1996. Observations and simulations of physical and biological processes at ocean weather station P, 1951-1980. Journal of Geophysical Research-Oceans 101 (C2), 3697-3713. Meyers, T., J. Sickles, R. Dennis, K. Russell, J. Galloway and T. Church. 2000. Atmospheric nitrogen deposition to coastal estuaries and their watersheds, in nitrogen loading in coastal water bodies: an atmospheric perspective, edited by Richard A. Valigura, American Geophysical Union. Miller, D. C., C. L. Muir and O. A. Hauser. 2002. Detrimental effects of sedimentation on marine benthos: what can be learned from natural processes and rates? Ecological Engineering 19 (3), 211-232. Moore, J. K., S. C. Doney, J. A. Kleypas, D. M. Glover and I. Y. Fung. 2002. An intermediate complexity marine ecosystem model for the global domain. Deep Sea Research II 49 (1-3), 403-462. Nagy, G. J., M. Gomez-Erache, C. H. Lopez and A. C. Perdomo, 2002. Distribution patters of nutrients and symptoms of eutrophication in the Rio de la Plata river estuary. Hydrobiologia 475 (1), 125-139. Newell, R. I. E., and E. W. Koch. 2004. Modeling seagrass density and distribution in response to changes in turbidity stemming from bivalve filtration and seagrass sediment stabilization. Estuaries, submitted. Orth R. J. and K. A. Moore. 1983. Chesapeake Bay-an unprecedented decline in submerged aquatic vegetation. Science 222 (4619), 51-53. 164 Oschlies, A. and V. Garcon, 1999. An eddy-permitting coupled physical-biological model of the North Atlantic- 1. Sensitivity to advection numerics and mixed layer physics. Global Biogeochemical Cycles 13 (1), 135-160. Pitkanen H., T. Tamminen, P. Kangas, T. Huttula, K. Kivi, H. Kuosa, J. Sarkkula, K. Eloheimo, P. Kauppila and B. Skakalsky, 1993. Late summer trophic conditions in the northeast gulf of Finland and the river Neva estuary, Baltic Sea. Estuarine Coastal and Shelf Science 37 (5), 453-474. Seitzinger, S. P. 1988. Denitrification in freshwater and coastal marine ecosystem: Ecological and geochemical significance. Limnology and Oceanography 33, 702-724. Seitzinger, S. P. and A. E. Giblin. 1996. Estimating denitrification in North Atlantic continental shelf sediments. Biogeochemistry 35, 235-260. Skogen, M. D., E. Svendsen, J. Berntsen, D. Aksnes and K. B. Ulvestad, 1995. Modeling the primary production in the North-Sea using a coupled 3- dimensional physical-chemical-biological ocean model. Estuarine Coastal and Shelf Science 41 (5), 545-565. Varela, R. A., A. Cruzado, J. Tintore and E. G. Ladona 1992 Modeling the deep- chlorophyll maximum-a coupled physical-biological approach. Journal of Marine Research 50 (3), 441-463. Venrick, E. L. 1991 Midocean ridges and their influence on the large-scale patterns of chlorophyll and production in the North Pacific. Deep-Sea Research Part A 38, S83-S102. 165 Wanninkhof, R.1992 Relatioship between wind speed and gas exchange over the ocean. Journal of Geophysical Research 97, 7373-7382. Wetzel, R. and H. Neckles. 1986. A model of Zostera marina L. photosynthesis and growth: simulated effects of selected physical-chemical variables and biological interactions. Aquatic Botany 26, 307-323. Xu, J., S.Y. Chao, R.R. Hood, H.V. Wang and W.C. Boicourt 2002 Assimilating high- resolution salinity data into a model of a partially mixed estuary, J. Geophys. Res. 107 (C7), 3074, doi, 10.1029/200JC000626. Xu, J., R. R. Hood and S. Y. Chao. 2004. A simple empirical optical model for simulating light attenuation variability in a partially mixed estuary. Estuaries. Submitted. 166 Figures 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 Chapter 5: Summary and Conclusions The main objective of this study is to explore if a simple biological model will work in a complex ecosystem such as Chesapeake Bay under proper physical forcing and what essential components are necessary in such a simple model to reproduce the observed variability in nutrient concentrations, oxygen and phytoplankton biomass. For this purpose we developed a relatively simple biogeochemical model that includes several paramterizations of key processes such as dentrification, P limitation, seasonal changes in ecosystem structure and temperature effects on remineralization, and a simple empirical light attenuation submodel. This biogeochemical model is coupled with a 3-D hydrodynamic model of Chesapeake Bay - WES-CH3D. WES-CH3D is 3-D primitive-equation model solving for salinity, temperature, water level elevation and velocities. The governing equations are recast in a boundary-fitted curvilinear coordinates system to cope with the irregular shoreline configuration and deep channel orientation of Chesapeake Bay. We decided to use WES-CH3D as our physical model because it has been tuned to Chesapeake Bay specifically and validated with a ten-year (1985-1994) time series of observations obtained by Chesapeake Bay Program and the Bay Program also uses its output to drive their water quality model. Therefore, the results from our simple biological model provide a direct reference to their more complicated water quality model. 191 Due to importance of physical forcing in biological processes we first undertook an effort to validate the physical model for our study period: 1995 and 1996, which was chosen because it encompasses a dry year (1995) and an exceptionally wet year (1996). The physical model reproduces the essential circulation features, such as the two-layered circulation in the main channel and major tributaries and the temperature and salinity structure compares reasonably well with observations. However, sizable discrepancies are found at depth and laterally in salinity profiles: the model tends to overestimate bottom water inflow and salinity at depth and the discrepancy increases toward the lower reaches of the bay so that the water column in the model tends to be more stratified, especially in the mid to lower bay. The enhanced bottom inflow with saltier bottom water in the model is a result of the coarse resolution of the model. The numerical representation of three- dimensional flow and density fields by a finite number of computation cells inevitably increases friction. To overcome the numerical damping, it is often necessary to enhance the bottom inflow of seawater from the estuarine mouth region in order to produce a realistic two-layered circulation well inside the estuary. To enhance further model realism without increasing grid resolution at the expense of computation speeds, an attempt was made to assimilate available high-resolution Scanfish salinity data in 1995 into the model. A Newtonian relaxation scheme was used in the salinity equation to receive the irregularly spaced time series of Scanfish data. After data assimilation the agreement between the model and the observations is generally improved and the improvement is generally more profound in the middle 192 reaches than in the upper and lower reaches of the bay. The improvement is maximal right after the data injections and decreases slowly thereafter. In the lower bay the improvement diminishes more quickly because the adjacent coastal ocean is excluded from data assimilation. The data assimilation also enhances the two-layered circulation. The speed enhancement ranges up to 4 cm s -1 . However, the improvement does not come without penalties. The data assimilation may trigger renewed gravitational circulation because the density structure is significantly altered. As mentioned above, the model tends to overestimate salinity at depths and towards the lower reaches of the bay and the intention of data assimilation is to reduce modeled salinity at depths and in the lower reaches of the bay. Because there are only two short periods of high-resolution data available the restoration rate has to be strong enough to make a lasting impact. When data are assimilated in the lower reaches of the bay, pressure is reduced near the mouth. The consequent increase in the seaward pressure gradient triggers the seaward expansion of the buoyant layer. The artificial seaward expansion is eliminated with continuous data injections. The strength and the longitudinal extent of gravitational readjustment are largely controlled by the restoration rate K. Precautions have to be taken in order to maximize the model improvement and minimize the volatile transients, which is necessary because of the limited availability of high-resolution data. Because only a few of narrow windows of high-resolution data are available, the nudging must be strong enough the make a difference. Even though the gravitational readjustment is transient it is problematic when running a 193 biogeochemical model because it can potentially have profound biogeochemical impacts. The seaward expansion of the buoyant layer may wash out the biological organisms and produce an unrealistic distribution of biological variables. Therefore, we made another attempt to improve the model performance by adjusting the salinity at the open-ocean boundaries. The adjustments were inferred from the observations at inner estuary by trial-and-error and they improved the salinity structure at depth, especially in the upper and mid reaches of the bay. After validating the physical model we embarked upon the development and implementation of the biogeochemical model. An immediate problem is how to reproduce the observed underwater light field variability. In aquatic environments, light tends to limit primary production in deep water. In Chesapeake Bay, light can be a limiting factor even in shallow waters under turbid conditions, and light attenuation is modulated by a variety of chromophoric substances, including dissolved organic matter and detritus. Because our goal is to keep the whole biogeochemical model as simple as possible, the complex spectral optical models, such as in Platt and Sathyendranath (1988), Smith et al. (1989) and Gellegos et al. (1990), were not considered. Moreover, these complex models are not easily applied in biogeochemical models where only a handful of the optically active constituents are represented. In our first attempt we used dual-wavelength light model as in Hood et al. (2001), but this model failed to capture the large range of light attenuation variability in Chesapeake Bay. Therefore, we derived an empirical light model specifically for Chesapeake Bay. Assuming the averaged light attenuation coefficient (K d (PAR) ) can be decomposed as a set of partial attenuation coefficients and each 194 partial attenuation coefficient is a linear function of concentration of each waterborne material, the values of each specific attenuation coefficient can be ?backed out? by multiple linear regression methods when data are available. A light model obtained by this approach will be simple and specific to the water body where the data are from. A total of 1348 data points from 1995 and 1996, where underwater light intensity (from which K d is derived) as well as water quality concentrations were measured simultaneously by Chesapeake Bay Program, were used to develop the empirical light model. In Chesapeake Bay, the main optically active constituents besides water itself are: phytoplankton, seston (non-phytoplankton particulate matter) and CDOM. Given the data availability and the constraints imposed by what we can represent in our biogeochemical model, chlorophyll, TSS and salinity were used in the linear regression relationship. Salinity is included as a proxy for CDOM because CDOM behaves conservatively like salinity in Chesapeake Bay (Rochelle-Newall and Fisher, 2002). The resulting model gives an R 2 of 0.72 between the calculated K d from this relation (K d _predicted) and the K d derived from direct light measurement (K d _observed). A stepwise statistical model showed that TSS is by far the most important factor in controlling light attenuation in Chesapeake Bay, explaining about 58% of the total variability in K d alone. Contrary to oceanic water, phytoplankton absorption only plays a minor role in controlling the light field in Chesapeake Bay. Chlorophyll only enters the model at 5% level and does not improve the R 2 . In addition, the specific light attenuation coefficient for chlorophyll turned out to be negative. We interpreted this to be a consequence of the strong effect of light 195 controlling phytoplankton growth in the system where K d is large and determined primarily by TSS and CDOM. Even though the regression relation for the entire bay gives a nice fit between K d _predicted and K d _observed it has one caveat: K d values can become negative at high salinity regions, which can be inconvenient in numerical modeling. Therefore, we divided the dataset into two groups by salinity (one for S ? 15 and one for S > 15) and for each group the same method was applied. This approach did not further improve the overall R 2 but it solved the negative K d problem. In addition, the changes in the specific light attenuation coefficients between the two salinity regimes demonstrated how the role of different optically active constituents can change over a wide range of salinity in an estuarine environment. Specifically, the intercept, the coefficients for salinity and TSS decreases from low to high salinity regions and the coefficient for chlorophyll changed sign (negative in low salinity, positive in high salinity). The smaller intercept and coefficient for salinity in high salinity regions show that the lower bay has less CDOM and CDOM has less influence on light attenuation. The decrease in coefficient for TSS indicates that TSS plays a less important role in attenuating light in lower bay. And we believe that the sign change in the coefficient for chlorophyll is results from two competing factors in determining the relationship between K d and chlorophyll. Namely, in turbid waters where constituents other than phytoplankton strongly influence K d , light controls phytoplankton growth and biomass, which will tend to give rise to a negative correlation between K d and chlorophyll; while in clearer water where phytoplankton growth and biomass are controlled by factors other than light, chlorophyll strongly 196 influences K d which will tend to give rise to a positive correlation between K d and chlorophyll. Because field data always contain information from both these two competing factors, an empirical model will tend to underestimate the effect of chlorophyll variability on K d . From the two relations for different salinity regime we estimated the contribution of each component in total K d . In the low salinity region, light attenuations by seston and water + CDOM are equally important while phytoplankton?s contribution is mostly below 10%. In the high salinity region, light attenuation due to water and CDOM dominates and phytoplankton plays a bigger role, while the contribution from seston becomes very small. However, one has to keep in mind that the contribution from phytoplankton estimated from empirical model tends to be underestimated, as we discussed above. Nevertheless, this analysis shows that in modeling the light field in estuaries, the first order importance is to reproduce both the mean distribution and variability of TSS and CDOM. In systems where CDOM behaves conservatively, CDOM concentration can be adequately represented by salinity. Therefore, in such systems TSS is the crucial component in modeling the underwater light field. In this part of study, we also tested the possibility of using Secchi depth (SD) derived K d in place of K d from direct light measurements. Three conversions from SD to K d were used: one fitted to the specific dataset, one obtained from Choptank River and another commonly used conversion. We then recalculated the regression equation relating these derived K d values to chlorophyll, TSS and salinity, which yielded three different relations (one for each SD to K d conversion equation). All three relations 197 describe much less of the observed K d variability than the model derived using direct light measurements. Moreover, they all tend to underestimate high K d and overestimate low K d . The severity of the biases at both extremes suggests that the SD derived K d values cannot be used in the derivation of an empirical light model in Chesapeake Bay. For the main body of the biogeochemical model, we adopted the simple five- compartment NPZD-type model (DIN, DON, phytoplankton, heterotrophs and detritus) in Hood et al. (2001). One DIP compartment was added to simulate phosphorus limitation in Chesapeake Bay. Two additional state variables were also included: dissolved oxygen (DO) and inorganic suspended solids (ISS). DO is included to simulate anoxia and hypoxia in the Bay and it also serves as a natural trigger to slow down the respiratory processes of heterotrophs under these conditions. ISS is included for dynamically modeling the underwater light field. The biological model is forced by nutrient input and suspended solids loads from major tributaries, atmospheric deposition of DIN and DON and air-water exchange of O 2 . In light of the important role played by nutrients in the biogeochemical activities in the system we assimilated available observations in the upper reaches of the tributaries in the model. This approach allows us to not only have a relatively correct estimate of the nutrient loading but also bypass the difficulty of accounting for some uncertain nutrient sinks/sources, such as nutrient retention/release by marshes and wetlands at the upper tributaries. The ecosystem model was then embedded in the hydrodynamic model of Chesapeake Bay and run for 1995 and 1996. In the process of diagnosing the model 198 performance we gradually added some simple parameterizations of important biochemical processes in the Bay: 1) P limitation is enforced by using a uniformly low value for DIP in winter and spring relative to a fixed half saturation constant for phosphorus uptake. 2) Under water light attenuation is calculated using our simple empirical optical model. 3) Denitrification loss of DIN is taken to be proportional to DIN concentrations in the bottom layer. 4) The temperature-dependent grazing rate of heterotrophs is represented by a step function. 5) The seasonal change in food-web structure is represented by changes of grazing preferences. 6) The sinking rate of phytoplankton is modified seasonally to represent the changes of species composition and account for the bottom accumulation of phytoplankton in spring. To show the overall performance of this coupled physical-biological model, both the surface distributions and vertical profiles (except for K d ) of DIN, phytoplankton, DO, TSS and K d were examined. With its simple configuration, the model successfully produced the general distribution of each field (both the mean and variability) in the mainstem of the bay and reproduced the observed seasonal and interannual patterns. However, some significant discrepancies were also observed. With DIN, the main problem resides in winter and spring: In winter, observed DIN is high at surface and decreases with depth. However, the modeled DIN is low at surface and high near bottom. In spring, the model generally overestimates DIN at depth. For phytoplankton, both the surface and subsurface distributions in the model generally compare favorably with observations. Discrepancies are mainly in the upper bay in 1996 and the model generally produces better chlorophyll distributions in the mid to lower bay where the variability is lower. Even though DO is simply cycled 199 with DIN in the model with air-sea exchange at the surface the model produces reasonable DO distributions all the time over the whole Bay. One pronounced problem is the overestimation of DO near bottom over the deep channel in summer, which may be related to too much reminerization near bottom and hence insufficient organic matter accumulation in spring. In Chesapeake Bay, inorganic suspended solids make up a major portion of TSS. In our model, ISS is simply modeled with loading, uniform sinking, deposition and resuspension rate. There is no size-structure and no aggregation effect. However, the modeled distribution of TSS is reasonable throughout the bay and all seasons except in lower bay where the model tends to underestimate both the mean and variability. The problem in the lower bay may be related to insufficient lateral loading and resuspension. Because K d becomes less dependent on TSS this does not affect the K d field very much. The modeled K d has less variability in the lower bay mainly because the empirical light model cannot explain much of the variability in the lower bay. The spatial patterns of DIN and chlorophyll in both years agree with our general understanding of the bay dynamics. In the low flow season/year (1995), high DIN concentrations are confined to the northernmost part of the Bay and the chlorophyll maximum occurs in the upper bay as well. During the high flow season/year (1996), due to the large amount of TSS loading the upper bay becomes light-limited so that chlorophyll maximum moves further south, and high DIN concentration also extends down the bay. However, the horizontal spatial plots of phytoplankton biomass for both years also revealed that the model produced unrealistically high chlorophyll concentrations in the shoal areas, especially along the 200 eastern shore, from upper to mid bay. We tested several hypotheses and suggest that this problem may be due to grazing control at the shoals and/or physical processes that transport phytoplankton laterally from the shoals to the deep channels. We test the former hypothesis by implementing bivalve filtration which succeeds in reducing the high chlorophyll concentrations in these shoal regions. However, the lateral transport hypothesis is difficult to test with our current model setting. Simply increasing horizontal mixing does not work. Taking advantage of the simple configuration we also carried out sensitivity studies on some key parameters. Varying the degree of P-limitation shows that P- limitation is important in both years but to a larger extend in 1996 when freshwater inflow was much higher and that P-limitation tends to be more important in the upper Bay. The model is quite sensitive to the light saturation parameter I k for phytoplankton growth in winter and spring. While in summer and fall, the effect is confined to the upper bay only, revealing the seasonal and spatial extent of light limitation. As the closure of the model, we expect the grazing parameters to be very important. Interestingly, the effects of adjusting the grazing parameters is mainly manifested in the upper bay except in spring when they extend much further south, which indicates that at most of time from mid to lower bay phytoplankton biomass is not grazing limited, but rather nutrient limited in the main stem. We have a sinking flux specified for three variables in our model: phytoplankton, detritus and ISS. Phytoplankton and ISS sinking is necessary to transport mass to deeper water and the rate greatly affects their vertical distributions. With phytoplankton sinking, the model becomes less sensitive to the detritus sinking rate. 201 The major conclusions that we derive from this study are: 1) Assimilating high-resolution data, even only in a short period, can improve the performance of the physical model; 2) a simple empirical optical model is capable of producing the light attenuation variability but it is necessary to carry TSS in a numerical model to dynamically model underwater light field in a complex estuarine environment such as Chesapeake Bay; 3) a relatively simple biochemical model is capable of reproducing the major features in phytoplankton biomass, nutrient, oxygen and TSS distributions and is robust enough to generate reasonable results in both wet and dry years. Potential future work includes: 1) more detailed examination of the physical model, especially with regard to the mixing parameters with efforts focused on generating more realistic stratification and mixing depth; 2) including an explicit sediment layer in the ecosystem model to allow a different zone for remineralization other than in the water column, and a more realistic interface for nutrient flux. An explicit sediment layer would also allow representation of sediment transport across the bed and a better representation of remineralization and denitrification; 3) implementation of phytoplankton and TSS biomass dependent sinking rate for phytoplankton that can capture the enhanced sinking rates and export flux that seem to occur in turbid, productive estuaries; 4) including realistic representation of grazing control in shoal regions in Chesapeake Bay to help reduce the high chlorophyll concentration in shoal regions at both the head of bay and along the shores. 202 References Hood, R.R., N.R. Bates, D.G. Capone and D.B. Olson. 2001. Modeling the effect of nitrogen fixation on carbon and nitrogen fluxes at BATS. Deep-Sea Research II, 48, 1609-1648. Gallegos, C. L., D. L. Correll, and J. W. Pierce. 1990. Modeling spectral diffuse attenuation, absorption, and scattering coefficients in a turbid estuary. Limnology and Oceanography 35: 1486-1502. Platt, T., and S. Sathyendranath. 1988. Ocean primary production: estimation by remote sensing at local and regional scales. Science 241: 1613-1620. Rochelle-Newall, E. J. and T. R. Fisher. 2002. Chromophoric dissolved organic matter and dissolved organic carbon in Chesapeake Bay. Marine Chemistry 77: 23-41. Smith, R.C., B.B. Prezelin, R. R. Bidigare, and K. S. Baker. 1989. Bio-optical modeling of photosynthetic production in coastal waters. Limnology and Oceanography 34: 1524-1544. 203 References Airoldi, L. 2003. The effects of sedimentation on rocky coast assemblages. Oceanography and Marine Biology 41, 161-236. Anning,T., H. L. MacIntyre, S. M. Pratt, P. J. Sammes, S. Gibb and R. J. Geider. 2000. Photoacclimation in the marine diatom Skeletonema costatum. Limnology and Oceanography 45(8): 1807-1817. Anthes, D. L. T., Data assimilation and initialization of hurricane-predicting models, J. Atmos. Sci., 31, 702-719, 1974. Arrigo, K. R. and C. W. Sullivan, 1994. A high resolution bio-optical model of microalgal growth: tests using sea-ice algal community time-series data. Limnology and Oceanography 39 (3): 609-631. Baker, E. T., D. A. Tennant, R. A. Feely, G. T. Lebon and S. L. Walker. 2001. Field and laboratory studies on the effect of particle size and composition on optical backscattering measurements in hydrothermal plumes. Deep-Sea Research Part I 48 (2): 593-604. Belval, D. L. and L. A. Sprague, 2000. Monitoring nutrients in the major rivers draining to Chesapeake Bay, USGS water resources investigations report 99- 4238. Benjamin, T. B., Gravity currents and related phenomena, J. Fluid Mech., 31, 209- 248, 1968. Bicknell, B., J. Imhoff, J. Kittle, A. Donigian Jr., R. Johanson and T. Barnwell. 1996. Hydrologic Simulation Program-Fortran user?s manual for release 11. U.S. 204 Environmental Protection Agency Environmental Research Laboratory, Athens, GA. Blumberg, A. F., and G. L. Mellor, A description of a three-dimensional coastal ocean circulation model, in Three-dimensional coastal ocean models, Coastal Estuary Sci., Vol. 4, edited by N.S. Heaps, 1-16, AGU, Washington, D. C., 1987. Bogden, P. S., and J. O?Donnell, Generalized inverse with shipboard current measurements: Tides and nontidal flows in Long Island Sound, J. Mar. Res., 56, 995-1027, 1998. Bogden, P. S., P. Malanotte-Rizzoli, and R. Signell, 1996. Open-ocean boundary conditions from interior data: Local and remote forcing of Massachusetts Bay, Journal of Geophysical Research 101: 6487-6500. Bowers, D. G., and D. Evans, D. N. Thomas, K. Ellis, P. J. le B. Williams. 2004. Interpreting the colour of an estuary. Estuarine, Coastal and Shelf Science 59: 13-20. Bowers, H. A., E. Tengs, H. B. Glasgow, J. M. Burkholder, P. A. Rublee and D. W. Oldach, 2000. Development of real-time PCR assays for rapid detection of Pfiesteria piscicida and related dinoflagellates. Applied and Environmental Microbiology 66 (11): 4641-4648. Boynton, W. R., J. H. Garber, R. Summers and W. M. Kemp. 1995. Inputs, transformations, and transport of nitrogen and phosphorus in Chesapeake Bay and selected tributaries. Estuaries 18 (1B): 285-314. 205 Boynton, W. R. and W. M. Kemp, 1985. Nutrient regeneration and oxygen consumption by sediments along an estuarine salinity gradient. Marine Ecology Progress Series 23, 45-55. Boynton, W. R., W. M. Kemp and C. G. Osborne, 1980. Nutrient fluxes across the sediment-water interface in the turbid zone of a coastal plain estuary. Estuarine Perspectives, Academic Press. Carpenter J. H., D. W. Pritchard and R. C. Whaley. 1969. Observations of eutrophication and nutrient cycles in some coastal plain estuaries. In: Eutrophication: causes, consequences, correctives-proceedings of a symposium. US national academy science publ 1700, Washington DC, pp210- 221. Cerco, C. F. and T.M. Cole. 1994 Three-dimensional eutrophication model of Chesapeake Bay, Volume I: Main Report. U.S. Army Corps of Engineers Waterways Experiment Station,Vicksburg, MS. Technical Report EL-94-4. Cerco, C. F. and M. R. Noel. 2004. The 2002 Chesapeake Bay eutrophication model. EPA 903-R-04-004. Prepared for Chesapeake Bay Program Office, U.S. Environmental Protection Agency. U.S. Army Corps of Engineers Waterways Experiment Station,Vicksburg, MS. Chang, J., P. Middleton, W. Stockwell, C. Walcek, J. Pleim, H. Lansford, S. Madronich, F. Binkowski, N. Seaman and D. Stauffer. 1990. The regional acid deposition model and engineering model, NAPAP SOS/T report 4. In National Acid Precipitation Assessment Program: State of science of 206 technology, Vol. 1, National Acid Precipitation Assessment Program, Washington, D.C. Ciotti, A. M., M. R. Lewis and J. J. Cullen. 2002. Assessment of the relationships between dominant cell size in natural phytoplankton communities and the spectral shape of the absorption coefficient. Limnology and Oceanography 47 (2): 404-417. Cooper S. R. and G. S. Brush, 1991. Long-term history of Chesapeake Bay anoxia. Science 254: 992-996. Cullen, J. J. and M. R. Lewis. 1988. The kinetics of algal photoadaptation in the context of vertical mixing. Journal Of Plankton Research 10 (5): 1039-1063. D?Elia, C. F., D. M. Nelson, W. R. Boynton, 1986. Nutrient enrichment studies in a coastal plain estuary: phytoplankton growth in large-scale, continuous cultures. Canadian Journal of Fishery and Aquatic Science 43: 397-406. Deegan, L.A. and R.H. Garritt, 1997. Evidence for spatial varability in estuarine food webs, Marine Ecology Progress Series 147: 31-47. Denman K. L. 2003 Modelling planktonic ecosystems: parameterizing complexity. Progress in Oceanography 57 (3-4), 429-452. Dennis, R. 1996. Using the Regional Acid Deposition Model to determine the nitrogen deposition airshed of the Chesapeake Bay watershed. In Atmospheric Deposition to the Great Lakes and Coastal Waters. Ed.: Joel Baker. Society of Environmental Toxicology and Chemistry. 207 Di Toro, D. and J. Fitzpatrick. 1993. Chesapeake Bay sediment flux model. U.S. Army Corps of Engineers Waterways Experiment Station, Vicksburg, MS. Contract report EL-93-2. Dieguez, M.C. and J. J. Gilbert. 2003. Predation by Buenoa Macrotibialis (Insecta, Hemiptera) on zooplankton: effect of light on selection and consumption of prey. Journal Of Plankton Research 25 (7): 759-769. Doney, S. C., D. M. Glover and R. B. Najjar, 1996. A new coupled, one-dimensional biological-physical model for the upper ocean: Applications to the JGOFS Bermuda Atlantic time-series study (BATS) site. Deep-Sea Research II 43 (2- 3), 591-624. Edinger, J. E., D. K. Brady, and J. C. Geyer, Heat exchange and transport in the environment, Report 14, EPRI Publication No. 74-049-00-3, prepared for electric power research institute, Palo Alto, CA, 1974. Elliott, A.J. and D.P. Wang, 1978. The effect of meteorological forcing on the Chesapeake Bay: the coupling between an estuarine system and its adjacent coastal waters, In Hydrodynamics of Estuaries and Fjords, edited by NCN. Nihoul, Elsevier, Amsterdam, The Netherlands. Ezer, T., and G, K, Mellor, 1994. Continuous assimilation of Geosat altimeter data into a three-dimensional primitive equation Gulf Stream model. Journal of Physical Oceanography 24: 832-847. Ferrari, G. M, and M.D. Dowell. 1998. CDOM absorption characteristics with relation to fluorescence and salinity in coastal areas of the southern Baltic Sea. Estuarine, Coastal and Shelf Science 47: 91-105. 208 Fisher, D.C. and M. Oppenheimer, 1991. Atmospheric nitrogen deposition and the Chesapeake Bay estuary, AMBIO, vol. 20, No.3-4. Fisher, T.R., P.R. Carlson and R.T. Barker, 1982. Sediment nutrient regeneration in three North Carolina estuaries. Estuarine, Coastal and Shelf Science 14: 101- 116. Fisher, T. R., A. B. Gustafson, G. M. Radcliffe, K. L. Sundberg and J. C. Stevenson, 2003. A long-term record of photosynthetically available radiation (PAR) and total solar energy at 38.60 degrees N, 78.2 degrees W. Estuaries 26 (6): 1450- 1460. Fisher, T. R., A. B. Gustafson, K. Sellner, R. Lacouture, L. W. Haas, R. L. Wetzel, R. Magnien, D. Everitt, B. Michaels and R. Karrh, 1999. Spatial and temproral variation of resource limitation in Chesapeake Bay. Marine Biology 133 (4): 763-778. Fisher, T. R., L. W. Harding, D. W. Stanely, L.G. Ward, 1988. Phytoplankton, nutrients, and turbidity in the Chesapeake, Delaware, and Hudson estuaries. Estuarine, Coastal and Shelf Science 27 (1): 61-93. Fisher, T. R., E. R. Peele, J. W. Ammernan, L. W. Harding, 1992. Nutrient limitation of phytoplankton in Chesapeake Bay. Marine Ecology Progress Series 82, 51- 63. Forbes, C., and O. Brown, 1996. Assimilation of sea surface height data into an isopycnic ocean model. Journal of Physical Oceanography 26: 1189-1213. 209 Fox, L.E., S.L. Sager and S.C. Wofsy, 1985. Factors controlling the concentrations of soluble phosphorus in the Mississippi estuary. Limnology and Oceanography 30: 826-832. Friedrichs, M. A. M., R. R. Hood and J. D. Wiggert, 2004. Ecosystem model complexity versus physical forcing, quantification of their relative impact with assimilated Arabian Sea data. Deep Sea Research II. Gal, G., E. R. Loew, L. G. Rudstam and A. M. Mohammadian. 1999. Light and diel vertical migration: spectral sensitivity and light avoidance by Mysis relicta. Canadian Journal Of Fisheries And Aquatic Sciences 56(2): 311-322. Gallegos, C. L. 2001. Calculating optical water quality targets to restore and protect submerged aquatic vegetation: overcoming problems in partitioning the diffuse attenuation coefficient for photosynthetically active radiation. Estuaries 24: 381-397. Gallegos, C. L., D. L. Correll, and J. W. Pierce. 1990. Modeling spectral diffuse attenuation, absorption, and scattering coefficients in a turbid estuary. Limnology and Oceanography 35: 1486-1502. Gallegos, C. L., and K. A. Moore. 2002. Factors contributing to water-column light attenuation, p. 16-27. In R. A. Vatiuk, P. Bergstrom, W. M. Kemp, E. Koch, L. Murray, J. C. Stevenson, R. Bartleson, V. Carter, N. B. Rybicki, J. M. Landwehr, C. Ballegos, L. Karrh, M. Naylor, D. Wilcox, K. A. Moore, S. Ailstock, and M. Teichberg (eds.), Chesapeake Bay Submerged Aquatic Vegetation Water Quality and Habitat-based Requirements and Restoration 210 Targets: A Second Technical Synthsis. U.S. Environmental Protection Agency, Chesapeake Bay Program, Annapolis, Maryland. Glibert, P. M., D. J. Conley, T. R. Fisher, L. W. Harding, T. C. Malone. 1995. Dynamics of the 1990 winter spring bloom in Chesapeake Bay. Marine Ecology Progress Series 122 (1-3), 27-43. Goodrich, D. M., W. C. Boicourt, P. Hamilton and D. W. Pritchard, Wind-induced destratificaiton in Chesapeake Bay, J. Phys. Oceanogr., 17, 2232-2240, 1987. Graham, W. M., F. Pages and W. M. Hammer. 2001. A physical context for gelatinous zooplankton aggregations: a review. Hydrobiologia 451: 199-212. Greene, K. and L. Linker. 1998. Chesapeake Bay Watershed Model application and calculation of nutrient and sediment loadings-phase IV Chesapeake Bay Watershed Model-appendix a: model hydrology calibration results. EPA 903- R-98-004, CBP/TRS 196/98, Chesapeake Bay Program Office, Annapolis, MD. Gundersen, J. S., W. D. Gardner, M. J. Richardson and I. D. Walsh 1998 Effects of monsoons on the seasonal and spatical distributions of POC and chlorophyll in the Arabian Sea. Deep-Sea Research 45, 2103-2132. Haltiner, G. J., and R. T. Williams, Numerical Prediction and Dynamic Meteorology, 477pp, John Wiley, New York, 1980. Harding, L. W. Jr. 1994 Long-term trends in the distribution of phytoplankton in Chesapeake Bay: role of light, nutrients and streamflow. Marine Ecology Progress Series 104, 267-291. 211 Harding, L. W. Jr., A. Magnuson, and M. E. Mallonee. 2004. SeaWiFS retrievals of chlorophyll in Chesapeake Bay and the mid-Atlantic bight. Estuarine Coastal and Shelf Science, in press. Harding, L. W. Jr., M. E. Mallonee and E. S. Perry 2002 Toward a predictive understanding of primary productivity in a temperate, partially stratified estuary. Estuarine Coastal and Shelf Science 55 (3), 437-463. Harding, L. W., B. W. Meeson, T. R. Fisher. 1986. Phytoplankton production in two east coast estuaries: photosynthesis-light functions and patterns of carbon assimilation in Chesapeake and Delaware Bays. Estuarine, Coastal and Shelf Science 23, 773-806. Heinle, D.R., C.F. D?Elia, J.L. Taft, J.S. Wilson, M. Cole-Jones, A.B. Caplins and L.E. Cronin, 1980. Historical review of water quality and climate data from Chesapeake Bay with emphasis on effects of enrichment, USEPA Chesapeake Bay Program final report, Grant No. R806189010, Chesapeake Research Consortium, Inc. Publ. No. 84, Annapolis, MD. Rep. No. TR-002E. UMCEES 80-15CBL. Hicks, S.D., 1964. Tidal wave characteristics of Chesapeake Bay. Chesapeake science 5, 103-113. Hood, R.R., N.R. Bates, D.G. Capone and D.B. Olson. 2001. Modeling the effect of nitrogen fixation on carbon and nitrogen fluxes at BATS. Deep-Sea Research II 48, 1609-1648. 212 Hood, R. R., K. E. Kohler, J. P. McCreary, S. L. Smith. 2003. A four-dimensional validation of a coupled physical-biological model of the Arabian Sea. Deep- Sea Research II 50 (22-26), 2917-2945. Hood, R. R., H. V. Wang, J. E. Purcell, E. D. Houde and L. W. Harding Jr.. 1999. Modeling particles and pelagic organisms in Chesapeake Bay: Convergent features control plankton distributions. Journal of Geophysical Research 104: 1223-1243. Houde, E. D. and E. S. Rutherford. 1993. Recent trends in estuarine fisheries ? predictions of fish production and yield. Estuaries 16 (2): 161-176. Huntley, M. E. and M. D. G. Lopez. 1992 Temperature-dependent production of marine copepods: a global synthesis. The American Naturalist 140 (2), 201- 242. Ippen, A. T., Tidal dynamics in Estuaries, part I: Estuaries of rectangular section, in Estuary and coastline hydrodynamics, edited by A. T. Ippen, 744pp, 1966. Johnson, B. H., K. W. Kim, R. H. Heath, H. L. Butler and B. B. Hsieh 1991 Development and verification of a three-dimensional numerical hydrodynamic, salinity, and temperature model of the Chesapeake Bay, Technical Report HL-91-7, US Army Engineer Waterways Experiment Station, Vicksburg, MS. Jones, K. J. and R. J. Gowen. 1990. Influence of stratification and irradiance regime on summer phytoplankton composition in coastal and shelf seas of the British- isles. Estuarine, Coastal and Shelf Science 30: 557-567. 213 Jordan, T. E, D.L. Correll, J. Miklas and D.E. Weller, 1991. Nutrients and chlorophyll at the interface of a watershed and an estuary. Limnology and Oceanography 36(2): 251-267. Jordan, T. E, D.L. Correll and D.F. Whigham, 1983. Nutrient flux in the Rhode River: tidal exchange of nutrients by brackish marshes. Estuarine, Coastal and Shelf Science 17: 651-667. Keith, D. J., J. A. Yoder, and S. A. Freeman. 2002. Spatial and temporal distribution of coloured dissolved organic matter (CDOM) in Narragansett Bay, Rhode Island: implications for phytoplankton in coastal waters. Estuarine, Coastal and Shelf Science 55: 705-717. Kemp, W.M., P. A. Sampou, J. Caffrey, M. Mayer, K. Henriksen and W.R. Boynton, 1990. Ammonium recycling versus denitrification in Chesapeake Bay sediments. Limnology and Oceanography 35(7): 1545-1563. Kemp, W.M., P. A. Sampou, J. Garber, J. Tuttle and W. R. Boynton. 1992. Seasonal depletion of oxygen from bottom waters of Chesapeake Bay: roles of benthic and planktonic respiration and physical exchange processes. Marine Ecology Progress Series 85, 137-152. Kirk, J.T.O. 1994. Light and Photosynthesis in Aquatic ecosystems. Cambridge University Press, Cambridge. Koenings J. P., and J. A. Edmundson. 1991. Secchi disk and photometer estimates of light regimes in Alaskan Lakes ? effects of yellow color and turbidity. Limnology and Oceanography 36: 91-105. 214 Krom, M.D and R.A Burner, 1980. The diffusion coefficients of sulfate, ammonium and phosphate ions in anoxic marine sediments. Limnology and Oceanography 25: 327-337. Kundu, P. K., A numerical investigation of mixed layer dynamics, J. Phys. Oceanogr.,10, 220-236,1980. Lapointe, B. E. and M. W. Clark, 1992. Nutrient input from the watershed and coastal eutrophication in the Florida Keys. Estuaries 15 (4), 465-476. Launder, B. E., and D. B. Spalding, The numerical calculation of turbulence flows, Computer methods in applied mechanics and engineering, vol. 3, 269-289, 1974. Lee, K. Y., T. R. Fisher and E. Rochelle-Newall, 2001. Modeling the hydrochemistry of the choptank river basin using GWLF and Arc/Info: 2. Model validation and application. Biogeochemistry 56 (3): 311-348. Levitus, S., Climatological atlas of the World ocean, NOAA Professional Paper, No. 13, US Government Printing Office, Washington, DC, 173pp, 1982. Levy, H. and W. J. Moxim 1987. Fate of U.S. and Canadian combustion nitrogen emissions. Nature 328, 414-416. Li, M., A. Gargett and K. Denman. 2000. What determines seasonal and interannual variability of phytoplankton and zooplankton in strongly estuarine systems? Application to the semi-enclosed estuary of Strait of Georgia and Juan De Fuca Strait. Estuarine Coastal and Shelf Science 50 (4): 467-488. 215 Lima, I. D. and S. C. Doney, 2004. A three-dimensional, multinutrient, and size- structured ecosystem model for the North Atlantic. Global Biogeochemical Cycles 18 (3), Art. No. GB3019. Linker L. C., G. W. Shenk, P. Wang, K. J. Hopkins and S. Pokharel. 2001. A short history of Chesapeake Bay modeling and the next generation of watershed and estuarine models. U.S. EPA Chesapeake Bay Program, Annapolis, MD. Logan, J. A. 1983. Nitrogen oxides in the troposphere: global and regional budgets. Journal of Geophysical Research 88, 10785-10807. Lohrenz, S. E., A. D. Weidemann and M. Tuel. 2003. Phytoplankton spectral absorption as influenced by community size structure and pigment composition. Journal of Plankton Research 25 (1): 35-61. Madden, C. and W. M. Kemp. 1996. Ecosystem model of an estuarine submersed plant community: calibration and simulation of eutrophication responses. Estuaries 19 (2B), 457-474. Magnuson, A., L. W. Harding Jr., M. E. Mallonee, and J. E. Adolf. 2004. Bio-optical model for Chesapeake Bay and the Middle Atlantic Bight. Estuarine Coastal and Shelf Science 61, 403-424. Malone, T.C. 1992. Effects of water column processes on dissolved oxygen, nutrients, phytoplankton and zooplankton, in oxygen dynamics in the Chesapeake Bay, a synthesis of recent research, edited by D.E. Smith, M. Leffler and G.E. Machiernan, Maryland Sea Grant College, Univ. of Maryland, College Park, Maryland, pp61-112. 216 Malone, T. C., D. J. Conley, T. R. Fisher, P. M. Glibert, L. W. Harding, K. G. Sellner, 1996. Scales of nutrient-limited phytoplankton productivity in Chesapeake Bay. Estuaries 19: 371-385. Malone, T. C., L.H. Crocker, S.E. Pike and B.W. Wendler, 1988. Influences of river flow on the dynamics of phytoplankton production in a partically stratified estuary. Marine Ecology Progress Series 48, 235-249. Malone T. C. and H. W. Ducklow. 1990. Microbial biomass in the coastal plume of Chesapeake Bay: Phytoplankton-bacterioplankton relationships. Limnology and Oceangraphy 35 (2), 296-312. Malone T. C., H. W. Ducklow, E. R. Peele and S. E. Pike. 1991. Picoplankton carbon flux in Chesapeake Bay. Maine Ecology Progress Series 78 (1): 11-22. Malone, T. C., W.M. Kemp, H.W. Ducklow, W.R. Boynton, J. H. Tuttle and R. B. Jonas, 1986. Lateral variation in the production and fate of phytoplankton in a partically stratified estuary. Marine Ecology Progress Series 32: 149-160. Marra, J. and C. Ho, 1993. Initiation of the spring bloom in the northeast Atlantic (47- degrees-N, 20 degrees-W)- A numerical-simulation. Deep-Sea Research II 40 (1-2): 55-73. McClain, C. R., K. Arrigo, K. S. Tai and D. Turk, 1996. Observations and simulations of physical and biological processes at ocean weather station P, 1951-1980. Journal of Geophysical Research-Oceans 101 (C2): 3697-3713. McMahon, T. G., R. C. T. Raine, T. Fast, and J. W. Patching. 1992. Phytoplankton biomass, light attenuation and mixing in the Shannon Estuary, Ireland. Journal of Marine Biological Association of the United Kingdom 72 (3): 709-720. 217 Merrill, J.Z., 1999. Tidal freshwater marshes as nutrient sinks: particulate nutrient burial and denitrification, Ph. D dissertation, Univ. of Maryland. Meyers T., J. Sickles, R. Dennis, K. Russell, J. Galloway and T. Church, 2000. Atmospheric nitrogen deposition to coastal estuaries and their watersheds, in nitrogen loading in coastal water bodies: an atmospheric perspective, edited by Richard A. Valigura, American Geophysical Union. Miller, D. C., C. L. Muir and O. A. Hauser. 2002. Detrimental effects of sedimentation on marine benthos: what can be learned from natural processes and rates? Ecological Engineering 19 (3), 211-232. Monahan, E.C., Pybus, M. J. 1978. Colour, ultraviolet absorbance and salinity of the surface waters off the west coast of Ireland. Nature 274: 782-784. Moore, J. K., S. C. Doney, J. A. Kleypas, D. M. Glover and I. Y. Fung. 2002. An intermediate complexity marine ecosystem model for the global domain. Deep Sea Research II 49 (1-3), 403-462. Nagy, G. J., M. Gomez-Erache, C. H. Lopez and A. C. Perdomo, 2002. Distribution patters of nutrients and symptoms of eutrophication in the Rio de la Plata river estuary. Hydrobiologia 475 (1): 125-139. Newell, R. I. E., and E. W. Koch. 2004. Modeling seagrass density and distribution in response to changes in turbidity stemming from bivalve filtration and seagrass sediment stabilization. Estuaries, submitted. Nixon, S.W., 1980. Between coastal marshes and coastal waters-a review of twenty years of speculation and research on the role of salt marshes in estuarine 218 productivity and water chemistry, in estuarine and watershed processes, edited by P. Hamilton and K.B. Macdonald, pp. 437-525, plenum press, New York. Nixon, S.W. and M.G. Pilson, 1983. Nitrogen in estuarine and coastal marine ecosystems, in nitrogen in the marine environment, edited by E.J. Carpenter and D.G. Capone, academic press, New York, pp565-648. Odum, E. P., 2000. Tidal marshes as outwelling/pulsing systems, in Concepts and controversies in tidal marsh ecology, edited by M.P. Weinstein and D.A. Kreeger, Kluwer Academic Publishers, The Netherlands. Orth R. J. and K. A. Moore. 1983. Chesapeake Bay-an unprecedented decline in submerged aquatic vegetation. Science 222 (4619), 51-53. Oschlies, A. and V. Garcon, 1999. An eddy-permitting coupled physical-biological model of the North Atlantic- 1. Sensitivity to advection numerics and mixed layer physics. Global Biogeochemical Cycles 13 (1): 135-160. Pitkanen H., T. Tamminen, P. Kangas, T. Huttula, K. Kivi, H. Kuosa, J. Sarkkula, K. Eloheimo, P. Kauppila and B. Skakalsky, 1993. Late summer trophic conditions in the northeast gulf of Finland and the river Neva estuary, Baltic Sea. Estuarine Coastal and Shelf Science 37 (5), 453-474. Platt, T., and S. Sathyendranath. 1988. Ocean primary production: estimation by remote sensing at local and regional scales. Science 241: 1613-1620. Price J. F., R. A. Weller and R. Pinkel, 1986. Diurnal cycling: Observations and models of the upper ocean response to diurnal heating, cooling, and wind mixing. Journal of Geophysical Research 91 (C7): 8411-8427. 219 Richardson, M. J. 1987. Particle-size, light-scattering and composition of suspended particulate matter in the north-atlantic. Deep-Sea Research Part A 34 (8): 1301-1329. Rijstenbil, J. W. 1987. Phytoplankton composition of stagnant and tidal ecosystems in relation to salinity, nutrients, light and turbulence. Netherlands Journal Of Sea Research 21: 113-123. Risovic, D. 2002. Effect of suspended particulate-size distribution on the backscattering ratio in the remote sensing of seawater. Applied Optics 41 (33): 7092-7101. Rochelle-Newall, E. J. and T. R. Fisher. 2002. Chromophoric dissolved organic matter and dissolved organic carbon in Chesapeake Bay. Marine Chemistry 77: 23-41. Rosner, B. 1995. Fundamentals of Biostatistics, 4th ed. Duxbury Press. Sarmiento, J. L., and K. Bryan, An ocean transport model for the North Atlantic, J. Geophys. Res., 87, 394-408, 1982. Sathyendranath, S., and T. Platt. 1989a. Remote sensing of oceanic primary production: computations using a spectral model. Deep-Sea Research 36: 431- 453. Sathyendranath, S., and T. Platt, 1989b. Computation of aquatic primary production: extended formalism to include effects of angular and spectral distribution of light. Limnology and Oceanography 34: 188-198. Schubel, J.R. and D.W. Pritchard, 1986. Responses of upper Chesapeake Bay to variations in discharge of the Susquehanna River. Estuaries 9: 236-249. 220 Schubel, J.R. and D.W. Pritchard, 1987. A brief physical description of the Chesapeake Bay, in contaminant problems and management of living Chesapeake Bay resources, edited by S.K. Majumdar, I.W. Hall and H.M. Austin, the Pennsylvania Academy of Science. Seitzinger, S. P. 1988. Denitrification in freshwater and coastal marine ecosystem: Ecological and geochemical significance. Limnology and Oceanography 33, 702-724. Seitzinger, S. P. and A. E. Giblin. 1996. Estimating denitrification in North Atlantic continental shelf sediments. Biogeochemistry 35, 235-260. Sheng, Y. P. 1986. A three-dimensional mathematical model of coastal,estuarine and lake currents using a boundary fitted grid, Rep. 585. ARAR Group, Titan System, Princeton, NJ. Siegel, D. A., S. Maritorena, N. B. Nelson, D. A. Hansell and M. Lorenzi-Kayser. 2002. Global distribution and dynamics of colored dissolved and detrital organic materials. Journal of Geophysical Research-Oceans 107 (C12), Art. No. 3228. Skogen, M. D., E. Svendsen, J. Berntsen, D. Aksnes and K. B. Ulvestad, 1995. Modeling the primary production in the North-Sea using a coupled 3- dimensional physical-chemical-biological ocean model. Estuarine Coastal and Shelf Science 41 (5): 545-565. Smedstad, O. M. and J. J. Obrien, 1991. Variational data assimilation and parameter- estimation in an equatorial Pacific-Ocean model. Progress in Oceanography 26 (2): 179-241. 221 Smith, Jr., W.O. 1982. The relative importance of chlorophyll, dissolved and particulate material, and seawater to the vertical extinction of light. Estuarine, Coastal and Shelf Science 15: 459-465. Smith, R.C., B.B. Prezelin, R. R. Bidigare, and K. S. Baker. 1989. Bio-optical modeling of photosynthetic production in coastal waters. Limnology and Oceanography 34: 1524-1544. Spitz, Y. H., and J. M. Klinck, Estimate of bottom and surface stress during a spring- neap tide cycle by dynamical assimilation of tide gauge observations in the Chesapeake Bay, J. Geophys. Res-Oceans, 103, C6, 12761-12782, 1998. Staver, K.W, and R.B Brinsfield, 1996. Seepage of groundwater nitrate from a riparian agroecosystem into the Wye River Estuary. Estuaries 19: 359-370. Stefan, H.G., J.J. Cardoni, F. R. Schiebe, and C.M. Cooper. 1983. Model of light penetration in a turbid lake. Water Resources Research 19: 109-120. Stuart, V., S. Sathyendranath, T. Platt, H. Maass and B. D. Irwin. 1998. Pigments and species composition of natural phytoplankton populations: effects on the absorption spectra. Journal of Plankton Research 20 (2): 187-217. Stumm W. and J.J. Morgan, 1981. Aquatic chemistry, 2nd edition, Wiley. Taft, J. L., A. J. Elliot and W. R. Taylor, 1978. Box model analysis of Chesapeake Bay ammonium and nitrate fluxes. In: Wiley, M. L. (ed.) Estuarine interactions. Academic Press, New York, pp115-130. Teal, J.M., 1962. Energy flow in the salt marsh ecosystem of Georgia. Ecology 43: 614-624. 222 Teal, J.M. and B.L. Howes, 2000. Salt marsh values: retrospection from the end of the century, in Concepts and controversies in tidal marsh ecology, edited by M.P. Weinstein and D.A. Kreeger, Kluwer Academic Publishers, The Netherlands. Turner, R.E. and N.N. Rabalais, 1994. Coastal eutrophication near the Mississippi river delta. Nature 368: 619-621. Ullman, D. S. and R. E. Wilson, 1998. Model parameter estimation from data assimilation modeling: temporal and spatial variability in bottom drag coefficient. Journal of Geophysical Research-Oceans 103 (C3): 5531-5549. Valiela, I., M.L. Cole, J. Mcclelland, J. Hauxwell, J. Cebrian and S.B. Joye, 2000. Role of salt marshes as part of coastal landscapes, in Concepts and controversies in tidal marsh ecology, edited by M.P. Weinstein and D.A. Kreeger, Kluwer Academic Publishers, The Netherlands. Valiela, I., J. M. Teal and N. Y. Persson, 1976. Production and dynamics of experimentally enriched salt-marsh vegetation-belowground biomass. Limnology and Oceanography 21 (2): 245-252. Varela, R. A., A. Cruzado, J. Tintore and E. G. Ladona 1992 Modeling the deep- chlorophyll maximum-a coupled physical-biological approach. Journal of Marine Research 50 (3), 441-463. Venrick, E. L. 1991 Midocean ridges and their influence on the large-scale patterns of chlorophyll and production in the North Pacific. Deep-Sea Research Part A 38, S83-S102. 223 Wang, D. P. 1979a. Subtidal sea level variations in the Chesapeake Bay and relation to atmospheric forcing, J. Phys. Oceanogr., 9, 413-421. Wang, D. P., 1979b. Wind-driven circulation in the Chesapeake Bay, winter 1975, J. Phys. Oceanogr., 9, 564-572. Wang D. P. and A. J. Elliott, 1978. Non-tidal variability in the Chesapeake Bay and Potomac River: evidence for non-local forcing, Journal of physical oceanography 2: 225-232. Wang, H. V., and R. S. Chapman, Application of vertical turbulence closure schemes in the Chesapeake Bay circulation model ? A comparative study, ASCE Estuarine and coastal modeling, 283-297, 1995. Wang, M., D. R. Lyzenga and V. V. Klemas. 1996. Measurement of optical properties in the Delaware Estuary. Journal of Coastal Research 12 (1): 211- 228. Wanninkhof, R.1992 Relatioship between wind speed and gas exchange over the ocean. Journal of Geophysical Research 97, 7373-7382. Webb, K. L., 1988. Comment on ?Nutrient limitation of phytoplankton growth in brackish coastal ponds? by Caraco, Tamse, Boutros and Valiela (1987). Canadian Journal of Fishery and Aquatic Science 45: 380. Weisberg, S. 1985. Applied Linear Regression, 2nd ed. John Wiley & Sons, Inc. Wetzel, R. and H. Neckles. 1986. A model of Zostera marina L. photosynthesis and growth: simulated effects of selected physical-chemical variables and biological interactions. Aquatic Botany 26, 307-323. 224 Wu, C. R., P. T. Shaw and S. Y. Chao, 1999. Assimilating altimetric data into a South China Sea model. Journal of Geophysical Research-Oceans 104, C12, 29987- 30005. Xu, J., S.Y. Chao, R.R. Hood, H.V. Wang and W.C. Boicourt 2002 Assimilating high- resolution salinity data into a model of a partially mixed estuary, J. Geophys. Res. 107 (C7), 3074, doi, 10.1029/200JC000626. Xu, J., R. R. Hood and S. Y. Chao. 2004. A simple empirical optical model for simulating light attenuation variability in a partially mixed estuary. Estuaries. Submitted. 225