The State of Washington Department of Ecology performed a total maximum daily load (TMDL) assessment for the Spokane River and Lake Spokane (Fig. 1) in eastern Washington. Lake Spokane, with a typical depth at the dam of ∼45 m, is a 38 km long reservoir located along the Spokane River (Fig. 2). The TMDL was established to limit pollutants affecting dissolved oxygen (DO) in Lake Spokane, including ammonia, total phosphorus (TP), and carbonaceous biochemical oxygen demand (CBOD; Cusimano 2004, Moore and Ross 2010). A water quality and hydrodynamic model was developed to evaluate strategies for meeting water quality standards. The model developed for both the Spokane River and Lake Spokane, CE-QUAL-W2 (Cole and Wells 2008), was applied to 2001, a critical low-flow year. A series of reports were produced documenting the modeling procedure and the analysis of management strategies to meet water quality objectives (Annear et al. 2001, 2005, Berger et al. 2002, 2003, 2007, 2009, Slominski et al. 2002, Wells, et al. 2003).

Modeling the response of dissolved oxygen to phosphorus loading in Lake Spokane

All authors

Scott A. Wells & Chris J. Berger

https://doi.org/10.1080/10402381.2016.1211910

Published online:

05 August 2016

Figure 1. Spokane River and Lake Spokane in eastern Washington, USA.

Display full size

Modeling the response of dissolved oxygen to phosphorus loading in Lake Spokane

All authors

Scott A. Wells & Chris J. Berger

https://doi.org/10.1080/10402381.2016.1211910

Published online:

05 August 2016

Figure 2. Water quality monitoring sites at Lake Spokane.

Display full size

Brett et al. (2016) evaluated the TMDL model and concluded that it does not represent Lake Spokane's hypolimnetic oxygen response to changes in phosphorus concentrations entering the lake. In addition, they claimed that the TMDL model was calibrated “for the wrong reasons,” had a “critical flaw,” did not respond “as expected to changes in external nutrient inputs,” was “not able to represent the TP-Chl-DO conditions in Lake Spokane sufficiently well to accurately predict the effect of the TMDL,” was unable to “reproduce the fundamental causal relationships underlying the eutrophication problem,” “should be recalibrated to a much more rigorous dataset,” has a “cryptic phosphorus source,” and “was poorly calibrated to field TP and Chl data, and the quality of the field phosphorus data used in the model was problematic.”

To assess the validity of these statements, we will explore the following main points:

Calibration dataset

Brett et al. (2016) suggested that 2001 was not a good calibration year because TP data were sparse and of questionable value. The model was originally calibrated to field data in 1991 and 2000 (Annear et al. 2001, Berger et al. 2002, 2003, Slominski et al. 2002, Wells et al. 2003) and only later applied to 2001, a critical low-flow year that would be used for the TMDL analysis. TP and DO data available from 1991 included 60 profiles from 5 sites, and 27 profiles from 7 sites were available from 2000. The 2001 profile data were only measured at 2 sites, including only 4 vertical profiles of TP and 10 DO profiles. Applying the model to 2001 was termed a “calibration check” rather than a full-scale change in the model structure and parameterization from the earlier calibration years. Instead of looking at the entire model calibration period (i.e., 1991, 2000, and 2001), Brett et al. (2016) focused only on 2001 and admitted “that the quality of the field phosphorus data used in the model was problematic.” It is not clear why Brett et al. (2016) did not focus their analysis more broadly on years with better quality field data. A thorough analysis of the appropriateness of the model should also have included the main model calibration years of 1991 and 2000.

The year 2001 was chosen as the TMDL year because the hydrology best represented low river flow conditions. Moore and Ross (2010) explained, “The low river flow period is expected to be the most critical period for pollutant loading effects in the river and Lake Spokane due to less dilution of nutrient concentrations and a longer retention time, both of which can exacerbate dissolved oxygen shortages. By using a representative critical low flow year, the water quality in Lake Spokane and the Spokane River should be adequately protected…” Therefore, the 2001 year is a unique low-flow year that would not necessarily represent the same water quality conditions in Lake Spokane as other hydrological years.

Total phosphorus in the model

The TP_in formulation used in Brett et al. (2016) neglected 7 modeled constituents that contain phosphorus, including labile dissolved organic matter phosphorus (LDOMP), refractory dissolved organic matter phosphorus (RDOMP), labile particulate organic matter phosphorus (LPOMP), refractory particulate organic matter (RPOMP), and the phosphorus contained in the 3 algal groups. Brett et al. (2016) incorrectly calculated TP_in using the following equation: (1)

where PO₄ is orthophosphate, and CBODP_1-12 is the phosphorus contained in organic matter originating from 12 different dischargers and tributaries. During the summer months when metabolic rates are higher and travel time in the system is longer, phosphorus is cycled through the PO_4, LDOMP, RDOMP, LPOMP, RPOMP, epiphyton, and phytoplankton compartments. For example, CBODP will be oxidized and release PO₄, which then can be utilized for growth by free-floating algae (phytoplankton) and attached algae (epiphyton). Phosphorus contained in phytoplankton and epiphyton can pass into LDOMP and/or LPOMP through mortality, respiration, and excretion. LDOMP and LPOMP will then oxidize, transforming phosphorus back into PO₄ (to be utilized by phytoplankton and epiphyton again) or into the slower decaying RDOM-P and RPOM-P compartments. The correct equation to calculate TP_in is: (2)

where ALGP_1-3 is the phosphorus contained in the 3 algal groups. As mentioned earlier, during the critical summer months the difference between the value calculated by Brett et al. (2016) and the correct TP_in was largest. An example is the difference in magnitude between the correctly and incorrectly calculated TP_in for the TMDL scenario (Berger et al. 2009).

Computation of TP is clearly shown in the CE-QUAL-W2 User Manual (Cole and Wells 2016) and can be written out directly from the model. Because TP is not a state variable of the model but is computed from other constituents, the following formula is shown in the User Manual:

Total phosphorus: (3)

where δ_P is the stoichiometric ratio of P to organic matter; δ_PISS is the stoichiometric ratio of P to inorganic matter; Φ_RDOM is the concentration of refractory dissolved organic matter; Φ_LDOM is the concentration of labile dissolved organic matter; Φ_algae is the concentration of algae biomass; Φ_zooplankton is the concentration of zooplankton biomass; Φ_LPOM is the concentration of labile particulate organic matter; Φ_RPOM is the concentration of refractory particulate organic matter Φ_PO4 is the concentration of dissolved orthophosphorus as P; Φ_CBOD is the concentration of CBOD both dissolved and particulate; and Φ_ISS is the concentration of inorganic suspended solids.

Brett et al. (2016) also identified a “cryptic source of phosphorus” within the Spokane River model. Setting initial phosphorus concentrations and phosphorus in the inflows to zero for the Lake Spokane model, predicted phosphorus concentrations in the lake were zero, and we found no evidence of a cryptic source. CE-QUAL-W2 also performs an internal mass balance check by comparing the spatially and temporally integrated change in phosphorus mass since the start of a simulation. The model will calculate the mass error between these values and then calculate the percent mass balance error based on the total mass change. For orthophosphate (PO₄), the percent error was −1.041 × 10⁻⁷%, which is close to machine accuracy. The CBOD groups 1–12 also had similar mass balance accuracy. Hence, CE-QUAL-W2 conserves mass, and we found no evidence of a “cryptic source of phosphorus.” Their alleged cryptic source may have been related to an error in their computation of TP.

Mean hypolimnetic dissolved oxygen

Brett et al. (2016) stated that “the minimum volume-weighted dissolved hypolimnetic DO concentration (DO_min) was … calculated for reservoir Segment 36 in the model output file SPR6 for depths below 15 m to be consistent with the observational record (Welch et al. 2015),” but they used only one model segment near the dam, Segment 36, to assess minimum volume-weighted DO <15 m. By contrast, Welch et al. (2015) used a hypolimnetic average DO at longitudinal sampling stations 15 m below the water surface. Approximately 5 field stations could be used to compute hypolimnetic average DO (Fig. 2). In the mathematical model, the segment lengths near the dam were about 1 km. Brett et al. (2016) used approximately the last kilometer of the reservoir below 15 m and assumed it was representative of the entire hypolimnetic average DO used by Welch et al. (2015). The minimum summer hypolimnetic DO concentration of each model segment for the Lake Spokane TMDL scenario (Fig. 4) shows that the entire hypolimnetic average DO at the downstream segment was significantly lower than other model segments during 2001. Hence, the Brett et al. (2016) use of the model segment near the dam biased their results to low DO values. In addition, the volume of this segment below 15 m at the segment near the dam was only 8% of the total hypolimnetic volume below 15 m in Lake Spokane. Brett et al. (2016) biased their analysis by using an average of only a small fraction of the total hypolimnetic volume in the reservoir, which skewed their results toward lower DO with little variability. Their choice to base their calculation exclusively at this segment is surprising because they acknowledged that “[s]egment 36 … tends to have the lowest DO in the historical record” and that “[s]egments 9 through 36 represent all of the model segments deep enough to have thermally stratified water columns in the summer.”

Modeling the response of dissolved oxygen to phosphorus loading in Lake Spokane

All authors

Scott A. Wells & Chris J. Berger

https://doi.org/10.1080/10402381.2016.1211910

Published online:

05 August 2016

Figure 3. Total phosphorus inflows into Lake Spokane for the TMDL scenario during June through October 2001. The correctly calculated TP_in compared with the incorrect TP_in used by Brett et al. (2016). The difference is due to Brett et al. (2016) not including 7 modeled constituents consisting of phosphorus.

Display full size

Modeling the response of dissolved oxygen to phosphorus loading in Lake Spokane

All authors

Scott A. Wells & Chris J. Berger

https://doi.org/10.1080/10402381.2016.1211910

Published online:

05 August 2016

Figure 4. Minimum summer hypolimnetic DO concentrations of model segments for Lake Spokane TMDL scenario. Model cells below 15 m depth were used to determine the volume-weighted hypolimnetic DO. Segment 36 was the only segment used to calculate minimum hypolimnetic DO in Brett et al. (2016).

Display full size

The temperature and circulation pattern in the lake predicted by the model in summer 2001 (^{Figure 4}Fig. 5) show that the deeper segments near the dam were more stagnant as a result of stratification and outflow dynamics. The model predictions of water age and velocity (Fig. 6) show that water age, or the length of time water has resided in the reservoir, is longer for this segment than any other in the reservoir. A greater water age provides more time for DO to be depleted. The segment near the dam was not representative of the entire hypolimnion as a result of dead zones.

Modeling the response of dissolved oxygen to phosphorus loading in Lake Spokane

All authors

Scott A. Wells & Chris J. Berger

https://doi.org/10.1080/10402381.2016.1211910

Published online:

05 August 2016

Figure 5. Temperature contours and velocity vectors in Lake Spokane on 24 July 2001 as predicted by the TMDL model. Welch et al. (2015) defined hypolimnetic DO as the volume-weighted DO below a depth of 15 m, whereas Brett et al. (2016) only considered the reservoir volume below 15 m in the last model segment.

Display full size

Modeling the response of dissolved oxygen to phosphorus loading in Lake Spokane

All authors

Scott A. Wells & Chris J. Berger

https://doi.org/10.1080/10402381.2016.1211910

Published online:

05 August 2016

Figure 6. Water age and velocity distribution in Lake Spokane on 22 September 2001 as predicted by the TMDL model.

Display full size

Correspondence between model calibration year and field data correlation

Brett et al. (2016) referred often to Welch et al. (2015), who used field data to explore the relationship between changes in TP coming into Lake Spokane and lake hypolimnetic DO. A figure in Welch et al. (2015) showing the relationship between minimum hypolimnetic DO and inflow TP was used in Brett et al. (2016). Interestingly, field data from the model calibration periods—1991, 2000, and 2001, the year of the model TMDL evaluation—were not included in figure in Welch et al. (2015), and although an r² value of the correlation curve between TP_in and DO_other was provided, there were no other measures of expected error using the correlation.

Brett et al. (2016) claimed that field data showed that “during the hypereutrophic era” DO_MIN averaged 1.4 ± 2.7 mg/L and “currently (2010–2014)” DO_MIN averaged 6.5 ± 0.8 mg/L, but that the “model's DO_MIN response for similar TP_in concentrations averaged 3.8 ± 0.4 µg/L and 4.7 ± 0.04 mg/L, respectively, for analogous hypereutrophic and contemporary conditions.” For the years of field data used during the hypereutrophic period (1972–1985), inflows between June and October averaged 126 m³/s, and for the current conditions (2010–2014), inflows averaged 141 m³/s (Welch et al. 2015). As mentioned earlier, 2001 was a low-flow year with an average flow between June and October 2001 of only 77 m³/s. Should there be correspondence between hypereutrophic (1972–1985) and contemporary conditions (2010–2014) and 2001? Correspondence between different years must account for variations in inflows, outflows, meteorological forcing affecting the temperature regime, inflow density, and the sediment oxygen and nutrient flux dynamics. There may be little correspondence between a model set-up for 2001 and these other years. Welch et al. (2015) showed that the residence time had an impact on hypolimnetic DO. The approach of merely varying TP_IN with a model set-up for 2001 and saying the results correspond to other years is inappropriate. For instance, a higher flow year will result in greater dilution of phosphorus loads originating from wastewater treatment plants, causing decreased nutrient concentrations and reducing productivity. Higher flows would also decrease hypolimnetic residence time, reducing the depletion of DO.

Use of minimum volume-weighted hypolimnetic dissolved oxygen as a water quality indicator

Brett et al. (2016) used a volume-weighted minimum hypolimnetic DO value to evaluate the model response to changes in TP_IN. Welch et al. (2015) also used this same parameter as an indication of lake water quality as it relates to TP inflow concentrations. Does using a single parameter, the minimum volume-weighted DO concentration, provide a complete understanding of lake water quality response to TP concentrations? Although many water quality standards for DO are written in terms of a minimum DO or no change from natural conditions (Moore and Ross 2010), the Environmental Protection Agency (EPA 1986) suggested the use of the mean and mean minimums of DO concentrations over 7 to 30 days.

We can postulate many examples in which a lake can have the same volume-weighted minimum DO but have very different mean 7- or 30-day mean minimum DO over a summer period and a very different impact on lake water quality. Consider for example 2 hypothetical curves of volume-weighted hypolimnetic DO (Fig. 7). Both curves have the same minimum volume-weighted hypolimnetic DO for the June–October period, but one has a summer average of 6.6 mg/L, which is above the 30-day mean DO for cold water criterion for fish (EPA 1986), whereas the other has an average of 4.4 mg/L, which is below the 30-day mean DO cold and warm water criteria for fish (EPA 1986). Although the hypolimnetic volume-weighted minimum DO can be one of several valuable indicators of lake or reservoir water quality, taken alone it may not be indicative of the distribution of DO that can affect fish and aquatic life.

Modeling the response of dissolved oxygen to phosphorus loading in Lake Spokane

All authors

Scott A. Wells & Chris J. Berger

https://doi.org/10.1080/10402381.2016.1211910

Published online:

05 August 2016

Figure 7. Comparison example of hypothetical June–October hypolimnetic volume-weighted DO in which both curves have the same minimum volume-weighted hypolimnetic DO but very different ecological impacts based on cold and warm water DO criteria (EPA 1986).

Display full size

Modeling the response of dissolved oxygen to phosphorus loading in Lake Spokane

All authors

Scott A. Wells & Chris J. Berger

https://doi.org/10.1080/10402381.2016.1211910

Published online:

05 August 2016

Figure 8. June–October time-averaged and volume-weighted distribution of DO concentrations in Lake Spokane for the 2001 calibration simulation and TMDL scenario. Only model cells at depths >15 m (hypolimnion) were included in the calculation. Concentrations were time-averaged using time intervals equal to the model time step, which varied between 2 and 300 s. Values were then weighted according to the time step length. The overall time- and volume-averaged DO was 5.6 mg/L for the 2001 calibration and 6.3 mg/L for the TMDL scenario.

Display full size

Modeling the response of dissolved oxygen to phosphorus loading in Lake Spokane

All authors

Scott A. Wells & Chris J. Berger

https://doi.org/10.1080/10402381.2016.1211910

Published online:

05 August 2016

Figure 9. Model-predicted DO_min and TP_in for the 2001 calibration simulation and TMDL scenario compared with equation (y = 220.45x^−1.295) developed in Welch et al. (2015). The 2001 calibration model used a maximum zero order SOD of 0.6 g/m²/d in Lake Spokane, whereas the TMDL scenario used a maximum of 0.25 g/m²/d.

Display full size

For example, using the TMDL model for summer 2001, the model shows details of the spatial and temporal dynamics of DO concentrations <15 m (Fig. 8). Comparing the 2001 calibration to the TMDL scenario, the percentage of hypolimnetic volume with very low DO concentrations (<2 mg/L) was reduced from 22% to 3.6%. These histograms, in contrast to a minimum volume-weighted DO concentration, have more information to assess model performance and response to upstream nutrient loading.

Algal group representation

Brett et al. (2016) criticized the TMDL model saying that the “…half-saturation constant for phosphorus is set to 3 µg/L for all 3 algal groups,” “…simply put, the model … treats cyanobacteria as (practically) equal competitors for phosphorus uptake,” and “the general designation of phytoplankton minimizes the likelihood of phosphorus limitation in the system.”

Setting up a water quality model of a waterbody usually has limitations that may not exist in natural systems. Natural freshwater systems usually have numerous algal species present, whereas the TMDL model used only 3 algae groups. The first 2 algae groups were not identified explicitly in our TMDL model documentation, but the 3 groups could represent broadly mixed diatoms, mixed green algae, and mixed cyanobacteria algae. As shown in Bowie et al. (1985), significant differences between diatoms, greens, and cyanobacteria exist with regard to temperature preferences. In modeling applications using just one algal group, the modeler often specifies a general broad temperature preference and one set of half-saturation constants for the entire mixed population of algae. The changes in algal species dominance in the Lake Spokane model was a result of temperature preferences rather than phosphorus competition. For example, the mixed diatoms maximum growth rate was at 10–18 C, mixed green algae 16–20 C, and mixed cyanobacteria >27 C. The half-saturation constants were set to the same value, implying no competition by P limitation, only by temperature preference. Ideally, many more than 3 mixed algal groups would be represented in the model with overlapping temperature preferences and varying half-saturation constants. It would be challenging, however, to model more than 3 algae groups without additional data such as algae speciation, algae biomass, and kinetic coefficient data. For most TMDL and modeling projects, these data are often not available due to time and cost restraints.

Were the half-saturation constants used in the model outside the range found in the literature? Bowie et al. (1985) showed that half-saturation constants for phosphorus for total phytoplankton can range from 0.0005 to 0.02 mg/L. The maximum values of the half-saturation coefficient can be even higher, up to 1.0 mg/L, for some species (Jeanjean 1969, Spencer and Lembi 1981). Bowie et al. (1985) showed that phosphorus half saturation constants can vary from 0.001 to 0.025 mg/L for diatoms, 0.002 to 0.03 mg/L for greens, and 0.0025 to 0.06 mg/L for cyanobacteria. Hence, phosphorus half-saturation constants between large mixed algal groups can indeed overlap. Half saturation constants in the Lake Spokane model were calibrated by comparing nutrient, chlorophyll a (Chl-a), and pH data with model predictions.

Using the TMDL model for predicting the lake recovery

Brett et al. (2016) stated that using the TMDL model was a “unique opportunity to test the ability of mechanistic models to simulate natural conditions.” The premise for this test was figure 3 from Welch et al. (2015), in which minimum volume-weighted hypolimnetic DO (<15 m) was compared to a volume-weighted riverine TP at station LL5 (see Fig. 2). It might seem a relatively easy task to use the 2001 TMDL model and vary the TP in the inflow and record the change in hypolimnetic DO to see whether the model reproduces figure 3 from Welch et al. (2015). Just changing TP in the inflow, however, merely reflects changes in TP for the conditions of 2001, using its hydrology, meteorology, operational conditions, and sediment conditions. Because the TMDL model year 2001 was an extreme year, we may not expect 2001 to follow a relationship similar to figure 3 in Welch et al. by just varying TP_IN.

In their simulations, Brett et al. (2016) applied the TMDL model by assuming that the model's zero order sediment oxygen demand (SOD) was constant and independent of upstream phosphorus loading. In the 2001 calibration simulation, SOD of model segments was set to values between 0.1 and 0.6 g/m²/d, with higher SOD in model segments closer to the dam. For the TMDL scenario, which simulated reduced upstream phosphorus loadings, SOD was reduced to a maximum of 0.25 g/m²/d (PSU 2010).

The correctly calculated model predicted DO_min and TP_in for the year 2001 calibration simulation and TMDL scenario were compared them with the regression equation developed in Welch et al. (2015; Fig. 9). Because the TMDL model used a fixed SOD for 2001, the model user would also need to assess how the zero order SOD rate would be affected by changes in TP loading for other years because the zero order model is not linked to processes in the water column. For example, in modeling 2010–2014, in addition to using hydrologic, operational, and meteorological data from each of those years, different values of SOD for those years would also need to be used rather than the SOD for 2001. Even using a predictive sediment diagenesis model (such as in CE-QUAL-W2 Version 4), the model user would still have to adjust the initial conditions of the sediment for a simulation if using the model for only one year. Brett et al. (2016) did note the model sensitivity to SOD but neglected to account for how SOD would change with TP loading.

Model calibration

Brett et al. (2016) claimed that “the original DO calibration may have been achieved by adjusting the particle settling velocity and organic matter degradation submodels within CE-QUAL-W2 until observed and modeled hypolimnetic DO concentrations matched.” Because they did not see the TP_in–DO_MIN relationship match with Welch et al. (2015), they concluded that the “calibration process achieved the correct answer for the wrong reasons.” The largest improvement in DO profiles came not from settling velocity or organic matter decay rate adjustments, but from adjustment of the wind field based on local information. Wind direction data originally used by the model was collected at the Spokane Airport and not at the reservoir. Using anecdotal observations of the wind field on the lake, adjustments were made to airport wind data to represent more accurately the wind field on the lake. This change affected lake circulation and improved model predictions. If Brett et al. (2016) had evaluated all the model calibration years (1991, 2000, and 2001), they would have noted that the model used the same values of organic matter settling velocities and organic matter degradation rates for all years. As shown in Berger et al. (2009), the modeled DO fluxes during 2001 for Lake Spokane were dominated by algal and epiphyton production, reaeration, and organic matter decay.

Peer review of model studies

Models are complex and require external peer review to assess if mistakes or differences in professional judgment need to be resolved during model set-up, calibration, application to management scenarios, or interpretation of model results. The National Research Council (2007) includes “peer review” as part of the model evaluation process and recommends evaluation of any modeling tool for the following items:

EPA (2015) also includes guidelines for peer review and states the reason for the peer review of the model:

Because of the complexity of water quality and hydrodynamic models used to set water quality standards, peer review should be a way to improve the technical competence of a model study or even to guide the interpretation of another's review. Successful peer review is not alone attained by a large number of authors or by journal manuscript review. To improve the quality of a modeling study, peer review benefits from collaboration between reviewers and developers.

Summary

Brett et al. (2016) analyzed the relationship between the mean summer TP loading into Lake Spokane (TP_in) and the minimum summer hypolimnetic DO (DO_min) using the TMDL model of Berger et al. (2009). They concluded that the model did not “reproduce the fundamental causal relationships underlying the eutrophication problem.” We disagree with their critique of the model behavior for the following reasons: