In the summer of 2010, the CDPHE conducted a High-Quality Water Supply study to assess the impact of algal growth in Colorado reservoirs on TOC concentrations and the potential to form DBPs. Most of the 30 reservoirs (Figure 1 and Table 1) have been classified for the Water Supply use, although only a few currently provide raw water for treatment facilities. The reservoirs varied in elevation (∼1,100 m a.s.l–∼3,000 m a.s.l), area (110–17,875 acres), capacity (1,330–829,600 acre ft), and secchi depth (0.3–5.1 m). The reservoirs were chosen to include a balance of lakes in the mountains and the plains, representing diverse climatic and watershed characteristics across Colorado. The lower elevation reservoirs tend to be shallow and have warmer temperatures compared to the high elevation reservoirs. In addition, most of the plains reservoirs are strongly influenced by agriculatural runoff and discharge of treated wastewater from large urban areas along the Front Range. In this study, the reseroivrs at and above 1700 m were considered mountain reservoirs (n = 13) and those below 1700 m a.s.l were included in the plains category (n = 17). Many of these reservoirs were located in state or county parks, wildlife and/or recreation areas where fishing, boating, hiking and camping are common actvities.

The 30 lakes were sampled during July and August (Table 1) at the peak of summer stratification. The Boulder and Westminster drinking water reservoirs were sampled biweekly from May through September 2010. The CDPHE staff collected 1,000 ml surface-water grab samples from the 30 lakes. The biweekly sampling was conducted by staff at the respective utility. The raw water for DOC, TOC and SUVA₂₅₄ was stored in pre-cleaned and pre-combusted amber glass bottles and stored on ice. Nutrient samples were collected in plastic HDPE bottles and stored on ice. Upon delivery to the University of Colorado by the following day, water samples were immediately refrigerated at 4°C. For 25 of the 30 lake samples which were analyzed for spectroscopic analysis and DOC, the samples were filtered through 0.7 µm pore size pre-combusted glass fiber filters (GFF). A 1000 ml sample of raw water was transferred to the Kiowa Laboratory at the Institute of Arctic and Alpine Research for TOC, Total Nitrogen (TN), inorganic nitrogen, and Total Phosphorous (TP) analysis. In addition, 1,000 to 1,500 ml of water was filtered on 0.7 µm GFF filters and frozen for Chl-a analysis. The analysis of phytoplankton for 14 of the 30 lakes was done for 500 ml surface grab samples collected in amber Nalgene bottles, preserved with Lugol’s solution, and stored at 4°C until analysis.

TOC, DOC and TN were measured on raw water using a Shimadzu TOC-V CSN at the Kiowa Laboratory. Total inorganic nitrogen species (nitrate, nitrite and ammonium), as well as TP, were analyzed at the Kiowa Laboratory on a Lachat 8,500. Chl-a concentrations were measured using a fluorescence method with acetone extraction (Wetzel and Likens, 2000). Accurate dilutions of a Chl-a standard from Fisher Scientific were analyzed to establish a linear calibration curve and a linear calibration factor between the UV-Vis Spectrophotometer and the Fluoromax 2 (F2) fluorometer. The fluorescence method was used due to higher sensitivity than the UV-Vis Spectrophotometer method. 3 ml aliquots of extracted Chl-a were measured on the F2 and were corrected for phaeophytin through acidification of 100 µL of HCL. Chl-a concentrations were determined using the formula established by Wetzel and Likens (2000).

Analysis of the phytoplankton using the FlowCAM^® was conducted to understand the relationships between the phytoplankton species, Chl-a, nutrient levels, physical characteristics of the lake, and TOC concentrations. For a subset of 14 reservoirs (eight mountain and six plains), phytoplankton identification and enumeration was conducted with a Fluid Imaging Technologies FlowCAM^®. Unlike traditional microscopy, the FlowCAM^® enables rapid monitoring of particles in fluid by combining flow cytometry with microscopy. The FlowCAM^® automatically counts and images each particle, while also evaluating characteristics of the digital image, such as shape and intensity. Such imaging microscopes are becoming more commonly used by water treatment utilities to monitor phytoplankton abundance in source water lakes and reservoirs, as well as to identify invasive species.

For these analyses, 150 ml of the 500 ml grab sample was transferred to a 250 ml centrifuge tube and settled for 24 h. Then, 130 ml of the sample was aspirated from the top of the sample in a manner that did not disturb the settled particles. The concentrated sample was then transferred to a 50 ml centrifuge tube. If the sample looked visibly cloudy, it was filtered with a 100 µm mesh net to avoid clogging in the flow cell. The 10X objective was used with a 100 µm flow cell. Acetone was run for five minutes to clean the flow cell and tubing. The FlowCAM^® was focused using a small volume of spare sample. A 2 ml of sub-sample was then run through the FlowCAM^®. After the sample finished running, image library files were made through the interactive data platform and sorted based on image characteristics associated with each of three dominant algal groups, Cyanobacteria, diatoms and chlorophytes (green algae). An automated filter was set up based on image analysis, and then files were verified by human-eye. Algal groups were then quantified by assessing the concentration as cells or colonies per mL, taking into account the amount of sample imaged. Area based diameter (ABD) calculated by the FlowCAM^® was then multiplied by the number of counts in the file for each algal group to get an approximation of biovolume. The number of cells or colonies for each group with respect to total count was multiplied by ABD in order to calculate % biovolume of each group. Further linear regressions in Excel and two-tailed t-tests were applied to interpret key physical relationships driving the DOM quality.

A subset of 25 filtered water samples from the 30 lakes sampled were analyzed on a Fluoromax 3 to obtain Excitation Emission Matrices (EEMs), following the methods of Cory and McKnight (2005). The full sample set was not analyzed due to logistical constraints. The FI (McKnight et al., 2001; Cory and McKnight, 2005) was calculated as an indicator of the DOM contribution from phototrophic microbes on a scale of 1.3 (terrestrial source) to 1.9 (microbial source). The Parallel Factor Analysis Model (PARAFAC) model developed by Cory and McKnight (2005) was used to quantify the loading and relative distribution of components. Although the samples span lakes and reservoirs from both mountains and plains across the state of Colorado, this individual regional dataset is not subjected to the accuracy and sensitivity limitations of PARAFAC analysis models of diverse samples that have been observed across large global datasets (Yu et al., 2015). The percent of the total fluorescence associated with protein-like components was calculated by summing the contributions of components 8 and 13. The distribution of quinone-like components was used to calculate the redox index (RI) as described by Miller et al. (2009). In surface water lakes, RI is a potential indicator for the extent of photobleaching of the humic fraction of the DOM, which is also related to the residence time of the DOM in the lakes. For example, Miller et al. (2009) showed that the RI n a high alpine lake in CO decreased throughout the summer as the residence time increased. Photobleaching also has the potential to reduce SUVA₂₅₄ independent of the DOM source.

The formation of DBPs by chlorination of reservoir water was evaluated using the uniform formation conditions method (Summers et al., 1996), which provides a representative simulation of typical U.S. drinking water distribution systems. Samples were chlorinated with a sodium hypochlorite solution, adjusted to pH 8 with a borate buffer and incubated in the dark for 24 h (±1 h) at room temperature (20 ± 1.0°C). Chlorine doses were determined from demand studies to obtain a chlorine residual of 1.0 mg/L (±0.4 mg/L) after 24 h. Chlorine residuals were measured with the N,N-diethyl-p-phenylenediamine, or DPD colorimetric method, Standard Method 4500-Cl G (Standard Methods 1998) and were quenched with ammonium chloride. The four regulated TTHMs (EPA Method 551.1) and the five regulated haloacetic acids (EPA Method 552.2) were measured by gas chromatography with an electron capture detector.

All analyses were based on data from the synoptic survey of lakes (Table 1) and were conducted using R v3.5.1 (R Core Team 2018). Pairwise relationships were assessed using simple linear regression based on ggplot2::smooth, and plots include a confidence interval that represents the standard error of the predicted value (Wickham 2016). The analysis of models for DBP efficiency that included the natural log of the DOC concentration and DOM quality indices is shown in Table 2. The models were evaluated based on the difference between the AIC_C for the model with the lowest AIC_C and the AIC_c for that model, referred to as deltaAIC (Burnham and Anderson, 2004).

We explored linear models explaining variation in DBP yield efficiency (ug DBP/mg DOC) each for TTHMs or HAA5s. We used the MuMIn package (Bartoń, 2018) for R to compare all possible candidate models using AIC_C scores to evaluate the plausibility of each model (Burnham and Anderson, 2004). The global model in the analysis was y ∼ FI + RI + SUVA₂₅₄ + ln (DOC +1), where y is the ln (x+1) transformed yield efficiency of the DBP of interest, FI is fluorescence index, RI is redox index, SUVA₂₅₄ is Specific UV Absorbance at 254 nm, and DOC is mg L⁻¹ dissolved organic carbon. All predictor variables were centered to a mean of 0 and standardized to unit variance. When global models were assessed, all predictor variables had a variance inflation factor <3, indicating the predictor variables included in the analysis were not multicolinear. During the model comparison analysis, we considered the global model as well as all possible combinations of predictor variables and a null model where the only predictor was the global mean of the data. We report unstandardized model coefficients because all data in the models have been centered and scaled. Only models with a delta AIC_C < 10 were considered in the final model comparison results. For each analysis, all lakes/sites were used for which the necessary data were available. In all cases, there was good representation (in terms of replicated observations) across the variables being assessed. Each model analysis compared candidate models against null models using AICc values, which provides a method to evaluate if/when a fixed effect is more likely, given the data, than the null hypothesis. Our sample sizes were sufficient to detect statistically significant effects given our results.

The nutrient concentrations in the reservoirs studied encompassed a broad range (0.16–1.98 mg/L TN) (Table 3). TP to TN molar ratios are compared to the 16:1 N:P Redfield ratio and the 23:1 Wetzel (2001) ratio that is considered to be an indicator of severe P limitation (Figure 2). Furthermore, higher total inorganic nitrogen (TIN) concentrations were associated with greater algal biomass as measured by Chl-a concentrations (Figure 3), with the strongest log-log relationship being between NH₄ ⁺ and Chl-a.

The characterization of the phytoplankton community composition was applied for a subset of 14 reservoirs. Filamentous cyanobacteria dominated the %Biovolume of phytoplankton species in five of the eight mountain reservoirs and in five of the six plains reservoirs (Figure 4). A representative FlowCam^® image of filamentous cyanobacteria from Prewitt Reservoir (Supplementary Figure 1) shows the presence of distinctive heterocysts. The greatest frequency of heterocysts for the filamentous cyanobacteria were observed in a mountain reservoir, Grand Lake (GR), with about 28% of the imaged filaments containing one or two heterocysts. The total nutrient concentrations in Grand Lake are at the low end of the range in this data set (Figure 2). In the three plains reservoirs for which heterocysts observed, Jumbo Reservoir (JB), Jackson Reservoir (JR), and Prewitt Reservoir (PW), the frequency ranged from 6-8% of the imaged filaments containing heterocysts. These three reservoirs are at the high end of the range for total nutrients (Figure 2).

The TOC concentrations in the 30 reservoirs ranged from a high value of 11.5 mg C/L to a low value of 2.1 mg C/L, and a strong positive log-linear relationship was found between Chl-a concentrations and TOC concentrations (Figure 5). To consider the details of this relationship in the context of DBP formation, the 30 reservoirs were grouped by TOC range, high (>8.1 mg/L, moderate (4.1–8.0 mg/L), and low (≤4.0 mg C/L) (Table 3).

As noted by Saunders et al. (2015), characterization of DOM chemical quality using spectroscopic indices provides complementary data for evaluating the underlying processes that yield the relationships found between nutrients, phytoplankton and TOC concentrations. The SUVA₂₅₄ values of the 30 reservoirs ranged from 1.1 to 3.3 L/mg-m (Table 3). The fluorescence index (FI) ranged from 1.28 to 1.63, encompassing the range generally found in surface waters, with higher values corresponding to DOM derived from phytoplankton, such as microbes and algae. Figure 6 shows that FI was strongly positively correlated in a log-linear relationship (r ² = 0.63) with Chl-a, and that the reservoirs with the high FI values (>1.50) generally had low SUVA₂₅₄ values (<2.0 L/mg-m: JM, AC, NS, HE, PW, JK, JB, BR) and all of these were plains reservoirs. Further, when the FI was used as a multiplicative factor for TOC concentration to increase the weight of algal-derived DOM, the strength of the log-linear relationship with Chl-a concentration increased slightly (Figure 5B). The relationship of SUVA₂₅₄ to the N:P ratio was somewhat weaker (r ² = 0.26) than the relationship for FI (r ² = 0.33). The relationship of FI to the log of the N:P ratio was also significant, but not as strong as the relationship with Chl-a. The plains reservoirs which had the highest FI values also had low N:P ratios, below the 23:1 threshold, and were associated with high TN. The high FI, low SUVA₂₅₄, high nutrients and relatively low N:P are consistent with wastewater influence common to all of these reservoirs.

Other characteristics of the fluorescence spectra that were quantified using the Cory and McKnight (2005) model did not show significant relationships with Chl-a or N:P ratios (Figures 6, 7). In particular, the percent of the fluorescence intensity accounted for by protein-like fluorescence had a broad range from 23% to very low values, with high values occurring for reservoirs along the full range of Chl-a concentrations and SUVA₂₅₄ values. Therefore, because %protein had a similar pattern to the RI (Figure 7), and represents a smaller portion of the EEM, which could be subject to more error than the RI, we included RI as a candidate predictor for DBP yield efficiency, but not %protein.

As expected, the highest formation of the two regulated organic DPBs, TTHMs and HAA5s occurred in the reservoirs with high TOC concentrations comprised predominantly of plains reservoirs (Table 3). However, it should be noted that increased HAA5 formation occurred in a wide range of TOC concentrations, overlapping across low, moderate and high TOC groups. Whereas, for TTHM formation, the range in the low TOC group of reservoirs (2–4 mgC/L) comprised predominantly of mountain reservoirs did not overlap with the range in the high TOC group (8–11 mgC/L).

The models for predicting TTHM and HAA5 yield efficiency differed in level for explaining variation and in the influence of spectral characteristics of the DOM (Table 2). None of the plausible models for TTHM or HAA5 based on a deltaAIC of less than five included only lnDOC; however, three of the seven plausible models for HAA5 did include a positive influence of ln DOC. For TTHM efficiency, there were seven plausible models and the null model was marginally plausible with a deltaAIC of 4.945. The model which included only SUVA₂₅₄ had the lowest AIC_C value, but had a low level for explaining variation (r ² = 0.271). Furthermore, all of the other six plausible models also included SUVA_254. as a predictor variable with a positive influence, and had a low level for explaining variation. Addition of FI as a predictor variable yielded the model with the lowest deltaAIC and had a positive influence of FI.

For HAA5 efficiency there were also seven plausible models and one that was marginally plausible. All models included a positive influence of SUVA₂₅₄ and had high levels of explaining variation, with r ² ranging from 0.79 to 0.815. The null model was not plausible at all. The model with the lowest AIC_c for HAA5 efficiency included a negative influence by FI.

Most the lakes and reservoirs studied had TN:TP ratios exceeding 23, indicating that P deficiency predominated during the mid-summer sampling campaign. In general, as TN and TP concentrations increased, the ratios decreased consistent with expectations for anthropogenic nutrient loading. For example, Barr Lake, which receives treated wastewater and urban runoff from the Denver area, was the only lake with a N:P ratio that indicated N deficiency (Figure 2); Barr Lake also had high TN (1.6 mg/L) and Chl-a (21.2 μg/L).

Road Canyon Reservoir was the only mountain reservoir among the eight reservoirs in the high TOC range (8.1 mgC/L to 11.5 mgC/L). The plains reservoirs in this group were Barr, Fruitgrowers, Henry, Jackson, Milton-Greeley Seaman, North Sterling, and Prewitt Reservoirs. This high TOC group also contained some of the highest concentrations of Chl-a (43.5 μg/L) and TN (2.0 mg/L), along with higher percentages of total inorganic nitrogen (%TIN). Road Canyon Reservoir had both high TN (2.0 mg/L) and TP (0.27 mg/L) values. Except for Barr Lake, the other reservoirs had N:P ratios close to or above the N:P ratio of 23 indicating P limitation. The three reservoirs with complete datasets in this group (Jackson, Prewitt and Road Canyon) all had phytoplankton communities dominated by filamentous cyanobacteria.

The 11 reservoirs with moderate TOC concentrations (4.1 mgC/L–7.99 mgC/L) included roughly the same number of mountain reservoirs (5) and plains reservoirs (6). Several of these reservoirs had low Chl-a concentrations, and all had TN:TP ratios indicating P limitation. Diatoms were dominant in three mountain reservoirs (Antero, Skaguay and Steamboat Reservoirs), and filamentous cyanobacteria were dominant in two mountain reservoirs (Granby Lake, and Stagecoach Reservoir) and one in the plains (Jumbo Reservoir). It is also notable that the phytoplankton community was dominated by filamentous cyanobacteria, which commonly thrive in warm lakes in the summer (Goldman and Horne, 1983), for 8 of the 11 reservoirs in the high and moderate TOC groups for which the phytoplankton community was characterized.

The group of 11 reservoirs with TOC concentrations between two and 4 mgC/L all still exceeded the lower limit for TOC removal prior to chlorination (2 mgC/L). This set included seven mountain reservoirs and 4 in the plains. These reservoirs generally had lower Chl-a concentrations (average 1.8 μg/L ± 1.3 μg/L), and also had lower TN (average 0.3 ± 0.2 mg/L) and lower %TIN. Filamentous cyanobacteria and diatoms were both abundant in the two mountain reservoirs for which the phytoplankton were characterized (Elevenmile Reservoir and Grand Lake), with the abundance of cyanobacterial filaments containing heterocysts in grand lake suggesting the possibility for nitrogen fixation occurring. In contrast, filamentous cyanobacteria dominated in one plains reservoir (Boulder Reservoir), and diatoms dominated in the other (Westminster Reservoir).

Given the evidence for high TOC concentrations being an outcome of TOC production by summer phytoplankton blooms, the findings of the highest formation of the two regulated organic DPBs, TTHMs and HAA5s, occurring in reservoirs with high TOC concentrations is generally consistent with findings of Saunders et al. (2015). In particular, the results in Table 3 indicate that controlling nutrient loadings that promote such blooms will limit the likelihood of excessive TTHM formation at drinking water treatment facilities.

A further question is whether or not the influence of algal production on DOM quality results in a change of the efficiency of DBP formation from DOM. All the plausible models for TTHM and HAA5 efficiency included a positive influence of SUVA₂₅₄, and the models for HAA5 efficiency had high levels of explaining variation compared to models for TTHM efficiency. The inclusion of FI contributed to the level of explaining variation for both TTHM and HAA5 efficiency, however with opposing influences. Because photobleaching and other processes occurring in the water column of the reservoir may have a greater influence on compared to FI, here may be value in including FI as an additional spectral property for some reservoirs used for water supply.

It is notable that three plausible models for HAA5 efficiency included a positive influence of lnDOC. Other influences on DOM quality beyond algal growth, such as wastewater inputs, may not be directly reflected in FI or SUVA₂₅₄. Thus, this finding may reflect the influences of other DOM sources in contributing to greater DBP formation efficiency. Nonetheless, the results of this analysis with a broad range of reservoirs strongly support the finding of Saunders et al. (2015) that SUVA₂₅₄ is a highly valuable predictor of DBP formation efficiency.

Overall, the strong relationships between Chl-a, N:P and FI provides additional support to the interpretation that nutrient loading is a driver for high TOC concentrations by increasing phytoplankton productivity in these reservoirs. Further, the stronger relationships for FI compared with SUVA₂₅₄ suggest that fluorescence measurements may be useful complementary measurements to SUVA₂₅₄ for establishing the contribution of phytoplankton growth to DOM production for some drinking water supply reservoirs.

In general, the reservoirs in the low TOC group were predomiantly located in the mountains and conversely the high TOC group was composed of plains reservoirs with one exception. Given the broad range of reservoirs and conditions represented in this study, we expect that these results will be useful for comparison to other regions that have both mountain and plains lakes, especially since monitoring and treating freshwater HABs is an increasingly important issue in the global drinking water supply (e.g. Backer et al., 2015; Schmale et al., 2019; Glibert, 2020). However, it should be noted that other regions may have different regional factors such as climate, altitude and precipitation patterns that may influence local water quality and treatment, warranting further investigation. As mentioned earlier, in Spring 2012, the Colorado Water Quality Control Commission adopted new regulations based on the results of this study (excluding phytoplankton analysis). As part of this plan, a Chl-a criterion is now available for protecting water supplies classified as Direct Use Water Supply Reservoirs (Colorado Water Quality Control Commission, 2013). The additional DOM optical characterization and phytoplankton identification suggests that monitoring algal species, and specifically filamentous Cyanobacteria, may provide useful ancillary information for understanding the drivers of high DOM concentrations during the summer.

There is a potential for watershed management to control the sources of excess nutrients to water supply reservoirs, which can include agricultural runoff, leakage from residential septic systems and upstream wastewater treatment effluent, augmented by atmospheric deposition of nitrogen. Therefore, prevention of elevated TOC levels through pre-treatment by management of nutrients in source waters could be more efficient and cost effective for drinking water utilities than removal of TOC that exceed standards. The results further suggest that watershed management could provide further benefits of limiting blooms of filamentous Cyanobacteria, which have the potential to develop into toxic blooms.

Furthermore, longer periods of ice-free conditions in spring and summer may cause greater light levels to support algal growth. Similarly, warmer lake temperatures could help promote the growing season for filamentous Cyanobacteria, which were commonly associated with DOM derived from microbial phototrophs. As a result, these findings may help future drinking water utilities plan and adapt to earlier lake cover ice out as a result of a changing climate (Sharma and Magnuson, 2014). Additionally, drought years, which result in a smaller snowpack, may increase the need for reservoir storage. Therefore, improved understanding of how dominant phytoplankton species and algal-derived DOM impact the formation of DBPs could help provide targeted approaches for water utilities in maintaining storage water quality in a changing climate.