Report to MD Dept. of Environment

Development of an integrated modeling system of watershed and estuarine ecology for management of the Patuxent basin

Thomas R. Fisher, lead PI, HPL
Michael Kemp, HPL
Raleigh Hood, HPL
Michael Williams, HPL
Gregory Radcliffe, HPL
Walter Boynton, CBL
Keith Eshleman, AL
Dan Fiscus, AL
Carl Cerco, ERDC
Sung Chan Kim, ERDC

Table of Contents
Introduction
Progress Report
Bathymetry of the Patuxent
Estuarine Transect Line
Estuarine Modeling Grid
Marshes in the Patuxent
Watershed Model Validation
Weather Data
Patuxent Nutrient Budgets
Diffuse Source Models
Miscellaneous
References

Introduction (back to table of contents)
The University of Maryland Center for Environmental Science (UMCES), jointly with the Waterways Experiment Station (WES), is developing a coupled modeling system for nutrient management of the Patuxent basin. Funds for this project became available from MD Department of the Environment (MDE) in January 2002, and this report covers the first year of progress on this project (Jan. 2002 - Dec. 2002). The main goal of this research effort is to deliver to MDE by December 2004 a watershed model (HSPF) spatially linked to a 3-dimensional estuarine circulation model (CH3D) and water quality model (CE-QUAL-ICM). This coupled modeling system of the Patuxent will be used by MDE to develop TMDLs in the first half of 2005 with assistance from UMCES and WES. All models are those in current use by the EPA Chesapeake Bay Program (CBP) using existing Bay Program, USGS, and other available data.

The original intent of this project was an independent calibration of the watershed and estuarine models using all available data. Following calibration to minimize errors in model output, validation with a reserved set of data was proposed to be done to generate independent estimates of model errors in segments of the terrestrial basin. However, at the request of MDE in the summer of 2002, six months after the project began, MDE expressed concern about the legal ramifications of more than one set of watershed model output (CBP and MDE), and they requested that we use exactly what the CBP uses for watershed model output. We have agreed to work only with the existing watershed model output, and to focus our efforts on validation of the watershed model and on the investigation of possible improvements to future versions of the watershed model (specifically Phase 5). More details on these changes are summarized in a recent letter to MDE, which initiated the development of an addendum to the MOU between MDE and UMCES that describes these changes requested by MDE.

The estuarine models will be run as described in the original proposal. We have improved the original modeling grid with a higher resolution version (see below) in order to improve the simulation of the effects of bathymetry and shoreline morphology. The new estuarine grid is embedded in the current grid of the entire Bay in order to simulate and estimate the role of exchange with the mainstem Bay as a net nutrient source or sink for the Patuxent. To support the modeling efforts, watershed and estuarine data are being compiled in GIS databases as needed for modeling and mapping of model output.

In addition to the main goal of delivering the coupled modeling system to MDE, we also have four research objectives. After we have validated the models in routine use by the Chesapeake Bay Program, we will do the following

• compare these models with alternative approaches (watershed: landuse-specific yield coefficients, the Patuxent Landscape Model, GWLF, SPARROW; estuary: box and statistical models).
• compare the importance of various sources of N and P in the basin (direct atmospheric deposition, diffuse and point sources, main stem Bay advective inputs).
• investigate the processing of the watershed inputs within the estuary (conversion of inputs to phytoplankton and sediments).
• explore net estuarine retention within the Patuxent vs. export to the mainstem Bay.

Progress Report (back to table of contents)
This progress report consists of a series of task descriptions. These have been conducted collectively and by individual principal investigators and associates and have been underway for the last year in order to achieve the main goal of this project. The descriptions are neither complete nor comprehensive, and the material shown is a sample of the activities which are underway to give examples of some of the project accomplishments.

Bathymetry of the Patuxent (back to table of contents)
One of the first tasks undertaken in the project was the construction of a bathymetric database for the Patuxent. Using data collected by MDE in August, 2001 (prior to the start of the project), as well as NOAA point soundings, we created a point coverage and grid file in ArcGIS. The point coverage was the actual sounding data with spatial coordinates, whereas the grid file was extrapolated from the point coverage to provide a spatially continuous and smoothed depth field for the development of the new estuarine voxel grid for estuarine modeling (see below). The point coverage of the original data is illustrated below in Fig. 1.

Estuarine Transect Line (back to table of contents)
An estuarine transect line is a convenient reference for location within the estuary in terms of distance from the mouth. We created an arc coverage which ran through the deepest parts of the estuary from the mouth to above Western Branch. The line was segmented at 1 km intervals from the mouth to km 105 in the little Patuxent River. The origin at the mouth of the estuary corresponds to the origin used by the estuarine circulation and water quality models, and also corresponds to the origin used by Cronin and Pritchard (1975). Fig. 1 shows the estuarine transect line, with the 1 km points along the line.

Estuarine Modeling Grid (back to table of contents)
One of the more important improvements to the estuarine modeling was the development of a more detailed voxel grid. The goal of this effort was to attempt to match the model voxels with the true bathymetry and shoreline as much as possible. The voxel grid in use for the Patuxent at the beginning of this project was developed from the perspective of Chesapeake Bay as a whole, and had relatively low resolution (136 by 96 by 19). The spatial fidelity relative to shorelines and bathymetry was low when viewed from within the Patuxent, particularly in the upper estuary (see old grid in upper estuary in Fig. 2). The bathymetry data described above and a shoreline arc file created from USGS topographic maps were used to develop a higher resolution voxel grid (152 by 166 by 19) which preserved more of the true shoreline and bottom bathymetry, and which also incorporated more of the tributary structure of the Patuxent, including Western Branch in the upper estuary. Note the additional tributaries and the more realistic rendition of Point Patience in the lower Patuxent (lower right corner of Fig. 2), as well as the more accurate representation of the upper Patuxent (upper right in Fig. 2). Not shown in Fig. 2 is the more realistic bathymetry which is rendered by subsurface voxels by extracting data from the NOAA NGDC Coastal Relief Model into the new grid. It is hoped that the more realistic bathymetry and shoreline in the voxel grid will improve the simulation of the estuarine circulation and mixing, and initial model runs have already been made successfully with the improved grid.

Marshes in the Patuxent (back to table of contents)
Marshes are an important landscape component that have not been included in previous versions of the Bay Program’s models. Intertidal and non-tidal marshes were not considered as landscape sinks for TSS, N, and P in either the watershed or estuarine modeling. The omission of intertidal marshes from the estuarine water quality model may be responsible for the significant overestimation of simulated N and P concentrations compared to those observed in the upper estuary where the marshes are very abundant (see also nutrient budget section below). In order to include these important landscape components in the coupled Patuxent modeling system, we have created a detailed polygon coverage of wetlands from USGS topographic maps. Three kinds of polygons were defined: (1) those which exchange laterally with the estuarine voxels, (2) those which occur at the end of a creek and exchange with the end of a voxel, and (3) those which are surrounded by land and interact with the land only (see Fig. 3). In addition to the attributes normally assigned by ArcGIS, each intertidal marsh polygon has attributes consisting of the adjacent estuarine voxel id with which it exchanges water and materials as well as the number of the nearest estuarine km from the estuarine transect line described above in Fig. 1. In 2003 we will complete the marsh polygon coverage, adding additional attributes of NWI class, dominant species, CNP plant biomass, CNP burial rates, and denitrification rates, using both observed and extrapolated data.

Watershed Model Validation (back to table of contents)
The modification of the objectives of watershed modeling have changed our primary focus to validation of the watershed model. Below we give some preliminary examples of the current validation errors for water discharge, total N, total P, and total suspended solids in the phase 4.3 watershed model.

Simulating water discharge accurately is essential for the coupled modeling system. Using data from Chesapeake Bay Program model version 4.3 obtained from Gary Shenk, and USGS data, we compared model output to observed flow on a daily basis. Scatter plots of model vs observed flow for Bowie and Laurel USGS stations are shown in Figs. 4 and 5. The model output for the Bowie segment shows an unbiased relationship to observed flow (Fig. 4), with rms errors of " 82% (ave. root-mean-square difference between predicted and observed relative to observed). This error reflects the scatter about the 1:1 line in Fig. 4) and indicates that on a daily basis most of the predicted flow values are within a factor of two of the observed flows.

The Laurel model output agrees less well with observations (Fig. 5). There is more scatter and many zero values for flow (not shown, but about 1400 out of 5100 daily data points). This difference between stations is probably due to the fact that Bowie was one of the calibration stations, whereas Laurel was not, which makes Fig. 5 a true validation test of the model output.

Daily data were also aggregated to annual runoff for 1984-1997. At Bowie, this comparison showed good agreement between the average observed and predicted runoff (36 vs 38 cm 1/y, respectively), with annual rms errors of + or - 18% (see Fig. 6). This indicates that model precision improved from daily predictions (rms errors = 82%) to annual predictions (rms errors = 18%), as model errors tend to cancel at longer time scales. For the entire 14 year period (1984-1997), the cumulative error declined still further to 6%. For Laurel, the average observed and predicted annual runoff values were 21 and 17 cm 1/y, respectively, with annual rms errors of 38% and a cumulative error of 19%. The comparison of annual runoff values suggests no substantial model bias for the Bowie station (approximately equal scatter about the 1:1 line in Fig. 6), with 18% rms errors. However, for Laurel there was a negative bias with larger rms errors of 38% (Fig. 6).

The Bowie segment is the only one with no predicted flow values of zero. For Bowie, spectral analysis was done on both the observed and predicted daily time series for the 1984-1997 period (Fig. 7). The results indicate that the phase 4.3 model does a reasonable job of matching the frequency domain and scaling behavior of the actual data; however, the biggest difference was in the high frequency portion of the spectrum (time scales <7 days, right vertical line in Fig. 6 at log f = -0.85), where the spectral slope for the observed data was 1.98 and the model was 1.47. This suggests less attenuation of higher frequency variance in the model relative to actual, which may be due to a hydrograph recession constant which is too small or a time resolution problem. The model output also does not have the same peak at the 1-year period as the actual data (left vertical line in Fig. 6 at log f = -2.56). We plan additional work to clarify the results of the frequency analysis, identify the causes, and propose corrections to the phase 5 watershed model.

Simulation of N concentrations in streams is another essential feature of watershed modeling. To validate the model output, we have compared the daily average predicted total N (TN, mg L-1) at Bowie with water chemistry observations by USGS in grab samples (Fig. 8). In this case, the data points scatter about the 1:1 line, indicating little apparent bias and rms errors of 43%. This can be seen more clearly in Fig. 9, which is a frequency distribution of the % difference between predicted and observed TN values relative to observed TN at Bowie. Errors are approximately log normally distributed and range from -70% (low predicted values) to +150% (high predicted values). The median error is only 0.7%, but the mode is -15%, indicating a slight tendency of the model to underestimate TN. Several measures of validation errors for TN are given in Fig. 9, but it is clear that most predicted values fall within "50% of the observed ones.

TP and TSS estimates by the phase 4.3 watershed model were weaker than for TN. In Fig. 10, the comparison of modeled and observed TP shows evidence of model bias (slope of 0.33) and considerable scatter (rms errors of 68%). As for TN (Fig. 9), the % errors for TP were log normally distributed and ranged from -90% to +360%, with most falling between "100% (not shown). The mode of the distribution was -25%, and median errors were -8%, both again suggesting a tendency to underestimate TP. For TSS (total suspended solids), even larger errors were found (Fig. 11). There was more scatter for TSS (rms errors of 240%), with evidence of model bias (slope of 0.50). In a large fraction of the paired values in Fig. 11, there is a large range of observed TSS values (2-500 mg L-1) for which modeled values range only over 7-30 mg L-1. This may be indicative of mismatches in the timing of discharge events (e.g., modeled water discharge lags or precedes observed discharge), or some other model inaccuracy. We plan to investigate the causes of the apparent biases in modeled values of both TP and TSS, and for phase 5 we hope to provide suggestions to reduce the scatter in all predicted watershed values compared to observations (Figs. 4-11).

Weather Data (back to table of contents)
Initially we were compiling weather and other input data for the watershed calibration. This activity occurred only in the first half of 2002 prior to the MOU change described above. While these data will not now be used, a considerable archive was acquired, particularly for weather data. A total of 25 stations for precipitation, temperature and other climate data were identified (Fig. 12). The 4 WBAN stations (excluding Sterling, VA) are the ones that have all the data required by the Chesapeake Bay HSPF model. The other co-op stations may also have useful data for precipitation and temperature.

Patuxent Nutrient Budgets (back to table of contents)
For the past several years efforts ave been underway to develop updated nutrient budgets for the Patuxent River estuary. These efforts follow completion of an N and P budget published by Boynton et al. (1995) based largely on information available for the late 1980’s. Results from that effort indicated the following:

• the Patuxent was moderately loaded with both N and P relative to other estuarine systems
• both point and diffuse sources were important
• direct atmospheric deposition of N was less important
• burial and denitrification were
  important N losses
• burial in sediments was the primary loss term for P

export to the Chesapeake Bay was small for N and almost zero for P

As part of the Patuxent TMDL project, these budgets are being updated. Using a nutrient budget framework, this re-evaluation uses the following new information: (1) the 1985-2000 time series of input data which contains both wet, normal and dry years; (2) measurements of long-term nutrient burial (based on 210Pb profiling) and denitrification (based on the new N2 /Ar method) in subtidal sediments; (3) an evaluation of N and P losses due to burial and denitrification in the tidal marshes developed by Merrill (1999); (4) a box model of the tidal Patuxent that enables estimates of nutrient transport at key portions of the system (Hagy et al. 2000); (5) availability of several watershed models estimating diffuse source nutrient inputs; (6) a system-wide “experiment” of attempted N reductions, concluded with the completion of biological nitrogen reduction (BNR) technology at all of the major sewage treatment plants in the basin, which allows us to see where and how much loads have been reduced.

This activity is still underway, and here we only present a summary of TN loads and a preliminary examination of diffuse source load estimates from three different approaches. Average annual TN inputs from all sources were organized for the period 1985-2000 and the years having the highest and lowest overall TN loading rates were extracted and results depicted in Fig. 13. The loads and estimated transport from the middle to lower estuary and from the estuary to Chesapeake Bay are shown for the years with the lowest (1991) and highest (1996) loads. The fall line load was compiled by USGS measurements of flow and nutrient concentration. This approach uses statistical modeling and includes point, diffuse and direct atmospheric deposition of TN to the river surface. The loads to the middle estuary included direct atmospheric deposition to estuarine surface waters, estimated septic inputs, point source inputs, and all other diffuse source inputs as estimated from the Chesapeake Bay Program land use model. Inputs to the lower basin were the same as for the middle basin except that there were no significant point source inputs. Transport from the middle to the lower basin and from the lower basin to Chesapeake Bay was estimated using the box model referenced above.

There was just over a factor of two difference between the lowest and highest load years on record. This represents a very substantial difference, considering that nutrient management goals for Chesapeake systems generally aim for smaller reductions than indicated here for inter-annual variability. The loads at the fall line, middle and lower basins were all higher during the wet year of 1996, but especially so in the middle basin. These increased loads occurred after full implementation of BNR at the major sewage treatment plants indicating that diffuse sources are of key importance, even in a basin typically thought to be point source dominated. This finding is further supported by the observation that the lowest loads occurred in 1991, prior to full implementation of BNR technology.

Losses within the estuary were significant in this new budget. Of the nitrogen that entered the tidal Patuxent, a relatively small proportion (21-23%) was transported out of the estuary into Chesapeake Bay, with little difference between wet and dry years. However, there were very significant N losses during transit from the head of tide (river km 75, see Fig. 1) to the beginning of the mesohaline estuary (river km 40). During the dry year about 34% of all N entering the estuary upstream of the mesohaline zone was removed, while in the wet year 45% of all inputs upstream of the mesohaline zone were removed. These losses appear to be related to denitrification and long-term burial of N in both sub-tidal and tidal marsh sediments. Direct measurements of these losses have yet to be summarized, but preliminary estimates suggest that measured rates are sufficient to account for these large estimated losses in the budget.

Diffuse Source Models (back to table of contents)
One of the key issues related to the preparation of a TMDL concerns accurate estimation of nutrient loading rates. For some of the sources typical of estuarine situations, estimates are either quite accurate (e.g. point source discharges) or are relatively small (e.g., direct atmospheric deposition to surface waters ) and therefore not centrally important. In cases where most of the drainage basin is above the head of tide, excellent estimates can be developed from high frequency flow and concentration measurements (e.g., USGS data sets). However, in those cases where a significant portion of the estuarine drainage basin is below the head of tide, other methods need to be utilized to estimate nutrient inputs. One of the most common methods is to use a watershed model to compute loads.

In the case of the Patuxent, about 60% of the drainage basin is located below the head of tide. Hence there is an important need for modeling nutrient inputs from these tidally influenced, primarily coastal plain portions of the basin. As is the case with numerous other variables of environmental interest, the Patuxent has a relatively long history of diffuse source modeling. Models have been developed by the Chesapeake Bay Program (HSPF model), University of Maryland scientists (Costanza et al. 2001, spatially explicit, mechanistic model), USGS (statistical models of load at the head of tide), and researchers from the Smithsonian Environmental Research Center (a large collection of small watershed monitoring with statistical modeling of loads). Several other models have been initiated in the past, but current status of those modeling activities or results are not known at this time.

To compare the load estimates from different investigations, TN loads for the head of tide (e.g., Bowie, MD) from three programs are summarized on a seasonal basis in Fig. 14. There appears to be general correspondence among all three estimates in terms of seasonal pattern with highest loads associated with winter or spring and smallest loads associated with summer or fall. The coherence between the empirical USGS load estimate and the CBP model results is tighter than with the Costanza et al. (2002) results. In general, the Costanza model estimates higher loads than either of the other models. While on the one hand the differences between these approaches raises a certain amount of concern, we plan to gain a more comprehensive understanding of processes regulating diffuse source inputs to this estuarine system by examining these approaches more closely in this project.

Miscellaneous (back to table of contents)
A research project with nine scientists at four locations generates additional research activities as well as administrative overhead. Brief summaries of some of the more pertinent ones are provided below.
An additional activity was created when a coalition of land owners from Island Creek in the Patuxent approached us concerning water quality problems. They are concerned about increased boat traffic (especially jet skis and ski boats associated with a waterfront bar) and the effects on shoreline erosion and water clarity in Island Creek. With some ideas from us, the landowners applied for and received WRAS funds from DNR to investigate this issue. This offshoot of our main project will get underway this summer.
There is also collaboration between this project and one headed by Dr. Tom Jordan at the Smithsonian Environmental Research Center (SERC). A letter of endorsement was given to Jordan for his proposal to NOAA, which addressed nutrient inputs to the Patuxent and which also included funds for workshops to compare various sources of data, as in Fig. 14 above. This effort was recently funded, and we anticipate considerable collaborative efforts between the two projects over the next two years.

Funds within our existing budget were internally rearranged to create a postdoctoral research position. After an extensive search, Dr. Michael Williams of the Ecosystem Center at Woods Hole was hired as of January, 2003. Michael commutes to Horn Point Lab for a week each month, and will be in residence full-time in June (after his daughter completes first grade). His primary activities are validation of the watershed and estuarine models.
To maintain communication between project participants, we have held quarterly project meetings. We have met both in person at meetings of opportunity and on the UM Interactive Video Network (IVN) on the following dates: Nov. 6, 2001 (ERF meeting); Dec. 3, 2001 (CBP Modeling Subcommittee meeting); Dec. 17, 2001 (IVN); Jan. 23, 2002 (IVN); May 15, 2002 (IVN); Aug. 29, 2002 (IVN); and Dec. 20, 2002 (IVN). The next project meeting is planned for Mar. 24, 2003 on IVN. Summaries of all meeting have been previously provided to MDE, along with brief monthly activity reports through the end of 2002.

References (back to table of contents)
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: 285-314

Costanza, R,. A. Voinov, R. Boumans, T. Maxwell, F. Villa, L.Wainger, and H. Voinov. 2002. Integrated ecological economic modeling of the Patuxent River watershed, Maryland. Ecological Monographs 72: 203-231

Cronin, W. B. and D. W. Pritchard. 1975. Additional statistics on the dimensions of the Chesapeake Bay and its tributaries: cross-section widths and segment volumes per meter depth. CBI Spec. Rep. 42 Johns Hopkins University, 475 pps.

Hagy, J. D., W. R. Boynton, and L. P. Sanford. 2000. Estimation of net physical transport and hydraulic residence times for a coastal plain estuary using box models. Estuaries 23: 328-340

Merrill, J. Z. 1999. Tidal freshwater marshes as nutrient sinks: particulate nutrient burial and denitrification. PhD thesis, University of Maryland, 342 pps.