Seasonal and longitudinal water quality dynamics in three effluent-dependent rivers in Arizona

Site description

We sampled water quality across six effluent-dependent river reaches in Arizona: three reaches in the Santa Cruz River basin, two reaches in the Rio de Flag basin, and one reach in the Salt River basin (Fig. 1). Groundwater levels are at least 40 m below the streambed surface for all six reaches, so each reach would dry completely during baseflow if effluent was not being added to the channel (Carlson et al., 2011; Arizona Department of Environmental Quality, n.d). Additionally, all reaches with the exception of the Rio de Flag, fall within Active Management Areas of the Arizona Department of Water Resources, and groundwater withdrawals are monitored to achieve safe yield. Study reach lengths ranged from 3 to >24 km. The reach in the Salt River and both reaches of the Rio de Flag were short (~3–5 km). In the Santa Cruz River, one reach was medium in length (~9 km) and two others were long (>24 km) (Table 1). Downstream of all six study reaches, the channels are dry except during periods of flooding.

The Santa Cruz River basin (22,000 km²) is located in southern Arizona and northern Mexico (Webb et al., 2014). The climate is characterized by hot summers (average high July temp: 38 °C) and moderately cool winters (average low January temp: 5 °C). Annual average precipitation is ~300 mm and rainfall is bimodal, with the months of January and August being the wettest and May and November being the driest (data provided by NOAA ESRL Global Monitoring Division, Boulder, Colorado, USA (http://esrl.noaa.gov/)). The uppermost study reach in the Santa Cruz Basin (~25 km long) is supported by the Nogales International Wastewater Treatment Plant. This plant was constructed in 1943, upgraded in 2009, and discharges approximately 57 million liters of tertiary-treated effluent into the river each day (Santa Cruz County, 2019). The second Santa Cruz reach is shorter (~9 km) and is supported by the Agua Nueva Water Reclamation Facility (WRF), a relatively new facility completed in 2014 that replaced the Roger Road Wastewater Treatment Plant constructed in 1951, and discharges approximately 30 million liters of tertiary-treated effluent into the river each day (Sonoran Institute, 2017). The third and longest (~31 km) reach is supported by the Tres Rios WRF, which was built in 1978, upgraded and modernized in 2014, and discharges approximately 115 million liters of tertiary-treated effluent each day (Sonoran Institute, 2017). In addition to baseflow from the treatment plants, all three study reaches in the Santa Cruz River experience seasonal flood events from precipitation runoff, which can surpass 280 m³/s. The three perennial study reaches are separated by dry reaches (Fig. 1) and are only briefly connected during large flood events.

The Salt River basin (35,000 km²) is located in central Arizona, with its headwaters in the White Mountains of eastern Arizona. The study reach is characterized by hot summers (average high July temp: 41 °C) and mild winters (average low December temp: 7 °C). Annual average precipitation is ~204 mm, and rainfall is bimodal. The months of March and July are the wettest, and June and October the driest (data provided by NOAA ESRL Global Monitoring Division, Boulder, Colorado, USA (http://esrl.noaa.gov/)). Our study reach (~6 km long) is supported by the 91st Avenue Wastewater Treatment Plant near Phoenix, which began discharging into the river in 1958, was upgraded in 2011, and processes an average of 530 million liters of tertiary treated effluent per day (US EPA, 2016).

The Rio de Flag basin (518 km²) is located in the conifer forest-dominated, high elevation mountains (>2,000 m) of north-central Arizona (Figs. 1 and 2). The climate is characterized by moderate to warm temperatures in summer (average high July temp: 27 °C) and cold winters (average low January temp: −12 °C). Annual average precipitation is ~588 mm, and rainfall is bimodal. The months of February and August are wettest, and June and November are driest (data provided by NOAA ESRL Global Monitoring Division, Boulder, Colorado, USA; http://esrl.noaa.gov/). However, because of underlying volcanic geology, only 0.25 cm per year becomes stream flow within the Rio de Flag, making it one of the lowest percentage runoff systems in Arizona (Bills & Enyedy, 2015). Our ‘Central’ reach (~3 km long) is supported by the Rio de Flag Water Reclamation Plant in central Flagstaff. This plant began discharging into the river in 1993, was upgraded in 2009, and discharges approximately 2.4 million liters of tertiary treated effluent per day into the reach (City of Flagstaff, 2020; personal correspondence: Erin Young, City of Flagstaff). Our other Rio de Flag study reach, ‘Wildcat’ (~4 km long), is located approximately 10 km downstream of the Central reach, and is supported by the Wildcat Hill Water Reclamation Plant. This plant began discharging into the river in 1971, was upgraded in 2008, and discharges approximately 12 million liters of tertiary treated effluent per day into the reach (City of Flagstaff, 2020; personal correspondence: Erin Young, City of Flagstaff).

Data collection

To understand how water quality changes downstream from effluent outfalls, at least three sites were sampled for water quality within each of the six study reaches. The exact location and number of sampling sites for each reach were selected based on the availability of public access points and the length of reach, resulting in a range of three to eight sampling sites per reach (Table 1). A minimum of 0.8 km separated each progressive downstream sampling site within a reach, with an average distance of 3.5 km between sites (Fig. 1). To examine seasonal variation in water quality, we took measurements (described below) at each site at least once per season (winter = January to March, spring = April to June, summer = July to September, fall = October to December) in the year of 2018. To provide additional temporal resolution, one reach was sampled bi-monthly (Santa Cruz River at Nogales: 6x per year) and one reach was sampled monthly (Santa Cruz River at Tres Rios: 12x per year) (Table 1).

During each sampling event, we measured basic water quality parameters (dissolved oxygen, temperature, pH, specific conductance, and alkalinity) and nutrients (total phosphorus, ammonia, and nitrate). Specific conductance is hereafter referred to as conductivity. Basic water quality parameters were measured using in situ water quality probes: dissolved oxygen (Apera Instruments AI480 DO850 probe), pH, temperature, and conductivity (Apera Instruments SX823-B multiprobe). For nutrient concentrations (total phosphorus, total ammonia nitrogen, and nitrate), we collected 500 mL water samples from each site, immediately transported them to the laboratory at University of Arizona, and stored at 4 °C until analyzed using a YSI EcoSense Model 9500 Photometer (within 48 h). Total phosphorus was analyzed using phosphate LR method, then the phosphate was converted into total phosphorus; the ammonia test is based on an indophenol method; and finally the nitrate was analyzed using a nitratest method (YSI Inc., 2010). The intensity of the color produced in the test is proportional to those parameter concentrations and was quantitatively optically measured. The detection limits of total phosphorus, ammonia, and nitrate were 0.001, 0.01, and 0.2 mg/L, respectively. When nutrient analysis results exceeded the test’s calibration range, samples were diluted and rerun, and the appropriate dilution factor correction was applied.

Data analysis

Each of the eight water quality parameters were modeled with linear mixed effects models (LMEs: Pinheiro & Bates 2000) using distance from effluent outfall (Distance), Season, and the interaction between distance and season (Distance*Season) as predictors, with reach location (Site) as a random factor. Models were run in the statistical program R (version 3.5.1: R Core Team, 2019) with the package ‘glmmTMB’ (version 1.0.2.1; Brooks et al., 2017). Histogram analysis indicated that there were no outliers in the dataset and all variables were approximately normally distributed, except for ammonia, which was log-transformed prior to analysis (data not shown). Finally, to visualize how overall water quality varied across the year at all sites, we included all eight water quality parameters in a principal component analysis (PCA). The resulting plot allowed us to identify clustering of water quality conditions across river basins, reaches, and seasons. PCA was run with PC-ORD version 7.08 (McCune & Mefford, 2018).

Distance effect

The effect of longitudinal distance from effluent outfalls on water quality varied widely among measured parameters. Effluent is generally warmer than natural baseflow temperatures, so heat exchange between the atmosphere and stream likely led to longitudinal cooling trends (Boyle & Fraleigh, 2003), especially during winter (Beganskas & Toran, 2021; also see “4.2 Seasonal effects” below). As temperatures cooled, dissolved oxygen levels rose. Increased oxygen levels were likely due to a combination of cooler temperatures, photosynthesis by macrophytes and algae, and atmospheric reaeration (O’Connor, 1967; Correa-González et al., 2014).

In contrast, longitudinal variation in pH was not significant and values remained circumneutral to slightly basic. These results are in line with other studies of streams receiving effluent, which reported relatively stable pH values downstream of effluent outfalls (Chen et al., 2009; Prat et al., 2013; Matamoros & Rodríguez, 2017). Relatedly, we also found no longitudinal patterns in alkalinity, which indicates the buffering capacity to neutralize acids and bases of a water body. Some previous studies have reported a decrease in alkalinity with distance downstream (Birge et al., 1989; Boyle & Fraleigh, 2003), a result likely due to dilution by other water sources. However, our effluent-dependent study reaches are not fed by any perennial tributaries that can dilute the effluent during baseflow.

Ammonia levels decreased significantly with increasing distances from effluent outfalls. High levels of ammonia have commonly been reported near effluent outfalls, with decreasing levels downstream (Sebenik, Cluff & DeCook, 1972; Hamdhani, Eppehimer & Bogan, 2020). Although trends across reaches were equivocal, in the long Tres Rios reach of the Santa Cruz River, the average rate of ammonia loss over ~31 km in the fall season was about 10 times greater than nitrate gain (−0.11 mg/L/km of N ammonia; 0.01 mg/L/km of N nitrate). In summer, the average rate of ammonia loss was double that of nitrate gain (−0.08 mg/L/km of N ammonia; 0.04 mg/L/km of N nitrate) (Fig. 5). We presume the seasonal differences in nutrient gain and loss are likely due to temperature. Thus, nitrification occurring in the stream channel (Sebenik, Cluff & DeCook, 1972) does not fully account for ammonia removal, suggesting mechanisms besides oxidation of ammonia to nitrate contributed to the observed pattern. These other mechanisms could include volatilization, adsorption by cation exchange complex, fixation by clays, and fixation by organic matter (Van Schreven, 1968; Weiler, 1979; Freney, Simpson & Denmead, 1981; Simon & Kennedy, 1987; Kowalchuk et al., 1998; Lin et al., 2019). Since pH was always <9 across all reaches and seasons (see Supplemental Material), total ammonia nitrogen was present mainly as ammonium ion (NH4⁺), thus volatilization is not expected to have contributed significantly to ammonia loss in any of the study reaches.

In contrast to our findings with ammonia, total phosphorus and nitrate did not change significantly with increasing distance from effluent outfalls. Previous studies in the other effluent-fed streams (Birge et al., 1989; Chen et al., 2009) indicated phosphorus concentration decreased with distance downstream. This reduction is likely due to biotic and abiotic nutrient uptake along the reaches (Davis & Minshall, 1999; Stutter, Demars & Langan, 2010). When we look at our longest reaches, there appears to be a trend toward decreasing total phosphorus, but a lack of response over shorter distances in other reaches may have precluded the global LME model from detecting an effect of distance.

In natural streams, conductivity tends to be low in headwater reaches and then increase in downstream reaches, especially when streams pass through urban areas (Chusov, Bondarenko & Andrianova, 2014). However, previous studies of effluent-fed streams have reported elevated conductivities near outfalls and then relatively stable values downstream (Chen et al., 2009; Prat et al., 2013; Matamoros & Rodríguez, 2017). In this study, we did not observe a significant effect of distance from effluent outfall on conductivity; instead values tended to be stable. However, patterns in one reach were noticeably different. The Salt River reach had both higher baseline levels of conductivity and also a significant increase in values moving downstream from the effluent outfall (Fig. 6). These patterns could be due to urban inputs (i.e., non-point sources) that were not mapped as part of our research, or due to underlying soil conditions. For example, the Salt River Valley has many areas that are high in soluble mineral salt (Harper, 1931).

Seasonal effects

Seasonal atmospheric and climatic factors, such as air temperature, solar radiation, relative humidity, cloud cover, and wind speed, work together to shape water temperatures through heat exchange between the atmosphere and the stream (Sinokrot & Stefan, 1994; Beganskas & Toran, 2021). This effect is usually stronger in unregulated stream systems (Bowles, Fread & Greeney, 1977) and relatively weak in effluent-dependent stream systems, especially near effluent outfalls (Hamdhani, Eppehimer & Bogan, 2020). However, we did observe a significant effect of season on water temperatures across all six study reaches, with warmer water observed in spring and summer than in fall and winter, regardless of reach length (Table 2; Fig. 3). This pattern demonstrates the strong influence of atmospheric and climatic conditions on water temperature despite the stable, warm temperature of effluent at its point of release (Beganskas & Toran, 2021). Trends in dissolved oxygen levels across seasons were generally opposite to those of temperature, with the lowest values observed in summer; however, this pattern was only significant in three of the six reaches (Table 2; Fig. 4). In the other three reaches, dissolved oxygen levels might be more strongly driven by non-temperature factors, such as variation in biochemical oxygen demand, water volume, and photosynthetic activity (Guasch & Subater, 1995; Mulholland, Houser & Maloney, 2005; Mandal, Upadhyay & Hasan, 2010).

Nutrient levels varied across seasons but not in a predictable way across reaches (Table 2), which may be related to seasonal variation in raw wastewater influent and performance of each wastewater treatment plant. The most substantial seasonal variation in nutrient levels occurred on the Tres Rios reach of the Santa Cruz River, which was also our longest study reach (Fig. 5). In summer and fall, ammonia levels were much higher and nitrate levels were at their lowest, relative to concentrations observed in winter and spring, especially nearer to the effluent outfall. This finding was unexpected because conversion of ammonia to nitrate during wastewater and potable water treatment is usually highest in warmer seasons (Siripong & Rittmann, 2007; Liu et al., 2017). We suspect that manipulation in treatment operations at the Tres Rios WRF contributed to lessened nitrification capacity in summer (Hegg, Rakness & Schultz, 1979; Liu et al., 2012). Finally, phosphorus levels in effluent from the Tres Rios WRF were higher during summer and spring, and lower during winter and fall. We speculate that this may be related to the type of Accumulibacter involved in phosphorus removal processes during wastewater treatment. Some types of Accumulibacter alter their activity patterns with seasonal-induced temperature changes (Flowers, Cadkin & McMahon, 2013).

Conductivity and pH differed significantly by season across five of the six study reaches, but we could not detect any consistent patterns (Table 2; Fig. 5). For example, conductivity was often highest in the summer, but in some cases it was lowest in the summer (Fig. 5), and similar inconsistencies were observed with pH (Supplemental Materials). Both factors can be influenced by precipitation and runoff, especially in urban areas where runoff may carry pollutants (Vega et al., 1998; Laudon, Westling & Bishop, 2000). However, we avoided sampling during periods of precipitation and subsequent runoff, thus all water quality measurements reflect conditions when flow was 100% effluent. Further study is needed to understand why these factors vary inconsistently by season in the six study reaches.

Geography and climate

Geographic and climatic settings strongly influence the structural and functional features of natural rivers and streams (Gasith & Resh, 1999; Shi et al., 2019), but it is unclear how strongly effluent-dependent systems might be influenced. For example, we expected that water temperature would be fairly consistent across all reaches, regardless of climate zone, given the warm temperature of effluent resulting from the treatment process (Brooks, Riley & Taylor, 2006; Hamdhani, Eppehimer & Bogan, 2020). However, when considering all measured factors across all reaches, sampling sites, and sampling dates together, distinct clusters were observed for high (>1,000 m), medium (~700 m), and low elevation sites (<300 m) (Fig. 7). As previously mentioned, the lowest elevation reach (91^st Ave, Salt River) occurs in a stream basin with natural salt deposits (Harper, 1931). The high mineral content in the water likely contributed to the distinctness of that reach (Fig. 7A). Higher nitrate levels also were observed in the higher elevation reaches. Because rates of biochemical reactions, including denitrification, are temperature-dependent (Dawson & Murphy, 1972; Xu, Dai & Chai, 2019), lower air temperatures may reduce nitrate removal capacity of the treatment plants at higher altitudes.

In contrast to the distinct geographic clusters in the PCA ordination, we observed broad overlap of water quality conditions by season (Fig. 7B). This lack of pattern is likely due to idiosyncratic seasonal changes across the factors and reaches, as described in the previous section and observed in previous studies (Hamdhani, Eppehimer & Bogan, 2020). Thus, although individual reaches may exhibit clear seasonal trends in water quality, the variation in the significance of those trends across reaches in different locations obscures detection of any overall seasonal patterns (Gardner & McGlynn, 2009; Pratt & Chang, 2012).

Potential for effluent to support native aquatic species

The ability of effluent to enhance or re-create habitat for native aquatic species is still actively debated, but some conservation successes have been noted in recent years (Halaburka et al., 2013; Luthy et al., 2015). In our six study reaches, effluent has restored flow to rivers that had been dry for decades due to extraction of groundwater and surface water (e.g., Webb et al., 2014). However, prior to wastewater treatment plant upgrades, water quality was generally poor and aquatic biodiversity was low (Cordy et al., 2000; Walker, Goforth & Rector, 2005; Sonoran Institute, 2017). In particular, high water temperatures, low dissolved oxygen levels, and elevated concentrations of ammonia have been cited as the cause of low aquatic biodiversity in effluent-dependent streams (e.g., Monda, Galat & Finger, 1995; Hamdhani, Eppehimer & Bogan, 2020; Eppehimer et al., 2020).

Many aquatic animals are adapted to specific temperature ranges (Carveth, Widmer & Bonar, 2006; Eliason et al., 2011; Nikel et al., 2021) and warmer water can affect their growth, behavior and survival (Crawshaw, 1977; Schneider & Connors, 1982; Marine & Cech, 2004). Among 11 native Arizona fish species tested by Carveth, Widmer & Bonar (2006), the most sensitive to high water temperatures was the Speckled dace (Rhinichthy osculus), which becomes disoriented at 34 °C and perishes at 36 °C. One of the most tolerant fishes identified was the Gila topminnow (Poeciliopsis occidentalis), which can withstand temperatures as high as 38–39 °C for short periods of time (Carveth, Widmer & Bonar, 2006). The highest temperature we recorded across all of our sites and seasons was 34.2 °C, but not all instantaneous measurements were taken at the hottest times of the day for each reach, so it is possible that temperatures reach 1 °C to 3 °C higher (Lowney, 2000). These findings suggest that effluent-dependent in Arizona are generally thermally suitable for native fish, but sensitive species could be negatively impacted in the warmest reaches. In these cases, a gradient in fish composition could occur, where thermally sensitive species are only found further downstream from the outfalls. A similar pattern has been reported for the thermally-sensitive Greenside darter (Etheostoma blennioides) in a Canadian effluent-fed stream (Brown et al., 2011).

Dissolved oxygen concentration also can be a direct indicator of the ability of a waterbody to support aquatic life. We occasionally observed lower dissolved oxygen levels near effluent outfalls in summer, with concentrations as low as 2.9 mg/L, as has been observed in other effluent-fed systems (Birge et al., 1989; Boyle & Fraleigh, 2003; Matamoros & Rodríguez, 2017). Mortality or loss of equilibrium for fishes and other aquatic organisms can occur at concentrations between 1 and 3 mg/L (US EPA, 1986), and chronic exposure can cause behavioral changes that make individuals vulnerable to predation or other risk factors (Dean & Richardson, 1999). For native Arizona fishes, dissolved oxygen levels ranging from 0.22 to 1.47 mg/L have been reported as being lethal (Lowe, Hinds & Halpern, 1967). These findings suggest that oxygen levels near effluent outfalls might be stressful for native fishes during warmer seasons.

Aquatic invertebrates can also be affected by low oxygen levels, with tolerant worms (Oligochaeta: Martins, Stephan & Alves, 2008) and fly larvae (e.g., Chironomus, Chironomidae: Lencioni et al., 2008) replacing sensitive mayflies (Ephemeroptera) and stoneflies (Plecoptera) near effluent outfalls (Hamdhani, Eppehimer & Bogan, 2020). In fact, a previous study from the Santa Cruz and Salt Rivers reported low oxygen levels and only a few tolerant invertebrate taxa near effluent outfalls (Cordy et al., 2000). However, wastewater treatment plants supplying these reaches have since been upgraded, and recent findings from the Santa Cruz River demonstrate a robust aquatic invertebrate community with few apparent dissolved oxygen limitations (Eppehimer et al., 2020). Although invertebrate studies are not available for the other two rivers (Salt and Rio de Flag), dissolved oxygen measurements suggest that diverse aquatic invertebrate communities could be found there too.

Ammonia has direct toxic effects on aquatic species (e.g., Richardson, 1997; Hickey & Vickers, 1994), with concentrations greater than 2 mg/L as N (pH 7.0 and temperature 20 °C) causing impairment of aquatic life (Constable et al., 2003; Yeom et al., 2007; US EPA, 2013). The vast majority of our measurements in Arizona were below this threshold, but at one site on the Santa Cruz River we did find concentrations as high as 3-5 mg/L in summer and fall. In a prior study of the Santa Cruz River, before treatment plants were upgraded, the absence of sensitive mayfly taxa was likely a result of elevated ammonia concentrations (Boyle & Fraleigh, 2003). However, recent work has shown that both mayflies and fish are now found near effluent outfalls that experience episodic, but not chronic, high ammonia levels (e.g., lower Santa Cruz River: Sonoran Institute, 2017; Eppehimer et al., 2020). Together, these findings suggest that ammonia is no longer a primary concern for aquatic species in effluent-dependent streams of Arizona.

We suggest that in order to improve water quality for native aquatic species, wastewater treatment plants should only release tertiary treated effluent. Furthermore, designing minimum discharge goals to maintain flow for at least several kilometers is recommended (Eppehimer et al., 2020). These longer lengths of flow allow additional natural remediation of tertiary treated effluent to occur, as shown in our improvements to dissolved oxygen and ammonia, for example. Low dissolved oxygen near effluent outfalls can also be improved through installation of physical structures that promote aeration of the water (Alp & Melching, 2011). Additionally, stable flow with significant volume will allow riparian gallery forests to develop, the shade from which helps moderate water temperatures and increase dissolved oxygen levels during summer months (Ghermandi et al., 2009; Dugdale et al., 2018). Finally, diurnal flow fluctuations are common in effluent dependent river systems, which can be problematic for aquatic biota (Eppehimer et al., 2021). Retention or detention structures can be engineered to help slow the flow of the water and reduce the magnitude of diurnal effluent fluctuations in recipient streams.