Environ. Sci. Technol. 2003, 37, 2043-2050
Assessing TMDL Effectiveness Using Flow-Adjusted Concentrations: A Case Study of the Neuse River, North Carolina C R A I G A . S T O W * ,† A N D MARK E. BORSUK‡ Environmental Science and Policy Division, Nicholas School of the Environment and Earth Sciences, Duke University, Durham, North Carolina 27708
Integrated control of both point and nonpoint source water pollution using Total Maximum Daily Load (TMDL) assignments will be a major regulatory focus over the next decade. We propose the use of “flow-adjusted” pollutant concentrations to evaluate the effectiveness of management actions taken to meet approved TMDLs. Pollutant concentrations are usually highly correlated with streamflow, and flow is strongly weather-dependent. Thus, pollutant loads, which are calculated as pollutant concentration multiplied by streamflow, have a large weather-dependent variance component. This natural variation can be removed by calculating flow-adjusted concentrations. While such values are not a direct measure of pollutant load, they make it easier to discern changes in streamwater quality. Additionally, they are likely to be a better predictor of pollutant concentrations in the receiving waterbody. We demonstrate the use of this technique using longterm nutrient data from the Neuse River in North Carolina. The Neuse River Estuary has suffered many eutrophication symptoms, and a program to reduce nutrient loading has been in place for several years. We show that, in addition to revealing recent reductions in nutrient inputs, annual flow-adjusted riverine nutrient concentrations show a more pronounced relationship with estuarine nutrient concentrations than do annual nutrient loads. Thus, we suggest that the calculation of flow-adjusted concentrations is a useful technique to aid in assessment of TMDL implementation.
Introduction Since the inception of the 1972 Clean Water Act, the emphasis of water quality regulation in the U.S. has been on technology-based, point-source controls. However, improvements associated with point-source controls have approached an asymptote that falls short of meeting water quality standards in many waterbodies. Thus, emphasis is shifting toward achievement of standards by establishing * Corresponding author phone: (919)684-8741; fax: 919-613-8105; e-mail:
[email protected]. † Address as of August 2003: Department of Environmental Health Sciences, Arnold School of Public Health, University of South Carolina, Columbia, SC 29208. ‡ Current address: Department of Systems Analysis, Integrated Assessment, and Modeling (SIAM), Swiss Federal Institute for Environmental Science and Technology (EAWAG), P.O. Box 611, 8600 Du ¨ bendorf, Switzerland. 10.1021/es020802p CCC: $25.00 Published on Web 04/12/2003
2003 American Chemical Society
Total Maximum Daily Loads (TMDLs) to reduce pollutant inputs from both point and nonpoint sources. TMDL development requires estimation of the maximum pollutant load (loading capacity) that an impaired waterbody can accept and still meet the applicable water quality standards. The Environmental Protection Agency (EPA) estimates that upward of 40 000 TMDLs must be completed over approximately the next 10 years to meet currently imposed deadlines (1), indicating that considerable resources will be directed toward this effort. Once a TMDL is established and watershed-based pollutant controls are initiated, there is likely to be general interest in evaluating the actual achievement of pollutant reductions. Indeed, adaptive implementation of TMDLs will require the periodic assessment of the systems’ response to management actions (1). In practice, the term Total Maximum Daily Load is a misnomer, and the effectiveness of pollutant management actions is likely to be evaluated using loads estimated on an annual, rather than a daily, basis. However, even on an annual scale, discerning pollutant reductions using load estimates may be difficult, as pollutant loads have a large weather-driven component of variation. This natural variation arises because load is calculated as the product of concentration and streamflow, and pollutant concentration is usually highly correlated with flow (2). To address this point, we propose the calculation of annual flowadjusted concentrations as useful indicators of temporal changes in pollutant inputs to a waterbody. In addition to making trends easier to discern, such values are likely to show a stronger relationship with pollutant concentrations in the receiving water. To demonstrate the calculation and interpretation of flowadjusted concentrations, we analyzed long-term nutrient data from the Neuse River and Estuary in North Carolina. The Neuse is an example of a coastal system that has been severely impaired by nutrient overenrichment, or eutrophication. Problems associated with eutrophication include low dissolved oxygen over large areas (3), fishkills, and concerns regarding toxic microorganisms (4). The issue of coastal eutrophication and its management using TMDLs has received increased national attention in recent years (5, 6), and the Neuse River has been prominent in this discussion. Concern about nutrient loads in the Neuse River watershed began in the 1980s and led to a phosphate detergent ban that became effective in 1988 (7), resulting in a pronounced step-decrease in phosphorus concentrations throughout the watershed (8). In the 1990s, a series of large fishkills in the estuary stimulated public and political interest in undertaking additional management actions to control watershed nutrient loading. A disagreement has arisen regarding the specific cause of the fishkills, with some favoring the formation of bottom-water hypoxia as the cause (9) and others promoting toxic microorganisms as the dominant antecedent (10). However, while this debate ensues, there is general agreement that the ultimate culprit is nitrogen loading, with nonpoint sources dominating. In 1997, the North Carolina Division of Water Quality developed the Neuse Nutrient Sensitive Waters Management Strategy (commonly referred to as the “Neuse Rules”) that sought to reduce total nitrogen loading to the Neuse Estuary by 30% by the year 2003. Because the Neuse Estuary is federally listed as an impaired waterbody under section 303(d) of the Clean Water Act, a TMDL assignment process for nitrogen has been required by the U.S. EPA. The final version, developed by the North Carolina Division of Water Quality and approved by VOL. 37, NO. 10, 2003 / ENVIRONMENTAL SCIENCE & TECHNOLOGY
9
2043
for loss in transit via denitrification (13, 14) make it desirable to estimate loading as close to the estuary as possible. Previously, this downstream estimation has been limited by a lack of streamflow data below the confluence of Contentnea Creek and the Neuse River mainstem (Figure 1). However, in 1996 the USGS instituted flow-monitoring at Fort Barnwell, below the confluence and approximately 20 km above Streets Ferry, where the river widens forming the estuary (Figure 1). Now, with over 4 years of daily flow-monitoring data at Fort Barnwell as well as more than 30 years of daily flow data at upstream stations on Contentnea Creek at Hookerton and on the Neuse main stem at Kinston (Figure 1), it is possible to back-estimate the long-term flow at Fort Barnwell and expand the historical record of load estimations. Load Estimation. To estimate daily flow at Fort Barnwell back through 1979, we used Kinston and Hookerton daily average flow data from October 1996 through December 2000 to fit the following model
log(flowfb) ) β0 + β1log(flowh) + β2log(flowk) +
FIGURE 1. Neuse River watershed map depicting the Hookerton, Kinston, and Fort Barnwell monitoring stations (top panel) and the estuary depicting sampling stations 1-5 (bottom panel). the EPA in 2001, also requires the achievement of a 30% reduction in nitrogen load relative to 1991-1995 levels. Despite the recent implementation of nutrient controls in the Neuse River watershed, corresponding declines in nitrogen loading to the Neuse estuary are not readily apparent, as we will show. This has been discouraging news for politicians, the public, and the regulated community who are eager to see the results of their actions. However, the inability to detect nitrogen load changes may be due, in part, to recent weather conditions, including two hurricanes and an extended drought, that have caused large fluctuations in streamflow. Therefore, we investigate the potential for discerning trends in flow-adjusted concentrations, which remove the component of variation related to streamflow. We also investigate the relationship between nutrient inputs, measured as both loads and flow-adjusted concentrations, and the corresponding nutrient concentrations in the Neuse Estuary. It is not our objective to find the best possible predictor of estuarine nutrient concentrations but rather to demonstrate that removing flow-variability has advantages for evaluating the ability of TMDLs to achieve reductions in both pollutant inputs and associated receiving water concentrations.
Study Site Description and Statistical Methods Site Description. The Neuse River originates northwest of Raleigh, NC, flows approximately 320 km through the central piedmont to the coastal plain and comprises a watershed area of 16 108 km2 (Figure 1). The major land types in the basin are agriculture (35%) and forests (34%) with the remainder consisting primarily of developed areas, wetlands, and open water. In 1983, the U.S. Army Corps of Engineers completed a dam that impounded the upper 35 km of the river to create Falls Lake, a multipurpose recreation and drinking water reservoir (11). Historical nitrogen and phosphorus loads have been previously calculated for several monitoring stations on the mainstem of the Neuse River in upper parts of the watershed (12), but the nonconservative nature of nitrogen and potential 2044
9
ENVIRONMENTAL SCIENCE & TECHNOLOGY / VOL. 37, NO. 10, 2003
(1)
where flowfb ) Fort Barnwell flow, flowh ) Hookerton flow, and flowk ) Kinston flow, all in m3/day. Due to the temporal density of the observations, residuals from the regression exhibited extreme temporal autocorrelation, as indicated by graphical examination, calculating correlations among lagged residuals, and fitting the residuals in an autoregressive model. The main effect of this autocorrelation was that it made the model standard error estimates inappropriately small. To remove the autocorrelation we eliminated observations at successively higher intervals until the correlation between lagged residuals consistently exceeded a p-value of 0.05. During this process, the model parameter estimates changed little and the R2 and mean squared error changed slightly. The primary effect was that the estimates of the model standard errors increased as the autocorrelation among residuals decreased. For the final model, only every 14th set of daily observations was included. Daily flow estimates in the natural metric were obtained from predicted logtransformed values via exponentiation, with the application of a retransformation bias correction (15). The NC DENR, or its predecessor agencies, began periodic monitoring of nutrients and other water quality characteristics in the Neuse watershed in the late 1960s. By 1979 this had developed into a regular ambient monitoring program with approximately monthly sampling at 16 stations in the Neuse River mainstem (Figure 1). Corresponding routine monthly sampling began at two sites on Contentnea Creek, the other main tributary in the watershed, in 1985 (Figure 1). In 1996, the frequency of sampling was increased to almost daily at both the Kinston and Fort Barnwell stations. To estimate nutrient loads at the Hookerton, Kinston, and Fort Barnwell stations, we used available concentration data and daily flow data to develop concentration:flow regression models for total Kjeldahl nitrogen (TKN), nitrate (nitrate + nitrite), total nitrogen (nitrate + TKN), and total phosphorus. We estimated models at each station as
log(C) ) βyear(year) + β(log flow) +
(2)
where C is nutrient concentration, βyear and β are model intercept and slope parameters, respectively, year is a categorical variable, and is the model disturbance term. This model, analogous to an analysis of covariance, allowed the intercepts to differ by year but maintained a constant slope across years. Almost daily concentration measurements at the Hookerton, Kinston, and Fort Barnwell stations from 1996 to 2000 resulted in considerable temporal autocorrelation among model residuals. The residuals from these models exhibited significant lags of up to 21 days during this
period. To minimize the effect of the autocorrelation, we used the procedure described above for the Fort Barnwell flow model, resulting in the use of data from every 14th day (from 1996 to 2000) to estimate the concentration:flow models. We experimented with models that included terms for seasonality and changes during rising or falling hydrographs, but these additional terms made negligible differences in the annual load estimates and were not included in the final models. Because a large proportion of the Fort Barnwell flows was estimated rather than measured, we considered the possibility that an “errors-in-variables” effect could result in attenuated slope estimates for the Fort Barnwell concentration:flow models (16). However, a comparison of the prediction variances from eq 1 to the marginal variance of Fort Barnwell flows indicated this effect was inconsequential. Daily concentrations in the natural metric were obtained via exponentiation with a bias correction applied (15). Daily concentration prediction variance in the natural metric was estimated as
(
(
σ2c ) exp[2log(C) + σµ2] exp[σµ2] 1 -
)
2σ2 m
TABLE 1. Summary Statistics for the Nutrient Concentration:Flow Models station
n mse R2 Model F slope (s.e.) Kinston n mse R2 Model F slope (s.e.) Fort Barnwell n mse R2 Model F slope (s.e.)
Hookerton
total total Kjeldahl total phosphorus nitrogen nitrate nitrogen 217 0.213 0.53 4393 -0.240 (0.029) 304 0.089 0.72 11645 -0.241 (0.019) 314 0.128 0.55 7991 -0.177 (0.022)
218 0.127 0.22 5738 -0.017 (0.022) 305 0.089 0.23 9014 0.040 (0.019) 312 0.096 0.36 8442 0.011 (0.020)
218 0.116 0.46 5725 -0.139 (0.021) 305 0.322 0.43 2204 -0.357 (0.036) 312 0.149 0.38 4714 -0.222 (0.024)
218 0.059 0.50 9493 -0.082 (0.015) 305 0.053 0.55 11403 -0.170 (0.014) 312 0.053 0.42 11542 -0.124 (0.014)
-m/2
(
-
1-
))
2 -m
σ m
where log(C) is the predicted concentration in the log metric, σµ2 is the prediction variance of the mean (log metric), σ2 is the variance of the error term, and m is the error degrees of freedom (15). We estimated daily loads by multiplying flow and concentration, using measured values when they were available and estimated values when there were no measured values. When estimated concentrations were used to calculate daily load, we estimated daily load variance as σ2c × (flow)2. The prediction variance for estimated Fort Barnwell flows was sufficiently low that it was not included in the load variance calculations. We estimated annual total load and annual load variance by summing the estimated daily loads and variances, respectively. Flow-Adjusted Concentrations. The concentration:flow models estimated with eq 2 provide the basis for deriving flow-adjusted concentrations. Conceptually, eq 2 represents log concentration vs log flow as a set of parallel lines, with each line representing a different year. If concentrations are changing with time, the yearly intercept estimates (βyear) will change. This result can be shown graphically by plotting the intercepts vs year, or, because the lines are parallel, by choosing any reference flow and graphing the associated yearly concentration estimates and their uncertainty. We chose to present these flow-adjusted concentrations at the median daily flow from 1979 to 2000 for each site. This was a choice made for statistical reasons and not because the value has any physical significance. The long-term median was chosen as a reference flow because it is invariant to log transformations and falls near the middle of the regression lines where parameter estimation uncertainty will be near the minimum. Relationship to Estuarine Concentrations. To represent nutrient conditions in the estuary, we selected five of the stations sampled by NC DENR that fully span the longitudinal salinity gradient (Figure 1). Median annual estuarine concentrations at these stations were compared to both annual riverine nutrient loads and flow-adjusted concentrations at Fort Barnwell to evaluate which measure of river input is most closely related to conditions in the estuary. This comparison was done both graphically and using simple linear regression. For the regression analysis, we calculated
p-values representing a test of H0, stating that the estimated regression slopes equal zero, versus H1, the optimum slope estimates. All statistical analyses were done using SAS software.
Results The model selected to estimate Fort Barnwell flow from flow at Hookerton and Kinston is
log(flowfb) ) 0.822 + 0.132 log(flowh) + 0.828 log(flowk) + with parameter standard errors ) (0.180) (0.041) (0.052), respectively, R2 ) 0.97, n ) 111, mean squared error ) 0.032, and model F ) 1569. Flow conditions from 1996 to 2000, the period for which this model was estimated, covered almost the full range of observed flow conditions from 1979 to 2000 so that only a small amount of extrapolation was needed at low flows. Average Kinston flow is approximately four times average Hookerton flow and the greater dependency of flow at Fort Barnwell to Kinston flow is reflected in the estimated slopes. Results of the concentration:flow model indicate that all nutrients except TKN exhibit an inverse flow relationship at all three river sites (Table 1). TKN shows no discernible flow relationship at any of the sites, implying an increasing ratio of reduced to oxidized nitrogen forms with increasing flow. The total phosphorus slope estimates are similar among sites, while the nitrate and total nitrogen slopes show distinct site differences. A comparison of the total phosphorus and total nitrogen slope estimates at each site indicates that phosphorus decreases more rapidly with flow than does total nitrogen. Together these results indicate that higher flows are associated with lower total N and total P concentrations, higher N:P ratios, and a higher proportion of reduced nitrogen. Total phosphorus loads over time show considerable yearto-year variability (Figure 2). Although significant reductions in phosphorus inputs have occurred since the imposition of the phosphate detergent ban in 1988 (7), declines in load are obscured by the heavy river flows and associated nutrient loads that occurred during several major hurricanes in 1996 and 1999, near the end of the time-series. Nitrogen loads over time are also highly variable and exhibit no obvious long-term trends in any of the components (Figure 3). VOL. 37, NO. 10, 2003 / ENVIRONMENTAL SCIENCE & TECHNOLOGY
9
2045
FIGURE 2. Estimated total phosphorus loads from 1979 to 2000. An interval of (2 standard error units is shown. Time trends are more apparent in flow-adjusted total phosphorus concentrations, which show substantial drops at all three river stations since the mid 1980s, consistent with the phosphate detergent ban (Figure 4). The Hookerton and Kinston stations have been fairly stable since then but continue to show slight overall declines through the 1990s. Concentrations at Fort Barnwell have shown not only more year-to-year variability but also an overall decline during the 1990s. Patterns in the flow-adjusted nitrogen fractions are not as consistent among the three stations as are those for phosphorus (Figure 5). At Hookerton, total nitrogen showed a steady decline from 1985 to approximately 1994 and has been stable from 1994 to 2000. The decline was principally the result of a decline in nitrate, though TKN showed a slight decrease as well. At Kinston, the flow-adjusted total nitrogen concentration was fairly stable from 1979 to 1995 but has shown a decline from 1995 through 2000. This recent decline appears principally as a nitrate decrease, though, as at Hookerton, TKN shows a slow steady decline over most of the period of record. Patterns at Fort Barnwell are similar to those at Kinston, with recent declines in total nitrogen and 2046
9
ENVIRONMENTAL SCIENCE & TECHNOLOGY / VOL. 37, NO. 10, 2003
FIGURE 3. Estimated nitrogen loads from 1979 to 2000 (total N, nitrate/nitrite, and TKN depicted by circles, triangles, and squares, respectively). An interval of (2 standard error units is shown. nitrate and a slow decline in TKN over most of the 22-year record. Figures 6 and 7 depict the response of the estuary to annual nutrient load change, across the transition from an advective riverine system to a lagoonal embayment. Though phosphorus and nitrogen loads have been highly variable from 1979 to 2000, this year-to-year load variability does not translate into a clear relationship with estuarine phosphorus (Figure 6) or nitrogen (Figure 7) concentrations. While consistent nutrient concentration decreases are apparent moving from the river out into the estuary, there is no statistically discernible relationship between annual phosphorus load and estuarine P concentration and only a slight relationship between annual total nitrogen load and estuarine N concentration (Table 2).
FIGURE 4. Flow-adjusted total phosphorus concentrations over time. The median long-term flow from each site is used as the reference flow. An interval of (1 standard error unit is shown There is, however, a clear positive relationship between the flow-adjusted P concentration at Fort Barnwell and estuarine P concentration (Figure 6) that is statistically discernible at the 0.05 level at all five estuarine sites (Table 2). The analogous total N relationship (Figure 7) is also consistently positive but only statistically discernible at the 0.05 level in the upper two estuarine stations (Table 2). We believe that the N relationship is less pronounced than the P relationship because Fort Barnwell flow-adjusted total N has shown relatively little variability over the period of record, whereas flow-adjusted P has shown considerable change. If the flow-adjusted nitrogen decrease that began in the mid 1990s continues, then we expect that the nitrogen relationship will become more pronounced. We point out that, because flow-adjusted phosphorus and nitrogen are both estimated with uncertainty, and the variance describing this uncertainty is large relative to the total variance of flow-adjusted P and N, there is an “errors-in-variables” or “observation error”
FIGURE 5. Flow-adjusted nitrogen concentrations over time (total N, nitrate/nitrite, and TKN depicted by circles, triangles, and squares, respectively). The median long-term flow from each site is used as the reference flow. An interval of (1 standard error unit is shown. effect that causes these relationships to be understated (16). In other words, the slope estimates (Table 2) for the flowadjusted relationships are biased toward zero. Based on the ratio of the estimation error variance to the total variance, approximate correction factors to adjust for this bias are 1.1 and 1.4 for flow-adjusted phosphorus and nitrogen, respectively. Flow, the other component of load, shows an inverse relationship to estuarine total P concentration in the upper stations, while the relationship between flow and total N is variable and closely mirrors the nitrogen load relationship (Figures 6 and 7). Both relationships with flow are only statistically discernible at the 0.05 level for site 1.
Discussion Recent historical flow at the downstream Fort Barnwell station was highly predictable from the flow at upstream Hookerton VOL. 37, NO. 10, 2003 / ENVIRONMENTAL SCIENCE & TECHNOLOGY
9
2047
FIGURE 6. Regression lines depicting median annual total phosphorus concentrations at Fort Barnwell (top dashed line, circles) and estuarine sites (1 ) plus, 2 ) star, 3 ) square, 4 ) diamond, 5 ) triangle, solid lines, top to bottom in each panel) vs Fort Barnwell annual total phosphorus load (top panel), Fort Barnwell flow-adjusted total phosphorus concentration (middle panel), and total annual flow at Fort Barnwell (bottom panel). and Kinston. This allowed us to significantly extend the historical flow record and, subsequently, the record of nutrient load estimates. Accurate nutrient load estimation required the development of a concentration:flow model to estimate concentrations for days in which measurements were not made. Concentration and streamflow were found to be highly correlated in the Neuse, and the relationship could be represented with the simple linear model:
log(C) ) β0 + β1log (flow) +
(3)
There is high temporal autocorrelation in the daily nutrient data, collected from 1996 to 2000, suggesting that daily data 2048
9
ENVIRONMENTAL SCIENCE & TECHNOLOGY / VOL. 37, NO. 10, 2003
FIGURE 7. Regression lines depicting median annual total nitrogen concentrations at Fort Barnwell (top dashed line, circles) and estuarine sites (1 ) plus, 2 ) star, 3 ) square, 4 ) diamond, 5 ) triangle, solid lines, top to bottom in each panel) vs Fort Barnwell annual total nitrogen load (top panel), Fort Barnwell flow-adjusted total nitrogen concentration (middle panel), and total annual flow at Fort Barnwell (bottom panel). are largely redundant and that lowering the sampling frequency would not seriously compromise the information content of the monitoring program. In fact, to appropriately estimate the uncertainty of our concentration:flow models, we used only every 14th nutrient measurement, without loss of information, at the annual scale. However, all of the measured values were used to calculate the daily load values and were supplemented by model estimates only when measurements were not available. In the relationship given in eq 3, the value of β1 is generally negative (Table 1), indicating that nutrient concentrations decrease with increased flow. Johnson (2) illustrated two types of general relationships between nutrient concentration and flow. In the first, concentration continually decreases as flow increases, via dilution. In the second, concentration initially
TABLE 2. Summary Statistics for Regressions of Median Annual Estuarine Nutrient Concentrations on Annual Riverine Total Load, Flow-Adjusted Concentration, and Stream Flow Phosphorus slope estimate standard error
p-value
site 1 site 2 site 3 site 4 site 5
Total P Load 0.000011 0.000033 0.000019 0.000033 0.000036 0.000033 0.000027 0.000034 0.000025 0.000036
0.737 0.564 0.278 0.432 0.493
site 1 site 2 site 3 site 4 site 5
0.712 0.519 0.277 0.222 0.170
site 1 site 2 site 3 site 4 site 5
Flow-Adjusted P 0.084 0.084 0.084 0.085 0.086
