Development and application of a WASP model on a large Texas reservoir to assess eutrophication control

Abstract

Ernst, M.R. and J. Owens. 2009. Development and application of aWASP model on a large Texas reservoir to assess eutrophication control.

The Tarrant Regional Water District developed and calibrated an 11-yr WASP-Eutro model for Cedar Creek Reservoir, a 13,350 ha (33,000 ac) drinking water supply and recreation reservoir located near Fort Worth, Texas, that is experiencing eutrophication. Cedar Creek was partitioned into 22 segments with up to three vertical segments in the main pool. Loading to the reservoir was divided into four categories: (1) watershed, (2) nine wastewater treatment plants (5.6 mgd), (3) atmospheric loading, and (4) internal NH₄ and OPO₄ sediment flux. Watershed loading was derived from a SWAT model of the 1007 mi² watershed. Wastewater treatment plant loading was based on 1 yr of weekly plant nutrient and flow data. Atmospheric loading was estimated from periodic rain sampling adjacent to the reservoir. Internal flux was estimated from hypolimnetic accumulation during stratification. Calibration was judged by conformity of median observed and predicted water quality data and longitudinal profiles in the reservoir. Sensitivity analysis of these four loadings revealed that Chla and total phosphorus were most influenced by the watershed nonpoint source load. Internal flux was the next most influential load; form, timing and location of loading are more important than actual magnitude of this load. Estimates of the nonpoint source reduction necessary for a significant reduction in Chla were about 30%, while estimates of reductions in OPO₄ flux to significantly reduce Chla were 75–100%. A combination of watershed and internal nutrient control are needed to control eutrophication in Cedar Creek Reservoir.

Key words:

• eutrophication
• flux
• models
• sensitivity analysis
• WASP

Eutrophication concerns in water supply reservoirs are a common concern because of taste, odor, and general increases in treatment costs (Taylor et al. 2006). Nutrient management is typically the course for water quality improvement (Smith et al. 2006). In Cedar Creek Reservoir near Fort Worth, Texas, eutrophication manifests as planktonic algae, often dominated by nitrogen-fixing blue-green algae. Given these conditions, phosphorus management is likely the best approach (Reckhow and Simpson 1980), but only if it can be sufficiently controlled to affect the algae population. Thomann and Mueller (1987) noted that reducing an abundant nutrient may not have the desired effect until concentration reaches a sufficiently low level to limit nutrient growth.

The U.S. Environmental Protection Agency's Water Analysis Simulation Program (WASP) Version 6.2 is a powerful water quality model used to predict and interpret water body responses to various nonpoint source loads and point-source pollution (Wool et al. 2003). This model was selected by the Tarrant Regional Water District (TRWD) to quantify the phosphorus balance within the reservoir and develop a mechanistic relationship between phosphorus and algae growth. Eight state variables were simulated, including dissolved oxygen, biochemical oxygen demand, inorganic phosphorus (P), organic P, ammonia, nitrate-nitrite, organic nitrogen (N), and algae as chlorophyll a (Chla). Similar work was done by Debele et al. (2006) for a 5-yr period on Cedar Creek Reservoir with CE-QUAL-W2, but this effort concentrated on the mechanics of integrating a watershed model with CE-QUAL-W2 more than the analysis of management implications for improving reservoir water quality. A number of other entities have used WASP for analysis of systems (Lung and Larson 1993, Jin et al. 1998), but the literature is void of an application on a large southwestern United States reservoir. This work provides a framework for modeling large water supply reservoirs common to the southwestern United States and assessing the impact of nutrient loads from various sources.

Methods

Physical depiction of Cedar Creek Reservoir in WASP

Cedar Creek is a 13,350 ha (33,000 ac), recreational and water supply reservoir located in north central Texas about 100 km southeast of Fort Worth, Texas. The reservoir, built by TRWD in the early 1960s, is a primary source of water for Fort Worth and the surrounding cities. It has an average hydraulic detention time (volume/outflow) of approximately 2.4 yr. The trend in Cedar Creek Chla over the past 18 yr has a significant positive slope with an annual rate of increase of 3.85% (Fig. 1). Median concentration for the reservoir is 19.6 μ g/L, and the third quarter median concentration is 27.4 μ g/L. Cedar Creek has been (303d) listed as impaired by the State of Texas for high pH violations that resulted from excessive algae photosynthesis and low alkalinities (lakewide median 54 mg/L CaCO₃). Cedar Creek currently ranks number 11 of 104 Texas reservoirs for highest Chla. Additional descriptions of Cedar Creek can by found in Debele et al. (2006) and Grover et al. (1999).

Figure 1 Annual median + 25th percentile Chla and the trend for 18 years of lakewide Cedar Creek Data.

The reservoir was divided into 22 segments to represent relatively homogeneous sections (Fig. 2). Much of the discussion in this paper will be focused at Segment 4, the location of the District's raw water pump station. Thermoclines, longitudinal gradients in water clarity, phosphorus and Chla concentrations, and the availability of observed data formed the basis for the segmentation scheme.

Figure 2 Cedar Creek Reservoir segmentation for the WASP model.

Segments are defined either as surface or subsurface segments. Surface segments include the main thalweg segments (1–6) and the eight cove segments (15–22). Depth of the main thalweg segments approximates 2 m, which is the typical photic zone for the reservoir. This segmentation scheme was set up to allow positive growth of algae in these segments. Depth of the cove segments was calculated as the volume divided by the surface area. Subsurface segments define the remaining volume of the Cedar Creek Reservoir. Subsurface segments 7, 8, 9, 10, and 11 characterize the Cedar Creek WASP model to the depth of the typical thermocline at approximately 7 m. These are called mid-segments. Subsurface segments 12, 13, and 14 define the three hypolimnetic segments in the Cedar Creek WASP model from the main pool area of the reservoir.

Dispersion

Due to large segment interfacial areas in the Cedar Creek WASP model, horizontal and vertical dispersion serves as an important transporter of mass. Horizontal dispersion was estimated from the 4/3 Power Law (Bowie et al. 1985). Horizontal dispersion ranged from 1 m²/sec to 10 m²/sec. Vertical dispersion between the surface segments and mid-segments (2–7, 3–8, 4–9, 5–10, and 6–11) were arbitrarily set at a high rate (0.001 m²/sec) to ensure uninhibited vertical mixing between these segments.

Hypolimnetic vertical dispersion has two components: the rate of dispersion and a time variable function. The duration of stratification was based on thermal plots from Cedar Creek field data at Segment 4 comparing 1 m below the surface to 1 m above the reservoir bottom to determine the period of distinct temperature differentials each summer and minimal mixing for hypolimnetic segments (9–12, 10–13, and 11–14). Based on these temperature plots, TRWD determined the duration of the temperature differentials for each year at Cedar Creek Reservoir and applied typical lake vertical dispersion rates listed by Chapra (1997). These rates varied from 0.0005 m²/sec for well-mixed time periods to 0.00001 m²/sec for summer stratification time periods. Dissolved oxygen modeling allowed for calibration of these dispersion rates. For example, 1996 data showed weak stratification and limited hypolimnetic anoxia; hence, vertical dispersion was increased in this summer period to allow more mixing and better simulate the observed data.

Hydrodynamics

To accurately model the transport and transformation of the nutrients in a water body, the hydrodynamics must be accurately represented within the model. For this effort, we used a hydrodynamic flow model developed by Alan Plummer and Associates, Inc. (Fort Worth, TX). This program uses the external flows to the system (precipitation, evaporation, pumpage, discharge, and tributary inflows) as recorded by TRWD and the corresponding geometry of each segment to solve for advective flows between adjacent segments. The program specifies a matrix solution employing the criteria of minimum kinetic energy, and the solution is input into the appropriate flow field for each segment in WASP. Important controls in the program include allowing only Seg 7, 8, 12, 13 and 14 to change in volume; all other segments maintain a constant volume and depth. Flows were input on a monthly basis.

The majority of tributary inflow comes from Kings and Cedar creeks into Segments 22 and 21, respectively. Most of the flow leaving the reservoir is through Segment 9, which has a water intake for the water supply pump station and the spillway for release of flood waters.

The 11-yr calibration period encompassed the near normal early 1990s and droughts of 1996 and 2000 (Fig. 3). During the start of the WASP simulation time-period, Cedar Creek Reservoir was at 92% of conservation pool volume. Because the initial volume is crucial in the expression of nutrient mass as a concentration, the initial volumes for the five segments (7, 8, 12, 13, and 14) that are allowed to vary in volume were adjusted downward to account for the difference in starting pool capacity.

Figure 3 Cedar Creek Reservoir hydrology surrounding the 11-yr WASP modeling time period 1991-2001. Conservation pool volume of the reservoir is 637,109 ac-ft at an elevation of 322 ft m.s.l.

Settling rates

The physical settling of particulate matter in any reservoir is an important transport phenomenon of particulate nutrients and often leads to a distinct longitudinal gradient in concentration. In WASP, settling is controlled by the fraction of the particulate nutrient and the settling velocity assigned. Only three of the eight state variables were assigned settling velocities: organic N, organic P and Chla. Laboratory analysis of filtered and unfiltered organic N and P samples allowed for estimation of the fraction dissolved for each species. Organic N averaged 64% dissolved and organic P 31% dissolved, but these were input on a segment-specific basis. We found it necessary to manipulate this term in boundary Segments 20–22 and set both organic N and P to 100% dissolved so that mass would be routed to the reservoir. In reality these segments often become shallow during dry periods and are short-circuited, but in the model they are a set volume with tremendous settling potential. We also had to increase the percent dissolved in Segment 1, from limited measured data suggesting 0.58 for organic N and 0.15 for organic P to 0.65 for both nutrient species, to allow this segment to conform more closely to the observed data.

Algae were assigned a settling rate of 5.0 × 10^{− 7} m/sec (0.04 m/day). The fraction of organic N and organic P determined to be in the particulate phase were given settling velocities of 8.0 × 10^{− 7}m/sec (0.07 m/day) and 1.6 × 10^{− 6} m/sec (0.14 m/day), respectively. Organic P was given a higher settling velocity than organic N to better fit observed data. Organic P has consistently shown longitudinal decrease in concentration from headwaters to dam in our reservoirs, a trend not as apparent with organic N. We speculate that organic P may bind with inorganic sediment more readily than N, but nitrogen fixation possibly plays a role in maintaining a more consistent concentration of nitrogen throughout the reservoir. The velocities we used conform to those in Chapra (1997) and are validated by the observed data longitudinal profiles. WASP does allow for time-variable settling rates, but we did not use that option in this model.

Environmental time functions

WASP requires time functions for important environmental functions such as temperature, incident light, light extinction, photoperiod, and wind. For this type of water-quality modeling, water temperature, light extinction, and incident light are critical. Water temperature and incident light are optimal in most Texas reservoirs and not a limiting factor for algae growth. Light can be limiting in the water column due to nonalgal turbidity. Four light curves were used to represent the longitudinal gradient from the turbid north end of Cedar Creek to the relatively clearer waters in the southern end near the dam. Curve 2 represents the north end of Cedar Creek (Segments 2, 7, and 20), while Curve 4 represents the relatively clearer southern end of Cedar Creek near mid-lake (Segments 4, 9, 12, 15, and 16; Fig. 4). For this effort, the nonalgal, light extinction coefficient (Ke') was calculated using (Ernst 1995):

where z is the Secchi depth in meters and Chla is in ug/L.

Figure 4 Nonalgal light extinction (1/m) curves for the upper end (dashed) and dam area of Cedar Creek Reservoir.

Nutrient loading

The Cedar Creek WASP model includes four types of nutrient loading systems:

1.	Lakeside point source loading
2.	Atmospheric deposition
3.	Benthic flux loading
4.	Watershed loading (point and nonpoint sources)

Lakeside point source discharges

This nutrient loading system includes two wastewater treatment plants (WWTPs) that discharge treated effluent directly into Cedar Creek. Plant 1 had weekly nutrient and discharge data available for May 2001 through April 2002 and was used to calculate the monthly load of nutrients to Segment 19 of Cedar Creek in kg/d. Plant 2 discharge data was available for October 2001 through September 2002 and was applied to Segment 4 of Cedar Creek. No flow data are associated with the WWTP load data input to the model because the contribution of flow was only 2310 m³/day (0.61 mgd). Because only 1 yr of nutrient discharge data was available, this year of data was input on a monthly time step and recycled for the 11-yr simulation time-period for the calibration model. The WWTPs modeled in SWAT were dealt with the same way.

Atmospheric loads

Nutrient loading from the atmosphere was calculated using rainfall amounts and nutrient data (NH₃, NO_X, TKN, total and ortho P) measured by TRWD on rainwater samples. The loads were then converted to a constant daily rate and applied to the model's surface segments (Table 1). The rates were applied as constants throughout the 11-yr simulation.

Table 1 Atmospheric loading rates to Cedar Creek Reservoir estimated from 26 rainfall samples over the period 1990–2004 at the reservoir and the average annual rainfall of 39 in/yr for this area.

Variable	Rate (mg/m^2-day)
Atmospheric Deposition of Nitrate	0.826
Atmospheric Deposition of Ammonia	0.75
Atmospheric Deposition of Orthophosphate	0.092
Atmospheric Deposition of BOD	0.0
Atmospheric Deposition of Organic Nitrogen	1.396
Atmospheric Deposition of Organic Phosphorus	0.06

Benthic flux

Benthic flux in the form of ammonia (NH₄) and orthophosphate (OPO₄) was added to the three hypolimnetic segments (12, 13, and 14). Analysis of intensive survey data (bimonthly to monthly sampling frequency) from two summers (1989 and 2000) allowed estimation of release rates from hypolimnetic increases in concentration. These rates of release and their duration seemed to be the best estimate of reality and were used in the model. Ammonia flux was estimated at 65 mg/m²/day, and phosphorus was estimated at 7 mg/m²/day for a duration of 43 days each summer beginning 20 June and ending 1 August. The rate and duration were constant for each year of the simulation. Data are insufficient to determine if some years have higher rates or longer periods of benthic flux.

Watershed loads

Nutrient loading from the watershed included point-source (PS) discharges from seven WWTPs located in the Cedar Creek watershed, and overland flow from approximately 1000 mi². These combined nutrient loads from the watershed were estimated by a separate consultant using the SWAT model and supplied to TRWD as an external nonpoint source (NPS) file for linkage to WASP (Debele et. al 2006). The loads were supplied at the point where tributaries met the reservoir. Loading from areas directly around the reservoir was assigned to the adjacent reservoir segment. Loads were given for all eight state variables on a daily basis.

Kinetics

Effort was made to keep all the kinetics within the ranges specified in the WASP model or literature (Bowie et al. 1985; Table 2). Kinetics that effect the nutrient cycles and algae growth deserved the most attention. We were fortunate to have an empirical basis for the nitrogen half saturation constant of 0.0485 mg/L from a number of Cedar Creek bioassays performed by Sterner and Grover (1998) and Grover et al. (1999). The phosphorus half saturation constant of 0.007 mg/L P was tuned during calibration. The half-saturation constants have an N:P ratio of 6.92. Similarly, the N:C ratio of 0.15 and P:C ratio of 0.022 was fit to the Cedar Creek model during calibration and suggests an N:P stoichiometry of 6.82. Both ratios show a bias toward P limitation, by making P relatively more necessary than under a ratio of 10:1 or higher. This was determined to best represent Cedar Creek and its large proportion of nitrogen-fixing blue-green algae. Kinetics that favor 10:1 or higher N:P ratios are most representative of green algae communities and underestimate the late summer Chla in Cedar Creek. Ideally, a WASP model that allowed simulation of two algal groups would circumvent this problem, but this feature was not available.

Table 2 Kinetics used in the 11-yr WASP model of Cedar Creek Reservoir. The “Avg Range” is based on guidance provided in WASP Version 6.2 and Bowie et al. (1985).

WASP Kinetic Constant Type	Avg. Range	Cedar Creek	Units
Nitrification Rate @ 20°C	0.001–0.2	0.15	day ^− 1
Nitrification Temp Coeff.	1.02–1.08	1.04	NA
Half Saturation: Nitrification Oxygen Limit	0.5–2.0	1.0	mg O2/L
Denitrification Rate @ 20°C	0–0.09	0.06	day ^− 1
Denitrification Temp Coeff.	1.02–1.09	1.06	NA
Half Saturation: Denitrification Oxygen Limit	0–2.0	2	mg O2/L
Phytoplankton Growth Rate @ 20°C	1.0–3.0	2.0	day ^− 1
Phytoplankton Growth Temp Coeff.	0–1.07	1.06	NA
Phytoplankton Light Formulation Switch (1 = DiToro)	NA	1 = DiToro	NA
Phytoplankton Max Quantum Yield Constant	NA	NA	NA
Phytoplankton Self Shading Extinction	NA	NA	NA
Phytoplankton Carbon::Chlorophyll Ratio	0–200	50	mg carbon/mg chla
Phytoplankton Optimal Light Saturation	0–350	200	Ly/day
Phytoplankton Half Saturation Constant: Nitrogen	0.01–0.06	0.0485	mg-N/L
Phytoplankton Half Saturation Constant: Phosphorus	0.0005–0.05	0.007	PO4-P/L
Phytoplankton Endogenous Respiration Rate @ 20°C	0–0.5	0.05	day ^− 1
Phytoplankton Respiration Temp Coeff.	1.0–1.08	1.045	NA
Phytoplankton Death Rate Non-Zooplankton Predation	0–0.25	0.05	day ^− 1
Phytoplankton Zooplankton Grazing Rate	0–5	NA	L/cell-day
Nutrient Limitation Option (0 = Min; 1 = Multiplicative)	0, 1	0	NA
Phytoplankton Decay Rate in Sediments @ 20°C	0–0.02	0.02	day ^− 1
Phytoplankton Decay Rate Temp Coeff.	1.0–1.08	1.08	NA
Phytoplankton Phosphorus::Carbon Ratio	0–0.24	0.022	mg P/mg C
Phytoplankton Nitrogen::Carbon Ratio	0–0.43	0.15	mg N/mg C
Phytoplankton Half Saturation for N and P	0–1.0	0	NA
BOD Decay Rate @ 20°C	0.05–0.4	0.1	day ^− 1
BOD Decay Rate Temp Correction	1.0–1.07	1.04	NA
BOD Decay Rate in Sediments @ 20°C	0.0004–1.0	1.0	day ^− 1
BOD Decay Rate in Sediments Temp Coeff.	1.0–1.08	1.08	NA
BOD Half Saturation Oxygen Limit	0.5–1.0	1.0	NA
Waterbody Type for Wind Driven Aeration	1.0–3.0	NA	NA
Oxygen::Carbon Stoichiometeric Ratio	0–2.67	2.67	mg O2/mg C
Reaeration Rate Constant @ 20°C	0.5–3.0	1	day ^− 1
Reaeration Rate Option (sums Wind and Hydraulic Ka)	0–1.0	NA	NA
Dissolved Organic N Mineralization Rate @ 20°C	0.02–0.075	0.02	day ^− 1
Dissolved Organic N Mineralization Temp Coeff.	1.0–1.08	1.045	NA
Organic N Decay in Sediments @ 20°C	0.0004–0.01	0	day ^− 1
Organic N Decay in Sediments Temp Coeff.	1.0–1.08	1.045	NA
Fraction of Phytoplankton Death Recycled to ON	0–1.0	1.0	NA
Dissolved Organic P Mineralization Rate @ 20°C	0–0.22	0.045	day ^− 1
Dissolved Organic P Mineralization Temp Coeff.	1.0–1.08	1.045	NA
Organic P Decay in Sediments @ 20°C	0.0004–0.01	0	day ^− 1
Organic P Decay in Sediments Temp Coeff.	1.0–1.08	1.08	NA
Fraction of Phytoplankton Death Recycled to OP	0–1.0	1.0	NA

Results

Cedar Creek WASP model calibration

The Cedar Creek WASP model was calibrated for an 11-yr period starting 1 January 1991 and ending 31 December 2001. Calibration involved comparison of median and percentile observed data for main channel segments to median predicted data for the entire 11-yr period. Attempts to calibrate on a more frequent basis, such as each observed sampling event or even quarterly, were too variable. Our approach to calibration of the WASP model generally begins with phosphorus modeling because it has fewer pathways and kinetics than nitrogen or algae. The two most critical components to phosphorus calibration are the settling rate and mineralization of the organic P. Longitudinal profiles in the reservoir and measurements of both organic and inorganic dissolved P facilitates setting these parameters (Fig. 5). We had a significant correlation (p < 0.05) between median concentrations both annually and seasonally for all three species of phosphorus and for P limitation (Table 3). Phosphorus limitation is the Monod equation with the half-saturation constant of 0.007 mg/L (Table 3). Insight into P limitation is often inhibited by high P detection limits. In this study, dissolved inorganic P had a minimum detection limit of 0.005 mg/L, yet some samples were reported less than this amount. Censoring those data to 1% of the detection limit allows for spreadsheet calculation of the P limitation and reveals just how low the P limit could actually be. Not censoring the data will bias the P limit to high values. Nitrogen modeling was the second step and performed reasonably for total N (TN; Fig. 5) but did not produce significant correlations between observed and predicted data for ammonia, nitrate-nitrite or organic N annually (Table 3). Excessive effort to calibrate the various species of N was not exercised because it is not the nutrient of concern in this study. The TN:TP ratio for the main lake segments (Segments 4, 5, and 6) were very close, providing for a good basis to model algae.

Figure 5 WASP model calibration statistics. Total phosphorus (TP), total nitrogen (TN), Chla, TN:TP, P-limitation and censored data P-limitation. Observed and predicted median and 25th percentiles for segments of Cedar Creek Reservoir ranging from the headwaters (Segment 1) to the dam (Segment 6). One mid-depth segment (9) was added to the graph because of a wealth of data at this site.

Table 3 Statistical comparison of median observed and predicted data from an 11-yr WASP model of Cedar Creek Reservoir. Seasonal analysis based on April–October. Bold r-square values significant at p = 0.05. RPD (Relative Percent Difference) calculated as average ABS(o-p)/avg(o + p)*100 (ABS = absolute value). Lab QC (quality control) is the RPD of 32 duplicate field samples for these parameters from 2000 to 2005 on Cedar Creek Reservoir.

			R-square Values		Relative Percent Difference
Parameter	Observed Data Annual	Seasonal	Annual	Seasonal	Annual	Seasonal	Lab QC
NH3 (mg/L)	0.053	0.051	0.1622	0.3372	15.2%	42.4%
NOX (mg/L)	0.108	0.033	0.5637	0.0121	63.3%	67.2%
Org N (mg/L)	0.744	0.697	0.2976	0.9340	15.5%	31.4%
TN (mg/L)	0.925	0.869	0.9050	0.9361	12.9%	20.9%	17.2%
OPO4 (mg/L)	0.020	0.019	0.9209	0.7601	13.1%	37.4%
Org P (mg/L)	0.073	0.074	0.8455	0.9430	22.7%	17.6%
TP (mg/L)	0.103	0.100	0.9345	0.9200	10.0%	15.2%	16.8%
TN:TP	10.2	10.1	0.6721	0.3183	18.2%	21.4%	28.3%
Chl'a' (ug/L)	19.950	19.460	0.2332	0.0130	17.7%	26.4%	21.0%
N-limit	0.78	0.70	0.3315	0.0841	9.0%	11.6%
P-limit	0.72	0.71	0.9704	0.8739	4.1%	14.1%

Algae, represented by Chla in WASP, is both too complicated and too simple: complicated by the numerous kinetics available, and simplified by only one algae group that does not fix nitrogen. After numerous attempts we felt we had a model that reasonably represented the reservoir and responded appropriately; however, the correlations of observed and predicted data were not significant (P > 0.05). The relative percent difference in the model results were very close to those observed in our Quality Control (QC) program where field samples are split and submitted to the lab. An approximate 20% variation in Chla splits is common, and model results within 20% of predicted data were considered acceptable calibration (Hutchinson et al. 1991).

Cedar Creek nutrient balance

An 11-yr TP mass balance of Cedar Creek showed that the majority (87%) of TP came from nonpoint sources, followed by WWTPs (7%), sediment flux (3%) and atmospheric deposition (3%). Annually, nonpoint source loads were the only variable loads. The WWTPs, benthic flux, and atmospheric loading were kept constant in the model because information was insufficient to vary these loadings. The phosphorus budget was simulated in WASP, broken down into organic and inorganic forms (Fig. 6). The values represent average loads for all inputs and outputs over the 11-yr simulation. Nonpoint sources contributed the highest load to the reservoir, with the majority being organic P. The majority of the atmospheric and WWTP loads were inorganic P, as well as the entire load released during the flux period. However, these loads combined were a much smaller contribution than the phosphorus brought in from the nonpoint source. The reservoir acts as a sink for phosphorus, retaining an average of 70% of the inorganic P and 80% of the organic P.

Figure 6 Average annual phosphorus budget for Cedar Creek Reservoir. All values are in kg/yr.

Sensitivity analysis

The Chla response in the calibrated WASP model of four nutrient-loading scenarios was evaluated by systematically shutting each off. When looking at Segment 4, Chla was most sensitive to the watershed load (P < 0.05) and not significantly affected by the other sources of loading (Fig. 7). Watershed loading also played the most significant role in Chla production at the dam site (Segment 6), but here benthic P flux also reduced Chla significantly (P < 0.05). While watershed loading is the most significant forcing function, its impact is relatively less in Segment 6 than in Segment 4, probably a result of the Cedar Creek hydrology routing much of the flow out through the pipeline and spillway in the vicinity of Segment 4 rather than the more conventional location of the dam.

Figure 7 Chla sensitivity in Segment 4 to systematic removal of sources of nutrients. The calibated model has all four sources of nutrients operating, while the No Watershed Loading simulation has that source shut off. Likewise, the No WWTP has the nine wastewater treatment plants turned off; No flux has the benthic P -flux turned off, and No ATM has the atmospheric loading turned off. Simulations with the same letter (A) are not significantly different (Kruskal Wallis test, P = 0.05).

Load reduction scenarios

Watershed reductions

Five systematic load reductions were simulated by scaling the NPS file to create reductions ranging from 15 to 65%. This would scale both N and P loads. A significant reduction (P < 0.05) in Chla concentration was not realized until a 35% reduction in the watershed loading was achieved (Fig. 8, top). Benthic flux reduction scenarios were done only with phosphorus since this is the most manageable nutrient to sequester (alum treatment). Reductions in P flux were done in a systematic fashion similar to that for the watershed loads. Phosphorus flux was reduced from 15 to 75%. No significant reduction in Chla could be achieved at Segment 4 (Fig. 7) even with up to a 100% reduction, but Segment 6 did show a significant reduction in Chla when flux was reduced by 75% (Fig. 8).

Figure 8 Chla sensitivity to systematic reductions in the water loading (Segment 4, top graph) and benthic P-flux (Segment 6, bottom graph). Plots are median and 25th percentiles for 11-year simulation. Simulations with the same letter are not significantly different (Kruskal Wallis test, P = 0.05).

Combination scenarios

Combination scenarios where an in-lake reduction in phosphorus was coupled with watershed management practices were explored. Five combinations were assembled:

1.	Reduce flux by 25% in the hypolimnion and NPS loading by 15%.
2.	Reduce flux by 50% in the hypolimnion and NPS loading by 25%.
3.	Reduce flux by 50% in the hypolimnion and ban the use of phosphorus fertilizer in a 2000-ft perimeter band of the reservoir.
4.	Reduce flux by 100% in the hypolimnion and NPS by 15%.
5.	Reduce flux by 75% and NPS by 30% in the two most prominent subwatersheds of Cedar Creek Reservoir that constitute 59% of the total watershed.

The results show that combination 2, 4, and 5 reduced Chla significantly (P < 0.05) in Segment 4 (Fig. 9). Segment 6 responded similarly.

Figure 9 Chla sensitivity in Segment 4 to potential combination of management options to reduce phosphorus loading to Cedar Creek Reservoir. See text for detailed descriptions of each scenario. Simulations with the same letter are not significantly different (Kruskal Wallis test, P = 0.05).

Discussion

The WASP model is a tool for understanding the nutrient budgets and algae growth of large reservoirs. This modeling effort is actually a culmination of 12 years of working with this model. Most of the previous models were built around one or two year years of data. Moosmann et al. (2006) suggested averaging over a 5–10-yr period for phosphorus budgets to improve accuracy. This 11-yr effort presented a great overall picture by encompassing some of the hydrologic variability and a large quantity of observed data. Nevertheless, the model is most valuable in assessing the relative impact of a change in loading (e.g., how much the removal of WWTPs changes the calibrated model). This approach differs from a more “absolute” approach where one attempts to hit a target Chla concentration. The fact that we have simplified Cedar Creek Reservoir to one algae group with no zooplankton or fish communities should be sufficient warning to any user that targeting a specific Chla concentration is very risky.

It does appear that WASP mimicked the phosphorus cycle and algal growth dynamics more closely than CE-QUAL-W2 (Debele et. al. 2006). When comparing statistics, Debele did not report statistically significant relationships between observed and predicted data for TP (r = 0.283) or TN (r = 0.279) or (r = 0.036). This paper reports excellent agreement for TP (r = 0.967) and TN (r = 0.951) and reasonable agreement for Chla (r = 0.483; Table 3). We strive for a visual fit of the observed to predicted data because few statistics can capture seasonality, long term trends, outliers, autocorrelation, and off-set phases as well as the modeler. After we have achieved a “visual” calibration, statistics such as correlation help substantiate our effort; however, the large variability of Chla makes correlation a challenge. Relative Percent Difference (RPD) estimates provide insight into the fit of the Chla. We strive for Chla predictions to be within 20% of observed median data.

Debele et al. (2006) attributed their poor results to propagation and compounding of errors. In CE-QUAL-W2 much effort was put into calibrating the hydrodynamic model to approximate Cedar Creek's surface elevation, while in WASP this was a given. Because we did not model the hydrodynamics, we described them via a monthly time step using the most precise data available for elevation, volume, inflow, outflow, evaporation. Similarly, while much effort was spent calibrating CE-QUAL-W2 temperature profiles, the WASP model has four hardwired temperature curves based on observed data and assigned to respective segments. In essence, WASP has fewer “degrees of freedom” on the front end that may make subsequent nutrient and algae modeling less influenced by error.

We chose the vertical segmentation to approximate the top segment 2 m in depth and the middle segments from 2 m to the depth of the thermocline at 7 m or the bottom, whichever came first. This approach was developed because WASP averages the light over the segment depth and uses that intensity for growth of algae. If the segments are deep, then the average light condition is too limiting for realistic algae growth. We opted to segment the model to get sufficient growth in the top segment and provide sufficient dispersion into the midsegments to better represent the observed data. However, some separate model runs suggest that this may not be necessary. It seems that combining the top and middle segments produces the same algal growth as separating the two segments and allowing for mixing; two vertical segments would greatly simplify the modeling effort.

The internal flow method used requires all segments except for a select group of deep, large volume segments to remain constant in volume throughout the model simulation. WASP will develop numerical instability if a segment volume gets too small. Typically, cove segments reduce in volume sometime each year, and the tributary inflow follows the old creek channel out into the main body of the reservoir. In essence this is short circuiting through the segment. WASP sees the segment at the original volume specified and operates as a completely stirred reactor (CSTR) whereby settling and dilution can be exaggerated. We circumvent this problem by increasing the percent dissolved for organic N and P to facilitate their routing to the main lake.

The use of observed 11-yr medians from main lake sites as a calibration basis simplified this modeling effort. Initial efforts were focused on time-series plots of all the observed data at monthly to quarterly intervals versus predicted data at monthly intervals. Cove segments showed large variability in predicted data because of their relatively small volume and large loadings of mass. Main lake segments have much less variability. While it is useful to check the cove segments for order of magnitude concurrence, the main lake sites on a gross time scale offer a more consistent framework to calibrate a model.

Because Cedar Creek is dominated by blue-green algae that can fix nitrogen, the model must be parameterized to be more influenced by phosphorus than nitrogen. This was done by keeping the ratio of the N and P half saturation constants around 7:1, rather than the 10:1 or greater ratios more common for green algae communities (Rhee 1978). While the N:P ratio is an accepted key controller of algae populations (Sterner et al. 1992), it may not elucidate the actual limiting nutrient because either nitrogen or phosphorus could be in sufficient quantities for algal growth. The half saturation model appears to be a better assessment if the nutrient concentration is indeed limiting algae growth, assuming the half saturation constant is representative of the algal community.

The significance of benthic flux in algal growth was surprising, but it becomes apparent when you consider the timing (growing season) and location (main pool with ample light). Loading from WWTPs was a greater percent of the P budget than benthic flux, yet the impact from WWTPs measured by a change in Chla was not as noticeable. This seems to be a result of the WWTP load coming into the coves, and much of the load is settled as organic particulates.

Benthic flux was originally estimated by sediment samples for TP and use of Nurnberg's (1988) regression equations. We estimated P flux for Cedar Creek at 7.33 mg/m²/day. Rates derived by analysis of hypolimnetic concentration from two summers were very similar at 7.18 mg/m²/day. We feel confident with the flux rates and with the duration of hypolimnetic release but have concerns regarding the extent of flux into other areas of the reservoir. We know that areas of rich organic sediments become anoxic for short periods of time, such as a few days, but their contribution to the P pool is unquantified. Because Cedar Creek's phosphorus concentration is so constant, even through variable hydrology and NPS loading, the sediments would logically play a major part in maintaining an “equilibrium concentration” (Reddy et al. 2007).

In addition to the uncertainty in the quantity of P flux, its relationship to loading complicates management strategies. Benthic P flux is more of a mechanism for phosphorus to enter the water column than a source. The source is the NPS, WWTPs and the atmosphere. The question is how long P flux will mask any attenuation of true phosphorus sources. We are a long way from linking P flux to the sources.

There are some unquantified sources of phosphorus in the watershed, namely loadings from onsite sewage facilities and resuspension of reservoir sediment. Cedar Creek Reservoir has approximately 75% of the shoreline developed with a combination of primary and weekend homes. A portion of the developments are serviced by an organized sewer system, but many homes are required to use individual on-site sewage facilities (OSSF). To estimate the potential phosphorus contribution to the reservoir from shoreline OSSF systems, a method from Reckhow and Chapra (1983) was applied to Cedar Creek using the number of systems, phosphorus concentration for residential sewage, and soil type. The number of systems was estimated to be 1587 by counting the number of docks in non-sewered areas from aerial images of the reservoir. The phosphorus concentration was set at 10.45 mg/L, the average from a dataset of 20 grab samples collected by TRWD from aerobic and anaerobic septic tanks. The soil types were determined from a digitized soils layer for the area that was a combination of Computer Based Mapping System (CBMS) and Soil Survey Geographic (SSURGO) data. Assuming an average wastewater usage per person of 90 gal/d and an average of 3 people/residence, the phosphorus load was calculated to be 3623 kg/yr to the reservoir. In the nutrient budget for the reservoir, this is < 1% of the TP contribution; therefore, it was determined that septic systems are not a significant contributor of nutrients to the reservoir and were not included as input to the WASP model.

Resuspension was investigated by some simple tests using filtered reservoir water, sediment samples, and a shaker table. We found P levels increased in the water after sediment was suspended and allowed to settle. This did not happen on two other reservoirs we operate, possibly because the lack of calcium in the Cedar Creek sediments makes phosphorus soluble upon resuspension. We do not have sufficient data such as sediment mapping, daily wind speed, or daily turbidity to incorporate an estimate of this loading source into the model.

Management implications

The Tarrant Regional Water District (TRWD) is faced with developing a watershed management plan that improves or minimizes the chance for any further degradation of Cedar Creek's water quality. While there are a plethora of studies cautioning about the slow response of lakes to external load reductions because of internal loading (Fisher et al. 2005, James et al. 1994), there is strong consensus that a reduction in the external load will improve the general water quality of the reservoir (James et al. 2005, Jeppesen et al. 2005). Cedar Creek does not have an obvious “smoking gun” such as a large WWTP whose discharge could be diverted to quickly reduce loading by 50%. The major source of phosphorus is diffuse, nonpoint source contributions from a 1007 mi² watershed. Our plan will focus on decreasing the external phosphorus load from this watershed. Agriculture best management practices (BMPs) are being assessed with SWAT to see what can be done to reduce loading from 15–30%. Wastewater treatment plants have been surveyed to allow an estimate of the cost to retrofit these plants for TP limits of 1 and 0.5 mg/L as well as estimates of WWTP's projected flow in 2050. Educational programs on lawn fertilizing, rainfall capture, and construction site management will be addressed. Discussions have begun with aeration and alum companies on techniques to sequester the benthic flux of phosphorus. Effort will be made to mold various watershed management practices together into a toolbox capable of reducing the P load by 35%.

Acknowledgments

The authors wish to thank all the members of the North Central Texas Water Quality Committee for making this work possible.