EFFECTS OF FLOW REDUCTION ON THERMAL DYNAMICS OF STREAMS: IMPROVING AN IMPORTANT LINK IN MICHIGAN’S WATER WITHDRAWAL ASSESSMENT TOOL (WWAT) By Ryan Andrews A THESIS Submitted to Michigan State University in partial fulfillment of the requirements for the degree of Fisheries and Wildlife—Master of Science 2018 ABSTRACT EFFECTS OF FLOW REDUCTION ON THERMAL DYNAMICS OF STREAMS: IMPROVING AN IMPORTANT LINK IN MICHIGAN’S WATER WITHDRAWAL ASSESSMENT TOOL (WWAT) By Ryan Andrews The response of fish to human alterations of habitat conditions is of critical management and policy importance. For example, withdrawal of groundwater from stream ecosystems can result in altered thermal regimes, and changes in fish populations. Many streams are fed groundwater inputs that help maintain in-stream hydraulic conditions by stabilizing discharge as well as stream temperature. However, groundwater withdrawal through high-capacity wells is also important to the agricultural industry. Withdrawal can cause reductions in streamflow and typically results in increased stream temperature, and can initiate shifts in aquatic community assemblage. In Michigan, the Water Withdrawal Assessment Tool (WWAT) is used to estimate ecological impacts associated with water withdrawal. The model uses flow-fish response relationships to estimate the effects of flow reduction on downstream warming for Michigan streams, and subsequently estimates the degree of impact of on stream fish communities. The chapters of this thesis investigate potential improvements to the mechanisms of Michigan’s WWAT used to set policy regarding water withdrawal throughout the state. In the first chapter, I investigate the use of several benchmark detection methods for setting thermal benchmarks for stream fishes. The second chapter includes the development and analysis of several models which attempt to predict downstream warming rates within 24 streams located throughout Michigan’s Upper and Lower Peninsulas. ACKNOWLEDGEMENTS Special thanks go to the Michigan Department of Natural Resources and the Dr. Howard A. Tanner Fisheries Excellence Fellowship who provided funding for this project. I would first like to thank my advisor, Dr. Daniel Hayes for providing guidance and lending his knowledge throughout this project. He provided me with my first opportunity as a graduate student after countless attempts, and for that I am incredibly grateful. He has been greatly influential in my development as a professional, and has prepared me for success in any event. I would also like to thank Dr. Christine Mayer (University of Toledo) and Dr. Mark DuFour for giving me my first position as a research assistant, and providing the springboard I needed to advance my career. I would like to thank Troy Zorn and Dr. David Lusch for their assistance and insight with interpreting results, and for being excellent committee members. I would also like to thank Andrew Pawloski, Katie Kierczynski, and Phillip Ankley for their assistance in the field and long hours in the car. I am also thankful for my friends and family for supporting me as I pursued my degree at MSU. Finally, I am very appreciative for my wife, Megan, who helped support me with the highest quality life a graduate student could ask for. iii TABLE OF CONTENTS LIST OF TABLES ......................................................................................................................... vi LIST OF FIGURES ........................................................................................................................ x KEY TO SYMBOLS .................................................................................................................... xii CHAPTER 1: Application of benchmark detection methods to identify thermal thresholds of stream fishes along a thermal gradient............................................................................................ 1 Introduction ................................................................................................................................. 1 Methods ....................................................................................................................................... 7 Data collection ......................................................................................................................... 7 WWAT model development and threshold determination ...................................................... 8 Threshold Indicator Taxa Analysis (TITAN) .......................................................................... 9 Classification and regression tree .......................................................................................... 10 LOESS regression.................................................................................................................. 11 Stratified dataset simulation .................................................................................................. 11 Results ....................................................................................................................................... 12 Regional comparison ............................................................................................................. 15 Regional subset simulation .................................................................................................... 17 Discussion ................................................................................................................................. 17 Overall comparison of detected benchmarks ......................................................................... 17 Detection of region-specific benchmarks .............................................................................. 22 Policy implications ................................................................................................................ 24 Conclusion ................................................................................................................................. 27 CHAPTER 2: Quantifying downstream warming rates of Michigan streams and assessing effects of baseflow reduction on thermal dynamics ................................................................................. 30 Introduction ............................................................................................................................... 30 Methods ..................................................................................................................................... 35 Study sites .............................................................................................................................. 35 Data collection ....................................................................................................................... 36 Model hypotheses .................................................................................................................. 37 Processes influencing stream warming.................................................................................. 38 Latent heat flux ...................................................................................................................... 39 Sensible heat flux................................................................................................................... 40 Shortwave radiation ............................................................................................................... 40 Longwave radiation ............................................................................................................... 40 Computing heat energy of various flows ............................................................................... 41 Statistical models ................................................................................................................... 43 Air temperature – water temperature differential (Model 1) ................................................. 43 Discharge ratio (Model 2)...................................................................................................... 43 Flow temperature change and upstream discharge (Model 3) ............................................... 44 iv Downstream – upstream discharge differential (Model 4) .................................................... 44 Day length (Model 5) ............................................................................................................ 45 Altitude angle (Model 6): ...................................................................................................... 45 Including day length and sun altitude angle (Model 7) ......................................................... 46 Differential effects of flow sources (Models 8, 9, and 10) .................................................... 46 Incorporating components of surface heat exchange into a statistical model (Model 11) .... 47 Forcing model (Model 12) ..................................................................................................... 48 Model fitting, selection, and prediction ................................................................................. 48 Partial influence ..................................................................................................................... 50 Model performance................................................................................................................ 50 Pooling data ........................................................................................................................... 50 Model comparison with WWAT ........................................................................................... 51 Results ....................................................................................................................................... 52 Site information ..................................................................................................................... 52 Observed flux rates ................................................................................................................ 53 Model selection...................................................................................................................... 54 Seasonal patterns ................................................................................................................... 57 Peak events ............................................................................................................................ 58 Quantifying influence on downstream temperature flux ....................................................... 59 Comparison of pooled datasets .............................................................................................. 59 Baseflow reduction scenarios ................................................................................................ 60 Comparison with the Zorn et al. (2008) downstream warming rate model. .......................... 62 Discussion ................................................................................................................................. 65 Model comparison ................................................................................................................. 66 Model performance................................................................................................................ 67 Factors influencing flux rates ................................................................................................ 70 Pooling data ........................................................................................................................... 74 Effects of baseflow reduction ................................................................................................ 75 Comparisons with the Zorn et al. (2008) model .................................................................... 77 Implications and recommendations of findings..................................................................... 78 Conclusion ................................................................................................................................. 81 APPENDICES .............................................................................................................................. 83 APPENDIX 1.0: Supporting tables and figures for Chapter 1: Application of benchmark detection methods to identify thermal thresholds of stream fishes along a thermal gradient ... 84 APPENDIX 2.0: Supporting tables and figures for Chapter 2: Quantifying downstream warming rates of Michigan streams: Improving an important link in Michigan’s Water Withdrawal Assessment Tool.................................................................................................. 103 APPENDIX 2.1: Derivation of net incoming shortwave radiation (Meeus, 1991)................. 147 REFERENCES ........................................................................................................................... 149 v LIST OF TABLES Table 1.1: Table of temperature thresholds (°F) of species with upper thermal thresholds identified by each of the three analytical threshold detection methods, as well as visually estimated threshold interpreted from LOESS line. ....................................................................... 84 Table 1.2: Table of differences between detected thresholds of each method. ............................ 85 Table 1.3: Information of the upper 20% of abundances of each decreaser species including, number of optimal sites (N), WWAT optimum temperature (mean), minimum (Lower) and maximum (Upper) July mean water temperatures for optimum temperature calculation, standard deviation, and predicted WWAT threshold (temperatures in °F). ................................................ 86 Table 1.4: Detected benchmarks of all 49 species of each of the three datasets (ALL, NLPUP, and SLP). NA indicates that no benchmark was estimated due to lack of minimum data requirements. ................................................................................................................................. 88 Table 1.5: Table of detected benchmarks of species which were classified as decreasers in both the NLPUP and SLP subregions. Average differences were calculated between the two subregions for each method. ......................................................................................................... 90 Table 1.6: Purity, reliability, and directionality of threshold response for those species identified as decreasers in the ALL dataset. P-values are significant at P  0.05. Purity (Pur) is the proportion of responses detected by TITAN among bootstrap replicates that agree with the observed response. Reliability (Rel) is an estimate using the proportion of responses using bootstrap replicates whose IndVal scores that correspond to a probability level of P  0.05. Directionality of response indicates whether species show a negative (decreaser; -) or positive (increaser; +) threshold response along the thermal gradient. ...................................................... 91 Table 1.7: Table of means of detected thresholds using logistic function to simulate abundances (10 simulations) with a known threshold at 68 °F. ....................................................................... 92 Table 1.8: Two-way ANOVA (α = 0.05) examining differences between region and method of detected benchmarks using logistic function to simulate abundances with a known threshold at 68 °F. ............................................................................................................................................. 93 Table 2.1: Site information specific to each stream segment. Streams marked with an asterisk are those for which data was collected over both the 2015 and 2016 field seasons. Up and down refer to upstream and downstream locations. ...................................................................................... 103 Table 2.2: Model numbers and the parameters included within each model denoted with an X. Ta = air temperature (°C); Tw = water temperature (°C); Qup = upstream discharge (cms); Qdown = downstream discharge; S = day length (hr); ΔTflow = cumulative heat energy (°C); ΔTup = upstream heat energy (°C); ΔTbase = baseflow heat energy (°C); ΔTover = overland flow heat energy (°C); α = sun altitude angle (°) ........................................................................................ 105 vi Table 2.3: Comparison of a priori thermal classes and a posteriori thermal classes based upon July mean water temperatures. For those streams with two years of data (denoted by *), July mean water temperatures and overall discharges were averaged across both years. C =  17.5 °C; CT = > 17.5 °C and  19.5 °C; WT = > 19.5 °C and  21.0 °C; W = > 21.0 °C. ...................... 106 Table 2.4: Goodness of fit components for individual stream reaches using the forcing model (Model 12). ................................................................................................................................. 107 Table 2.5: Models and their goodness of fit components. Results are averages over each stream reach. K = number of parameters; SSE = sums of squared errors; r = correlation between predicted and observed temperature change; L = log likelihood component of AIC; AIC = Akaike Information Criteria; ωi = Akaike weight. Count is the total number of streams for which each model was identified as providing the best fit. Models with the smallest SSE, L, or AIC are best fitting. .................................................................................................................................. 108 Table 2.6: Model weights for each stream and the total count for which each model provided the best fit.......................................................................................................................................... 109 Table 2.7: Models and average r for streams with one and two years of data. ........................... 110 Table 2.8: Overall root-mean-square error (RMSE) (°C) between observed and predicted values of downstream temperature flux for each of the four best models (M8, M9, M10, and M11) across each of the 21 streams over 2015 and 2016. .................................................................... 111 Table 2.9: Compared performances of the four best models using root-mean-square error (RMSE) (°C) between observed and predicted values of downstream temperature flux across each of the 21 streams. ................................................................................................................ 112 Table 2.10: Overall root-mean-square error (RMSE) (°C) between observed and predicted values of downstream temperature flux for each of the four best models (M8, M9, M10 and M11) across each of the four thermal classes. ................................................................................................. 113 Table 2.11: Maximum and minimum values of downstream temperature flux rates (°C/km) for observed and predicted values of each of the four best models. ................................................. 114 Table 2.12: Residuals of observed downstream temperature fluxes and model predictions of the 1st and 99th percentiles. ................................................................................................................ 115 Table 2.13: Partial R2 values of each variable included in Model 10. Values reflect the influence of each variable on downstream temperature flux rates. Ta = air temperature (°C); Tw = water temperature (°C); Qup = upstream discharge (cms); Qdown = downstream discharge; S = day length (hr); ΔTup = upstream heat energy (°C); ΔTbase = baseflow heat energy (°C); ΔTover = overland flow heat energy (°C); α = sun altitude angle (°). ........................................................ 116 vii Table 2.14: Partial R2 values averaged over thermal class. Ta = air temperature (°C); Tw = water temperature (°C); Qup = upstream discharge (cms); Qdown = downstream discharge; S = day length (hr); ΔTup = upstream heat energy (°C); ΔTbase = baseflow heat energy (°C); ΔTover = overland flow heat energy (°C); α = sun altitude angle (°). ........................................................ 116 Table 2.15: Type III sums of squared errors for a model testing only for thermal class as a fixed effect on downstream temperature flux rates using GLM. ......................................................... 117 Table 2.16: Type III sums of squared errors for a model testing for both thermal class as a fixed effect and individual stream reach as a random effect on downstream temperature flux rates using GLM. ................................................................................................................................. 118 Table 2.17: Variance accounted for by each model from Tables 2.15 and 2.16. ....................... 119 Table 2.18: Diagnostics comparing observed values of downstream temperature flux (°C/km) with predicted values using Model 10 optimized specifically across thermal classes and across the stream network as a whole. Observed and Predicted refer to mean values of downstream temperature flux across the entire study period for each stream reach. ...................................... 120 Table 2.19: Diagnostics comparing observed values of downstream temperature flux (°C/km) with predicted values using Model 10 optimized specifically across thermal classes and across the stream network as a whole. Results are averaged among streams within thermal classes. .. 121 Table 2.20: Average monthly downstream temperature flux rates (°C/km) for baseflow reduction scenarios using Model 10............................................................................................................ 122 Table 2.21: Average monthly downstream temperature flux rates (°C/km) of thermal classes for baseflow reduction scenarios using Model 10. ........................................................................... 123 Table 2.22: August mean values used for predicting downstream temperature flux using Model 10. Ta = air temperature (°C); Tw = water temperature (°C); Qup = upstream discharge (cms); Qdown = downstream discharge; S = day length (hr); ΔTup = upstream heat energy (°C); ΔTbase = baseflow heat energy (°C); ΔTover = overland flow heat energy (°C); α = sun altitude angle (°). ..................................................................................................................................................... 124 Table 2.23: August mean values used for predicting downstream temperature flux using the Zorn et al. (2008) model. ..................................................................................................................... 125 Table 2.24: Observed August downstream temperature flux rates (°C/km) for each of the stream reaches and the predicted rates of Model 10 and the Zorn et al. (2008) model under a 0% reduction scenario. ...................................................................................................................... 126 Table 2.25: Average August downstream temperature flux rates (°C/km) experienced by each of the four stream thermal classifications and the estimated rates using Model 10 and the Zorn et al. (2008) model (King Creek and Squaw Creek excluded from Zorn et al. (2008) model estimates). ..................................................................................................................................................... 127 viii Table 2.26: Predicted downstream temperature flux rates (°C/km) following scenario analysis of baseflow reductions of 10% increments using Model 10. .......................................................... 128 Table 2.27: Predicted downstream temperature flux rates (°C/km) following scenario analysis of baseflow reductions of 10% increments using the Zorn et al. (2008) model. ............................ 129 Table 2.28: Parameterized beta coefficients associated with each stream reach to estimate downstream temperature flux using Model 10. Ta = air temperature (°C); Tw = water temperature (°C); Qup = upstream discharge (cms); Qdown = downstream discharge; S = day length (hr); ΔTup = upstream heat energy (°C); ΔTbase = baseflow heat energy (°C); ΔTover = overland flow heat energy (°C); α = sun altitude angle (°). ....................................................................................... 130 Table 2.29: Parameterized beta coefficients specific to each thermal class and across all streams to estimate downstream temperature flux using Model 10. Ta = air temperature (°C); Tw = water temperature (°C); Qup = upstream discharge (cms); Qdown = downstream discharge; S = day length (hr); ΔTup = upstream heat energy (°C); ΔTbase = baseflow heat energy (°C); ΔTover = overland flow heat energy (°C); α = sun altitude angle (°). ...................................................... 1312 Table 2.30: Partial R2 values of each variable included in Model 10 for summer months (July and August). Values reflect the influence of each variable on downstream temperature flux rates. Ta = air temperature (°C); Tw = water temperature (°C); Qup = upstream discharge (cms); Qdown = downstream discharge; S = day length (hr); ΔTup = upstream heat energy (°C); ΔTbase = baseflow heat energy (°C); ΔTover = overland flow heat energy (°C); α = sun altitude angle (°)............... 133 Table 2.31: Partial R2 values averaged over thermal class for summer months (July and August). Ta = air temperature (°C); Tw = water temperature (°C); Qup = upstream discharge (cms); Qdown = downstream discharge; S = day length (hr); ΔTup = upstream heat energy (°C); ΔTbase = baseflow heat energy (°C); ΔTover = overland flow heat energy (°C); α = sun altitude angle (°)............... 133 ix LIST OF FIGURES Figure 1.1: LOESS regression of logistically simulated data with a known threshold of 68 °F. . 94 Figure 1.2: Distribution of Brook Trout abundances along a July mean water temperature thermal gradient. Vertical lines correspond to identified benchmarks of TITAN, CART, and WWAT benchmark detection methods. The LOESS regression is depicted by the solid red regression line and is used in visual identification of inflection points. ............................................................... 95 Figure 1.3: Linear regression of TITAN and CART identified thresholds. Dashed line represents 1:1 relationship. Axis units are July mean water temperature (°F)............................................... 96 Figure 1.4: Linear regression of TITAN and WWAT identified thresholds. Dashed line represents 1:1 relationship. Axis units are July mean water temperature (°F). ............................ 97 Figure 1.5: Linear regression of CART and WWAT identified thresholds. Dashed line represents 1:1 relationship. Axis units are July mean water temperature (°F)............................................... 98 Figure 1.6: Linear regression of TITAN identified thresholds for NLPUP and SLP regions. Dashed line represents 1:1 relationship. Axis units are July mean water temperature (°F). ........ 99 Figure 1.7: Linear regression of CART identified thresholds for NLPUP and SLP regions. Dashed line represents 1:1 relationship. Axis units are July mean water temperature (°F). ...... 100 Figure 1.8: Linear regression of WWAT identified thresholds for NLPUP and SLP regions. Dashed line represents 1:1 relationship. Axis units are July mean water temperature (°F). ...... 101 Figure 1.9: Box plot showing results of threshold detection methods applied to the July mean temperature gradients of the NLPUP and SLP regions with abundance data simulated using a logistic function. ......................................................................................................................... 102 Figure 2.1: Stream gauge locations..............................................................................................134 Figure 2.2: Downstream temperature flux rates (°C/km) compared to upstream discharge (m3/s) for individual stream reaches among each of the four thermal classes. ...................................... 135 Figure 2.3: Downstream temperature flux rates (°C/km) compared to upstream discharge (m3/s) for all stream reaches within of each of the four thermal classes combined. ............................. 136 Figure 2.4: Hourly downstream temperature flux (°C) (Tdown – Tup) represented with individual points and overlain by a LOESS regression to provide a smoothed representation of downstream temperature flux over time. ......................................................................................................... 137 Figure 2.5: Model 10 predictions of downstream temperature flux rates (°C/km) compared with LOESS smoothed observed downstream temperature change for four different streams, one in x each of the four thermal classes: Cedar River (C), Fish Creek (CT), Honeyoey Creek (WT), North Branch Thunder Bay River (W). ...................................................................................... 138 Figure 2.6: Model 11 predictions of downstream temperature flux rates (°C/km) compared with LOESS smoothed observed downstream temperature change for four different streams, one in each of the four thermal classes: Cedar River (C), Fish Creek (CT), Honeyoey Creek (WT), North Branch Thunder Bay River (W). ...................................................................................... 139 Figure 2.7: Model 9 predictions of downstream temperature flux rates (°C/km) compared with LOESS smoothed observed downstream temperature change for four different streams, one in each of the four thermal classes: Cedar River (C), Fish Creek (CT), Honeyoey Creek (WT), North Branch Thunder Bay River (W). ...................................................................................... 140 Figure 2.8: Model 8 predictions of downstream temperature flux rates (°C/km) compared with LOESS smoothed observed downstream temperature change for four different streams, one in each of the four thermal classes: Cedar River (C), Fish Creek (CT), Honeyoey Creek (WT), North Branch Thunder Bay River (W). ...................................................................................... 141 Figure 2.9: Seasonal structures of error (Obs – Sim) of downstream temperature flux represented using LOESS regression for each of the four best models identified using model selection criteria. One stream from each of the four thermal classifications is represented: Cedar River (C), Fish Creek (CT), Honeyoey Creek (WT), North Branch Thunder Bay River (W). ................... 142 Figure 2.10: Comparison of model fit using a LOESS regression. Using nonlinear optimization, Model 10 was parameterized using a composite dataset of the entire stream network, and also specific to each thermal class. Model fit was then compared to observed values. ..................... 143 Figure 2.11: Baseflow reduction scenarios for example streams. Scenarios of 0%, 50%, and 90% baseflow reduction are plotted using LOESS regressions, along with observed values. LOESS regressions track seasonal and yearly trends of downstream temperature flux. ......................... 144 Figure 2.12: Downstream temperature flux rates in response to baseflow reductions. .............. 145 Figure 2.13: Downstream temperature flux rates in response to baseflow reductions. Rates were omitted for King Creek and Squaw Creek to maintain an appropriate scale on the y-axis. ....... 146 xi KEY TO SYMBOLS α altitude angle (°) Δ declination angle of the sun (°) ΔT Downstream temperature change (°C) ΔTbase baseflow heat energy (°C) ΔTdown upstream water temperature (°C) ΔTflow cumulative heat energy (°C) ΔTover overland flow heat energy (°C) ΔTr residual temperature change (°C) ΔTup upstream water temperature (°C) Δz elevation change (m) εa effective atmospheric emissivity εclr clear-sky emissivity σ Stefan-Boltzmann constant (5.67 x 10-8 W m-2 K-4) τ residence time (s) A surface area (m2) ARI Adverse Resource Impacts CART Classification and Regression Trees c Specific heat capacity of water (J kg-1 K-1) ccloud cloud cover (%) d stream depth (m) ea air vapor pressure (mb) xii es surface vapor pressure (mb) g gravitational acceleration constant (m s-2) G gauge reading (inches) H sensible heat (W m-2) IndVal Indicator Value JMT July Mean Temperature L reach length (m) Lclr clear-sky longwave radiation LE latent heat of evaporation (W m-2) LWnet net longwave radiation (W m-2) NLP Northern Lower Peninsula NLPUP Northern Lower Peninsula-Upper Peninsula p density of water (kg m-3) P atmospheric pressure (mb) q Heat flux (W m-2) Q discharge (m3 s-1) Qbase baseflow discharge (m3 s-1) Qdown downstream discharge (m3 s-1) Qup upstream discharge (m3 s-1) s wind velocity (m s-1) S day length (hr) sf shading factor SLP Southern Lower Peninsula xiii SWin incoming solar radiation (W m-2) SWnet net shortwave radiation (W m-2) TITAN Threshold Indicator Taxa ANalysis T temperature (°C) Ta air temperature (°C) moving air temperature average (°C) Tbase baseflow temperature (°C) Tup upstream water temperature (°C) Tw water temperature (°C) UP Upper Peninsula v volume (m3) V velocity (m s-1) w stream width (m) WWAT Water Withdrawal Assessment Tool xiv Application of benchmark detection methods to identify thermal thresholds of stream fishes CHAPTER 1 along a thermal gradient Introduction Stream thermal regimes are critical determinants of fish migration, growth, and survival. Stream temperature has been found to be a limiting factor in fish distribution and production (Wehrly et al., 2003), and for highly-valued stream fishes, such as trout, high water temperatures in summer can limit the ability to persist and thrive (Wehrly et al., 2007). Excessive water temperature limits migration, health, and performance of salmonids and can create competitive disadvantages (Mantua et al., 2010). The thermal regime of a stream also influences processes underlying lower trophic levels in aquatic ecosystems, and is subject to influences of local environmental characteristics such as riparian shading, groundwater input, and streamflow. Consequently, warming waters due to changing environmental conditions can create conditions unsuitable for coldwater stream fishes. Characterizing the response of stream fish communities to stream thermal characteristics is valuable for both watershed and fisheries resource management. Distinguishing patterns of population or community response to environmental conditions can help to establish action points or benchmarks associated with significant changes in species abundance or community composition. Ecological thresholds are used to describe transitions between alternative stable states once breakpoints along an environmental gradient are breached (Andersen et al., 2009). Alternative states can represent shifts in community assemblage or crossing of an extinction threshold where effects of a changing environmental gradient reduce reproduction or survival rates beyond a population’s capacity (Huggett, 2005). Given the current state of anthropogenic 1 stressors, threshold analysis has value in its ability to quantitatively assess species requirements and identify targets for conservation efforts (Groffman et al., 2006), and leverage this information as an early detection system to prevent crossing thresholds. There remains a need to define specific ecological benchmarks to aid in setting regulations for resource management. Accurate detection methods are not only necessary to avoid exceeding thresholds, but also to reduce scenarios that unnecessarily impose limitations when actual threshold responses may not exist. A particular challenge for water resources and fisheries managers in defining action points is the difficulty in maintaining a proper balance between natural streamflow regimes and artificial water withdrawals. Natural flow regimes are important for maintaining in-stream hydraulic conditions that support native aquatic communities, and disruption of flow variations can have immediate and long-term consequences for aquatic organisms. When flows are reduced, stream temperatures may more quickly equilibrate with ambient air temperatures. Decreased flows cause reductions in water velocity, water depth, and wetted width of the river, all of which can increase residence time of water and exposure to solar radiation which subsequently raises stream temperature. Changes in stream conditions associated with flow reduction in streams may result in loss of habitat volume, reduced connectivity, and water quality impairment (Labbe and Fausch, 2000; Richter et al., 2003). Under conditions of extreme flow restriction, summer stream warming rates can increase to the extent that they will no longer support coldwater fishes. Artificial water withdrawal can exacerbate flow reduction by removing groundwater that helps maintain the natural flow and thermal characteristics that coldwater species require. Water withdrawal can also create less favorable conditions for fluvial specialists, and in turn lead to conditions more favorable for macrohabitat generalists (Kanno and Vokoun, 2 2010). As such, altering natural streamflow processes can change thermal regimes resulting in a shift in aquatic community composition. A major challenge for managers is the regulation of withdrawals to minimize impacts on stream quality while allowing water usage for important human needs. As concerns over the impacts of climate and land use change on aquatic environments grows, it is important to focus efforts towards conserving cold- and coolwater species most vulnerable to temperature increases. Previous research has led to classifying species into thermal guilds based on temperature requirements (i.e., July mean water temperature) for efficient management practices (Lyons et al., 2009; Wehrly et al., 2003) as needs may vary between thermal regimes. There is notable concern over the management of coolwater streams because the response of fish assemblages occupying coolwater streams differs from that of coldwater or warmwater streams. In addition, coolwater streams provide an overlap in thermal habitat for both cold- and warmwater species. In conjunction with thermal guild association, it is important to estimate the upper thermal limit where cold- and coolwater fish assemblages begin to dissociate. Doing so will further support the efforts of Lyons et al. (2009) in providing managers with a clearer definition of how fish assemblages relate to stream thermal regimes. When regulating water withdrawal, successful water resources and fisheries management must consider the needs of coldwater game fish that provide significant economic benefits, but cannot overlook the value groundwater provides to other stakeholders. Groundwater is an important resource for irrigation, supporting the agricultural industry in Michigan, and is a major source of drinking water for the state of Michigan (Groundwater Conservation Advisory Council, 2006). In 2010, the United States Geological Survey (Maupin et al., 2010) estimated 693 Mgal/d (~1,072 cfs) of groundwater were withdrawn in Michigan. 3 Groundwater pumping can reduce the contributions of baseflow which would typically enter the stream by lowering the water table and consequently the hydraulic head. Recent efforts in Michigan have been aimed to manage rivers and streams to maintain urban and agricultural production, while also providing sufficient baseflow to preserve ecological flow requirements (Zorn et al, 2008). Baseflow inputs stabilize minimum stream flows and temperatures that define coldwater streams and the coldwater fish communities that occupy these waters. Hamilton and Seelbach (2011) reported that flow depletions between 2-4% in cold-transitional streams (July mean water temperatures > 63.5 °F and  67 °F) and small rivers in the Kalamazoo River Watershed met the threshold for potential ecosystem alteration. Management of coldwater fisheries has gained attention following the predicted decline in coldwater fish species and the high sensitivity of cold-transitional streams to reductions in baseflow (Zorn et al. 2008). Under anticipation of increased thermal stressors due to climate change and increased reliance upon groundwater resources, and the potential for ecosystem degradation, identification of optimal habitat conditions and benchmarks along thermal gradients is critical for fisheries biologists in order to protect aquatic communities. In order to fulfill requirements of the 2001 Annex to the Great Lakes Charter (Anonymous, 2001), which committed Great Lakes states and provinces to protection of water resources, Zorn et al. (2008) developed Michigan’s Water Withdrawal Assessment Tool (WWAT) to determine the potential for high-volume (>100,000 gal/d) water withdrawals to create adverse resource impacts to characteristic fish populations of Michigan streams. The development of Michigan’s WWAT required an assessment of species optimum for several characteristics (July mean water temperature, drainage area, and August 50% exceedance flow (or index flow) of in-stream habitat and an assessment of ecological responses to flow depletion. 4 This led to the creation of a dose-response model to relate the dynamics of fish populations with in-stream habitat conditions. The model is useful in predicting degradation patterns of individual species as well as community assemblages representative of the various thermal regimes of Michigan rivers. The model is used to set regional flow standards, but Zorn et al. (2008) recognize limitations in the determination of species optimum for the variables of interest. Since the development of the WWAT, new threshold detection methods have emerged that have been applied to numerous environmental and disturbance gradients (Baker and King, 2010; Brenden et al., 2008; Qian et al., 2003), and have been used to propose environmental regulations and legislation (Adams and Greeley, 2000; Richardson and Qian, 1999). Threshold detection methods operate on various assumptions regarding the distribution and other properties of empirical data, and are trained to identify disturbance thresholds. In an effort to identify thresholds of individual taxa along an environmental gradient, Baker and King (2010) developed Threshold Indicator Taxa ANalysis (TITAN), an analytical method to identify abrupt changes in abundance. TITAN incorporates occurrence frequency, abundance, and directionality of response to best capture the strength-of-association of a taxon to a particular location along an environmental gradient, and is particularly useful in identifying thresholds of rarer taxa. Classification and regression trees (CART) are another recent method for exploring ecological relationships by modeling relationships between response and explanatory variables (Qian and Anderson, 1999; Guisan and Zimmermann, 2000). CART uses splitting rules such that homogeneity is maximized, and enables researchers to identify the most influential predictor variable(s) on distribution (Prasad et al., 2006). Although these tools are intended to aide in estimating benchmarks to predict potentially detrimental effects associated with changes in environment, there also exists the possibility that detected benchmarks are statistical artifacts 5 rather than abrupt changes in species abundance. As an alternative, LOESS regressions are useful to represent the mean response of species abundance along an environmental gradient without the presumption of a threshold response. It is important to validate the presence of a threshold response using visual cues because threshold detection methods are designed to detect abrupt ecological changes along an environmental gradient, and in some cases a threshold may not be well defined and change is gradual (Qian et al., 2003). Additionally, the WWAT geographically stratifies data to reduce variation in describing optimal habitat conditions. Partitioning of data into subsets can reduce variation caused by intrinsic characteristics of distinct geographic regions. Individual taxa may show particular associations with distinct sub-regions that may be ignored when pooling data for analysis. However, data becomes sparser with an increasing degree of stratification, and may ignore rarer species that are often sensitive to changes and occupy a narrow range of environmental conditions. Further research is necessary to explore evidence supporting justification of using regional subsets in model development. For example, when estimating thermal niches of aquatic vertebrates Huff et al (2011) found that more local factors (ecoregion) are better predictors of aquatic assemblage than broad scale regional factors, although the appropriate scale and extent should be evaluated on a case-by-case basis. In 2006, the state of Michigan enacted a water law (2006 PA 33; Michigan Legislature 2006) to monitor large scale water withdrawals such that no adverse resource impacts (ARIs) are created. ARIs are characterized as the point beyond which a stream’s ability to support characteristic fish species becomes impaired. Fish assemblage structure and characteristic fish assemblages are predicted using habitat suitability criteria such as July mean water temperature. Fish response curves were developed to describe how fish assemblages respond to temperature 6 changes caused by incremental streamflow reductions and identify flows resulting in ARIs, but there is no clear establishment of a threshold beyond which population structure begins to deteriorate. As such, an investigation into thermal thresholds of Michigan stream fishes using WWAT and recently developed methods will improve the management of Michigan rivers. Thermal classifications of stream reaches provide a framework for understanding stream temperature as a limiting factor for stream fish community assemblages (Wehrly et al., 2003), as well as a common language for communication among managers, researchers, and stakeholders (Hudson et al., 1992). The purposes of this study were to (1) compare benchmark detection methods between WWAT, TITAN, and CART, (2) explore evidence supporting the need for regional stratification of data, and (3) discuss policy implications associated with incorporating alternative benchmark detection methods into the WWAT. A challenge for policy makers and managers is the high degree of variability in fish population data and varying responses of fish to habitat conditions. Data collection Methods Field data used to evaluate methods for threshold detection values were collected by the Michigan Department of Natural Resources Fisheries Division, the United States Forestry Service Hiawatha National Forest, and University of Michigan (Zorn et al., 2008). This large dataset consisted of samples from 1,389 fish assemblage surveys and 331 salmonid surveys at 1,119 unique stream reaches between 1980 and 2006. Most fish surveys were conducted using a backpack or towbarge electrofisher, but included rotenone surveys on large rivers. Catch data were reported as number of fish per lineal foot of stream sampled during electrofishing surveys. In order to make data comparable between survey types, fish catch data from non-single-pass 7 electrofishing (rotenone, multiple-pass depletion, and mark-recapture) surveys were corrected for sampling efficiency or effort. Zorn et al. (2008) did not explicitly identify a threshold beyond which species distribution and abundance or fish assemblages are reduced. Therefore, it was necessary to establish an appropriate threshold value based on the method for determining optimal habitat in the development of Michigan’s WWAT. Given the abundance scoring system developed by Zorn et al. (2008) where habitat conditions supporting characteristic species are within 1.5 standard deviations of species’ optimum, a value of 1.75 standard deviations from optimal habitat conditions was determined as the threshold (personal communication with Zorn). This study focused on July mean water temperature (°F) as the primary habitat variable. July mean water temperature has previously been identified as a critical determinant of stream fish distribution and production in Michigan (Wehrly et al., 2003; Zorn et al., 2004). July mean water temperature was estimated for each stream reach using a predictive model that combines both regression modeling and geostatistical kriging of water temperature data from 830 streams throughout Michigan (Zorn et al., 2008). WWAT model development and threshold determination Zorn et al. (2008) used the top 20% of sites based on catch per unit effort for each species in order to calculate optimal and upper thermal threshold values of July mean water temperature for species occurring in 50 or more sites. Therefore, threshold values were calculated for those species occurring at 50 or more sites. The mean of July mean water temperature for sites in the upper 20% of catch rates was used to define the optimal thermal condition for each species. Zorn et al. (2008) regionally stratified the dataset used to assign species-specific optima to the Southern Lower Peninsula (SLP) and Northern Lower Peninsula-Upper Peninsula 8 (NLPUP) to represent distinct ecoregions as determined by Albert et al. (1986). The resulting region-specific optima were viewed as more accurate descriptions of optimal conditions for each species by reducing the coefficient of variation of summarized habitat variables (i.e. July mean water temperature). As part of my analysis, species-specific optima were calculated for all sites (ALL) in order to determine whether regional stratification limits benchmark detection methods and whether benchmarks vary significantly compared with those determined by examination of the entire dataset. Threshold Indicator Taxa Analysis (TITAN) Baker and King (2010) developed TITAN in order to identify and detect change points in both species occurrence frequency and relative abundance. TITAN uses indicator value analysis (Dufrêne and Legendre, 1997) to create indicator value (IndVal) scores specific to each taxon. IndVal scores are used to associate taxa with a positive or negative response along an environmental gradient. IndVal scores are computed as the product of cross-group relative abundance (proportion of abundance among all sample units belonging to group i) and within- group occurrence frequency (proportion of sample units in group i with a positive abundance value). IndVal scores are computed using equation (1) as follows: (1) where Aij = Number of individualsij/Number of individualsi, and Bij = Number of sitesij/Number of sitesi. Taxa are partitioned into negative (decreaser) or positive (increaser) response groups after comparing IndVal scores above and below each candidate change point. The magnitude of IndVal scores on each side of a candidate change point determines whether a taxon shows greater association with a negative or positive response (Baker and King, 2010). 9 TITAN’s ability to estimate thresholds for rarer taxa is based on the methods of Dufrêne and Legendre (1997) which assigns higher IndVal scores to locations along an environmental gradient where species show strong presence at the majority of sites at one point along the gradient. TITAN is able to mitigate bias by group size by integrating occurrence frequency and relative abundance. For this reason, TITAN requires only a minimum of three occurrences for threshold detection. In this analysis, all taxa with less than five occurrences were removed from the dataset. TITAN quantifies uncertainty in change point identification using bootstrap resampling in order to estimate synchrony in TITAN’s ability to replicate species-specific change points. Purity and reliability are two measures of quality assurance and are measured using 500 bootstrap replicates. Using randomly distributed data, indicator purity measures the proportion of response directions (positive or negative) associated with a change point for each taxon in relation to the observed change point (i.e., high purity results in 95% of bootstrap runs resulting in same response direction). Indicator reliability provides an estimate of how substantial the probability of obtaining an equal or larger IndVal score differs when comparing the observed dataset to a randomly distributed dataset. For example, if 95% or more of bootstrap replicates result in a P-value less than a predetermined probability level of 0.05, then that indicator is reliable. Classification and regression tree Classification and regression trees (CART) are often used in ecology to describe the relationship of an outcome using a set of explanatory variables (Qian and Anderson 1999; Prasad et al., 2006; Steen et al., 2008). Both the response and explanatory variables may consist of a continuous gradient or discrete compartments. Beginning with a single node or root, 10 classification and regression trees expand downward based on splitting rules that partition explanatory variables based on the likelihood of outcome. Variables and splits are chosen such that the impurity of the outcome is minimized, and the tree will continue to grow until a terminal node is reached, or recursive partitioning creates no additional nodes. In order to avoid such errors, trees are pruned to meet tests of independence and/or cross-validation, as well as meeting a minimum error rate. Splitting rules are used to minimize impurity of each node of the regression tree (Moisen, 2008). Change points are used to specify impurities below which a node will not be split to avoid situations of overfitting. In this model, the explanatory variable used was July mean water temperature (°F). LOESS regression Loess, or locally weighted regression, fits a regression line to data through univariate smoothing (Cleveland, 1979) based on weighted averages. A LOESS regression was fit to equation (2), (2) where yi = species abundance, g = smooth function, xi = July mean water temperature, and ϵi = random errors with mean 0 and constant scale. I used a LOESS regression to obtain a smoothed representation of the data to visually identify apparent thresholds, if any exist. For example, Figure 1.1 is a graphical representation of a simulated dataset following a logistic function with a known threshold of 68 °F. The vertical line is located at the inflection point representing a true threshold response. Stratified dataset simulation The current WWAT procedure calls for stratification of fish distribution and abundance data into NLPUP and SLP sub-regions to reduce the coefficient of variation used to describe 11 species-specific optimal thermal habitat. Regional stratification can also capture variation in regional factors that can dictate fish assemblage structure (Huff et al., 2011) which is potentially overlooked when pooling data over large scales. However, stratification potentially limits the WWAT in its ability to detect benchmarks for rarer species which do not meet the required number of occurrences to estimate habitat optima and subsequent benchmarks. There also remains a possibility where region-specific benchmarks are biased due to differing thermal conditions in the different regions without any difference in the species-gradient function. I used a logistic function to simulate fish abundance with a known threshold to test for significant differences in benchmark detection between the NLPUP and SLP regions when compared with the ALL dataset. Simulated abundance across a thermal gradient was simulated using the logistic function in equation (3), (3) where L = the curve’s maximum value, e = the natural logarithm base, k = the steepness of the curve, and x0 = the sigmoid’s midpoint. For each simulation, the values of x (i.e., stream temperature) were taken from the observed data in each region. A random error term was introduced to create random abundances associated with each of the stream reaches included in the dataset. Ten datasets were simulated where a known threshold was defined as x0, each with a different set of randomly assigned error terms, and benchmark detection methods were then applied to examine the variance within and between regions. Results Data was subset representing three geographic regions of Michigan (ALL, NLPUP, and SLP) with 1119 unique stream reaches total, 576 located in the NLPUP, and 543 located in the 12 SLP, and ALL is a statewide pooled dataset. July mean water temperatures ranged from 50.9 °F to 77.2 °F across all sites. Temperatures tended to be warmer in the SLP with a median July mean water temperature of 69.3 °F, whereas the median temperature was 63.5 °F in the NLPUP. Across all sites, the median July mean water temperature was 65.3 °F. When examining the pooled ALL dataset which combines all stream reaches throughout Michigan, benchmarks were identified using WWAT, TITAN, and CART methods for 49 species in total. Of these 49 species, only the 12 species for which TITAN identified a decreasing association (z-) with the thermal gradient were examined to focus on cold or coolwater fishes as they are of primary regulatory interest (although Brook Stickleback was identified as a decreaser by TITAN, but the identified threshold (56.7 °F) was far below the WWAT optima (65.9 °F), and was therefore not included in further analysis). Table 1.1 presents the predicted thresholds for decreaser species of each of the three methods as well as the apparent visual threshold detected using the LOESS regression. Mean differences between detected thresholds of each method were calculated and are shown in Table 1.2. The mean difference between TITAN and CART benchmarks was 2.0 °F with more conservative estimates for CART, although CART predicted higher benchmarks for Rainbow Trout and Chinook Salmon. A minimum difference of 0.0 °F was calculated for Northern Redbelly Dace and a maximum difference of 8.2 F ° for Brook Trout (Table 1.2). Mean differences between TITAN and WWAT were 2.0 °F with TITAN providing more conservative estimates. A minimum difference of 0.2 °F was calculated for Blacknose Dace and a maximum difference of 4.4 F ° for Brook Trout (Table 1.2). The mean difference between WWAT and CART was 4.0 °F as CART estimated more conservative benchmarks. A minimum difference of 0.3 °F was calculated for Rainbow Trout and a maximum difference of 12.6 °F for Brook Trout (Table 1.2). The WWAT nearly always identified a higher 13 threshold when compared with TITAN and CART (both TITAN and CART identified higher thresholds for Chinook Salmon at 0.70 °F and 2.80 °F above WWAT, respectively). Thresholds were visually apparent via LOESS regression for eight of the 13 species with upper thermal thresholds (see supplemental files). Table 1.1 presents results of the identified thermal thresholds of each of the three methods along with the visually identified thresholds, where apparent. Visual thresholds were apparent for Brook Trout, Blacknose Dace, Brown Trout, Chinook Salmon, Coho Salmon, Longnose Dace, Mottled Sculpin, and Slimy Sculpin. Visual thresholds were not apparent for Burbot, Northern Brook Lamprey, and Northern Redbelly Dace (see supplemental files). The LOESS regression for Burbot showed a slight incline where CART identified a threshold, but gained little separation from the x-axis and showed a very gradual downward trend, with no clear downward inflection. Northern Brook Lamprey achieved a slight separation from the x-axis, but showed a gradual downward trend. Northern Redbelly Dace also showed a slight separation from the x-axis, and also displayed a gradual downward trend as opposed to a clear threshold response. Predicted thresholds of each method were plotted against each other with linear regressions (Figures 1.3 – 1.5). TITAN and CART (Figure 1.3) thresholds showed similarities of a near 1:1 relationship (slope = 1.08; R2 = 0.37), but CART thresholds tended to be more conservative as the regression was shifted below the 1:1 regression line towards TITAN. When compared with WWAT (Figure 1.4), TITAN thresholds tended to be more conservative at the lower range of temperatures, but thresholds began to converge at higher temperatures as the fitted regression converged with the 1:1 regression line (slope 0.50; R2 = 0.48). Thresholds of CART and WWAT (Figure 1.5) were the least comparable between each of the three methods as 14 even at the upper range of temperatures WWAT predicted thresholds were consistently higher (slope = 0.65; R2 = 0.07). Regression equations are displayed in Figures 1.3 – 1.5. Species abundances for sites were plotted along the gradient of July mean temperatures. Figure 1.2 is an example of how Brook Trout abundance varied along the thermal gradient with benchmarks representing detected thresholds for each of the three methods overlain on the plots along with the LOESS regression. Plots for all species are available in the supplemental files. Table 1.3 shows information on the optimum temperature calculated by WWAT along with the lower and upper temperatures of the top 20% of optimal sites, and the predicted thresholds of all 50 species that met the minimum criteria required by WWAT. WWAT thresholds were consistently at the upper range of the thermal gradient of observed species densities used to calculate optimum temperature. WWAT predicted thresholds for 28 species above their respective maximum temperatures used in optimum temperature calculation (i.e., WWAT predicted threshold for Coho Salmon was 1.6 °F above the upper temperature of 67.5 °F as seen in Table 1.3), whereas TITAN and CART had 5 and 1 occurrences, respectively. Regional comparison Following regional subsetting, of the 49 species for which benchmarks were detected by the three methods in the ALL dataset, detection ability was limited to 45 (90%), 46 (94%), and 32 (65%) species in the NLPUP for TITAN, CART, and WWAT, respectively. In the SLP, detection ability was limited to 45 (90%), 45 (92%), and 40 (82%) species for TITAN, CART, and WWAT, respectively. For 8 (16%), 7 (14%), and 24 (49%) species, abundance data did not meet the minimum number of sites required to detect benchmarks for one or both regions for TITAN, CART, and WWAT, respectively (Table 1.4). 15 Species-specific benchmarks were analyzed to compare benchmark detections between regions. Regional differences in benchmark estimates for each method were compared using the 13 species from the overall comparison if species were identified as decreasers by TITAN in each of the NLPUP and SLP regions. Table 1.5 contains region-specific benchmarks of the three species (Brook Trout, Brown Trout, and Mottled Sculpin) which met the criteria for inclusion in the regional comparison, while the remaining 10 species included in the overall analysis either did not show a decreasing response (z-) for both regions or benchmarks were only detected in one or no sub-regions. On average, temperature thresholds differed by 1.5 °F, 0.7 °F, and 1.1 °F for TITAN, CART, and WWAT, respectively (Table 1.5). Estimated benchmarks were higher in the SLP subregion for both TITAN (Figure 1.6) and WWAT (Figure 1.8), while CART (Figure 1.7) had higher estimated benchmarks for the NLPUP subregion. Estimated benchmarks using TITAN were compared between each of the three regional datasets and are presented in Table 1.8. Results varied between region for the number of benchmarks detected for each region, directionality of response, purity, and reliability (Table 1.6). In ALL, TITAN detected significant benchmarks for 12 species, but only ten and seven species for NLPUP and SLP, respectively. TITAN failed to detect significant benchmarks for five species in either one or both subregions (Burbot, Chinook Salmon, Northern Brook Lamprey, Northern Red Dace, and Slimy Sculpin). Each of the 12 species in ALL showed a negative response with the thermal gradient. In the NLPUP, of the ten species with significant benchmarks, four species (Brook Trout, Brown Trout, Coho Salmon, and Mottled Sculpin) showed negative responses along the thermal gradient, and the remaining six (Blacknose Dace, Burbot, Longnose Dace, Northern Brook Lamprey, Northern Red Dace, and Rainbow Trout) had a positive threshold response. In the SLP, all seven species for which TITAN identified a 16 significant threshold showed a negative threshold response (Burbot, Chinook Salmon, Northern Brook Lamprey, Northern Red Dace, and Slimy Sculpin did not have a significant threshold response). Regional subset simulation A logistic function was used to simulate abundance data with a known threshold (68 °F) in order to examine the differences in detected thresholds using TITAN, CART, and WWAT when data were regionally divided. Detected thresholds were 1.4 °F, 0.6 °F, and 1.0 °F higher in the SLP compared to the NLPUP for TITAN, CART, and WWAT, respectively (Table 1.7). Maximum thresholds detected were 70.1 °F, 68.3 °F, and 68.2 °F for TITAN, CART, and WWAT, respectively, and occurred within the SLP for each method (Figure 1.9). For each of the three datasets, TITAN identified average thresholds 1.2 °F and 1.5 °F higher than CART and WWAT, respectively (Table 1.7). Additionally, TITAN identified an average threshold response at or beyond the known threshold for each region (ALL = 69.3 °F; NLPUP = 68.0 °F; SLP = 69.4 °F). Table 1.8 consists of results of a two-way ANOVA comparing the mean response between region and method. Results indicate highly significant differences exist between detected benchmarks of each region, method, and region:method combination. Overall comparison of detected benchmarks Discussion Accurate determination of temperature benchmarks and identification of habitat available to fish is critical to the improvement of water and fisheries resource management. Using a robust fish abundance and stream temperature dataset, I was able to apply three benchmark detection methods that use different features of abundance and distribution data to detect and compare benchmarks specific to Michigan stream fish species. LOESS regressions also were useful in 17 providing visual cues to validate the presence of an estimated threshold. Applying the benchmark detection methods to the cumulative dataset allowed for comparisons over a broader range of species due to the minimum amount of occurrences required to provide unbiased benchmarks. Employing various methods allowed for comparing estimated benchmarks from an approach which assumes normally distributed (WWAT) abundances along a habitat gradient with approaches that do not (TITAN and CART). The WWAT always assigned upper thermal benchmarks as it was developed to relate the effects of baseflow reduction to increased water temperatures during summer conditions, and subsequently, the effects of increased water temperatures on fish assemblage. However, TITAN and CART identified lower thermal benchmarks for some species, and upper thermal benchmarks for other species. This is due in part to their abilities to incorporate directionality of response along the environmental gradient. Moreover, the LOESS smoother indicated that many species did not respond to temperature with a clear threshold response suggesting that each method may assign statistical change points which do not necessarily indicate the presence of an actual threshold. Estimates from Michigan’s WWAT, which is currently in use, generally were higher than estimates from TITAN and CART, which showed greater similarity than when compared with the WWAT. Benchmark estimation using the WWAT is the result of modeling fish responses to flow reduction. The WWAT assigns benchmarks as a function of species-specific optimal thermal habitat and assumptions of a normal distribution of abundance along a habitat gradient. For this reason, WWAT benchmarks can be assigned independent of an actual break point or inflection in the relationship between abundance and the thermal gradient. WWAT benchmark detection is useful in its ability to detect sub-optimal habitat because it predicts reduced abundance with changes in habitat conditions associated with water withdrawal. Thus, WWAT 18 can identify the point at which species are more likely to display abundance characteristics related to changes in habitat quality. However, this is not the same concept as an abrupt change, or threshold, in species abundance as related to the thermal gradient, which is better described by methods which distinguish linear patterns from nonlinear patterns (Ficetola and Denoël, 2009). The WWAT model is based on the assumption that water withdrawal will shift fish abundances to more closely represent that of another fish assemblage at a location with increased water temperature. The WWAT further assumes the breakpoint between statistical outliers and those within the general population is a fixed distance (mean + 1.75 std). Issues arise when anomalously high values of abundance occur at extreme temperatures relative to the species in question since the WWAT uses the thermal conditions of the top 20% of sites where abundance was the highest. These extreme values, often considered outliers, lie in the tails of the normal distribution, and theoretically would lie at or beyond the proposed threshold. However, statistical outliers typically belong to a different population because they originate from another process or source (Hampel et al., 1986). Incorporating extreme temperatures into the determination of the optimal temperature inherently biases the estimate because the statistical outliers are derived from a second (contaminating) distribution (Reimann et al., 2005) and are potentially influenced more so by other environmental variables. Currently, the WWAT operates under an assumption such that a threshold exists at the mid-point between 1.5 and 2.0 standard deviations from the optimal temperature which distinguishes thriving and/or characteristic populations (i.e. high abundance, multiple age classes, and good reproduction) from those at poorer, less suitable habitat conditions. As such, a threshold of the mean + 1.75 std was investigated for this analysis. However, as Characteristic populations are described as having an upper range of mean + 1.5 standard deviations the true 19 threshold should lie at the upper range of characteristic populations. Beyond this threshold, species should immediately begin to show declining trends in abundance similar to those at sites with sub-optimal habitat characteristics. Application of benchmark detection methodology to the statewide dataset allowed me to identify notable variations among thermal benchmarks of various stream fish. Of the 54 total species for which each method was able to detect a significant benchmark, only 13 of those species were found to have an upper thermal benchmark among all three methods. As TITAN and CART are able to distinctly identify directionality of response along an environmental gradient, it becomes difficult to apply these tools specifically for the purpose of identifying upper thermal benchmarks associated with the overall goals of Michigan’s WWAT. However, large differences in benchmarks existed among those 13 species, most notable of which was the nearly 2 °F and 4 °F average difference between WWAT and TITAN and WWAT and CART, respectively. WWAT estimates were often beyond the inflection points of the LOESS regression and sometimes occurred at or even beyond habitat conditions suitable for cold- or cool-water species. For example, WWAT identified an upper thermal benchmark of 69.3 °F for Brook Trout using the statewide dataset. Brook Trout acclimation temperatures, beyond which reduced feeding and growth of salmonids begin to occur, have typically been observed between 46 – 68 °F (Selong et al., 2001). Results are similar for Slimy Sculpin. With previously identified acclimation temperatures in the range of 41 – 68 °F (Otto and Rice, 1977), WWAT’s estimate of 68.7 °F is potentially beyond the range of preferred thermal habitat for some cold- and coolwater species. Although the goal of the WWAT is to assess the effects of increased stream temperatures on fish abundance for all stream fishes in Michigan, it is important to recognize that warmwater 20 species may not encounter temperatures near their upper thermal limit, and an upper threshold may not exist for these species in Michigan. A major advantage of TITAN and CART is the capability of assigning either a lower or upper thermal threshold. It is critical for fisheries management to identify upper thermal limits for those cold- and coolwater species which are particularly susceptible to increasing thermal stress. Lyons et al. (2009) developed water temperature criteria for coolwater streams to help differentiate them from cold- and warmwater streams and reassigned widespread species to thermal classes based on indicator species analysis. When comparing thermal guild classifications of Lyons et al. (2009) with estimated thermal thresholds, TITAN was able to detect upper thermal thresholds for all species showing significant associations with coldwater streams. TITAN also detected upper thermal thresholds for species associated with cold-transitional and warm-transitional streams. Additionally, there were no upper thermal thresholds detected for any species which Lyons et al. (2009) classified as a warmwater species. These findings highlight TITAN’s potential for identifying thresholds relevant to the goals of Michigan’s WWAT. Effective water resources and fisheries management not only relies on protecting against fish assemblage degradation, but also avoiding unnecessary restrictions on allowable water withdrawal which many stakeholders rely upon for economic purposes. In the cases of coldwater fishes such as Brook Trout and Coho Salmon, CART identified thermal benchmarks well below the optimal thermal habitat estimated by WWAT. Furthermore, LOESS regressions of these species provide visual evidence of multiple inflection points indicative of thermal thresholds, suggesting that more than one threshold or transition point can exist. Presence of multiple thresholds throughout a broad environmental gradient, such as the thermal gradient investigated in this study, can indicate the presence of multiple stable states in 21 which changes in environmental parameters influence the behavior of state variables (population abundance) (Beisner et al., 2003). Thus, population dynamics can fluctuate variably along an environmental gradient. For managers, it is critical to understand the consequences associated with crossing each threshold when multiple thresholds exist. Consideration must be given to biological processes such as growth, reproduction, and survival, each of which are susceptible to change at various points and ranges along an environmental gradient. Detection of region-specific benchmarks Stratification is potentially useful to minimize variability often attributed to local and regional factors such as competition, climate, land use, and biogeographical history (Huff et al., 2011). The field data used in my analysis was limited to only three species for which each method identified an upper thermal benchmark within both regions. Regional stratification affected each method differently in the detection of region-specific benchmarks. Subsetting the data decreased the number of sites available with which to calculate optimal habitat values according to the criteria required by the WWAT. In the case of TITAN, stratification reduced measures of purity and reliability which TITAN relies upon to ensure the quality of indicator response for each taxon. Variable directionality of response between regions was also a major factor contributing to the regional differences for each species. When considering decreaser species, stratification of the dataset resulted in reduced benchmark detection for each method. In comparison to the cumulative ALL dataset, benchmark detection was reduced by 15.4 % and 38.5% for TITAN, 0% and 30.8% for CART, and 0% and 61.6% for WWAT, in the NLPUP and SLP, respectively. Reduced detection of benchmarks by TITAN resulted from an inability to detect a significant threshold (P  0.05). Significant thresholds result when an IndVal score for a particular change point (location along the 22 environmental gradient, or in this case thermal gradient) is distinguishable as the maximum IndVal score when confidence intervals for taxon-specific change points do not overlap. Although significant thresholds are identified for some species, measures of purity and reliability can be used as diagnostic measures to measure the quality of indicator response for taxa (Baker and King, 2010). Table 1.6 contains measures of purity and reliability for each of the decreaser species among each of the three data subsets. Following regional stratification, the number of significant thresholds detected by TITAN decreased from nine species in ALL which met all three criteria of a significant P-value, purity, and reliability ( 0.05), to five in the NLPUP, and six in the SLP. Inability to detect thresholds by CART is the result of a failure to branch or split from a single node, and thus, no significant threshold is estimated. As stated previously, WWAT requires a minimum of 50 sites with abundance data in order to calculate optimal thermal habitat and, therefore, a thermal threshold. As salmonids and other species typically associated with coldwater streams such as slimy sculpin comprise a majority of the species included in our analysis, it is not surprising that the SLP region is most affected by data stratification. The SLP region (69.3 °F) had a median July mean water temperature 5.8 °F higher when compared to the NLPUP region (63.5 °F). TITAN relies upon the relative magnitude of each IndVal score on each side of potential change points along the environmental gradient to determine whether a specific taxon shows greater association with reductions or increases in the gradient variable. In our analysis, TITAN identified an upper thermal benchmark more often in the SLP, perhaps suggesting fish species in the SLP more often encounter streams near an upper thermal threshold and exhibit a response indicative of an actual threshold. Temperature is a critical factor in determining the life history of fish, but distribution is often influenced by other environmental variables (Shrode et al., 1982) and thus causes 23 intraspecific variation in thermal optima among regions, thereby complicating our ability to detect the true thermal threshold. Species-specific thermal optima and associated benchmarks are related to regional stream temperature (Christie and Smol, 1993). Huff et al. (2011) found such phenomenon in their comparison of taxon-specific thermal niches of aquatic vertebrates between basin and ecoregion where some species showed differences in thermal niches among regions where median temperatures varied in a similar way. Evidence from the simulation analysis presents similar trends when applying benchmark detection methods to thermal gradients characteristic of each region and using a known thermal threshold. Estimated thresholds for each method were highest for the SLP (median July mean water temperature = 69.3 °F) and lowest for the NLPUP (median July mean water temperature = 63.5 °F), while thresholds for ALL were consistently in between (median July mean water temperature = 65.3 °F). As the data for this analysis were simulated, the differences are purely artifacts of each method. Greater variation in temperatures of streams between regions allows for increased likelihood that regional temperature will play an important role in determining available habitat for taxa. The similarity in trends of increasing threshold temperature with median July mean water temperature are perhaps suggestive of a situation where the greater presence of warmer streams results in a shift in the estimated threshold for each method. Policy implications Our results suggest that current WWAT methods may be overestimating upper thermal benchmarks for coldwater species and assigning potentially spurious thresholds for warmwater species. A more conservative estimate of the upper limit of optimal thermal habitat would likely require further restriction on total allowable withdrawals, but there remains an issue of how to implement restrictions since the effects of flow reduction vary among each of Michigan’s four 24 thermal stream classes (Zorn et al., 2008). It is important to recognize that trout species typically associated with coldwater streams (Lyons et al., 2009) had upper thermal thresholds in our analysis near the equivalent of warm-transitional streams. This is also true of other species with upper thermal thresholds, however Blacknose Dace, Burbot, and Northern Brook Lamprey had upper thermal thresholds within warm streams, but these species either show significant association with warmer thermal guilds or have no strong association with one particular thermal guild. In the case that further restrictions are required on cold-transitional streams which are most sensitive to flow reductions, perhaps lighter regulations on warm streams will mitigate the loss in water availability to commercial and recreational water withdrawals. Although further research is necessary to understand the upper thermal limits of warmwater fish species, further refinement of thermal habitat requirements of stream fish can lead to more efficient water resources management. Further withdrawal restrictions placed on cold-transitional streams will likely handicap stakeholders who rely upon these streams because of their sensitivity to reductions in index flow. Under current legislation, just 2-4% reductions in index flow are predicted to cause adverse resource impacts to the stream ecosystem. However, estimates of the upper thermal thresholds for Brook Trout, Brown Trout, and Slimy Sculpin using TITAN or CART are within cold- transitional streams as opposed to estimates from WWAT which are in the range of warm- transitional streams. Using these more conservative estimates could restrict water withdrawal within cold-transitional basins completely, or require case-by-case investigations for all cold transitional streams to assess a stream reach’s resilience to withstanding large quantity withdrawals and avoiding adverse resource impacts. It is worth mentioning that this could also affect total allowable withdrawal in circumstances where cold transitional streams are 25 downstream of cold stream withdrawals. In any case, it is crucial to the development of water withdrawal policy implementation to distinguish the point where coldwater and warmwater species overlap to properly enforce restrictions that maximize optimal thermal habitat and avoid unnecessarily restricting withdrawals to streams able to withstand a substantial amount of water withdrawal. Another confounding factor in setting water withdrawal regulations is the issue arising when species do not exhibit a true threshold response, but rather a more gradual response to increasing stream temperatures. The magnitude of change may not be indicative of a clear threshold due to insufficient data, or threshold behavior may be the response to an interaction of more than one causal agent (Huggett, 2005), some of which are not easily controlled by management measures. Other investigations caution against implementing regulation to prevent crossing the threshold of a single species (Huggett 2005; Lindenmayer and Luck, 2005) because not all species conform to the same landscape and habitat patterns. It may be more effective to set thresholds for specific threats such as habitat deterioration or connectivity for a guild of species. The WWAT currently employs a similar strategy by constructing zones corresponding to the magnitude of impact that flow reduction will have on fish populations (Steinman et al., 2011) for each of the 11 stream size-thermal classification combinations throughout Michigan. By applying thresholds based on thermal guilds, management can avoid the use of species-based thresholds, and instead group species based on similarities in optimal habitat preferences to preserve and create specific habitat types. As many of the warmwater species in this analysis did not show evidence of a true threshold, but instead showed a more gradual decline in abundance beyond a certain point on the thermal gradient, it may be reasonable to manage warm streams without adhering to the same assumptions regarding thresholds in cold transitional streams. 26 Doing so could potentially create water availability for those stakeholders who may be affected by further restrictions protecting against adverse resource impacts in coldwater streams. The implementation of a water resources management strategy based around threshold analysis requires an understanding of species and thermal guild responses to flow reduction to most efficiently and effectively regulate water withdrawal such that ecological effects are minimized and allowable withdrawal is properly allocated to all stakeholders. Adjustments to the current impact zones used by the WWAT will require careful consideration for stakeholders who will be affected, and mitigating factors must be given thought. If possible, actions which are most cost effective and which have strong feedbacks to the rest of the system should be prioritized as long as protection goals are being met (Suding et al., 2004). For example, as stream temperature can be regulated in this case by regulating flow reduction, management can control for species whose distributions are limited more so by flow requirements. Effective water resources management requires a balance of science-based decision-making and broad stakeholder involvement to encourage buy-in and eventual implementation of new regulations aimed towards proper allocation of resources in the event of increased demand (Steinman et al., 2011). Michigan’s WWAT was created to aid in the protection of aquatic ecosystems such that Conclusion artificial water withdrawal does not cause ecological functional impairment or harm a stream’s ability to support characteristic fish populations (Steinman et al., 2011). Such changes in stream health and fish assemblage composition can be characterized by an ecological regime shift caused by breaching of ecological thresholds. As such, applying threshold detection methods to estimate important benchmarks related to fish abundance distribution along a thermal gradient 27 can be used to set action points to avoid degradation of optimal thermal habitat. By comparing WWAT’s previously determined species-specific thresholds with those estimated using methods which do not assume a normal distribution, we estimated more conservative benchmarks for fish species typically associated with coldwater streams and were able to distinguish between true threshold responses and spurious upper thermal thresholds detected for warmwater species. Current water withdrawal legislation may be overestimating the upper thermal limit of coldwater fish and assigning false upper thermal limits to many warmwater species due to assumptions that fish abundance and distribution data is normally distributed and that upper thermal limits exist for all species. Alternative threshold detection methods like TITAN and CART can distinguish between lower and upper thresholds and also provide diagnostic measures to estimate uncertainty and protect against overfitting. However, CART requires subjective interpretation of directionality of response and may also provide multiple splits or branches which can complicate the detection of a single threshold. While these methods can estimate thermal thresholds, many models will identify a single threshold independent of the nature of the relationship or the amount of thresholds which actually exist (Clements et al., 2010). LOESS regressions are a valuable tool for examining species-specific relationships between fish abundance and environmental gradients. For the purposes of determining available thermal habitat in regards to the relationship between increasing thermal stressors and fish populations, I recommend using TITAN as it provides an enhanced ability to assign more accurate ecological benchmarks that are critical to the successful management of Michigan’s natural resources. With the use of IndVal scores which incorporate aspects of distribution such as occurrence and abundance, TITAN is able to detect 28 thresholds of rarer species. TITAN also provides diagnostic measures such as purity and reliability that measure the quality of indicator response for any taxon. TITAN provides an ability to further distinguish thermal guild associations of some highly-desired coldwater species which may occupy coolwater streams that are difficult to manage due to the degree in overlap of cold- and warmwater fishes. I also suggest the use of LOESS regressions to supplement threshold detection as a visual aid to verify the existence of a true threshold response. 29 Quantifying downstream warming rates of Michigan streams and assessing effects of baseflow CHAPTER 2 reduction on thermal dynamics Introduction Water temperature plays a critical role in biotic and abiotic processes in stream ecosystems. Temperature influences organismal physiology (Fry, 1971), is a limiting factor for species distributions throughout river reaches (Caissie, 2006), and can also affect rates of biological processes (Mantua et al., 2010). Stream temperature is susceptible to change due to alterations in streamflow through both direct, surface withdrawals and groundwater pumping, potentially leading to changes in ecosystem processes and community composition. Investigation of flow reduction on the thermal dynamics of streams is necessary in order for natural resource managers to successfully regulate water uses while maintaining the ecological services of streams. Groundwater discharge sustains streamflows during low-flow periods (Kendy and Bredehoeft, 2006) and helps stabilize stream temperature during summer by providing consistent inputs of cool groundwater reserves (average groundwater temperatures vary from about 10.0o C in southern Lower Michigan to about 5.5o C in northern Lower Michigan). While groundwater serves as a critical component in maintaining aquatic habitat, it is also a valuable resource in areas with limited or fully allocated surface-water supplies. Withdrawal from wells can affect the rate of baseflow input to streams by creating a cone of depression which can lower the water table, and subsequently reduce baseflows to the stream reach, potentially leading to lower streamflow (Leake et al., 2008). Flow modification through water withdrawal has been shown to impact the spatial and temporal variability of water temperature (Sinokrot and Gulliver, 2000). 30 Sustainability of coldwater streams in Michigan is dependent upon a consistent input of baseflows from aquifers that provide the thermal habitat required by valuable coldwater fish species. Effects of water withdrawal can vary on a seasonal basis. Summertime streamflows are of particular concern to managers as reductions in flows can cause water temperatures to more quickly equilibrate with ambient air temperatures, reducing habitat quality and availability for stream fish (Zorn et al., 2002; Wehrly et al., 2006). Further, demand for water for irrigation and cooling uses is highest during summer months, and is expected to increase with climate change. As a response to the development of Instream Flow Incremental Methodology (Bovee, 1982), which examines the impacts of flow reductions on fish and habitat, Nuhfer and Baker (2004) examined the effects of reduced summer streamflows on brook trout population levels in Hunt Creek, Michigan. As the study reach was relatively short (602 m), and the coldwater stream was near the thermal optimum for brook trout, effects of downstream warming had little impact on brook trout populations. However, through additional work, Nuhfer et al. (2017) reported that maximum daily temperatures would be predicted to increase to temperatures potentially uninhabitable by many trout species throughout Midwestern streams based on the results observed in Hunt Creek. Prompted by the findings of Nuhfer and Baker (2004), Zorn et al. (2008) developed a single physical process model to estimate the impact of reduced summer baseflows (identified as 50% exceedance flow of August, typically the lowest flow month of the year in Michigan) on the downstream warming rates of Michigan rivers. Michigan streams are classified into four types according to July mean water temperature (JMT) (cold (C) = JMT  17.5 °C (63.5 °F), cold- transitional (CT) = 17.5 °C (63.5 °F) < JMT  19.5 °C (67 °F), warm-transitional (WT) = 19.5 31 °C (67 °F) < JMT  21.0 °C (70 °F), warm (W) = JMT > 21.0 °C (70 °F)). The stream classes provide useful predictions of fish assemblages, and differ in their predicted response to incremental reductions in baseflow. Michigan streams are also classified based on catchment area, but this variable was not a focus of this analysis. Through development of dose-response relationships, fish assemblage responses to decreased baseflows were predicted to aid in determining percentages of baseflow reduction leading to adverse resource impacts to characteristic fish populations. However, the rate of heating model was unable to replicate the magnitude of warming observed in the Hunt Creek experimental stream reach. The uncertainty associated with the physical process model and lack of data for other stream types led to the recommendation of further study of downstream temperature flux (Zorn et al., 2008). There are many models capable of predicting stream temperature (Sinokrot and Stefan, 1993; Mohseni et al., 1998; Caissie et al., 2001; Caissie et al., 2007; Cheng and Wiley, 2016); however, few studies exist for estimating factors that influence downstream temperature change rates (Rutherford et al., 2004; Magnusson et al., 2012; Davis et al., 2016). As demands for groundwater increase due to effects of climate and land use change, there is a need to understand the impacts of withdrawal on thermal dynamics of streams. Stream temperatures are driven by a suite of both natural and anthropogenic variables. Meteorological features such as components of heat flux including solar radiation, longwave radiation, evaporation, and conduction/convection, and air temperature have an established influence on stream temperature change (Webb et al., 2008). Hydrological features including discharge, depth, and width, (Webb and Nobilis, 2007; Magnusson et al., 2012) as well as groundwater inputs (Mohseni and Stefan, 1999) and connections with hyporheic zones (Hannah et al., 2009) are also known drivers of stream temperature change. There remains a gap in the quantification and modeling of relationships 32 between meteorological variables, stream discharge, and downstream warming rates; however, achieving predictability in the physical processes governing reach-scale downstream temperature response to baseflow reduction is an important step in understanding thermal dynamics of streams (Bustillo et al., 2014). Methods for predicting downstream temperature flux vary depending on the types of available data. Deterministic or physical models use meteorological forcing variables and hydrological processes to solve heat budget equations between the river and surrounding environment (Webb et al., 2008). Magnusson et al. (2012) implemented an energy balance model to describe downstream temperature change in proglacial stream reaches to better understand influential processes of longitudinal temperature increases. However, forcing data can be expensive to acquire and requires sophisticated technology to observe in situ, and can be difficult to manage and maintain. Additionally, weather stations and other sources of data can be unreliable and produce variability in estimates dependent upon location and time period, and are often unavailable near the study location (Caissie et al., 2001). Methods for deriving components of heat flux are often challenging, and are dependent upon input parameters such as relative humidity, wind speed, and cloud cover that can be highly localized. While a deterministic modeling approach is useful when conducting scenario analysis and understanding cause-and- effect responses, it can be highly complex and require input parameters that may not be available for use by models at large geographical scales (e.g., statewide models). Alternatively, empirical or regression models rely on observed relations between weather parameters and catchment characteristics to make predictions on associated temperature flux (Neumann et al., 2003; Benyahya et al., 2007). Using regression-based analysis, stream temperature change can be quantified and predicted at various spatial and temporal scales 33 (Mohseni et al., 1998; Caldwell et al., 2013). Empirical models can be less computationally intensive than process-based models, and are able to be easily implemented and validated. However, there remains a need to collect hydrological data to describe empirical responses to flow reductions and air temperatures, and the models developed to date typically are site specific Although the examination of Hunt Creek by Nuhfer et al. (2017) did not reveal considerable effects of water withdrawal on brook trout populations, flow reduction did have significant downstream effects on stream temperature. Reductions in flow increased the warming rate to extremes which would likely have negatively impacted growth or survival of salmonids in downstream reaches. Although Zorn et al. (2008) were able to approximate observed downstream warming rates of the Hunt Creek study, further examination of processes influencing longitudinal temperature change is warranted given that no other warming rate data were available for assessing the model’s applicability to other stream types. As the Hunt Creek stream reach is groundwater dominated, temperature change in other stream thermal classes may be difficult to simulate without in situ data to accurately parameterize models. Improved modeling of temperature change rates is critical to improving the predictive power of the flow- fish response model used to establish flow reduction thresholds, an essential component in managing water resources for both human purposes and fish habitat. The impetus for this research is defined by a general lack of data available for understanding stream flow and temperature relationships. Further understanding of whether warming rates vary by stream type will improve habitat and fisheries management, as managers and ecologists have already begun identifying significant relationships between stream types and fish community composition (Zorn et al., 2002, Wehrly et al., 2003, Zorn et al., 2008). Furthermore, whether warming rates are predictable for differing stream types using a single 34 model will help identify potential improvements to current management tools. The objectives of this study are to (1) collect streamflow and temperature data to quantify and model relationships between meteorological variables, streamflow variables, and downstream temperature flux rates among different stream thermal classes; (2) evaluate process-based and regression-based models to predict downstream warming rates; (3) use model selection criteria to identify models that best fit data; (4) assess differences in model accuracy and determine whether downstream flux rates differ significantly between stream thermal classes; (5) compare estimates from developed models to the physical process model used by the Water Withdrawal Assessment Tool (WWAT); and (6) evaluate implications for improving the downstream warming module. Study sites Methods Six streams within each of the four WWAT thermal classes (C, CT, WT, and W) were strategically chosen to provide geographical representation of stream types across Michigan and, if possible, to obtain data for streams potentially affected by existing water withdrawal activities. As such, stream reach selection was guided by input from local fish managers and water withdrawal staff with the Michigan Department of Environmental Quality. Table 2.1 displays site information for each stream chosen for the study site including a priori thermal class (from the WWAT database), latitude and longitude, reach length, and elevations of upstream and downstream locations. The a priori thermal classifications were based on predicted July mean temperatures using regression modeling and geostatistical kriging methods described in Zorn et al. (2008). Stream reaches were located in one of three regions throughout Michigan: southern Lower Peninsula (SLP), northern Lower Peninsula (NLP), and Upper Peninsula (UP). Of the 24 streams chosen, eleven were located in the SLP, six in the NLP, and seven in the UP. Figure 2.1 35 displays the locations of each of the 24 streams at which stream gauges were installed to collect water temperature (and pressure) and discharge measurements, as well as air temperature and barometric pressure. Data collection Thermal, hydrological, and weather data were recorded over the summers of 2015 and 2016. In 2015, data were collected from 15 streams beginning in July and extending to early November. In 2016, data were collected from 21 streams beginning in May and extending through October. Paired stream gauges were installed at upstream and downstream locations within each stream. Stream gauges were built using PVC piping attached to a fence post driven into the stream bed. Staff rulers were attached to gauges to record water level for use in developing stage-discharge curves. HOBO® U20 Water Level Loggers measuring temperature and pressure at 15-minute intervals were housed in stream gauges. Numerous holes were drilled into the bottom of the stream gauges to allow water to flow through. Mesh lining was placed around the drilled holes to prevent sediment build up within the housing. Data loggers were calibrated against each other by placing loggers in an ice bath that slowly warmed to room temperature. All sensors agreed within  0.18 °C and were corrected for observed constant offsets. Air temperature and barometric pressure were also collected at 15-minute intervals using Monarch® Track-It data loggers. Measurements were ultimately averaged to hourly intervals. Discharge measurements were taken at each of the upstream and downstream locations for each of the 24 streams using a SonTek® Flowtracker. Discharge was determined three to five times each year to collect a wide range of river stages. Stage-discharge curves were developed using stage readings from staff rulers and discharge measurements. Barometric pressure readings were used to adjust the pressure readings from water level loggers which measured total pressure 36 from air and water. The true water pressure was obtained by subtracting barometric pressure from the total pressure recorded by the water level loggers. The true water pressure was then converted to inches of water, which was then aligned with staff gauge measurements when discharge measurements were conducted. Gauge adjustments were then used to convert inches of water to reflect an estimated gauge reading for each pressure observation recorded by the water level loggers. Discharges were then calculated for each observation using equation (1): (1) where Q = discharge (m3 s-1), and G = gauge reading (in), while a and b are parameters estimated by a power function relating discharge to stream stage. Power functions (equations 2-4) were also developed to estimate widths, depths, and velocities as a function of discharge for use in water temperature modeling: (2) (3) (4) where w = stream width (m), d = stream depth (m), V = velocity (m3 s-1), and c, e, f, h, i, and j are parameters estimated using a power function. Model hypotheses A series of physical process (termed here a forcing model) and statistical models based on hypothesized mechanisms that affect stream temperature dynamics were used to predict downstream temperature change. Table 2.2 lists the models and parameters within each model that were included in the model selection procedure. Models were parameterized using nonlinear optimization (R Core Team, 2013) for each process contributing to temperature change, and were individually fit to each stream. The forcing model used processes influencing stream 37 warming described below. All parameters of the models were considered unknown coefficients, including an intercept and the various components contributing to temperature change, thus nonlinear optimization was used to estimate values for each, maximizing the fit between observed and predicted water temperature flux. Processes influencing stream warming An energy balance equation similar to the methods of Magnusson et al. (2012) was used to calculate downstream temperature change using meteorological forcing variables as well as hydrologic data. Downstream temperature change ΔT (°C) over a stream reach of length L (m) and average width w (m) can be predicted using equation (5): (5) Here, q (W m-2) is the heat flux across the stream surface, c (J kg-1 K-1) is the specific heat capacity of water, ρ (kg m-3) is the density of water, Δz (m) is the change in elevation between the up- and downstream ends of the reach, g (m s-2) is the gravitational acceleration constant, and ΔTr (°C) is residual temperature change. The first term of equation (5) can also be calculated using surface area A (m2), volume V (m3), and mean residence time of water τ (s) by (Q A τ)/(V c p). The second term represents frictional heating. Equation (5) predicts downstream temperature change as a result of heat flux across the stream surface and frictional heating due to dissipation of gravitational energy. Surface heat transfer accounted for 88 – 91% of the variance in temperature increase in most stream reaches studied by Magnusson et al. (2012). The exchange rate of heat at the stream surface resulting in cooling or warming is proportional to heat flux, length of the stream reach, and average width of the stream reach, while also being inversely proportional to the discharge. Equation (5) also assumes that available potential gravitational energy dissipates into heat. The residual 38 temperature change term represents heat that cannot be accounted for due to surface heat flux or frictional heating. Additional heat sources of the residual temperature change include groundwater advection and overland flow advection. Surface heat flux is accounted for through the following four components included in equation (6): (6) where SWnet is net shortwave radiation (W m-2), LWnet is net longwave radiation (W m-2), LE is latent heat of evaporation (W m-2), and H is sensible heat (W m-2). Each of the other components of q is described below. Latent heat flux Latent heat flux was calculated using the methods of Brocard and Harleman (1976) and is shown in equation (7): (7) where V is wind velocity (m s-1), es is surface vapor pressure (mb), and ea is air vapor pressure (mb). Vapor pressure for both the water surface and air were calculated using water temperature and the saturation pressure function in equation (8): (8) where T is the water temperature (°C) for es, and dew point temperature (°C) for ea. The vapor pressure of air was determined by multiplying an optimal value for relative humidity obtained through nonlinear optimization methods (discussed below) with the measured air temperatures. A value for wind speed was also obtained through nonlinear optimization. 39 Sensible heat flux Sensible heat flux was also calculated using the methods of Brocard and Harleman (1976) and is shown in equation (9): (9) where P is atmospheric pressure (mb), Tw is water temperature (°C), and Ta is air temperature (°C). Shortwave radiation Net shortwave radiation was calculated as follows in equation (10): (10) where α is the surface albedo set at a fixed value of 0.03 (dimensionless) (Sinokrot and Stefan, 1993), sf is a shading factor, and SWin is the incoming solar radiation (W m-2). The shading factor was another term for which an optimal value was obtained for each site using nonlinear optimization. Appendix 2.1 describes in full detail the methods used for deriving the net shortwave radiation. Longwave radiation Longwave radiation was calculated using methods described in Magnusson et al (2012) where clear-sky emissivity εclr (dimensionless) was calculated (see equation 12) using the clear- sky longwave radiation Lclr (W m-2) calculated in equation (11): (11) 40 where w is precipitable water (cm) given by . Clear-sky emissivity was then estimated using equation (12): (12) where σ is the Stefan-Boltzmann constant (5.67 x 10-8 W m-2 K-4). Effective atmospheric emissivity εa (dimensionless) was then estimated using cloud cover ccloud (%) in equation (13): (13) Cloud cover was estimated as a function of the air vapor pressure in equation (14): (14) Net longwave radiation was calculated using equation 14 as follows (Morin and Couillard, 1990): (15) Computing heat energy of various flows In order to account for the heat energy from flows contributing to total discharge, a method was used to weight the heat energy contributions based on flow temperature and flow discharge. Water temperature is a measure of heat energy concentration in a stream. Water temperature is proportional to the heat energy divided by water volume, or in other words, the heat load divided by the flow rate. Increases or decreases in heat load affect stream temperature by altering the amount of heat energy in the system (Poole and Berman, 2000). In this analysis it was necessary to create parameters for both groundwater and overland flow because they have different thermal patterns that affect stream temperature differently. Change in water temperature 41 as a function of upstream water temperature ΔTup (°C) and the difference in upstream and downstream discharge was calculated using equation (16): (16) where Qup is the upstream discharge (m3 s-1), Tup is the upstream water temperature (°C), and Qdown is the downstream discharge (m3 s-1). Baseflow gain within each stream section (denoted by Qbase (m3 s-1)) was estimated for each stream by determining the gain in discharge from the upstream gauge to the downstream gauge for the lowest 7-day flow period for each year. Contributions of baseflow heat energy ΔTbase (°C) to downstream temperature flux were calculated using equation (17) as follows: (17) where Qbase is the baseflow discharge (m3 s-1) and Tbase is the baseflow temperature (°C). Baseflow temperature was set as a constant particular to each of the three regions where streams are located throughout Michigan (UP = 5.6 °C; NLP = 8.3 °C; SLP = 11.1 °C (USEPA, 2016)). Finally, heat energy of overland flow ΔTover (°C) was determined using equation (18), (18) where overland flow discharge was estimated as the remaining flow after subtracting upstream and baseflow discharges from the downstream discharge. Temperature of the overland flow was determined using a moving air temperature average (°C) where air temperature was averaged over a 12-hour period for each observation. As water temperature is proportional to the heat load divided by discharge, the downstream change in water temperature caused by heat energy inputs can be accounted for through the ΔTflow variable. Equation (19) combines equations 42 (16 – 18) and creates a variable ΔTflow (°C) which describes the cumulative heat energy due to each of the flow components: Statistical models (19) Models were built in a hierarchical fashion sequentially incorporating parameters for processes judged to be dominant in influencing downstream temperature flux. As such, the base model (Model 1) was the simplest and subsequent models became progressively complex. As water temperature is influenced by heat energy and water volume, parameters included in the models were based on hydrologic and meteorological principles known to have an influence on water temperature heating dynamics as described above. Air temperature – water temperature differential (Model 1) A major heat exchange process occurs at the air-water interface. As the rate of heating depends on the magnitude of temperature difference between ambient air temperature and water temperature, incorporating the air-water temperature differential into the model was a first step to model downstream temperature change. The following model was used as a base model and consisted of an intercept and beta coefficient associated with the differential between the air temperature and upstream water temperature: [1] Discharge ratio (Model 2) Discharge ratio was incorporated into the base model to reflect the amount of discharge gained between the upstream and downstream gauges. A baseline value of 1:1 can be used to evaluate the thermal dynamics of streams in situations where there is no loss or gain in 43 discharge. However, each additional source of discharge may alter the thermal inertia of the stream. Sources of additional discharge may include baseflow gain, overland flow, or precipitation directly into the stream channel. Model 2 incorporated the upstream-downstream discharge ratio into the base model as follows: [2] Flow temperature change and upstream discharge (Model 3) Variations in stream temperature have been shown to be inversely proportional to stream discharge (Magnusson et al., 2012). Incorporating a variable accounting for the upstream discharge may explain some of the variation in downstream temperature flux since reductions in stream discharge can lead to increased residence times of water within stream reaches, leading to increased exposure to elements of surface heat flux. A parameter was also added to describe the cumulative heat energy of the stream reach using the ΔTflow variable. Fluctuations in stream temperature are the result of altering the amount of heat energy added to the stream (Poole and Berman, 2001). Model 3 included a parameter for upstream discharge, as well as a parameter for the heat load of the three sources of discharge: [3] Downstream – upstream discharge differential (Model 4) As another measure of the gain in discharge, the downstream – upstream discharge differential was added to account for the potential effects of different sources of discharge on downstream temperature change. The discharge differential between the downstream and upstream locations was incorporated into model 4 as follows: [4] 44 Day length (Model 5) To account for the effects of solar radiation without having direct measurements of solar heat flux, day length was used as a proxy. Day lengths specific to each stream reach and day of year were calculated using equations 20 – 22 below (Meeus, 1991), and incorporated into model 5: [5] Day length S (hr) was derived using equation (20): (20) where τ is the local hour angle of the sun and derived by equation (21), (21) where lat is the geographical latitude. Δ is the declination angle of the sun calculated using equation (22) (22) where x is the number of days since the vernal equinox (March 21). Altitude angle (Model 6): Solar radiation is a driver of heat flux occurring at the stream surface. The sun altitude angle was used as a measure of the sun’s intensity at any given hour throughout the day. This allowed for differentiating between diurnal variation in downstream temperature change, as well as the sun’s intensity throughout the day. Altitude angle α was used as a measure of the sun’s intensity based on the hour of day and day of year. Using equations 23 – 26 below (Meeus, 1991), sun altitude angle was implemented to create model 6: [6] 45 Altitude angle was derived using equation (23) as follows: (23) where h is the hour angle derived by equation (24) (24) Equation (25) calculates the apparent solar time AST : (25) where LST is the local standard time (hr; adjusted for daylight savings time DST, if necessary, such that LST = DST – 1 hr), long is the geographical longitude, LSTM is the local longitude of standard time meridian where , and ET is the equation of time in minutes approximated by equation (26) (26) where , and N is the day of year. Including day length and sun altitude angle (Model 7) An additional model was created to determine whether inclusion of both day length and sun altitude angle improved model fit. [7] Differential effects of flow sources (Models 8, 9, and 10) Additional models were developed by creating separate parameters for the heat load of upstream flow, baseflow, and overland flow. Groundwater and overland flow have different thermal patterns, so individual parameters were included for each. Although Model 3 included a parameter evaluating the effects of cumulative heat energy of the three components of flow, 46 Models 8, 9, and 10 include coefficients to describe the individual contributions of upstream, baseflow, and overland flow from equations (14 – 16), respectively, to downstream temperature change to evaluate whether reductions in sums of squared errors outweighs increased parameterization. These models include either one or both of the day length and altitude angle variables to predict downstream temperature change. [8] [9] [10] Incorporating components of surface heat exchange into a statistical model (Model 11) A hybrid model was developed that incorporated aspects of the forcing model described above, into a statistical model with coefficients to describe the influence of drivers of heat flux on downstream temperature flux. The components of heat flux were modified to remove values considered as constants (ρ, c, L, W, and frictional heating) or approximations obtained through nonlinear optimization (s, dew point temperature, sf, α). As sensible heat flux H is dependent upon latent heat flux, the sensible heat flux term was not included in the hybrid model. In addition to the drivers of heat flux, the air – water temperature differential, baseflow heat load, and overland flow heat load were also included. When considering the shortwave radiation, altitude angle was used as a proxy for shortwave heat flux at the stream surface. 47 [11] Forcing model (Model 12) The deterministic forcing model was developed based on the methods of Magnusson et al (2012) using the principles of heat flux and processes influencing downstream warming described above. As mentioned previously, downstream temperature flux is dependent upon changes in discharge. Changes in stream temperature are directly related to the mean residence time of water in a stream reach, meaning that increased residence times lead to increases in downstream temperature flux. Furthermore, additional heat input to the stream channel through other sources of discharge can alter the thermal dynamics of the stream. In this case, it is important to account for the heat energy of baseflow and overland flow. Surface heat flux components were derived using in situ and site information data (i.e., elevation, latitude), as well as parameters obtained through nonlinear optimization for required variables which were not collected (wind velocity, relative humidity, shading). Using equation (5) which describes the components of heat flux influencing longitudinal stream temperature change, and equations (16) and (17) to account for the temperature change due heat loads of baseflow and overland flow, the forcing model in equation (27) was developed: (27) Model fitting, selection, and prediction Nonlinear optimization methods were used to parameterize unknown components of each of the aforementioned models in order to minimize the residual sums of squares between modeled and observed ΔT (equation (28)): 48 (28) As all parameters in the statistical and hybrid models were considered free fitting coefficients, no constraints were placed on the possible values. However, when using nonlinear optimization to parameterize unknown variables in the forcing model, constraints were placed on the variables such that lower and upper bounds limited parameterization to realistic values. For example, values for shading factor and relative humidity are proportions between 0 – 100, thus, these were the constraining values used for these variables. In order to assess each model’s accuracy in predicting ΔT, model selection was applied using Akaike’s Information Criterion (Akaike, 1973) (AIC), and models were prioritized using the weighting method of Buckland et al. (1997). AIC is one of many model selection methods that prioritize models based on fitting ability and level of parsimony (e.g., Thayer et al., 2007). Equations (29) and (30) display the AIC and likelihood L (Seber and Wild, 1989), respectively, for finite sample size: (29) (30) where k is the number of unknown parameters, and n is the sample size (number of hourly observations; varied for each stream). Using the above criteria, the models with the lowest AIC score fit the observed data best. Models were then weighted following equation (31) below: (31) 49 where M is the total number of models and Δi is the difference in AICc of an individual model compared to the best model. Model weights determined one model’s fitting ability relative to another. Partial influence Partial regression analysis was used to assess the influence of each of the variables included in the best statistical model on downstream temperature flux. Partial R2 values are useful for understanding the residual variation accounted for by a predictor variable that cannot be explained by a constrained model. In order to determine the partial R2 of a particular predictor variable, the proportion of residual variance can be accounted for by , where SSR is the sums of squared errors of the reduced model, and SSE is the sums of squared errors of the full model. Model performance Root mean square error (RMSE) which is given by: with Pi and Oi being the predicted and observed values, respectively, was used to compare relative performance among the three best models. The RMSE was calculated on a monthly, annual, and overall basis. Pooling data All of the above models were developed on a site-by-site basis, and as such, provide estimates specific to the sites studied. In order to determine if flux rates could be successfully modeled at a more aggregated level, I used a mixed general linear model (e.g., Littell et al. 2006), to determine whether streams among the a priori thermal classes exhibited significantly different flux rates and should be modeled separately. The structure of the mixed model was 50 based on model 10 above with the exception that data from each stream was treated as a random effect, and thermal class was treated as a fixed effect. The results of the mixed model were evaluated by comparing observed and predicted values for each stream reach. LOESS regressions (Cleveland, 1979) were then used to provide a smoothed graphical depiction of observed downstream temperature flux with predictions made using the best model parameterized specific to each pooled dataset. Summary statistics were also calculated for each stream and averaged over thermal classes. Model comparison with WWAT A critical component of this analysis was to compare the downstream warming rates predicted by the statistical models developed in this study, and those predicted by the downstream warming rate model developed by Zorn et al. (2008) that informed the fish assemblage response curves which provided the four management zones (A, B, C and D) used in the WWAT. The Zorn et al. (2008) downstream warming rate model uses components of a physical process model. Without the use of site-specific inputs such as shading factors, groundwater inputs, and land use/land cover, the formula is based on aspects of the energy balance equation and Newton’s Law of Cooling/Heating such that Tw can be derived using equation (32) as follows: (32) where TE is the equilibrium temperature (°C), Ti is the initial water temperature (°C), and t is the travel time of water (hours) in each stream segment. The heat exchange coefficient (1/h) k is determined using equation (33): (33) 51 where d is the depth (m). For the purposes of scenario analysis of incremental baseflow reductions, hydraulic geometry relationships were used to estimate depths and velocities at selected levels of baseflow reduction to determine downstream warming rates. Equation (34) was used to calculate depth using Qbase as the input variable: (34) Equation (35) was then used as a function to adjust velocity at the same rate as depth. The predicted velocity is scaled to the initial velocity Vi using the initial depth di as follows: (35) Downstream warming rates between the two models were predicted for the month of August since the WWAT uses the August 50% exceedance flow (referred to as the Index Flow) as one of the key habitat variables important to fish metabolism, survival, reproductive success, distribution, and abundance. Average values for the month of August were used. Downstream warming rates were compared between the overall best model based upon model weight, and the Zorn et al. (2008) downstream warming module. Scenario analyses were conducted to project the impact of reduced baseflow on downstream temperature change. Impacts of baseflow reduction on downstream warming rates used the best model, as determined by model weight, by dynamically reducing the amount of baseflow input such that reductions in baseflow ultimately reduced the downstream discharge. Downstream warming rates were scaled to reflect downstream temperature flux per kilometer of stream reach ΔT/L (°C/km). Site information Results Hydrological data were successfully collected from 21 of the 24 streams from Table 2.1. For each of the East Branch Waiska River, Hemingway Lake Outlet, and Unnamed Gun River 52 Tributary reaches, stream gauges were unable to be recovered or became dislodged during deployments. In 2015, the East Branch Waiska River was heavily inundated and discharge measurements were unable to be taken as water levels overtopped the stream gauge, which was ultimately lost. The Hemingway Lake Outlet stream reach was dammed by beavers leading to inaccurate flow measurements for the 2016 field season. Additionally, during the 2016 field season, the Unnamed Gun River Tributary experienced a flooding event that dislodged the stream gauges at both upstream and downstream locations, and the gauges were eventually reinstalled. After reinstallation, the downstream gauge was eventually lost, possibly due to another flooding event, or human removal. Over the 2015 and 2016 field seasons, two years of data were successfully collected for 10 streams, while the remaining 11 streams had one year of data. When selecting stream reaches for this study, streams were chosen based upon the a priori thermal classification predicted using the regression model developed by Zorn et al. (2008). Streams were selected to provide an even distribution from each of the four thermal classifications (C, CT, WT, and W). Table 2.3 displays the a priori and a posteriori thermal classes of each of the 21 streams in this study. July mean water temperatures (°C) were averaged across years if data was collected across two years. The a priori thermal classifications were correct for 11 of the 21 streams. Based upon the a posteriori thermal classifications, the final number of stream thermal classifications was 3, 10, 5, and 3 for cold, cold-transitional, warm- transitional, and warm streams, respectively. Observed flux rates Analyzing downstream temperature flux in comparison with stream discharges provides information about flow conditions under which streams display greater temperature fluctuation 53 rates. Figure 2.2 shows downstream temperature flux rates compared to the upstream discharge (m3 s-1) for representative stream reaches amongst each of the four thermal classes. Observed temperature fluctuations displayed a strong response to discharge. Downstream temperature flux rates appeared to be inversely related to discharge (Figure 2.2-2.3). Variation in flux rates among each stream were typically reduced as discharges increased (see supplemental files). These trends were similar when comparing flux rates with discharges among thermal classes as a whole (Figure 2.3). Model selection Models were fit to predict downstream temperature change at hourly intervals for the 21 streams with at least one year of data. Predicting downstream temperature change using the forcing model (model 12) provided a poor fit to the data (Table 2.4), and so was not considered in the model selection process. Results using the forcing model led to extremely high sums of squared errors and low correlation for nearly all stream reaches. However, the SSE using Model 12 was the lowest among the models for the East Branch Black River, although the correlation was not particularly notable. Parameterization of each of the shading, wind speed, and relative humidity variables resulted in values that were outside of the range of realistic values. For example, the shading coefficient of each stream was optimized to approximately 90%. While some smaller streams may show 90% shading, higher order streams were unlikely to reflect this amount of stream cover. Hourly observations of downstream temperature flux rates are displayed in Figure 2.4 for example streams (one stream from each thermal class). The vertical and horizontal scatter of hourly observations is indicative of diurnal and seasonal patterns in temperature change. A LOESS regression was used to provide a smoothed representation of the hourly data, and display 54 general seasonal patterns. Downstream warming rates typically peaked in summer, while values were often lowest during spring and late fall. The hourly observations indicate that streams experience warming and cooling events diurnally, and the extent of each varies by stream. The LOESS regressions were also useful for visually comparing model fit with observed data (Figures 2.5 – 2.8). The models ultimately selected for ranking were statistical models developed using hydrological and meteorological processes contributing to heat flux. Of all the models tested, only four models (Models 8, 9, 10, and 11) received consideration for best fit by model weight (Table 2.5). However, two models (Models 10 and 11) provided the best overall fit by model weight to any of the streams. The top ranked models contained separate parameters controlling for each of the individual components of stream heat load (upstream flow, baseflow, and overland flow). Additionally, these models controlled for at least one of the components of heat flux at the stream surface (e.g. day length and/or sun altitude angle). The best model (Model 10, Table 2.5) included both day length and sun altitude angle to account for the effects of solar heat flux at the stream surface. This model tracked overall trends well (r = 0.62), but would occasionally underestimate peak periods of downstream warming and cooling (Figure 2.5). For some streams, this model predicted peak downstream temperature change events earlier or later than observed (see supplemental files). This model provided the best fit for 76% (16/21) of streams (Table 2.6). Model 10 tracked trends the best for streams with both one (r = 0.66) and two (r = 0.58) years of data (Table 2.7). Average model weight for Model 10 was 0.74 (Table 2.5), far greater than the second best model (Model 11) with an average weight of 0.24. 55 The second best model (Model 11, Table 2.5) was a statistical model developed to incorporate the effects of heat flux at the stream surface. Instead of including a parameter to control for the effects of the upstream heat load, this model included parameters to emulate evaporative, solar, and longwave heat flux along the downstream gradient. Additionally, this model included parameters for the overland flow and base flow components of heat load. Model 11 did not match trends as well (r = 0.50) as the best model (Figure 2.6). Model 11 also failed to match peak downstream heat flux timing and values for many streams (see supplemental files). This model did not perform better for any particular stream thermal class, but instead was identified as the best model for at least one stream in each of CT, WT, and W thermal classes. Model 9 provided the best fit by model weight for 5 of the 21 streams with an average model weight of 0.24 (Table 2.5). The third best model (Model 9, Table 2.5) included sun altitude angle as a substitute for solar heat flux. This model captured trends nearly as well as the best model (r = 0.56), but more often failed to accurately predict peak downstream temperature fluxes (Figure 2.7). This model also failed to capture the full range of downstream temperature fluxes of some streams (see supplemental files). Model 9 did not provide the best overall fit for any of streams, but received consideration for the best fit by model weight in four streams. Model 9 had the third highest average weight at 0.024. The fourth best model (Model 8, Table 2.5) included day length as a substitute for solar heat flux. This model captured trends nearly as well as the best model (r = 0.59), but failed to fully capture the range of peak downstream temperature fluxes (Figure 2.8; see supplemental files). Model 8 did not provide the best overall fit for any of the streams, and had a very low 56 model weight of 0.001. Model 8 received support by model weight in just one stream (Table 2.6). Seasonal patterns Model fit varied seasonally and yearly in each of the best models (Figures 2.5 – 2.8). Predictions of downstream temperature change were less accurate at the beginning and ends of each year. Using correlation as a measure of how well models tracked trends, models tended to track better in streams for which there was only one year of data collected compared to streams with two years of data (Table 2.7). Each of the three models presented issues in the timing of peak downstream temperature flux for both downstream warming and cooling. The seasonality of errors between observed downstream temperature fluxes and those predicted using the best models identified by model selection criteria are depicted in Figure 2.9. The differences that occur between the four different models are small (see supplemental files). Summertime downstream temperature fluxes tend to be underestimated by each of the four best models. Estimations of fall downstream temperature flux are more variable; models tended to underestimate the downstream cooling effects experienced by streams beginning in early fall. When using the RMSE of the four best models to evaluate performance across all streams, all models had similar overall RMSE (Table 2.8). Models 8, 9, and 10 performed best with an RMSE of 0.12 °C, while Model 11 performed similarly with an RMSE of 0.13 °C. Each of the four models performed best in October (M8 = 0.10 °C; M9 = 0.11 °C; M10 = 0.10 °C; and M11 = 0.11 °C) and worst in July (M8 = 0.13 °C; M9 = 0.13 °C; M10 = 0.13 °C; and M11 = 0.14 °C) (Table 2.9). Table 2.10 displays the performance based on RMSE of each of the models among the four thermal classifications. Models performed best in cold streams (M8 = 0.06 °C; 57 M9 = 0.06 °C; M10 = 0.06 °C; and M11 = 0. 065 °C) and worst in warm-transitional streams (M8 = 0.16 °C; M9 = 0.16 °C; M10 = 0.16 °C; and M11 = 0.17 °C). Peak events The maximum and minimum downstream temperature fluxes of both the observed and predicted values of each of the four best models are displayed in Table 2.11. The highest downstream temperature warming rate of 1.74 °C/km was observed in Squaw Creek, while the highest rate of downstream cooling of -1.43 °C/km occurred in Morgan Creek. The average maximum and minimum downstream temperature flux rates across all streams were 0.17 °C/km and -0.14 °C/km, respectively. Model 11 had both the highest average maximum value at 0.14 °C/km, and lowest average minimum value at -0.09 °C/km. Models were occasionally unable to capture the directionality of downstream temperature flux. For example, Nottawa Creek experienced a maximum downstream temperature flux of 0.03 °C/km, however, none of the models predicted a downstream temperature flux greater than -0.01 °C/km (Model 10; Table 2.11). A similar situation occurred for the East Branch Black River which had a minimum downstream temperature flux of -0.10 °C/km, but only one model (Model 9; Table 2.11) predicted a negative downstream temperature flux at any point. In order to examine each model’s ability to capture peaks and troughs of downstream temperature flux, the residuals between observed and model predictions of the 1st and 99th percentiles were compared (Table 2.12). Models performed similarly in their abilities to estimate peaks and troughs of downstream temperature flux. Models tended to underestimate peak downstream warming and cooling events, such that models were unable to fully capture the range of downstream temperature changes experienced by the streams. Model 10 most closely captured downstream temperature fluxes in both the 1st and 99th percentiles. 58 Quantifying influence on downstream temperature flux The influence of each variable included in Model 10 on the rate of downstream temperature flux was examined through partial regression analysis. On an overall basis, day length had the most influence on the flux rate with a partial R2 of 0.09, while the discharge differential between upstream and downstream gauges had the least at 0.02 (Table 2.13). When averaging partial R2 values over thermal classes, the variables with the greatest influence over flux rates was dependent upon stream thermal class (Table 2.14). Day length remained the most influential variable for both cold and cold-transitional streams, while baseflow heat load and altitude angle provide the greatest influence for warm-transitional and warm streams, respectively. Comparison of pooled datasets The GLM used to compare the differences in temperature flux rates revealed that incorporating the effects of individual stream reaches (i.e., thermal dynamics unique to each stream reach) as random effects improved model performance as indicated by the model R2. When accounting for the variance between individual stream reaches, R2 increases from 0.30 to 0.57 (Table 2.15). Incorporating the random effects shows evidence for stream to stream variation in thermal flux rates. Although significant differences exist in the downstream thermal flux rates of thermal classes (Table 2.15; 2.17), flux rates are also significantly different between individual stream reaches (Table 2.16 – 2.17). Random effects of individual streams should be taken into consideration when modeling downstream flux rates. The results from the GLM show that individual streams should be assessed relative to other streams within the same thermal class. 59 Pooling data allowed for comparing model accuracy across streams of similar thermal classes and across the entire set of streams sampled. Using LOESS regressions to visually compare fit among the example streams in Figure 2.10, the pooled data set which provides the best fit based on model optimization varies between thermal classes. Results are similar among the rest of the streams within the study streams (see supplemental files). Summary statistics including r, SSE, and RMSE tended to improve for individual stream reaches (Table 2.18) and overall thermal classes (Table 2.19) when optimizing Model 10 for each thermal class as compared to the pooled data set. Using the pooled data set, cold streams had the best RMSE at 0.16 °C, while warm streams had the worst RMSE at 0.35 °C. Trends were tracked best based upon r in warm streams (r = 0.52), while trends in warm-transitional streams tracked worst (r = 0.12) When parameterizing Model 10 specific to each thermal class, RMSE improved across all thermal classes. Based on RMSE, the model performed best in cold streams (RMSE = 0.12 °C), while performing worst in warm-transitional streams (RMSE = 0.34 °C). Based on r, the model performed similarly in warm (r = 0.45) and cold (r = 0.45), and most poorly in cold-transitional streams (r = 0.19). Baseflow reduction scenarios Downstream temperature flux rates show a dynamic response to baseflow reduction where effects vary by season and reduction scenario (Table 2.20; Figure 2.11; see supplemental files). In the model, as baseflow is reduced, so are each of the downstream discharge and downstream – upstream discharge differential variables. The dynamic effects lead to interesting results where some streams actually show downstream cooling in the 90% reduction scenario, and downstream warming in the 0% reduction scenario (e.g.., Fish Creek; Figure 2.11). More interestingly, Fish Creek also appears to show a greater degree of downstream cooling later in 60 the year under the 90% reduction scenario when compared with the 0% reduction scenario (i.e., October; Figure 2.11; Table 2.20). This is potentially due to the lesser volume of water in the stream channel, which becomes more susceptible to a greater cooling due to air temperature. Cedar Creek is estimated to experience a similar phenomenon under the 90% reduction scenario, where the flux rate is projected to cool at 0.46 °C/km. The reduced volume of water combined with the shaded study reach of Cedar Creek may provide a unique scenario allowing for increased cooling, although this may require additional inputs of cool groundwater somewhere further along the reach. Occasionally, streams displayed an increase in downstream warming rates following baseflow reduction. Cedar River, for example, was predicted to show an increase in monthly downstream warming for each of the reduction scenarios listed in Table 2.20. This trend also held true for several other streams within the study. Butterfield Creek was one stream particularly susceptible to reductions in baseflow. Under baseflow conditions, Butterfield Creek was predicted and observed to experience downstream cooling; however, under the most extreme reduction scenario, the flux rate shifted to downstream warming. Given these results, it is clear that individual streams show complex responses to baseflow reductions throughout the year. Thermal class responses to baseflow reduction varied by month and among each class. The warm thermal class was predicted to experience downstream cooling under a 0% reduction scenario for each month, while the other three classes were predicted to warm with downstream flow (Table 2.21). These results suggest that thermal sensitivity varies by thermal class. Cold streams showed the least thermal sensitivity between the 0% and 90% reduction scenarios, with the greatest difference in downstream flux rates of 0.12 °C occurring in October. The greatest rates of downstream warming under each flow reduction scenario examined were projected for 61 October within warm-transitional streams. Interestingly, warming rates were projected to decrease as baseflow reduction increased. When focusing on summer months, the greatest impact of baseflow reduction on downstream flux rates occurred within the warm thermal class. The downstream flux rates in warm streams were projected to increase by 0.59 °C and 0.30 °C in July and August, respectively. Comparison with the Zorn et al. (2008) downstream warming rate model. As part of this analysis, longitudinal downstream temperature flux rates (°C/km) for August were estimated using Model 10 and were compared to rates estimated by the Zorn et al. (2008) downstream warming rate model. August is often the lowest flow month of the year, and water temperature can peak during this time having a dominant effect on fish physiology, growth, and survival. Mean values used to parameterize Model 10 (Table 2.22) and the Zorn et al. (2008) model (Table 2.23) were calculated for August. The greatest downstream temperature flux rate was observed in Squaw Creek at 0.77 °C/km (Table 2.24), while Hasler Creek experienced the greatest cooling rate at -0.26 °C/km. On average, cold streams experienced the greatest gain in downstream temperature change of 0.14 °C, while warm streams experienced an average downstream temperature change of -0.11 °C (Table 2.25). During August baseflow conditions, Model 10 most accurately reflected the average downstream temperature flux of each thermal class (Table 2.24). The Zorn et al. (2008) model severely overestimated downstream temperature flux for both downstream warming and cooling events during baseline scenarios with no baseflow reduction for individual stream reaches and for stream thermal classes as a whole. Potential effects of incremental baseflow reduction on August downstream temperature flux rates were compared using both Model 10 and the Zorn et al. (2008) model. Incremental 62 reductions in baseflow using Model 10 resulted in concurrent reductions in downstream discharge. Using Model 10, the impact of simulated baseflow reduction was variable between individual stream reaches (Table 2.26). The flux rates of many streams were predicted to remain relatively unchanged following reductions in baseflow on average throughout the study period (Figure 2.12). For example, 15 streams were estimated to experience  0.10 °C/km. Although the flux rates were predicted to change following baseflow reduction for most streams, the response for nearly all streams was predicted to be less than 0.50 °C/km in the most extreme reduction scenario (Table 2.26). Cedar Creek was predicted to cool by an additional 0.57 °C/km at 90% baseflow reduction. The flux rates of eight streams were predicted to decrease following baseflow reduction, while nine streams were predicted to experience an increase in downstream flux rates. The greatest increase in downstream temperature flux rate was predicted in Squaw Creek which was predicted to gain 0.80 °C/km following a 10% reduction, ultimately resulting in a gain of 1.20 °C/km in the most extreme scenario of 90% baseflow reduction. The predicted flux rates in response to baseflow reduction for Butterfield Creek are notable, as it is predicted to experience downstream cooling, but following baseflow reduction it is predicted to experience downstream warming beginning at an 80% baseflow reduction scenario; however, this is the only stream predicted to experience this phenomenon. When using the Zorn et al. (2008) model to estimate downstream temperature flux under baseflow conditions, the magnitude of change was much greater than observed fluxes (Table 2.24). In addition, the Zorn et al. (2008) model predicted much greater rates of downstream temperature flux than those predicted using Model 10 under baseflow reduction scenarios. The rates of longitudinal temperature flux were heavily influenced by the air-water temperature differential such that the rate of heating or cooling was directly related to the difference between 63 air temperature and water temperature. For example, each stream for which the air temperature was warmer than water temperature was predicted to experience downstream warming, and the opposite was true for those streams where the water temperature was greater than air temperature (Table 2.23 – 2.24). There were also some streams for which the downstream temperature flux rate was beyond the range of the air-water temperature differential. For example, the East Branch Black River was predicted to gain 11.67 °C/km although the air-water temperature differential was 5.07 °C. This trend held true for 13 of 21 (62%) stream reaches in the study. Estimations of downstream temperature flux rates for the North Branch Thunder Bay River and King Creek were likely affected by the large differential between air and water temperatures of -9.31 °C and -8.15 °C, respectively. The large differential, as well as low depths and slow velocities likely contributed to the extreme flux rates for these two streams. Predicting the response of downstream temperature flux rates to baseflow reductions using the Zorn et al. (2008) model was a dynamic process whereby stream depth, velocity, and discharge were reduced with each incremental reduction in baseflow. The downstream temperature flux rate of streams displayed an exponential response to baseflow reduction (Figure 2.13). The rates for King Creek and Squaw Creek were omitted to maintain an appropriate scale on the y-axis. The flux rate of the North Branch Thunder Bay River experienced the greatest response of -0.73 °C/km to an initial 10% reduction in baseflow, while the flux rate of the Carp River was initially resistant to the 10% baseflow reduction. The stream most affected by the 90% reduction scenario was Hasler Creek, which had a predicted flux rate of -28.67 °C/km, while the Carp River was the least affected with a flux rate of 0.18 °C/km. 64 As water use throughout the United States steadily increases (Kenny et al., 2009), the Discussion resulting alteration of hydrologic conditions pose threats to aquatic ecosystems (Arthington et al., 2006). Demands for groundwater for agricultural, commercial, and public usage are expected to rise due to land use change and increasing thermal stressors through climate change. In response to increasing demands, it is important to understand stream thermal behavior in order to project the response of streams to hydrologic alteration. Michigan has devoted significant conservation efforts in an attempt to balance the needs between ecosystem function and economic benefit by developing the WWAT to estimate the potential effects of baseflow depletion on stream ecology. The primary objective of this analysis was to fill in a major data gap in the understanding of effects of flow reduction on downstream temperature flux rates. Data collection on small streams allowed for the design and comparison of models used to estimate downstream temperature flux rates and subsequently project stream thermal response to reductions in baseflow. While it was clear that thermal flux rates showed an inverse relationship with discharge (Figure 2.2; 2.3), it was important to determine whether flux rates differed between Michigan stream thermal classifications. Although results indicate that the thermal dynamics of different stream classes differed significantly, the greatest variability lies within the specific characteristics of individual stream reaches (Table 2.17). As such, the similarity in thermal dynamics within a stream class provides some level of predictive capability, however accounting for the random effects of individual reaches implies that site-specific measurements may be needed for accurate and precise predictions. 65 Model comparison The objectives of this analysis were to design and compare models that incorporate hydrological and meteorological processes to estimate downstream temperature flux rates, and to project the influence of baseflow reduction on flux rates. I also sought to compare these results with the Zorn et al. (2008) model that informed the fish assemblage response curves that underpin Michigan’s WWAT. Using a network of stream gauges throughout Michigan, a suite of statistical models were developed and compared using AIC model selection to identify the best model based on parsimony and goodness of fit. AIC was useful in identifying processes most influential in fitting observed data, while also considering the least amount of parameterization. Interestingly, the model identified as providing the best overall fit also contained the greatest number of parameters (Table 2.5); however, the reduction in sums of squared errors outweighed the penalties incurred through overparameterization. The multi-model selection approach was useful in this analysis given the hierarchical approach in model development where models became increasingly more complex, and often contained many of the same or similar parameters. Given the simultaneous incorporation of goodness of fit and parsimony in model selection, AIC was able to distinguish the best model from other similar models that may have gone unnoticed simply relying on other performance- based metrics such as RMSE, and visual comparisons using LOESS regressions. This method also proved to be useful in highlighting the importance of processes which may have been underestimated (Thayer et al., 2007). For example, although day length and sun altitude angle are similar variables, the inclusion of both in the best model was ultimately important in providing the best fit to the observed data at the fine level of temporal granularity (i.e., 1 hr intervals) used in modeling. 66 Model performance After identifying the four models that best predicted downstream temperature change, model accuracy and precision were examined through a series of model performance criteria. The models performed similarly in their accuracy of overall and peak estimations (Tables 2.8 – 2.12). The measure which likely separated the best model (Model 10) from the others in the suite was the ability to better match daily and seasonal trends due to the inclusion of variables representing day length and sun altitude angle. Comparing sums of squared errors of each of the models, the reduction using Model 10 indicates a greater fit to hourly observations. Based on LOESS regressions, Model 10 appears to provide a greater fit than other models throughout the study periods indicating a greater ability to match seasonal trends (Figures 2.5 – 2.8). Including both variables into the model helped account for the variability associated with fluctuations of solar radiation at the sub-daily time scale. An important consideration beyond the scope of this investigation is the temporal granularity of data collection and modeling, and the implications of the resolution of data collection. In this work, I collected environmental measurements at a 15-minute interval, but averaged data on an hourly interval to smooth out measurement variability. Thus, the results provided represent temperature flux at a fine temporal scale. Other investigators have commonly used daily means to represent stream temperature, particularly when relating temperature to biological responses. While neither temporal scale is inherently “correct”, interpretation of modelling results are somewhat dependent on the scale. For example, sun angle was found to be an important contributor to temperature flux at the hourly scale, but would not be a factor that could be included at a daily time resolution. As such, I recommend caution in interpreting and 67 comparing data and modeling results that have different temporal granularity or that cover different time periods (e.g., summer only, versus spring through fall data collection). One main issue associated with Model 10 and its ability to estimate downstream temperature change rates was the failure to fully capture downstream warming and cooling peaks. Using residuals of the 1st and 99th percentiles, Model 10 underestimated peak downstream warming and cooling by 0.11 °C/km and 0.09 °C/km, respectively (Table 2.12). The failure to capture peaks of downstream temperature flux is noteworthy when attempting to estimate the effects of baseflow reduction on downstream temperature flux rates. Underestimations of the effects of baseflow reduction on flux rates can be a complicating factor when establishing withdrawal limits. Failure to accurately estimate the impacts of water withdrawal on downstream warming rates may lead to underestimating the consequences on fish habitat loss. The range of predictability would likely be improved by direct readings of solar heat flux. Solar heating is a primary influencer of stream warming, while wind speed and evaporative heat flux could influence downstream cooling (Keith et al., 1998). There also remains the possibility that greater influxes of baseflow than were predicted could cool streams more so than was predicted using the best model. Furthermore, the spatial distribution of the baseflow would likely influence downstream temperature flux. If baseflow entered the stream primarily at the beginning of the study reach, the segment would be more susceptible to warming than if baseflow inputs occurred throughout the reach. Hyporheic exchange has also been hypothesized as a source of stream cooling (Moore et al., 2003) and an influence on creek temperature overall (Poole and Berman, 2001). Currently, common models used for representing stream temperature dynamics are based on concepts related to physical processes that incorporate energy and water balance equations to 68 estimate heat flux at the stream bed and stream surface (Sinokrot and Stefan, 1993; Caissie et al., 2005; Caissie et al., 2007). One of the main drawbacks to physical process models is that they rely upon sophisticated technology and intensive data processing to estimate surface/streambed heat exchange. Issues arise when estimating sub-daily downstream temperature flux across multiple seasons using optimized parameters that result from lack of data on wind speed, shading factors, and other components used to derive heat flux. Without these data, relying upon empirical formulas and nonlinear optimization methods to derive the components of heat flux led to inaccurate estimates of stream temperature dynamics. For example, constant values of shading – estimated through calibration – fail to accurately reflect seasonal trends in vegetative cover. Furthermore, using constant values for relative humidity and wind speed can under- and overestimate the amount of heat flux occurring at the stream surface, causing inaccurate estimates of downstream temperature change. Relative humidity is required to derive all components of heat flux but shortwave radiation. As values fluctuate diurnally, and since this study used hourly observations throughout different seasons, constant values fail to capture this variation. The value of the shading parameter used to buffer the amount of solar radiation reaching the stream surface in streams with significant vegetative cover is likely to fluctuate throughout the year due to seasonal transitions. When examining the effects of patchy shade on the rate of change of daily maximum temperatures among second order streams, Rutherford et al. (2004) highlighted a strong linear relationship between heating/cooling rates and change of shade. Riparian shading has been shown to have an effect on the microclimate of the stream corridor such as increases in air temperature, and decreases in relative humidity (Chen et al., 1995). Additionally, values for wind speed are influential in estimating evaporative heat flux. As with 69 estimating solar radiation at the air-water interface, wind speed, and thus evaporative heat flux, is buffered by vegetative cover and is dependent upon seasonality. Finally, estimating longitudinal temperature change using physical process models becomes difficult during low flow periods as surface heat exchange warms or cools streams at a rate proportional to the mean residence time of water in the stream reach, extending exposure time of the stream surface to elements of heat flux. Assuming a constant rate of heating without mitigation from shading or inputs from cooler subsurface flows can lead to inaccurate predictions of downstream temperature flux. Factors influencing flux rates The present analysis showed that it was possible to use site location, air temperature, stream temperature, and stream discharge data to estimate longitudinal rates of downstream temperature flux. Although flux rates displayed evidence of diurnal, seasonal, and yearly variation, the best statistical model was able to accurately capture trends in sub-daily observations across a wide range of stream thermal regimes. Downstream temperature flux rates are dependent upon various factors across differing temporal scales. While the inclusion of both day length and sun altitude angle was important in estimating heat flux at the sub-daily time scale, air temperature has been shown to be an adequate predictor of stream temperature at the weekly time scale (Stefan and Preud’homme, 1993; Mohseni et al., 1998). The results of this analysis, however, indicate that the effects of baseflow reduction can have an impact on stream temperature at the hourly time scale. The analysis of factors influencing downstream temperature flux on a finer scale was important in predicting the response of stream thermal dynamics to baseflow reduction. Although many of the factors used to estimate the effects of baseflow reduction on sub-daily downstream flux rates were critical to the empirical model developed in this analysis, some may be less important when looking at effects of flow reductions on changes 70 in July mean temperatures as the WWAT does. At longer time scales, simpler models may be able to adequately capture enough variability to be useful to natural resource managers. Although many models have been developed attempting to relate air temperatures (Erickson and Stefan, 2000; Caissie et al., 2001; Stefan and Preud’homme, 2003) and stream discharge (Sinokrot and Gulliver, 2000) to stream temperature dynamics, few have analyzed the impacts of groundwater withdrawal on longitudinal fluxes in stream temperature (Risley et al., 2010). In order to do so it was important to understand the influences of baseflow heat energy on the thermal regimes of streams. For this analysis, a simple mass balance method was used to predict the baseflow contributions to each stream reach where the difference in discharge between upstream and downstream gauges during the lowest seven-day flow period corresponded to the average baseflow input to each stream reach. This value was assumed to be representative of constant baseflow input to the stream reach between upstream and downstream locations. The location(s) at which baseflow enters the stream reach can have an important effect on stream thermal dynamics. Peak daytime energy inputs from solar radiation often have the potential to reduce the cooling effects of groundwater inputs (Story et al., 2003). Groundwater inputs would have the greatest impact on reducing downstream warming in shaded stream reaches, or when solar radiation is minimal (i.e., nighttime). The parameters included in Model 10 were significant predictors of longitudinal stream temperature change across many of the streams. Overall, the most influential variable as determined by partial R2 was day length followed closely by sun altitude angle (Table 2.13). Sun altitude angle and day length were used as proxies for estimating the intensity of solar radiation, or sunlight, acting on the stream surface. When focusing only on summer months (July and August; often the warmest months of the year in Michigan) results were similar. Altitude angle 71 accounted for the most variability (Table 2.30) overall (partial R2 = 0.08), and among three of the four thermal classes (baseflow heat energy accounted for the most variability within warm- transitional streams; partial R2 = 0.13; Table 2.31). Similar findings have been reported when examining the influences on energy input into streams to estimate both longitudinal stream temperature changes (LeBlanc et al., 1997; Magnusson et al., 2012) and stream temperatures (Caissie et al., 2007) using deterministic models. Day length was identified as a strong predictor of stream temperatures in a previous study by Risley et al. (2010) to assess the impacts of groundwater pumping on stream temperatures. When examining the factors influencing river heat budgets, Evans et al. (1998) estimated that net shortwave radiation dominated total energy gains, and on average, over 82% of total energy transfers occurred at the air-water interface. Sinokrot and Gulliver (2000) posited that shallow streams are more sensitive to solar heat flux. Thus, reductions in baseflow and overall discharge can have a stronger effect on the downstream temperature flux of small streams. Overland flow, or surface runoff, and its thermal components ranked as having the third most influence on downstream temperature flux in the study streams. Many factors, including increased impervious area (Leopold, 1968; Tong and Chen, 2002) and agricultural/forestry practices (Hidayat et al., 2012), have been shown to influence the amount of surface runoff entering stream channels. Impervious surfaces can store significant amounts of thermal energy that may transfer to streams during runoff events. Overland flow has the potential to warm more quickly to surrounding air temperature when compared to stream discharge. These characteristics of urbanization and surface runoff have been shown to have a direct response on average stream temperatures (Galli, 1990). As many of the streams in this study were located in proximity to agricultural fields and road crossings with the potential to allow surface runoff into the stream 72 channel, it is not surprising that overland flow was a highly ranked influencer of downstream temperature change. Heat energy contributions from baseflow inputs ranked as a moderately influential variable on downstream temperature change (Table 2.13). In regards to daily maximum stream temperatures, groundwater has been suggested to be a prerequisite for daytime cooling in forested stream reaches downstream of forest clearings (Story et al., 2003). However, Story et al. (2003) were able to account for only 40% of downstream cooling due to groundwater. Hyporheic exchange and bed heat conduction were responsible for the remaining 60% of variation in downstream cooling events. These results suggest that, although groundwater inputs are important in offsetting the effects of solar radiation and other warming factors, bed heat conduction and hyporheic flow have been speculated as significant sources of energy exchanges in the moderation of daily temperature extremes (Sinokrot and Stefan, 1993; Poole and Berman, 2001). Air temperature has been used extensively as a predictor for water temperatures (Stefan and Preud’homme, 1993; Pilgrim and Stefan, 1995; Mohseni et al., 1998; Ozaki et al., 2003). In this analysis, the air-water temperature differential proved to be another variable which explained moderate variation in downstream temperature fluxes. The rate of downstream temperature change is influenced by the magnitude of difference between air temperature and water temperature. Ozaki et al. (2003) showed that, while stream temperature was dependent on air temperature, other parameters such as solar radiation are also driving factors, especially in summer months. Upstream discharge and the heat energy associated with it ranked equally among variables responsible for variations in downstream temperature flux. The impact of discharge and 73 thermal inertia on water temperatures occurs primarily through increases in depth of the river (Sinokrot and Gulliver, 2000). This means that at higher discharges, water temperature variation is decreased, and also increases the frequency of lower water temperatures as solar radiation has less of an influence in deeper streams. When stream discharge is reduced, the decreased depth of the river is associated with more influence from solar radiation, and therefore more rapid rates of heating and cooling due to high surface area to volume ratio. Pooling data Model 10 was parameterized using different pooled data sets in order to determine the potential of using hydrologic and meteorological data to predict downstream temperature flux in streams without site-specific data. Results indicate that pooling data based on thermal classes (Table 2.19) has the potential to track seasonal trends in downstream temperature flux (Figure 2.10), although the accuracy of sub-daily observations is highly variable (Table 2.18). The model’s accuracy and ability to track trends is variable within and between thermal classes, but given the improvement over the composite data set, there exists the potential for pooled data to provide reasonable predictions of downstream temperature flux within ungaged streams. Data collection across additional years and streams could improve the model’s predictive power in streams without site-specific data. Compared with the current warming rate module used within the WWAT, Model 10 would provide a more accurate predictive model of downstream temperature change rates. Calibrating empirical models with statewide observational data captures widespread variability in trends in longitudinal stream temperature fluctuation rates, particularly so for the small streams examined in this analysis. The results in Table 2.24 highlight the limitations of the Zorn et al. (2008) model’s ability in estimating thermal flux rates 74 of small Michigan streams. However, future work should compare temperature flux rates across larger streams within these thermal classes. Although pooling data within thermal classes generally provided a better fit than the GLM that used data from all streams combined, some of the thermal classes showed variable responses to driving parameters included in the model. In particular, poor correlation between predicted and observed thermal flux was observed for the pooled model for the cold-transitional thermal class. This may be the result of an overall negative relationship between downstream temperature flux and baseflow heat energy, even though many of the individual cold-transitional streams show a positive relationship between the two variables. Since calibrating models specific to stream thermal class provided an improved fit, collecting data over a broader range of streams within each thermal class would likely parse out the true overall response of thermal class flux rates to atmospheric and hydrologic changes. Additionally, future investigations could group streams based on similarities in basin characteristics known to influence hydrologic patterns which may be able to capture greater variance when estimating downstream warming rates. Similar methods have been developed to estimate exceedance flows of Oregon streams (Cooper, 2002). The same factors which drive hydrologic patterns in basins having similar catchment characteristics may also control the response in downstream temperature flux rates to fluctuations in flow. Effects of baseflow reduction The development of the statistical models served two purposes in the present analysis. First, the robust hydrological dataset and stream gauge network allowed for a determination of the dominant factors affecting downstream temperature flux among each stream thermal regime. When examining parameter influence on downstream temperature flux within thermal classes, 75 results were similar with one exception; in warm-transitional streams baseflow heat load was identified as the most influential variable (Table 2.14). Second, by identifying the model which fits observed data and follows trends best, it was possible to simulate future behavior of stream thermal dynamics following reductions in baseflow. This becomes evident when examining the effects of baseflow reduction on thermal dynamics among the four thermal classifications. Although the downstream temperature flux rates of warm-transitional streams are nearly negligible during baseflow conditions, they become the most affected by each 10% reduction in baseflow. During the most extreme reduction scenario of 90% reduction, warm-transitional streams were predicted to gain 0.69 °C/km on average, much greater than the 0.22 °C/km predicted for cold streams which was the next highest rate. The influence of baseflow heat loads on temperature flux within warm-transitional streams suggests that the thermal inertia provided by baseflow within this thermal class may be an important balancing factor of the heat energy provided from upstream and overland flow. As stream temperatures increase, groundwater- related cooling processes would tend to operate at higher rates (Story et al., 2003). Effects of baseflow reduction on downstream temperature flux rates varied between seasons (Figure 2.11; Table 2.20). As incremental reduction in baseflow was projected to peak downstream warming events during summer months, streams were occasionally projected to experience downstream cooling during the fall, whereby effects became more pronounced at the more extreme reduction scenarios (Figure 2.11; see supplemental files). These projections are reasonable given that surface heat exchange warms or cools streams at a rate inversely proportional to discharge. Similar results were found by Magnusson et al. (2012) when investigating processes influencing stream warming. When investigating the effects of baseflow reduction across thermal classes based on July mean water temperatures, colder streams had a 76 natural tendency to experience greater downstream warming during the 0% reduction simulation when compared to warmer streams. Comparisons with the Zorn et al. (2008) model As mentioned previously, using physical process models to estimate stream temperature flux becomes increasingly difficult during low-flow periods in small stream systems. The Zorn et al. (2008) model used a rate of heating model based on physical relationships and a heat transfer coefficient. The stream temperature change from upstream to downstream is largely dependent on the travel time of water within each stream segment. Additionally, the heating coefficient is inversely proportional to stream depth, so as stream depth decreases, the heating coefficient increases as does estimated stream temperature. Estimating downstream temperature flux rates given average August conditions (often the lowest flow summer month in Michigan) using the Zorn et al. (2008) model resulted in unrealistic values for several streams (Table 2.24). These predictions largely occur from a lack of boundary conditions whereby heating or cooling ceases at a pre-defined condition such as maximum or minimum air temperature. Additionally, streams with very low velocities and discharges are highly susceptible to extreme rates of downstream temperature flux given the heating module implemented in the Zorn et al. (2008) model. An issue with implementing physical process models in the absence of shading factors is that the modules do not consider cooling effects from riparian shading, or hyporheic and groundwater flows that may provide an offset to warming effects of air temperature and solar radiation. When using the Zorn et al. (2008) model to estimate downstream temperature flux in this study, predictions of warming or cooling were directly related to the air – water temperature gradient. Although this was shown to be an important factor in determining downstream temperature flux, it was not the most significant factor (Table 2.13). The difference in which 77 variables primarily contribute to fluctuations in downstream warming rates could have played a major role in the apparent differences predicted by each model in stream responses to baseflow reduction based on July mean temperature (Figures 2.12 – 2.13). Based on my results, streams have been shown to experience downstream warming even when air temperatures were cooler than water temperatures (Table 22; Table 24). This could potentially be the result of two factors: either the temperature of baseflow entering the stream reach was greater than air temperature, or the heat input via solar radiation was greater than the cooling force of the air – water temperature differential and baseflow inputs. Given that scenario analyses for comparisons of Model 10 and the Zorn et al. (2008) model were conducted for August, the lowest mean August air temperatures were 10.34 °C for King Creek and North Branch Thunder Bay River (Table 2.23) located in the NLPUP (Table 2.1), where groundwater temperatures were set at 8.3 °C (Methods; Computed heat energy of various flows). When air temperature is greater than stream temperature, sensible heat fluxes towards the stream, and the opposite is true when air temperature is cooler than water temperature. This means that it is more likely that solar radiation heat flux was the primary cause of downstream warming, as evidenced by partial regression analysis (Table 2.13). Implications and recommendations of findings Resource managers employ a variety of modeling techniques to simulate stream temperature dynamics. Models are useful to investigate influences on thermal dynamics and to examine effects of varying parameters on changes to stream thermal regimes. As the effects of downstream temperature flux vary on a seasonal and yearly basis, sufficient data must be collected on a wide range of streamflow and weather conditions to account for within- and between-year variation in order to properly calibrate models to ensure accurate predictions. 78 Accurate predictions are necessary when conducting scenario analyses so that managers and policy makers can confidently rely on projections when implementing restrictions on water usage rates. As shown in this analysis, model correlation decreased in streams with more than one year of data. Although this was expected, it highlights the importance of multi-year data collection periods as stream thermal dynamics do not respond precisely the same from year to year. When conducting baseflow reduction scenarios using the Zorn et al. (2008) model, it was clear that this physical process model had difficulties in predicting downstream warming rates of small streams, even under 0% reduction scenarios. Many of the streams selected within this analysis had discharges (min = 2e-3 m3s-1; max = 1.75 m3s-1; mean = 0.49 m3s-1; sd = 0.49 m3s-1) similar to those used to develop the Zorn et al. (2008) warming model (min = 0.01 m3s-1; max = 1.93 m3s-1; mean = 0.47 m3s-1; sd = 0.45 m3s-1). Predicted and projected flux rates were often extreme, and could severely restrict water withdrawal in management units in close proximity to these types of stream basins. For such small streams, baseflow reduction itself may be the critical limiting factor before increased temperatures becomes an issue. The downstream warming module of the Zorn et al. (2008) model is used to relate the effects of baseflow reduction on temperature increases in order to predict corresponding changes to the stream fish community on streams and rivers of greater size. Since the findings of this analysis show that the Zorn et al. (2008) model may be overestimating downstream warming rates for Michigan streams of comparable baseflow yield in this study, it is recommended to conduct this type of study on a set of larger streams or rivers that would more clearly reflect the thermal properties of rivers originally chosen to calibrate the physical process model used in the WWAT. There remains a question of how to assess the impacts of baseflow reduction in the absence of site-specific data. The findings of previous research regarding the impacts of shading 79 on temperature flux rates within streams shows that riparian corridors can mitigate the effects of reduced flows and heat flux at the stream surface (Moore et al., 2003; Story et al., 2003). The statistical models developed in this analysis partly account for the effects of shading along the stream reach. This could explain the much lower observed downstream warming rates, and those predicted using Model 10, in comparison with those predicted by the Zorn et al. (2008) model, absent baseflow reduction. Although measures of model accuracy and correlation with observed temperatures were reduced following pooling of data by thermal classes, predictions were much better than those using the Zorn et al. (2008) model; however this is likely due to differences in the mechanics of each model. The accuracy of the model should be further examined when information on overland flow and downstream – upstream discharge differential are not readily available. In addition, the model should be calibrated with data across additional years and streams, although the possibility of estimating temperature flux in the absence of site-specific data remains promising. In the case that streams are not warming at a constant rate, as assumed by the physical process model, implementation of riparian shade corridors may significantly mitigate the warming effects associated with water withdrawal. Clear-cut stream reaches are more vulnerable to climatic impacts such as increased solar heat flux at the stream surface, increased wind speed, and advection of warm air from forest clearings. Although riparian shade buffers can reduce the magnitude of downstream warming, there are still gaps in our knowledge of the extent of stream cooling after flowing through shaded environments (Moore et al., 2005). It has been speculated that downstream cooling under shaded stream reaches still requires inputs of groundwater or hyporheic flows (Beschta et al., 1987; Story et al., 2003). In any case, riparian shade buffers can 80 provide relief from the previously mentioned drivers of stream temperature increases, and could potentially contribute to downstream cooling. This study used model selection criteria to compare a suite of statistical models in order Conclusion to determine the most important processes influencing downstream temperature flux in Michigan streams. The development and comparison of progressively complex models with many of the same parameters highlighted the importance of including variables related to solar heat flux, and the heat energy inputs of different sources of discharge. These variables were important in accounting for sub-daily variation of downstream temperature fluxes across 21 Michigan streams with differing thermal properties. This study provides new data on effects of discharge and air and water temperatures on downstream temperature flux rates in small Michigan streams. Comparisons of the best statistical model developed in this analysis with those of the Zorn et al. (2008) model showed substantial differences in the estimations and influencing factors of downstream temperature flux. The Zorn et al. (2008) physical process model had difficulties predicting downstream temperature flux rates in small streams, which likely led to inaccurate projections of flux rates following baseflow reduction. This was likely the result of this model’s failure to account for influences from riparian shading and intermittent groundwater inputs along the stream reach, which can act as an offset to heating factors. Using empirical relationships, the statistical models developed in this study were able to inherently capture these mitigating factors which physical process models may not be able to fully capture. Furthermore, by including a parameter which estimated the influence of baseflow inputs on downstream temperature flux rates, it was possible to conduct scenario analyses within individual stream reaches. Baseflow reduction scenarios allowed for comparisons of the impacts 81 of groundwater withdrawal on the thermal dynamics of a wide range of streams. Results of groundwater abstraction simulations presented a wide range of projections including increased downstream warming and cooling under extreme reduction scenarios. These results indicate that shallow streams such as those included in this analysis are vulnerable to sources of heating and cooling such that flux rates of streams with shallow depths show a rapid response to elements of heating and cooling. Research presented here shows the importance of understanding driving factors of stream temperature change, and stream thermal response to changes in flow regime. Although the present study showed the strong influence of several factors on stream temperature change, previous research shows that there are other variables critical in driving changes in stream temperature. Further examination of drivers of stream temperature change would be useful in improving the predictive range of temperature flux, which would subsequently provide a stronger projection of stream thermal response to baseflow reduction. However, increased collection of hydrological data on small streams would refine the best thermal flux model developed within this analysis through calibration, and improve predictive power across a robust range of stream types. 82 APPENDICES 83 APPENDIX 1.0 Supporting tables and figures for Chapter 1: Application of benchmark detection methods to identify thermal thresholds of Table 1.1: Table of temperature thresholds (°F) of species with upper thermal thresholds identified by each of the three analytical threshold detection methods, as well as visually estimated threshold interpreted from LOESS line. stream fishes along a thermal gradient. WWAT upper threshold 69.3 69.1 68.7 67.8 69.4 71.8 69.0 71.2 72.4 70.8 71.2 73.0 LOESS 65.0 68.0 65.0 68.0 70.0 Not apparent 65.0 69.0 Not apparent 68.0 69.0 Not apparent Species Brook Trout Coho Salmon Slimy Sculpin Chinook Salmon Brown Trout Northern Redbelly Dace Rainbow Trout Mottled Sculpin Northern Brook Lamprey Longnose Dace Blacknose Dace Burbot WWAT Optimum TITAN CART 62.3 62.9 63.1 63.6 64.1 64.1 64.2 64.5 65.3 65.4 65.9 66.2 56.8 60.5 65.4 70.6 64.7 67.7 68.8 68.8 66.7 68.8 68.8 68.3 64.9 67.7 65.4 68.5 65.9 67.7 68.3 68.8 70.0 68.9 71.0 72.5 84 Table 1.2: Table of differences between detected thresholds of each method. Species TITAN – CART TITAN – WWAT CART – WWAT Brook Trout Blacknose Dace Brown Trout Burbot Chinook Salmon Coho Salmon Longnose Dace Mottled Sculpin Northern Brook Lamprey Northern Redbelly Dace Rainbow Trout Slimy Sculpin Average difference 8.2 2.3 1.2 4.3 -2.1 7.2 0.1 0.1 3.7 0.0 -0.5 0.1 1.9 -4.4 -0.2 -3.5 -0.5 0.7 -1.4 -1.9 -2.4 -2.1 -4.1 -0.7 -3.3 -2.0 -12.6 -2.5 -4.7 -4.7 2.8 -8.6 -2.0 -2.5 -5.7 -4.1 -0.3 -3.3 -3.8 85 Table 1.3: Information of the upper 20% of abundances of each decreaser species including, number of optimal sites (N), WWAT optimum temperature (mean), minimum (Lower) and maximum (Upper) July mean water temperatures for optimum temperature calculation, standard deviation, and predicted WWAT threshold (temperatures in °F). Species Brook Trout Black Bullhead Black Crappie Blacknose Dace Blackside Darter Bluegill Bluntnose Minnow Brown Trout Brook Stickleback Burbot Central Mudminnow Central Stoneroller Channel Catfish Chinook Salmon Coho Salmon Common Carp Common Shiner Creek Chub Fathead Minnow Golden Redhorse Golden Shiner Grass Pickerel Green Sunfish Greenside Darter N 128 21 24 118 74 69 69 152 37 32 128 29 10 12 17 38 101 174 10 23 14 37 88 12 Optimum Temperature Lower Upper Std. WWAT TITAN CART 69.1 73.9 76.4 71.0 75.2 76.4 77.2 70.9 72.9 72.4 74.5 76.6 75.1 67.5 67.5 76.6 76.1 74.8 74.5 76.1 73.7 73.5 77.2 77.2 4.0 2.6 2.5 3.1 2.5 2.9 2.8 3.0 4.2 3.9 3.4 3.0 0.9 2.4 3.5 2.0 2.8 3.2 3.1 2.1 5.0 3.0 2.9 2.8 69.3 74.1 76.5 71.2 74.7 75.3 76.8 69.4 73.2 73.0 72.4 75.0 74.7 67.7 69.1 76.6 74.5 73.0 74.7 77.1 76.1 74.0 75.6 78.8 64.9 67.6 71.5 71.0 66.4 66.4 69.2 65.9 56.7 72.5 62.6 67.4 71.9 68.5 67.7 71.0 67.3 63.8 68.6 70.7 71.8 65.5 68.2 73.9 56.8 71.4 75.0 68.8 68.8 68.8 76.3 64.7 71.2 68.3 70.6 68.8 71.9 70.6 60.5 72.4 69.4 67.7 74.5 75.9 71.9 72.4 76.3 75.9 62.3 69.6 72.2 65.9 70.4 70.2 71.9 64.1 65.9 66.2 66.5 69.7 73.2 63.6 62.9 73.2 69.7 67.4 69.3 73.4 67.4 68.8 70.5 73.9 50.9 63.7 67.1 59.5 63.4 61.1 61.6 53.8 50.9 58.7 56.6 61.7 71.9 59.8 56.9 67.3 63.5 56.7 62.7 67.2 58.3 61.1 63.7 67.8 86 Table 1.3: (cont’d) Hornyhead Chub Johnny Darter Largemouth Bass Log Perch Longnose Dace Mottled Sculpin Northern Brook Lamprey Northern Hog Sucker Northern Redbelly Dace Northern Pike Pumpkinseed Rainbow Darter Rainbow Trout Redhorse Sucker Rock Bass Rosyface Shiner Shorthead Redhorse Slimy Sculpin Smallmouth Bass Spotfin Shiner Stonecat Walleye White Sucker Yellow Bullhead Yellow Perch 43 130 67 31 42 105 10 64 18 49 64 48 87 13 84 18 12 10 48 15 29 12 190 34 37 68.5 68.6 69.3 69.3 65.4 64.5 65.3 71.6 64.1 71.3 71.5 69.2 64.2 74.0 72.1 71.7 71.8 63.1 73.2 73.5 72.9 72.8 68.8 71.5 67.6 76.6 76.4 76.4 75.7 72.4 71.1 69.4 77.2 72.9 76.6 76.6 74.7 71.5 75.8 77.2 76.1 75.1 69.9 77.2 76.6 77.2 75.8 77.2 76.4 73.7 3.3 2.9 3.3 4.0 3.1 3.9 4.0 3.0 4.4 3.2 3.5 2.6 2.8 1.1 3.0 2.2 2.3 3.2 2.1 2.0 2.6 2.5 3.2 3.0 4.4 74.3 73.7 75.1 76.2 70.8 71.2 72.4 76.8 71.8 76.8 77.6 73.8 69.0 75.9 77.3 75.5 75.9 68.7 76.9 77.0 77.4 77.1 74.4 76.7 75.3 69.2 66.5 68.0 66.1 68.9 68.8 70.3 69.6 67.7 70.4 71.4 67.7 68.3 73.2 68.3 69.3 71.7 65.4 70.5 71.5 70.2 72.4 64.5 70.4 68.6 66.0 68.0 69.4 72.8 68.8 68.8 66.7 76.3 67.7 71.4 71.6 68.1 68.8 74.3 73.4 75.9 71.5 65.4 76.3 76.3 73.8 74.8 69.9 71.3 71.5 60.3 61.1 63.4 60.3 58.7 55.8 59.7 64.2 56.6 61.1 61.0 61.3 58.8 71.9 62.3 68.4 67.2 59.1 69.1 68.5 67.6 67.7 60.4 63.2 58.3 87 Table 1.4: Detected benchmarks of all 49 species of each of the three datasets (ALL, NLPUP, and SLP). NA indicates that no benchmark was estimated due to lack of minimum data requirements. Species Brook Trout Black Bullhead Black Crappie Blacknose Dace Blackside Darter Bluegill Bluntnose Minnow Brown Trout Brook Stickleback Burbot Central Mudminnow Central Stoneroller Channel Catfish Chinook Salmon Coho Salmon Common Carp Common Shiner Creek Chub Fathead Minnow Golden Redhorse Golden Shiner Grass Pickerel Green Sunfish Greenside Darter Hornyhead Chub TITAN ALL NLPUP 64.9 67.6 71.5 71.0 66.4 66.4 69.2 65.9 56.7 72.5 62.6 67.4 71.9 68.5 67.7 71.0 67.3 63.8 68.6 70.7 71.8 65.5 68.2 73.9 69.2 62.4 68.9 70.7 62.9 65.3 70.4 68.3 69.6 56.6 70.3 62.7 66.7 NA NA 67.5 66.3 64.7 63.9 65.2 70.5 72.4 NA 67.3 NA 67.0 CART SLP ALL NLPUP 66.3 71.9 74.9 71.0 69.4 69.4 69.5 69.3 71.6 NA 72.5 68.7 72.8 NA 61.1 70.5 69.2 71.4 68.6 71.5 70.5 73.6 69.4 74.3 69.4 56.8 69.2 71.2 66.4 69.1 69.1 70.5 69.4 72.8 71.8 64.8 66.8 NA 65.4 60.5 67.3 70.7 62.9 62.7 71.2 67.6 NA 71.8 NA 65.9 56.8 71.4 75.0 68.8 68.8 68.8 76.3 64.7 71.2 68.3 70.6 68.8 71.9 70.6 60.5 72.4 69.4 67.7 74.5 75.9 71.9 72.4 76.3 75.9 66.0 88 WWAT SLP ALL NLPUP SLP 59.4 67.3 74.6 71.4 76.6 75.0 71.1 72.2 75.0 68.8 75.1 69.4 76.3 76.6 70.9 64.4 72.6 71.2 69.6 NA 73.2 70.6 74.5 68.8 75.8 71.9 NA NA NA 60.5 76.4 76.3 75.9 74.8 72.5 71.3 NA 74.5 75.9 76.9 75.8 71.9 73.7 72.4 76.3 75.7 79.0 75.9 68.4 74.8 69.4 NA NA 69.0 73.5 72.6 71.7 68.3 71.6 73.4 70.1 NA NA 67.7 69.1 NA 73.0 71.3 NA NA NA NA 73.6 NA 71.0 69.3 74.1 76.5 71.2 74.7 75.3 76.8 69.4 73.2 73.0 72.4 75.0 74.7 67.7 69.1 76.6 74.5 73.0 74.7 77.1 76.1 74.0 75.6 78.8 74.3 Table 1.4: (cont’d) Johnny Darter Largemouth Bass Log Perch Longnose Dace Mottled Sculpin Northern Brook Lamprey Northern Hog Sucker Norther Redbelly Dace Northern Pike Pumpkinseed Rainbow Darter Rainbow Trout Redhorse Sucker Rock Bass Rosyface Shiner Shorthead Redhorse Slimy Sculpin Smallmouth Bass Spotfin Shiner Stonecat Walleye White Sucker Yellow Bullhead Yellow Perch 66.5 68.0 66.1 68.9 68.8 70.3 69.6 67.7 70.4 71.4 67.7 68.3 73.2 68.3 69.3 71.7 65.4 70.5 71.5 70.2 72.4 64.5 70.4 68.6 64.7 70.4 66.9 64.6 70.3 67.0 66.6 70.5 70.7 70.7 67.4 58.7 67.3 67.2 70.4 70.5 NA 67.0 68.1 67.5 70.5 64.2 69.5 68.1 67.8 68.4 72.2 68.4 71.0 NA 71.1 NA 71.1 69.4 68.1 63.2 73.1 70.5 69.4 72.2 NA 71.1 71.4 70.7 74.8 65.7 71.1 70.4 68.0 69.4 72.8 68.8 68.8 66.7 76.3 67.7 71.4 71.6 68.1 68.8 74.3 73.4 75.9 71.5 65.4 76.3 76.3 73.8 74.8 69.9 71.3 71.5 89 65.8 70.4 68.5 64.6 68.6 67.1 69.6 64.6 71.8 65.1 69.6 68.2 66.5 69.1 70.7 71.2 66.1 70.4 68.2 68.2 72.4 59.4 70.4 68.6 68.0 69.4 72.8 68.4 68.8 NA 76.3 NA 71.4 71.6 68.1 63.2 74.3 73.4 75.9 72.8 NA 76.3 76.3 73.8 74.8 69.9 71.3 71.5 73.7 75.1 76.2 70.8 71.2 72.4 76.8 71.8 76.8 77.6 73.8 69.0 75.9 77.3 75.5 75.9 68.7 76.9 77.0 77.4 77.1 74.4 76.7 75.3 71.9 71.1 72.9 70.7 69.2 70.7 72.2 71.7 75.7 72.4 72.6 68.4 NA 74.8 NA NA 68.3 73.0 NA NA NA 71.2 NA 72.7 73.3 75.3 76.3 NA 72.0 NA 77.1 NA 76.1 76.5 74.4 71.2 NA 77.2 76.1 75.1 NA 77.4 76.5 77.4 77.3 74.4 76.5 75.2 Table 1.5: Table of detected benchmarks of species which were classified as decreasers in both the NLPUP and SLP subregions. Average differences were calculated between the two subregions for each method. Species/Region Brook Trout Brown Trout Mottled Sculpin TITAN CART WWAT NLPUP SLP NLPUP SLP NLPUP SLP 67.3 70.9 72.0 59.4 64.4 68.8 62.4 69.6 70.3 66.3 69.3 71.0 56.8 69.4 68.6 69.4 68.3 69.2 Average difference -1.5 0.7 -1.1 90 Table 1.6: Purity, reliability, and directionality of threshold response for those species identified as decreasers in the ALL dataset. P- values are significant at P  0.05. Purity (Pur) is the proportion of responses detected by TITAN among bootstrap replicates that agree with the observed response. Reliability (Rel) is an estimate using the proportion of responses using bootstrap replicates whose IndVal scores that correspond to a probability level of P  0.05. Directionality of response indicates whether species show a negative (decreaser; -) or positive (increaser; +) threshold response along the thermal gradient. ALL NLPUP Pur P P Pur +/- Rel - 0.004 1.000 1.000 0.004 0.998 0.998 + 0.004 1.000 1.000 - 0.024 0.578 0.436 + 0.004 0.948 0.904 NA NA 0.004 1.000 1.000 - 0.004 0.998 0.998 + 0.004 1.000 1.000 - 0.016 0.894 0.888 + 0.008 0.986 0.962 + 0.004 1.000 1.000 + 0.004 1.000 1.000 NA NA P Rel 0.004 1.000 1.000 0.004 1.000 1.000 0.008 0.828 0.804 0.004 1.000 1.000 NA NA NA NA NA +/- - - - NA - - - 0.028 0.936 0.904 0.004 1.000 1.000 0.048 0.702 0.64 0.004 0.966 0.952 NA NA 0.024 0.922 0.828 NA NA 0.016 0.658 0.658 - NA NA NA NA SLP Pur Rel 0.004 1.000 1.000 0.004 1.000 1.000 0.004 1.000 1.000 NA NA NA NA 0.004 1.000 0.98 0.008 0.994 0.916 0.004 1.000 1.000 NA NA 0.004 1.000 1.000 NA NA NA NA Species Brook Trout Blacknose Dace Brown Trout Burbot Chinook Salmon Coho Salmon Longnose Dace Mottled Sculpin Northern Brook Lamprey Northern Red Dace Rainbow Trout Slimy Sculpin +/- - - - - - - - - - - - - 91 Table 1.7: Table of means of detected thresholds using logistic function to simulate abundances (10 simulations) with a known threshold at 68 °F. Method/Region TITAN CART WWAT Average ALL 69.30 67.80 67.3.0 68.13 NLPUP 68.00 67.40 67.00 67.50 SLP 69.40 68.00 68.00 68.40 Average 68.92 67.71 67.41 68.01 92 Table 1.8: Two-way ANOVA (α = 0.05) examining differences between region and method of detected benchmarks using logistic function to simulate abundances with a known threshold at 68 °F. Method Region Method:Region P-value 2.00e-16 3.70e-16 1.49e-4 Significant? Yes Yes Yes 93 Figure 1.1: LOESS regression of logistically simulated data with a known threshold of 68 °F. 94 Figure 1.2: Distribution of Brook Trout abundances along a July mean water temperature thermal gradient. Vertical lines correspond to identified benchmarks of TITAN, CART, and WWAT benchmark detection methods. The LOESS regression is depicted by the solid red regression line and is used in visual identification of inflection points. 95 Figure 1.3: Linear regression of TITAN and CART identified thresholds. Dashed line represents 1:1 relationship. Axis units are July mean water temperature (°F). 96 Figure 1.4: Linear regression of TITAN and WWAT identified thresholds. Dashed line represents 1:1 relationship. Axis units are July mean water temperature (°F). 97 Figure 1.5: Linear regression of CART and WWAT identified thresholds. Dashed line represents 1:1 relationship. Axis units are July mean water temperature (°F). 98 Figure 1.6: Linear regression of TITAN identified thresholds for NLPUP and SLP regions. Dashed line represents 1:1 relationship. Axis units are July mean water temperature (°F). 99 Figure 1.7: Linear regression of CART identified thresholds for NLPUP and SLP regions. Dashed line represents 1:1 relationship. Axis units are July mean water temperature (°F). 100 Figure 1.8: Linear regression of WWAT identified thresholds for NLPUP and SLP regions. Dashed line represents 1:1 relationship. Axis units are July mean water temperature (°F). 101 Figure 1.9: Box plot showing results of threshold detection methods applied to the July mean temperature gradients of the NLPUP and SLP regions with abundance data simulated using a logistic function. 102 APPENDIX 2.0 Supporting tables and figures for Chapter 2: Quantifying downstream warming rates of Michigan streams: Improving an important link in Michigan’s Water Withdrawal Assessment Tool. Table 2.1: Site information specific to each stream segment. Streams marked with an asterisk are those for which data was collected over both the 2015 and 2016 field seasons. Up and down refer to upstream and downstream locations. Stream Region Thermal Up Up Class Latitude Longitude Down Latitude Down Longitude Reach Length (m) Pokagon Creek* SLP Fish Creek* SLP C C 41.89517 -86.162632 41.915803 -86.175679 4050 43.245992 -84.964747 43.242022 -84.915223 5186 SLP CT 42.932887 -86.081828 42.91636 -86.146075 6550 Up Down Elevation Elevation (m) 224 240 186 (m) 219 234 180 Pigeon River* Unnamed Gun River Tributary Nottawa Creek* Hemingway Lake Outlet Honeyoey Creek Middle Branch Tobacco River Hasler Creek Prairie River* SLP CT 42.537894 -85.593867 42.530494 -85.562968 3131 231 226 SLP WT 42.192564 -85.060415 42.195998 -85.104618 3758 SLP WT 43.32678 -85.124515 43.330136 -85.154513 3197 SLP WT 43.433623 -84.701648 43.379136 -84.705982 6638 279 277 230 276 273 224 SLP WT 43.909194 -84.697312 43.929905 -84.666327 4091 258 249 SLP SLP W W 43.042332 -83.423206 43.083594 -83.442947 7586 41.801832 -85.116614 41.832568 -85.165065 5863 250 287 233 284 103 W C C C 41.90477 -85.297885 41.921249 -85.312047 2539 44.375846 -85.972647 44.369588 -85.999598 2551 44.956875 -85.132748 44.968664 -85.138993 1454 45.070651 -84.283728 45.089439 -84.284929 2879 263 232 232 277 359 223 262 225 214 272 352 211 Table 2.1: (cont’d) Swan Creek SLP NLP NLP NLP Cedar Creek Cedar River* East Branch Black River Butterfield Creek* NLP CT 44.273249 -85.094087 44.256377 -85.03362 5978 King Creek NLP WT 45.018848 -83.650705 45.047993 -83.634655 5822 North Branch Thunder Bay River Morgan Creek* Slapneck Creek Spring Creek* Carp River* Middle Branch Escanaba River East Branch Waiska River Squaw Creek NLP W 45.179007 -83.923148 45.191635 -83.891476 4630 241 237 UP UP UP UP C 46.519698 -87.504502 46.521351 -87.494782 1106 CT 46.354843 -86.946771 46.350637 -86.928918 1564 CT 46.512909 -90.156133 46.513418 -90.177011 1681 CT 46.509131 -87.418924 46.510534 -87.388497 2614 368 252 360 237 366 249 358 189 UP WT 46.420206 -87.797962 46.398398 -87.770883 6131 432 426 UP W 46.418818 -84.474418 46.406065 -84.499266 3623 186 185 UP W 46.057035 -87.18974 45.985396 -87.140559 1676 283 249 104 Table 2.2: Model numbers and the parameters included within each model denoted with an X. Ta = air temperature (°C); Tw = water temperature (°C); Qup = upstream discharge (cms); Qdown = downstream discharge; S = day length (hr); ΔTflow = cumulative heat energy (°C); ΔTup = upstream heat energy (°C); ΔTbase = baseflow heat energy (°C); ΔTover = overland flow heat energy (°C); α = sun altitude angle (°) Parameter Intercept Ta - Tw S α Qup Qdown - Qup Qdown/Qup ΔTflow ΔTup ΔTbase ΔTover (1/Qup)(Tw + 273.16)4 (1/Qup)[(Tw + 273.16)4 - (Ta + 273.16)4] (1/Qup)(α) X X X X M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X 105 Table 2.3: Comparison of a priori thermal classes and a posteriori thermal classes based upon July mean water temperatures. For those streams with two years of data (denoted by *), July mean water temperatures and overall discharges were averaged across both years. C =  17.5 °C; CT = > 17.5 °C and  19.5 °C; WT = > 19.5 °C and  21.0 °C; W = > 21.0 °C. Stream a priori a posteriori July Mean Thermal Class Thermal Class Temperature (°C) Discharge (m3 s-1) Cedar Creek Cedar River* East Branch Black River Pokagon Creek* Fish Creek* Morgan Creek* Pigeon River* Butterfield Creek* Slapneck Creek Carp River* Spring Creek* Middle Branch Tobacco River Middle Branch Escanaba River Honeyoey Creek King Creek Nottawa Creek* Prairie River* Squaw Creek Swan Creek Hasler Creek North Branch Thunder Bay River C C C C C C CT CT CT CT CT WT WT WT WT WT W W W W W C C C CT CT WT CT CT CT CT WT CT CT WT WT W CT CT WT W W 15.7 14.5 16.5 19.1 19.5 20.8 18.7 18.8 18.4 19.2 20.1 18.5 18.6 20.4 20.1 23.1 18.4 18.7 20.9 23.0 21.8 0.12 1.34 0.83 0.42 0.86 0.11 0.48 0.09 0.33 1.79 0.12 0.52 0.99 0.06 0.01 0.55 0.21 0.04 0.38 0.04 0.43 106 SSE 9473 8375 112 16510 2948 42096 588 110782 11788 56711 4.59e+6 -0.60 0.02 0.35 0.00 -0.06 -0.02 3.30e+11 -0.16 0.63 0.26 -0.07 0.07 0.03 0.01 0.00 0.01 -0.11 -0.24 0.09 0.09 0.11 0.50 1.43e+18 7.32e+8 30168 215297 85116 758770 2790 5820 3390 Table 2.4: Goodness of fit components for individual stream reaches using the forcing model (Model 12). r Name Cedar Creek Cedar River East Branch Black River Pigeon River Pokagon Creek Prairie River Fish Creek Middle Branch Tobacco River Butterfield Creek Slapneck Creek Squaw Creek Middle Branch Escanaba River Carp River Swan Creek Honeyoey Creek King Creek Morgan Creek Spring Creek Nottawa Creek Hasler Creek North Branch Thunder Bay River 107 Table 2.5: Models and their goodness of fit components. Results are averages over each stream reach. K = number of parameters; SSE = sums of squared errors; r = correlation between predicted and observed temperature change; L = log likelihood component of AIC; AIC = Akaike Information Criteria; ωi = Akaike weight. Count is the total number of streams for which each model was identified as providing the best fit. Models with the smallest SSE, L, or AIC are best fitting. Model no. K Ave. SSE Ave. r Ave. L Ave. AIC Ave. ω Count 1 2 3 4 5 6 7 8 9 10 11 2 3 4 5 6 6 7 8 8 9 7 225 202 196 183 173 173 168 161 164 156 184 0.21 0.32 0.40 0.45 0.54 0.50 0.58 0.59 0.56 0.62 0.50 -10875 -10661 -10499 -10368 -10128 -10224 -9971 -9911 -9990 -9745 -10130 21753 21328 21006 20746 20264 20461 19957 19838 19996 19508 20274 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.001 0.024 0.737 0.238 0 0 0 0 0 0 0 0 0 16 5 108 Table 2.6: Model weights for each stream and the total count for which each model provided the best fit. Name Cedar Creek Cedar River East Branch Black River Pigeon River Pokagon Creek Prairie River Fish Creek Middle Branch Tobacco River Butterfield Creek Slapneck Creek Squaw Creek Middle Branch Escanaba River Thermal Class M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 C C C CT CT CT CT CT CT CT CT CT 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 1.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 1.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 1.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 1.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 1.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 1.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.026 0.974 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.002 0.998 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 1.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 1.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 1.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 1.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 1.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 1.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 1.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 1.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.016 0.479 0.505 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 1.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.002 0.998 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 1.000 Carp River Swan Creek King Creek CT WT Honeyoey Creek WT WT WT WT W W Morgan Creek Spring Creek Nottawa Creek Hasler Creek North Branch Thunder Bay W 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 1.000 0.000 River Count 0 0 0 0 0 0 0 0 0 16 5 109 Table 2.7: Models and average r for streams with one and two years of data. Years of data Model 10 Model 11 1 2 0.66 0.58 0.55 0.44 Model 9 0.61 0.52 Model 8 0.62 0.55 110 Table 2.8: Overall root-mean-square error (RMSE) (°C) between observed and predicted values of downstream temperature flux for each of the four best models (M8, M9, M10, and M11) across each of the 21 streams over 2015 and 2016. RMSE Stream Cedar Creek Cedar River East Branch Black River Middle Branch Tobacco River Prairie River Fish Creek Pigeon River Pokagon Creek Butterfield Creek Slapneck Creek M10 M11 M9 M8 0.07 0.07 0.07 0.07 0.19 0.22 0.21 0.20 0.09 0.09 0.09 0.09 0.07 0.07 0.07 0.07 0.08 0.10 0.10 0.09 0.15 0.17 0.15 0.15 0.04 0.04 0.04 0.05 0.07 0.07 0.07 0.10 0.12 0.13 0.13 0.12 0.22 0.25 0.22 0.22 0.46 0.47 0.46 0.46 0.07 0.13 0.11 0.08 0.06 0.07 0.06 0.06 0.38 0.42 0.39 0.40 0.12 0.13 0.12 0.12 0.10 0.11 0.10 0.10 0.34 0.37 0.34 0.34 0.30 0.30 0.30 0.31 0.08 0.09 0.08 0.08 0.07 0.07 0.07 0.08 North Branch Thunder Bay River 0.08 0.07 0.08 0.08 Morgan Creek Spring Creek Nottawa Creek Hasler Creek Honeyoey Creek King Creek Middle Branch Escanaba River Squaw Creek Carp River Swan Creek Average 0.12 0.13 0.12 0.12 111 Table 2.9: Compared performances of the four best models using root-mean-square error (RMSE) (°C) between observed and predicted values of downstream temperature flux across each of the 21 streams. Model May June July August September October November 2015 2016 0.10 0.12 M10 0.11 0.13 M11 M9 0.11 0.12 0.10 0.12 M8 0.14 0.12 0.13 0.15 0.14 0.14 0.15 0.13 0.13 0.15 0.13 0.13 0.10 0.11 0.11 0.10 0.11 0.12 0.12 0.12 0.11 0.11 0.11 0.11 0.09 0.11 0.10 0.09 112 Table 2.10: Overall root-mean-square error (RMSE) (°C) between observed and predicted values of downstream temperature flux for each of the four best models (M8, M9, M10 and M11) across each of the four thermal classes. RMSE Thermal Class M10 M11 M9 M8 0.06 0.12 0.16 0.10 0.06 0.12 0.16 0.10 0.06 0.13 0.17 0.10 C CT WT W 0.06 0.11 0.16 0.10 113 Table 2.11: Maximum and minimum values of downstream temperature flux rates (°C/km) for observed and predicted values of each of the four best models. Stream Cedar Creek Cedar River East Branch Black River Pigeon River Pokagon Creek Prairie River Fish Creek Middle Branch Tobacco River Butterfield Creek Slapneck Creek Squaw Creek Middle Branch Escanaba River Carp River Swan Creek Honeyoey Creek King Creek Morgan Creek Spring Creek Nottawa Creek Hasler Creek North Branch Thunder Bay River Observed Model 10 Model 11 Model 9 Model 8 Max Min 0.14 -0.19 -0.41 0.66 Max 0.09 0.47 Min -0.07 -0.15 0.14 -0.10 0.09 0.00 0.04 0.19 0.08 0.06 -0.06 -0.12 -0.08 -0.03 0.01 0.05 0.06 0.03 -0.03 -0.01 -0.05 -0.01 Max 0.11 0.30 0.09 0.01 0.04 0.04 0.07 Min -0.07 0.03 0.01 -0.03 -0.02 -0.07 -0.01 0.14 -0.05 0.11 -0.03 0.13 -0.01 0.04 0.65 1.74 -0.13 -0.51 -0.53 0.07 0.56 0.81 -0.10 -0.19 -0.29 0.05 1.50 2.08 -0.07 -0.11 -0.98 0.06 -0.06 0.04 -0.05 0.01 -0.04 0.14 0.52 0.09 0.07 1.51 0.52 0.03 0.02 -0.09 -0.35 -0.04 -0.07 -1.43 -0.86 -0.13 -0.08 0.06 0.30 0.09 0.04 0.41 0.17 -0.01 0.01 -0.02 -0.22 -0.03 -0.02 -0.68 -0.20 -0.10 -0.06 0.08 0.21 0.10 0.11 0.39 0.23 -0.03 0.02 -0.03 -0.10 -0.02 -0.10 -1.14 -0.21 -0.11 -0.06 Max 0.08 0.32 0.09 0.01 0.04 0.05 0.03 0.11 0.05 0.51 0.80 0.02 0.06 0.25 0.09 0.04 0.42 0.18 -0.02 0.01 Min -0.09 -0.15 Max 0.10 0.38 Min -0.07 -0.12 -0.01 0.07 0.00 -0.03 -0.02 -0.04 -0.01 0.01 0.05 0.06 0.03 -0.02 -0.01 -0.05 -0.02 -0.03 0.10 -0.04 -0.08 -0.17 -0.29 0.07 0.58 0.80 -0.09 0.09 -0.27 -0.05 0.03 -0.04 -0.02 -0.17 -0.03 -0.03 -0.70 -0.20 -0.10 -0.06 0.06 0.28 0.09 0.04 0.39 0.17 -0.01 0.02 -0.02 -0.16 -0.03 -0.03 -0.68 -0.20 -0.09 -0.05 0.09 -0.05 0.06 -0.03 0.09 -0.02 0.05 -0.03 0.07 -0.03 Average 0.17 -0.14 0.10 -0.07 0.14 -0.08 0.09 -0.06 0.09 -0.06 114 Table 2.12: Residuals of observed downstream temperature fluxes and model predictions of the 1st and 99th percentiles. Stream Cedar Creek Cedar River East Branch Black River Pigeon River Pokagon Creek Prairie River Fish Creek Middle Branch Tobacco River Butterfield Creek Slapneck Creek Squaw Creek Middle Branch Escanaba River Carp River Swan Creek Honeyoey Creek King Creek Morgan Creek Spring Creek Nottawa Creek Hasler Creek North Branch Thunder Bay River Average Model 10 Model 11 Model 9 Model 8 Observed - Predicted P01 -0.07 -0.39 -0.09 -0.03 -0.07 -0.05 -0.02 -0.04 -0.06 -0.41 -0.65 -0.01 -0.06 -0.26 -0.04 -0.06 -0.77 -0.61 -0.36 -0.02 -0.05 -0.11 P99 0.05 0.30 0.07 0.03 0.09 0.05 0.02 0.04 0.01 0.28 0.93 0.03 0.06 0.35 0.01 0.03 0.87 0.47 0.11 0.01 0.04 0.09 P01 -0.07 -0.40 -0.09 -0.03 -0.07 -0.06 -0.02 -0.04 -0.06 -0.49 -0.73 -0.04 -0.06 -0.31 -0.04 -0.05 -0.86 -0.60 -0.36 -0.02 -0.04 -0.12 P01 0.05 0.38 0.07 0.03 0.09 0.07 0.01 0.03 0.02 0.28 0.77 0.05 0.07 0.38 0.02 0.04 0.91 0.44 0.12 0.01 0.03 0.10 P01 -0.07 -0.42 -0.09 -0.03 -0.07 -0.05 -0.02 -0.04 -0.06 -0.44 -0.65 -0.03 -0.06 -0.28 -0.04 -0.06 -0.77 -0.61 -0.36 -0.02 -0.04 -0.11 P99 0.05 0.34 0.07 0.03 0.09 0.05 0.02 0.04 0.01 0.29 0.93 0.04 0.06 0.36 0.01 0.04 0.87 0.47 0.11 0.02 0.04 0.10 P01 -0.07 -0.40 -0.10 -0.03 -0.07 -0.05 -0.02 -0.04 -0.10 -0.40 -0.65 -0.02 -0.06 -0.28 -0.04 -0.06 -0.78 -0.60 -0.36 -0.02 -0.05 -0.11 P01 0.05 0.33 0.07 0.03 0.09 0.05 0.02 0.05 0.03 0.29 0.93 0.03 0.07 0.37 0.02 0.03 0.86 0.48 0.11 0.01 0.04 0.10 115 Table 2.13: Partial R2 values of each variable included in Model 10. Values reflect the influence of each variable on downstream temperature flux rates. Ta = air temperature (°C); Tw = water temperature (°C); Qup = upstream discharge (cms); Qdown = downstream discharge; S = day length (hr); ΔTup = upstream heat energy (°C); ΔTbase = baseflow heat energy (°C); ΔTover = overland flow heat energy (°C); α = sun altitude angle (°). Variable Ta - Tw Qup Qdown - Qup Overall 0.04 0.02 0.03 Partial R2 S 0.09 ΔTup ΔTbase ΔToverland 0.03 0.05 0.04 α 0.07 Table 2.14: Partial R2 values averaged over thermal class. Ta = air temperature (°C); Tw = water temperature (°C); Qup = upstream discharge (cms); Qdown = downstream discharge; S = day length (hr); ΔTup = upstream heat energy (°C); ΔTbase = baseflow heat energy (°C); ΔTover = overland flow heat energy (°C); α = sun altitude angle (°). Partial Regression S 0.14 0.12 0.04 0.04 ΔTup ΔTbase ΔToverland 0.01 0.04 0.02 0.00 0.02 0.05 0.07 0.02 0.01 0.04 0.09 0.02 α 0.06 0.10 0.02 0.08 Thermal Class Ta - Tw Qup Qdown - Qup 0.02 0.04 0.08 0.01 0.01 0.05 0.01 0.02 0.00 0.02 0.05 0.01 C CT WT W 116 Table 2.15: Type III sums of squared errors for a model testing only for thermal class as a fixed effect on downstream temperature flux rates using GLM. Source DF Type III SS Mean Square F Value Pr > F 1004.89 <.0001 3 1 2259.62 <.0001 54.37 1 <.0001 287.89 <.0001 1 1 1377.65 <.0001 800.70 <.0001 1 347.83 <.0001 1 189.79 <.0001 1 1 1480.57 <.0001 3.79 20.06 96.01 55.80 24.24 13.23 103.18 3.79 20.06 96.01 55.80 24.24 13.23 103.18 Thermal Class Ta - Tw Qup Qdown - Qup S α ΔTup ΔTbase ΔToverland 210.10 157.48 70.03 157.48 117 Table 2.16: Type III sums of squared errors for a model testing for both thermal class as a fixed effect and individual stream reach as a random effect on downstream temperature flux rates using GLM. Source . . . 0.00 1744.81 DF Type III SS Mean Square F Value Pr > F 0 16 1 1 1 1 1 1 1 1 2102.80 <.0001 1089.75 <.0001 0.4668 432.23 <.0001 379.78 <.0001 375.25 <.0001 1231.38 <.0001 939.66 <.0001 611.04 <.0001 109.05 56.51 0.03 22.42 19.70 19.46 63.86 48.73 31.69 56.51 0.03 22.42 19.70 19.46 63.86 48.73 31.69 0.53 Thermal Class Stream Ta - Tw Qup Qdown - Qup S α ΔTup ΔTbase ΔToverland 118 Table 2.17: Variance accounted for by each model from Tables 2.15 and 2.16. r Thermal Class 0.30 0.57 Stream Coeff Var Root MSE ΔT Mean 1129.90 974.69 0.02 0.02 0.26 0.23 119 Table 2.18: Diagnostics comparing observed values of downstream temperature flux (°C/km) with predicted values using Model 10 optimized specifically across thermal classes and across the stream network as a whole. Observed and Predicted refer to mean values of downstream temperature flux across the entire study period for each stream reach. Composite Thermal Class Predicted r SSE RMSE 19 277 0.07 0.20 0.74 0.47 0.14 26 0.10 Stream Cedar Creek Cedar River East Branch Black River Pigeon River Pokagon Creek Prairie River Fish Creek Middle Branch Tobacco River Butterfield Creek Slapneck Creek Squaw Creek Middle Branch Escanaba River Carp River Swan Creek Honeyoey Creek King Creek Morgan Creek Spring Creek Nottawa Creek Hasler Creek North Branch Thunder Bay River Observed 0.02 0.22 0.09 -0.06 0.07 0.04 0.04 0.13 -0.11 0.17 0.57 -0.01 0.02 0.10 0.18 0.07 -0.17 -0.04 -0.71 -0.17 0.08 Predicted r -0.03 0.09 0.06 0.49 SSE RMSE 55 405 0.12 0.24 0.15 0.17 38 0.12 0.01 0.06 0.00 -0.01 0.42 0.19 0.23 -0.11 72 74 182 53 0.09 0.09 88 -0.02 0.14 0.00 0.63 0.17 0.12 257 260 1620 -0.02 0.34 49 0.00 0.04 0.01 0.01 0.02 -0.09 0.00 -0.05 0.02 -0.23 0.31 0.27 0.24 -0.36 0.13 0.30 0.72 0.55 108 598 192 53 1973 625 3899 80 37 120 0.10 0.11 0.17 0.09 0.16 0.19 0.30 0.77 0.12 0.13 0.43 0.24 0.13 0.58 0.34 0.79 0.16 0.11 0.02 0.22 0.09 0.01 0.06 0.00 -0.02 0.37 0.29 0.19 -0.24 74 65 191 61 0.09 0.16 83 -0.03 0.16 0.07 -0.02 0.02 0.00 0.10 0.02 0.13 0.13 -0.17 -0.16 0.03 0.55 0.21 0.20 0.36 -0.24 0.27 0.65 0.39 0.00 0.08 0.25 0.76 0.33 258 245 1411 49 86 621 86 44 1894 906 2632 23 43 C o l d C o l d - t r a n s i t i o n a l t r a n s i t i o n a l W a r m - W a r m 0.11 0.10 0.18 0.10 0.15 0.19 0.29 0.72 0.12 0.12 0.43 0.16 0.12 0.56 0.42 0.65 0.09 0.12 Table 2.19: Diagnostics comparing observed values of downstream temperature flux (°C/km) with predicted values using Model 10 optimized specifically across thermal classes and across the stream network as a whole. Results are averaged among streams within thermal classes. Thermal Class C CT WT W Composite Thermal Class Observed Predicted r SSE RMSE Predicted r SSE RMSE 0.11 0.09 0.03 -0.27 0.07 0.03 0.00 -0.01 0.24 0.19 0.12 0.52 166 276 688 1339 0.16 0.22 0.34 0.35 0.11 0.04 0.08 -0.10 0.45 0.19 0.28 0.45 107 252 711 900 0.12 0.21 0.34 0.28 121 Table 2.20: Average monthly downstream temperature flux rates (°C/km) for baseflow reduction scenarios using Model 10. Stream Cedar Creek Cedar River East Branch Black River Pigeon River Pokagon Creek Prairie River Fish Creek Middle Branch Tobacco River Butterfield Creek Slapneck Creek Squaw Creek Middle Branch Escanaba River Carp River Swan Creek Honeyoey Creek King Creek Morgan Creek Spring Creek Nottawa Creek Hasler Creek North Branch Thunder Bay River July 50% -0.10 0.29 0% 0.08 0.25 August September October 90% -0.46 0.40 0% 0.07 0.20 50% -0.11 0.24 90% -0.51 0.33 0% -0.04 0.14 50% -0.19 0.21 90% -0.47 0.25 0% -0.10 0.07 50% -0.21 0.12 90% -0.36 0.13 0.12 0.12 0.12 0.10 0.10 0.10 0.07 0.07 0.07 0.05 0.05 0.05 -0.04 0.13 0.11 0.08 -0.06 0.12 0.12 0.03 -0.07 0.12 0.13 -0.02 -0.05 0.10 0.06 0.06 -0.06 0.10 0.09 0.00 -0.06 0.10 0.12 -0.02 -0.08 0.06 0.00 0.01 -0.08 0.06 0.00 -0.02 -0.06 0.05 0.06 -0.03 -0.06 0.01 0.01 0.01 -0.06 0.01 0.01 -0.01 -0.06 0.00 0.03 -0.03 0.16 0.16 0.15 0.15 0.14 0.13 0.12 0.11 0.11 0.08 0.07 0.06 -0.16 0.31 0.85 -0.05 0.29 1.09 0.08 0.28 1.16 -0.18 0.14 0.76 -0.08 0.13 0.98 0.06 0.12 1.20 -0.13 0.06 0.62 -0.07 0.04 0.71 0.04 0.02 0.80 -0.01 -0.01 0.42 0.04 -0.04 0.45 0.10 -0.05 0.47 0.09 0.16 0.17 0.01 0.07 0.09 -0.12 -0.08 -0.05 -0.20 -0.15 -0.12 0.01 0.27 0.26 0.15 -0.31 -0.10 -0.23 -0.26 0.00 0.25 0.24 0.14 -0.49 -0.10 -0.21 -0.25 -0.01 0.24 0.23 0.13 -0.62 -0.09 -0.19 -0.16 0.02 0.16 0.22 0.11 -0.20 -0.12 -0.19 -0.26 0.01 0.16 0.22 0.10 -0.44 -0.10 -0.18 -0.26 0.02 0.16 0.25 0.09 -0.63 -0.03 -0.17 -0.11 0.02 0.08 0.16 0.07 -0.05 -0.05 -0.20 -0.19 0.03 0.06 0.18 0.07 -0.25 -0.01 -0.19 -0.18 0.04 0.05 0.42 0.07 -0.54 0.01 -0.18 -0.03 0.05 0.00 0.00 0.02 -0.15 0.07 -0.14 -0.08 0.06 -0.01 -0.06 0.02 -0.35 0.09 -0.13 -0.08 0.07 0.00 -0.14 0.02 -0.62 0.10 -0.13 0.04 0.15 0.15 0.15 0.13 0.13 0.13 0.06 0.06 0.06 0.05 0.05 0.04 122 Table 2.21: Average monthly downstream temperature flux rates (°C/km) of thermal classes for baseflow reduction scenarios using Model 10. Thermal Class Cold Cold-transitional Warm-transitional 0% 0.50 0.22 0.33 July 50% 0.33 0.33 0.32 August September October 90% 0.55 0.41 0.40 0% 0.48 0.05 0.47 50% 90% 0.32 0.16 0.52 0.53 0.25 0.59 0% 0.45 -0.06 0.67 50% 90% 0.35 0.00 0.68 0.41 0.07 0.81 0% 0.42 50% 0.34 90% 0.30 -0.11 -0.08 -0.04 1.85 1.74 1.68 Warm -0.84 -0.57 -0.25 -0.65 -0.53 -0.35 -0.80 -0.69 -0.43 -0.77 -0.69 -0.58 123 Table 2.22: August mean values used for predicting downstream temperature flux using Model 10. Ta = air temperature (°C); Tw = water temperature (°C); Qup = upstream discharge (cms); Qdown = downstream discharge; S = day length (hr); ΔTup = upstream heat energy (°C); ΔTbase = baseflow heat energy (°C); ΔTover = overland flow heat energy (°C); α = sun altitude angle (°). Site Cedar Creek Cedar River East Branch Black River Pigeon River Pokagon Creek Prairie River Fish Creek Middle Branch Tobacco River Butterfield Creek Slapneck Creek Squaw Creek Middle Branch Escanaba River Carp River Swan Creek Honeyoey Creek King Creek Morgan Creek Spring Creek Nottawa Creek Hasler Creek North Branch Thunder Bay River Reach length (km) 2.55 1.93 2.88 6.55 4.05 5.86 5.19 4.09 6.01 1.56 1.44 6.13 2.61 2.54 6.64 5.82 0.91 1.53 3.76 7.59 4.63 Ta – Tw (C) Qup Qdown - Qup S α ΔTbase ΔTover ΔTup 5.07 6.50 5.07 3.14 2.95 3.25 1.86 4.76 2.59 1.54 1.92 -2.12 0.83 0.54 2.28 -8.15 -0.14 -0.46 -0.84 -0.91 -9.31 0.36 0.42 -0.04 0.15 0.00 0.06 0.11 0.02 0.10 1.03 0.03 0.30 -0.23 0.00 0.03 0.00 0.03 0.09 -0.05 0.05 -0.06 0.09 1.33 0.68 0.44 0.45 0.17 0.95 0.51 0.06 0.27 0.02 0.70 1.62 0.40 0.04 0.00 0.10 0.06 0.68 0.02 0.34 124 14.32 14.38 14.38 14.20 14.13 14.12 14.23 20.28 19.98 19.86 20.26 20.66 20.58 20.33 -9.81 -10.96 -1.51 -11.69 -4.14 -11.31 -11.65 -13.92 -13.19 -1.50 -12.11 -3.54 -13.79 -15.74 -12.56 -3.38 -0.01 -4.25 -0.73 -6.28 -4.95 14.29 20.17 -10.14 -9.47 -0.78 14.31 14.27 14.26 20.00 19.08 18.91 -13.02 -16.50 -16.83 -15.04 -8.74 -11.47 -11.57 -13.43 -8.51 14.50 19.59 -14.56 -19.05 -6.47 14.50 13.98 14.25 14.37 14.51 14.48 14.13 13.99 19.48 19.82 20.24 19.88 19.49 19.50 20.54 20.98 -5.90 -15.07 -6.84 -7.53 -15.57 -14.61 -5.77 -13.34 -6.20 -14.90 -19.32 -8.03 -16.54 -17.73 -8.96 -22.79 -0.86 -0.37 -10.52 -2.01 -4.56 -12.76 -3.94 -16.45 14.39 19.85 -0.07 -0.07 0.00 Table 2.23: August mean values used for predicting downstream temperature flux using the Zorn et al. (2008) model. Site Cedar Creek Cedar River East Branch Black River Pigeon River Pokagon Creek Prairie River Fish Creek Middle Branch Tobacco River Butterfield Creek Slapneck Creek Squaw Creek Middle Branch Escanaba River Carp River Swan Creek Honeyoey Creek King Creek Morgan Creek Spring Creek Nottawa Creek Hasler Creek North Branch Thunder Bay River Reach length (km) Qup (cms) Ta (°C) Tw (°C) Depth (m) Velocity (m/s) 20.86 20.43 21.64 21.47 21.45 21.25 20.99 15.79 13.93 16.57 18.34 18.50 18.00 19.13 22.72 17.96 20.33 19.68 19.63 17.74 18.14 17.71 16.93 19.05 18.98 21.93 22.41 10.34 18.98 18.76 21.55 22.14 18.14 21.39 20.13 18.49 19.12 19.21 22.39 23.05 10.34 19.65 0.17 0.42 0.60 0.48 0.44 0.20 0.46 0.29 0.31 0.29 0.28 0.42 0.33 0.39 0.33 0.13 0.21 0.29 0.28 0.22 0.37 2.55 1.93 2.88 6.55 4.05 5.86 5.19 4.09 6.01 1.56 1.44 6.13 2.61 2.54 6.64 5.82 0.91 1.53 3.76 7.59 4.63 0.09 1.33 0.68 0.44 0.45 0.17 0.95 0.51 0.06 0.27 0.02 0.70 1.62 0.40 0.04 3.95 e-03 0.10 0.06 0.68 0.02 0.34 125 0.16 0.36 0.15 0.10 0.15 0.09 0.13 0.33 0.04 0.17 0.02 0.13 0.53 0.11 0.04 0.01 0.11 0.03 0.21 0.02 0.10 Table 2.24: Observed August downstream temperature flux rates (°C/km) for each of the stream reaches and the predicted rates of Model 10 and the Zorn et al. (2008) model under a 0% reduction scenario. Stream Cedar Creek Cedar River East Branch Black River Pigeon River Pokagon Creek Prairie River Fish Creek Middle Branch Tobacco River Butterfield Creek Slapneck Creek Squaw Creek Middle Branch Escanaba River Carp River Swan Creek Honeyoey Creek King Creek Morgan Creek Spring Creek Nottawa Creek Hasler Creek North Branch Thunder Bay River °C/km Observed M10 0.07 0.24 0.11 -0.04 0.11 0.06 0.06 0.16 -0.15 0.16 0.77 -0.01 0.03 0.21 0.24 0.12 -0.19 -0.07 -0.22 -0.26 0.16 0.07 0.20 0.10 -0.05 0.10 0.06 0.06 0.15 -0.18 0.14 0.76 0.01 0.02 0.16 0.22 0.11 -0.20 -0.12 -0.19 -0.26 0.13 Zorn et al. (2008) model 4.21 1.81 11.67 5.48 3.59 2.53 1.65 0.65 17.74 2.45 106.64 -1.77 0.05 1.52 13.48 -221.41 -0.79 -12.06 -0.22 -8.93 -15.85 126 Table 2.25: Average August downstream temperature flux rates (°C/km) experienced by each of the four stream thermal classifications and the estimated rates using Model 10 and the Zorn et al. (2008) model (King Creek and Squaw Creek excluded from Zorn et al. (2008) model estimates). Thermal Class Observed M10 0.12 0.11 0.03 -0.11 0.14 0.11 0.06 -0.11 C CT WT W Zorn et al. (2008) model 5.90 3.60 0.54 -8.33 127 Table 2.26: Predicted downstream temperature flux rates (°C/km) following scenario analysis of baseflow reductions of 10% increments using Model 10. Stream Cedar Creek Cedar River East Branch Black River Pigeon River Pokagon Creek Prairie River Fish Creek Middle Branch Tobacco River Butterfield Creek Slapneck Creek Squaw Creek Middle Branch Escanaba River Carp River Swan Creek Honeyoey Creek King Creek Morgan Creek Spring Creek Nottawa Creek Hasler Creek North Branch Thunder Bay River 0% 0.07 0.20 0.10 -0.05 0.10 0.06 0.06 0.15 -0.18 0.14 0.76 0.01 0.02 0.16 0.22 0.11 -0.20 -0.12 -0.19 -0.26 0.13 10% 0.04 0.21 0.10 -0.05 0.10 0.06 0.05 0.15 -0.17 0.14 0.80 0.02 0.02 0.16 0.24 0.10 -0.24 -0.12 -0.19 -0.26 0.13 20% 0.01 0.21 0.10 -0.05 0.10 0.07 0.05 0.14 -0.15 0.14 0.84 0.03 0.02 0.16 0.26 0.10 -0.28 -0.12 -0.19 -0.26 0.13 30% -0.02 0.22 0.10 -0.06 0.10 0.07 0.03 0.14 -0.13 0.13 0.88 0.04 0.02 0.16 0.26 0.10 -0.33 -0.11 -0.19 -0.26 0.13 °C/km 40% -0.06 0.23 0.10 -0.06 0.10 0.08 0.01 0.14 -0.11 0.13 0.93 0.06 0.01 0.16 0.25 0.10 -0.38 -0.11 -0.19 -0.26 0.13 50% -0.11 0.24 0.10 -0.06 0.10 0.09 0.00 0.14 -0.08 0.13 0.98 0.07 0.01 0.16 0.22 0.10 -0.44 -0.10 -0.18 -0.26 0.13 60% -0.17 0.26 0.10 -0.07 0.10 0.09 0.00 0.14 -0.05 0.13 1.03 0.08 0.01 0.15 0.22 0.10 -0.52 -0.08 -0.18 -0.24 0.13 70% -0.25 0.28 0.10 -0.06 0.10 0.10 0.00 0.13 -0.03 0.13 1.08 0.08 0.01 0.16 0.24 0.10 -0.55 -0.07 -0.18 -0.22 0.13 80% -0.36 0.30 0.10 -0.06 0.10 0.11 -0.01 0.13 0.01 0.12 1.14 0.09 0.02 0.16 0.24 0.09 -0.59 -0.05 -0.18 -0.18 0.13 90% -0.51 0.33 0.10 -0.06 0.10 0.12 -0.02 0.13 0.06 0.12 1.20 0.09 0.02 0.16 0.25 0.09 -0.63 -0.03 -0.17 -0.11 0.13 128 Table 2.27: Predicted downstream temperature flux rates (°C/km) following scenario analysis of baseflow reductions of 10% increments using the Zorn et al. (2008) model. 0% 4.21 1.81 10% 4.48 1.91 20% 4.81 2.03 30% 5.21 2.18 °C/km 50% 6.39 2.61 40% 5.72 2.37 60% 7.32 2.95 70% 8.74 3.45 80% 11.25 4.33 90% 17.45 6.43 11.67 12.21 12.85 13.60 14.51 15.66 17.17 19.28 22.59 29.15 Stream Cedar Creek Cedar River East Branch Black River Pigeon River Pokagon Creek Prairie River Fish Creek Middle Branch Tobacco River 5.48 3.59 2.53 1.65 0.65 5.71 3.77 2.67 1.74 0.69 5.97 3.98 2.83 1.83 0.73 6.27 4.23 3.02 1.94 0.79 6.64 4.54 3.26 2.08 0.86 7.09 4.94 3.57 2.26 0.96 7.66 5.46 3.98 2.49 1.09 8.45 6.20 4.57 2.82 1.30 9.64 7.39 5.54 3.35 1.66 29.92 5.74 11.96 9.85 7.64 4.45 2.56 42.38 8.26 Butterfield Creek Slapneck Creek 17.74 2.45 18.24 2.59 18.82 2.75 19.52 2.95 20.40 3.21 21.54 3.53 23.12 3.98 25.53 4.63 Squaw Creek Middle Branch Escanaba River Carp River Swan Creek Honeyoey Creek King Creek Morgan Creek Spring Creek Nottawa Creek Hasler Creek North Branch Thunder Bay River 106.64 110.40 114.92 120.52 127.70 137.38 151.43 174.44 221.96 415.44 -1.77 -1.86 -1.96 -2.08 -2.23 -2.42 -2.67 -3.02 -3.58 -4.72 0.05 1.52 13.48 -221.41 -0.79 -12.06 -0.22 -8.93 0.05 1.59 13.89 -240.91 -0.83 -12.40 -0.23 -9.25 0.05 1.67 14.39 -267.21 -0.88 -12.80 -0.25 -9.64 0.06 1.77 14.98 -304.95 -0.95 -13.28 -0.26 -10.13 0.06 1.89 15.74 -364.37 -1.02 -13.88 -0.29 -10.75 0.07 2.04 16.73 -473.47 -1.12 -14.66 -0.32 -11.59 0.08 2.24 18.13 -747.38 -1.25 -15.74 -0.36 -12.82 0.09 2.51 20.31 0.12 2.94 24.47 0.18 3.80 37.64 -2940.45 -2.E+148 -6.E+56 -1.44 -17.38 -0.42 -14.85 -1.75 -20.39 -0.53 -19.12 -2.46 -28.97 -0.78 -37.60 -15.85 -16.58 -17.43 -18.44 -19.65 -21.18 -23.17 -25.94 -30.24 -38.84 129 Table 2.28: Parameterized beta coefficients associated with each stream reach to estimate downstream temperature flux using Model 10. Ta = air temperature (°C); Tw = water temperature (°C); Qup = upstream discharge (cms); Qdown = downstream discharge; S = day length (hr); ΔTup = upstream heat energy (°C); ΔTbase = baseflow heat energy (°C); ΔTover = overland flow heat energy (°C); α = sun altitude angle (°). Site Thermal Class intercept Ta – Tw Qup Qdown – Qup S ΔTup ΔTbase ΔTover α - - 0.001 0.128 0.053 -0.059 0.286 -0.278 0.024 0.112 0.015 -0.010 - 0.169 0.006 0.003 0.020 - 0.003 - 0.010 0.201 -0.235 0.119 0.021 0.001 0.010 - - 0.001 0.056 - - 0.016 0.397 0.028 -0.026 0.384 -0.359 Swan Creek WT -2.471 0.155 -0.179 -1.988 0.249 130 Cedar Creek Cedar River East Branch Black River Pigeon River Pokagon Creek Prairie River Fish Creek Middle Branch Tobacco River Butterfield Creek Slapneck Creek Squaw Creek Middle Branch Escanaba River Carp River C C C CT CT CT CT CT CT CT CT CT CT -2.732 -0.008 4.463 0.107 0.107 0.128 -0.045 -0.672 0.014 -0.277 0.076 0.128 - 0.005 0.025 -0.003 0.029 -0.105 -0.011 -0.241 -0.361 0.033 0.006 0.823 -0.912 0.905 -1.248 0.042 0.127 0.145 0.089 -1.747 0.005 0.000 -0.345 0.155 -2.263 -0.033 -2.910 1.156 0.309 - - 0.006 0.048 - - 0.005 0.084 0.106 -0.051 0.021 -0.019 - 0.013 0.080 0.002 0.029 -0.335 0.009 0.059 0.038 0.002 0.005 0.681 0.032 -1.350 0.405 0.006 0.019 - 0.022 - 0.334 -0.440 0.048 -3.320 -0.286 0.244 -0.661 -0.017 -2.451 0.195 0.122 0.918 0.011 -14.190 0.576 -0.025 -7.020 -0.011 -0.312 -0.301 0.545 0.283 0.014 -0.029 -0.112 -0.017 WT WT WT WT W W W Table 2.28: (cont’d) Honeyoey Creek King Creek Morgan Creek Spring Creek Nottawa Creek Hasler Creek North Branch Thunder Bay River -0.056 0.101 -4.613 1.750 0.135 - 0.007 0.071 0.159 -0.078 -1.905 -0.012 0.883 -3.095 0.178 -0.910 0.019 1.666 8.136 0.433 -0.032 -0.040 0.597 0.018 - - 0.004 0.015 - - 0.008 -0.024 0.003 0.009 0.002 0.148 -0.127 0.077 -0.075 0.006 0.011 0.447 0.002 -0.037 -0.288 -0.073 - 0.005 0.006 0.038 -0.024 0.663 0.000 2.294 5.003 -0.038 0.015 0.035 0.063 0.058 -0.784 0.026 -0.369 0.131 0.115 0.005 0.085 0.231 -0.232 131 Table 2.29: Parameterized beta coefficients specific to each thermal class and across all streams to estimate downstream temperature flux using Model 10. Ta = air temperature (°C); Tw = water temperature (°C); Qup = upstream discharge (cms); Qdown = downstream discharge; S = day length (hr); ΔTup = upstream heat energy (°C); ΔTbase = baseflow heat energy (°C); ΔTover = overland flow heat energy (°C); α = sun altitude angle (°). Thermal Class C CT WT W Composite intercept -1.163 -0.748 -1.514 0.044 -0.677 Ta – Tw 0.005 0.016 0.019 -0.042 0.020 Qup 0.112 -0.140 -0.084 -0.279 -0.141 Qdown – Qup S 0.117 0.111 0.795 0.880 0.211 0.098 0.090 0.159 -0.014 0.088 α 0.002 -0.002 -0.002 0.008 -0.003 ΔTup 0.018 0.049 0.014 0.091 0.045 ΔTbase 0.046 -0.017 0.070 0.016 0.003 ΔTover -0.039 0.027 -0.034 0.007 0.016 132 Table 2.30: Partial R2 values of each variable included in Model 10 for summer months (July and August). Values reflect the influence of each variable on downstream temperature flux rates. Ta = air temperature (°C); Tw = water temperature (°C); Qup = upstream discharge (cms); Qdown = downstream discharge; S = day length (hr); ΔTup = upstream heat energy (°C); ΔTbase = baseflow heat energy (°C); ΔTover = overland flow heat energy (°C); α = sun altitude angle (°). Variable Ta - Tw Qup Qdown - Qup Overall 0.04 0.01 0.04 Partial R2 S 0.04 ΔTup ΔTbase ΔToverland 0.03 0.05 0.05 α 0.08 Table 2.31: Partial R2 values averaged over thermal class for summer months (July and August). Ta = air temperature (°C); Tw = water temperature (°C); Qup = upstream discharge (cms); Qdown = downstream discharge; S = day length (hr); ΔTup = upstream heat energy (°C); ΔTbase = baseflow heat energy (°C); ΔTover = overland flow heat energy (°C); α = sun altitude angle (°). Partial Regression S 0.04 0.06 0.02 0.04 ΔTup ΔTbase ΔToverland 0.02 0.03 0.03 0.00 0.02 0.04 0.13 0.00 0.01 0.06 0.09 0.00 α 0.07 0.11 0.02 0.09 Thermal Class Ta - Tw Qup Qdown - Qup C CT WT W 0.02 0.05 0.04 0.02 0.06 0.04 0.02 0.08 0.00 0.01 0.01 0.00 133 Figure 2.1: Stream gauge locations. 134 Figure 2.2: Downstream temperature flux rates (°C/km) compared to upstream discharge (m3/s) for representative stream reaches among each of the four thermal classes. 135 Figure 2.3: Downstream temperature flux rates (°C/km) compared to upstream discharge (m3/s) for all stream reaches within of each of the four thermal classes combined. 136 Figure 2.4: Hourly downstream temperature flux (°C) (Tdown – Tup) represented with individual points and overlain by a LOESS regression to provide a smoothed representation of downstream temperature flux over time. 137 Figure 2.5: Model 10 predictions of downstream temperature flux rates (°C/km) compared with LOESS smoothed observed downstream temperature change for four different streams, one in each of the four thermal classes: Cedar River (C), Fish Creek (CT), Honeyoey Creek (WT), North Branch Thunder Bay River (W). 138 Figure 2.6: Model 11 predictions of downstream temperature flux rates (°C/km) compared with LOESS smoothed observed downstream temperature change for four different streams, one in each of the four thermal classes: Cedar River (C), Fish Creek (CT), Honeyoey Creek (WT), North Branch Thunder Bay River (W). 139 Figure 2.7: Model 9 predictions of downstream temperature flux rates (°C/km) compared with LOESS smoothed observed downstream temperature change for four different streams, one in each of the four thermal classes: Cedar River (C), Fish Creek (CT), Honeyoey Creek (WT), North Branch Thunder Bay River (W). 140 Figure 2.8: Model 8 predictions of downstream temperature flux rates (°C/km) compared with LOESS smoothed observed downstream temperature change for four different streams, one in each of the four thermal classes: Cedar River (C), Fish Creek (CT), Honeyoey Creek (WT), North Branch Thunder Bay River (W). 141 Figure 2.9: Seasonal structures of error (Obs – Sim) of downstream temperature flux represented using LOESS regression for each of the four best models identified using model selection criteria. One stream from each of the four thermal classifications is represented: Cedar River (C), Fish Creek (CT), Honeyoey Creek (WT), North Branch Thunder Bay River (W). 142 Figure 2.10: Comparison of model fit using a LOESS regression. Using nonlinear optimization, Model 10 was parameterized using a composite dataset of the entire stream network, and also specific to each thermal class. Model fit was then compared to observed values. 143 Figure 2.11: Baseflow reduction scenarios for example streams. Scenarios of 0%, 50%, and 90% baseflow reduction are plotted using LOESS regressions, along with observed values. LOESS regressions track seasonal and yearly trends of downstream temperature flux. 144 August downstream temperature flux response to baseflow reduction using Model 10 1.50 1.00 m k / C ° 0.50 0.00 -0.50 -1.00 0% 20% 40% 60% 80% 100% Cedar Creek Cedar River East Branch Black River Pigeon River Pokagon Creek Prairie River Fish Creek Middle Branch Tobacco River Butterfield Creek Slapneck Creek Squaw Creek Middle Branch Escanaba River Carp river Swan Creek Honeyoey Creek King Creek Morgan Creek Spring Creek Nottawa Creek Hasler Creek Figure 2.12: Downstream temperature flux rates in response to baseflow reductions. % Reduction North Branch Thunder Bay River 145 August downstream temperature flux rate response to baseflow reduction using Zorn et al. (2008) m k / C ° 50.00 40.00 30.00 20.00 10.00 0.00 -10.00 -20.00 -30.00 -40.00 -50.00 Cedar Creek Cedar River East Branch Black River Pigeon River Prairie River Pokagon Creek Fish Creek Middle Branch Tobacco River Butterfield Creek Slapneck Creek Middle Branch Escanaba River Carp River Swan Creek Honeyoey Creek Morgan Creek Spring Creek Nottawa Creek Hasler Creek 0% 20% 40% 60% 80% 100% North Branch Thunder Bay River % Reduction Figure 2.13: Downstream temperature flux rates in response to baseflow reductions. Rates were omitted for King Creek and Squaw Creek to maintain an appropriate scale on the y-axis. 146 Derivation of net incoming shortwave radiation (Meeus, 1991). APPENDIX 2.1 I0 = solar irradiation values throughout the year (W/m2) ISC = solar constant (1353 W/m2) N = day of year δ = declination angle LST = local standard time for time zone (adjusted for daylight savings time; LST = DST - 1) Long = longitude LSTM = local longitude of standard time meridian ET = equation of time in minutes H = hour angle (azimuth angle of sun’s rays caused by earth’s rotation) θz = zenith angle β1 = altitude angle 147 L = latitude α1 = solar azimuth angle IDN = direct normal irradiance to the ground (W/m2) A = apparent extraterrestrial solar intensity (W/m2) B = atmospheric extinction coefficient p/p0 = pressure at location of interest (atm) z = altitude above sea level (ft) Hs = net solar heat flux (W/m2) SF = shading factor (%) α = surface albedo (.1) 148 REFERENCES 149 REFERENCES Adams, S. M. & Greeley, M. S. 2000. Ecotoxicological indicators of water quality: using multi- response indicators to assess the health of aquatic ecosystems. Water, Air, and Soil Pollution 123:103-115. Akaike, H. 1973. Maximum likelihood identification of Gaussian autoregressive moving average models. Biometrika 60:255-265. Albert, D. A., Denton, S. R., & Barnes, B. V. 1986. Regional landscape ecosystems of Michigan. School of Natural Resources, The University of Michigan, Ann Arbor. Andersen, T., Carstensen, J., Hernández-Garcia, E., & Duarte, C. M. 2009. Ecological thresholds and regime shifts: approaches to identification. Trends in Ecology and Evolution 24:49- 57. Anonymous. 2001. The Great Lakes Charter Annex: a supplementary agreement to the Great Lakes Charter. Council of Great Lakes Governors, Chicago. Arthington, A. H., Bunn, S. E., Poff, N. L., & Naiman, R. J. 2006. The challenge of providing environmental flow rules to sustain river ecosystems. Ecological Applications 16:1311- 1318. Baker, M. E. & King, R. S. 2010. A new method for detecting and interpreting biodiversity and ecological community thresholds. Methods in Ecology and Evolution 1:25-37. Beisner, B. E., Haydon, D. T., & Cuddington, K. 2003. Alternative stable states in ecology. Frontiers in Ecology and the Environment 1:376-382. Benyahya, L., Caissie, D., St-Hilaire, A., Ouarda, T. B., & Bobée, B. 2007. A review of statistical water temperature models. Canadian Water Resources Journal 32:179-192. Beschta, R. L., Bilby, R. E., Brown, G. W., Holtby, L. B., & Hofstra, T. D. 1987. Stream temperature and aquatic habitat: fisheries and forest interactions. In Salo, E. O. and Cundy, T. W. (ed.): Streamside Management: Forestry and Fishery Interactions, Contribution No. 57, p. 191-232. University of Washington, Institute of Forest Resources, Seattle, Washington. Bovee, K. D. 1982. A Guide to Stream Habitat Analysis Using the Instream Flow Incremental Methodology. Instream Flow Information Paper 21. FWS/OBS-82/26, U.S. Department of Fisheries and Wildlife Services, Washington, D.C., U.S.A. Brenden, T. O., Wang, L., Su, Z. 2008. Quantitative identification of disturbance thresholds in support of aquatic resource management. Environmental Management 42:821-832. 150 Brocard, D. N. & Harleman, D. R. 1976. One-dimensional temperature predictions in unsteady flows. Journal of the Hydraulics Division 102:227-240. Buckland, S. T., Burnham, K. P., & Augustin, N. H. 1997. Model selection: an integral part of inference. Biometrics 53:603-618. Bustillo, V., Moatar, F., Ducharne, A., Thiéry, D., & Poirel, A. 2014. A multimodel comparison for assessing water temperatures under changing climate conditions via the equilibrium temperature concept: a case study of the Middle Loire River, France. Hydrological Processes: An International Journal 28:1507-1524. Caissie, D. 2001. Modelling of maximum daily water temperatures in a small stream using air temperatures. Journal of Hydrology 251:14-28. Caissie, D., Satish, M. G., & El-Jabi, N. 2005. Predicting river water temperatures using the equilibrium temperature concept with application on Miramichi River catchments (New Brunswick, Canada). Hydrological Processes: An International Journal 19:2137-2159. Caissie, D. 2006. The thermal regime of rivers: a review. Freshwater Biology 51:1389-1406. Caisse, D., Satish, M. G., & El-Jabi, N. 2007. Predicting water temperatures using a deterministic model: application on Miramichi River catchments (New Brunswick, Canada). Journal of Hydrology 336:303-315. Caldwell, R. J., Gangopadhyay, S., Bountry, J., Lai, Y., & Elsner, M. M. 2013. Statistical modeling of daily and subdaily stream temperatures: application to the Methow River Basin, Washington. Water Resources Research 49:4346-4361. Chen, J., Franklin, J. F., & Spies, T. A. 1995. Growing-season microclimatic gradients from clearcut edges into old-growth Douglas-fir forests. Ecological Applications 5:74-86. Cheng, S. -T. & Wiley, M. J. 2016. A Reduced Parameter Stream Temperature Model (RPSTM) for basin-wide simulations. Environmental Modelling & Software 82:295-307. Clements, W. H., Vieira, N. K. M., Sonderegger, D. L. 2010. Use of ecological thresholds to assess recovery in lotic ecosystems. Journal of the North American Benthological Society 29:1017-1023. Cleveland, W. S. 1979. Robust locally weighted regression and smoothing scatterplots. Journal of the American Statistical Association 74:829-836. Cooper, R. M. 2002. Determining surface water availability in Oregon. Oregon Water Resources Department Open File Report SW 02-002. 151 Davis, L. J., Reiter, M., & Groom, J. D. 2016. Modelling temperature change downstream of forest harvest using Newton’s law of cooling. Hydrological Processes: An International Journal, 30:959-971. Dufrene, M. & Legendre, P. 1997. Species assemblages and indicator species: the need for a flexible asymmetrical approach. Ecological Monographs 67:345-366. Erickson, T. R. & Stefan, H. G. 2000. Linear air/water temperature correlations for streams during open water periods. Journal of Hydrologic Engineering 5:317-321. Evans, E. C., McGregor, G. R., & Petts, G. E. 1998. River energy budgets with special reference to river bed processes. Hydrological Processes: An International Journal 12:575-595. Ficetola, G. F. & Denoël, M. 2009. Ecological thresholds: an assessment of methods to identify abrupt changes in species-habitat relationships. Ecography 32:1075-1084. Fry, F. E. J. 1971. The effect of environmental factors on the physiology of fish. In Hoar, W. S. & Randall, D. J. (ed.): Fish physiology, Vol. 6, Environmental relations and behavior, p. 1-98. Academic Press, New York. Galli, J. 1990. Thermal impacts associated with urbanization and stormwater best management practices. Metropolitan Washington Council of Governments, Washington, D.C., 188 p. Groffman, P. M., Baron, J. S., Blett, T., Gold, A. J., Goodman, I., Gunderson, L. H., Levinson, B. M., Palmer, M. A., Paerl, H. W., Peterson, G. D., Poff, N. L., Rejeski, D. W., Reynolds, J. F., Turner, M. G., Weathers, K. C., & Wiens, J. 2006. Ecological thresholds: the key to successful environmental management or an important concept with no practical application?. Ecosystems 9:1-13. Groundwater Conservation Advisory Council. 2006. Final report to the Michigan Legislature in response to Public Act 148 of 2003. Michigan Department of Environmental Quality. Guisan, A. & Zimmermann, N. E. 2000. Predictive habitat distribution models in ecology. Ecological Modelling 135:147-186. Hamilton, D. A. & Seelbach, P. W. 2011. Michigan’s water withdrawal assessment process and screening tool. Michigan Department of Natural Resources, Fisheries Special Report 55, Lansing. Hampel, F. R., Ronchetti, E. M., Rousseeuw, P. J., & Stahel, W. A. 1986. Robust statistics: the approach based on influence functions. John Wiley & Sons, New York, New York. Hannah, D. M., Malcolm, I. A., & Bradley, C. 2009. Seasonal hyporheic temperature dynamics over riffle beforms. Hydrological Processes: An International Journal 23:2178-2194. 152 Hidayat, Y., Murtilaksono, K., & Sinukaban, N. 2012. Characterization of surface runoff, soil erosion and nutrient loss on forest-agriculture landscape. Journal of Tropical Soils 17:259-266. Hudson, P. L., Griffiths, R. W., & Wheaton, T. J. 1992. Review of habitat classification schemes appropriate to streams, rivers, and connecting channels in the Great Lakes drainage basin. In Busch, W.D.N. & Sly, P.G. (ed.): The development of an aquatic habitat classification system for lakes, p. 73-107. CRC Press, Boca Raton, Florida. Huff, D. D., Hubler, S. L., & Borisenko, A. N. 2005. Using field data to estimate the realized thermal niche of aquatic vertebrates. North American Journal of Fisheries Management 25:346-360. Huggett, A. J. 2005. The concept and utility of ‘ecological thresholds’ in biodiversity conservation. Biological Conservation 124:301-310. Kanno, Y. & Vokoun, J. C. 2010. Evaluating effects of water withdrawals and impoundments on fish assemblages in southern New England streams, USA. Fisheries Management and Ecology 17:272-283. Keith, R. M., Bjornn, T., Meehan, W. R., Hetrick, N. J., & Brusven, M. A. 1998. Response of juvenile salmonids to riparian and instream cover modifications in small streams flowing through second-growth forests of southeast Alaska. Transactions of the American Fisheries Society 127:889-907. Kendy, E. & Bredehoeft, J. D. 2006. Transient effects of groundwater pumping and surface- water-irrigation returns on streamflow. Water Resources Research 42. Kenny, J. F., Barber, N. L., Hutson, S. S., Linsey, K. S., Lovelace, J. K., & Maupin, M. A. 2009. Estimated use of water in the United States in 2005. U.S. Geological Survey Circular 1344, 52 p. Labbe, T. R. & Fausch, K. D. 2000. Dynamics of intermittent stream habitat regulate persistence of threatened fish at multiple scales. Ecological Applications 10:1774-1791. Leake, S. A., Pool, D. R., & Leenhouts, J. M. 2008. Simulated effects of ground-water withdrawals and artificial recharge on discharge to streams, springs, and riparian vegetation in the Sierra Vista Subwatershed of the Upper San Pedro Basin, southeastern Arizona. U.S. Geological Survey Scientific Investigations Report 2008-5207, 14p. LeBlanc, R. T., Brown, R. D., & FitzGibbon, J. E. 1997. Modeling the effects of land use change on the water temperature in unregulated urban streams. Journal of Environmental Management, 49(4), 445-469. Leopold, L. B. 1968. Hydrology for urban land planning: A guidebook on the hydrologic effects of urban land use. 153 Lindenmayer, D. B. & Luck, G. 2005. Synthesis: thresholds in conservation and management. Biological Conservation 124:351-354. Littell, R. C., Milliken, G. A., Stroup, W. W., Wolfinger, R. D., & Schabenberger, O. 2006. SAS for mixed models 2nd ed. SAS Inst., Cary, NC. Lyons, J., Zorn, T., Stewart, J., Seelbach, P., Wehrly, K., & Wang, L. 2009. Defining and characterizing coolwater streams and their fish assemblages in Michigan and Wisconsin, USA. North American Journal of Fisheries Management 29:1130-1151. Magnusson, J., Jonas, T., & Kirchner, J. W. 2012. Temperature dynamics of a proglacial stream: Identifying dominant energy balance components and inferring spatially integrated hydraulic geometry. Water Resources Research 48. Mantua, N., Tohver, I., & Hamlet, A. 2010. Climate change impacts on streamflow extremes and summertime stream temperature and their possible consequences for freshwater salmon habitat in Washington State. Climatic Change 102:187-223. Maupin, M. A., Kenny, J. F., Hutson, S. S., Lovelace, J. K., Barber, N. L., & Linsey, K. S. 2014. Estimated use of water in the United States in 2010. U.S. Geological Survey Circular 1405, 56 p. Meeus, J. H. 1991. Astronomical Algorithms. Willman-Bell, Incorporated. Michigan Legislature. 2006. Public Acts of 2006, Act No. 33: 93rd Legislature, Regular Session of 2006. Mohseni, O., Stefan, H. G., & Erickson, T. R. 1998. A nonlinear regression model for weekly stream temperatures. Water Resources Research 34:2685-2692. Mohseni, O. & Stefan, H. G. 1999. Stream temperature/air temperature relationship: a physical interpretation. Journal of Hydrology 218:128-141. Moisen, G. G. 2008. Classification and regression trees. In Jørgensen, S. E. & Fath, B. D.: Encyclopedia of Ecology, Vol. 1, p. 582-588. Elsevier, Oxford, UK. Moore, R. D., Macdonald, J. S., & Herunter, H. 2003. Downstream thermal recovery of headwater streams below cutblocks and logging roads. Forestry Impacts on Fish Habitat in the Northern Interior of British Columbia: A Compendium of Research from the Stuart-Takla Fish-Forestry Interaction Study, 179. Morin, G. & Couillard, D. 1990. Predicting river temperatures with a hydrological model. Encyclopedia of fluid mechanic, surface and groundwater flow phenomena 10:171-209. 154 Neumann, D. W., Rajagopalan, B., & Zagona, E. A. 2003. Regression model for daily maximum stream temperature. Journal of Environmental Engineering 129:667-674. Nuhfer, A. J. & Baker, E. A. 2004. A long-term field test of habitat change predicted by PHABSIM in relation to brook trout population dynamics during controlled flow reduction experiments. Michigan Department of Natural Resources, Fisheries Research Report 2068, Ann Arbor. Nuhfer, A. J., Zorn, T. G., & Willis, T. C. 2017. Effects of reduced summer flows on the brook trout population and temperatures of a groundwater-influenced stream. Ecology of Freshwater Fish 26:108-119. Otto, R. G. & Rice, J. O. 1977. Responses of a freshwater sculpin (Cottus cognatus gracilis) to temperature. Transactions of the American Fisheries Society. 106:89-94. Ozaki, N., Fukushima, T., Harasawa, H., Kojiri, T., Kawashima, K., & Ono, M. 2003. Statistical analysis on the effects of air temperature fluctuations on river water qualities. Hydrological Processes: An International Journal 17:2837-2853. Pilgrim, J. M. & Stefan, H. G. 1995. Correlations of Minnesota stream water temperatures with air temperatures. Project Report 382, St. Anthony Falls Laboratory, University of Minnesota, Minneapolis. Poole, G. C. & Berman, C. H. 2000. Pathways of human influence on water temperature dynamics in stream channels. Environmental Management, 20 p. Poole, G. C. & Berman, C. H. 2001. An ecological perspective on in-stream temperature: natural heat dynamics and mechanisms of human-caused thermal degradation. Environmental Management 27:787-802. Prasad, A. M., Iverson, L. R., Liaw, A. 2006. Newer classification and regression tree techniques: bagging and random forests for ecological prediction. Ecosystems 9:181-199. Qian, S. S. & Anderson, C. W. 1999. Exploring factors controlling the variability of pesticide concentrations in the Willamette River basin using tree-based models. Environmental Science & Technology 33:3332-3340. Qian, S. S., King, R. S., & Richardson, C. J. 2003. Two statistical methods for the detection of environmental thresholds. Ecological Modelling 166:87-97. Reimann, C., Filzmoser, P., & Garrett, R. G. 2005. Background and threshold: critical comparison of methods of determination. Science of the Total Environment 346:1-16. Richardson, C. J. & Qian, S. S. 1999. Long-term phosphorus assimilative capacity in freshwater wetlands: a new paradigm for sustaining ecosystem structure and function. Environmental Science & Technology 33:1545-1551. 155 Richter, B. D., Mathews, R., Harrison, D. L., & Wigington, R. 2003. Ecologically sustainable water management: managing river flows for ecological integrity. Ecological Applications 13:206-224. Risley, J. C., Constantz, J., Essaid, H., & Rounds, S. 2010. Effects of upstream dams versus groundwater pumping on stream temperature under varying climate conditions. Water Resources Research 46. Rutherford, J. C., Marsh, N. A., Davies, P. M., & Bunn, S. E. 2004. Effects of patchy shade on stream water temperature: how quickly do small streams heat and cool?. Marine and Freshwater Research 55:737-748. Seber, G. A. F. & Wild, C. J. 1989. Nonlinear regression: Wiley series in probability and mathematical statistics. Selong, J. H., McMahon, T. E., Zale, A. V., & Barrows, F. T. 2001. Effect of temperature on growth and survival of bull trout, with application of an improved method for determining thermal tolerances in fishes. Transactions of the American Fisheries Society 130:1026-1037. Shrode, J. B., Zerba, K. E., & Stephens Jr, J. S. 1982. Ecological significance of temperature tolerance and preference of some inshore California fishes. Transactions of the American Fisheries Society 111:45-51. Sinokrot, B. A. & Stefan, H. G. 1993. Stream temperature dynamics: measurements and modeling. Water Resources Research 29:2299-2312. Sinokrot, B. A. & Gulliver, J. S. 2000. In-stream flow impact on river water temperatures. Journal of Hydraulic Research 38:339-349. Steen, P. J., Zorn, T. G., Seelbach, P. W., & Schaeffer, J. S. 2008. Classification tree models for predicting distributions of Michigan stream fish from landscape variables. Transactions of the American Fisheries Society 137:976-996. Stefan, H. G. & Preud’homme, E. B. 1993. Stream temperature estimation from air temperature 1. Journal of the American Water Resources Association 29:27-45. Steinman, A. D., Nicholas, J. R., Seelbach, P. W., Allan, J. W., & Ruswick, F. 2011. Science as a fundamental framework for shaping policy discussions regarding the use of groundwater in the state of Michigan: a case study. Water Policy 13:69-86. Story, A., Moore, R. D., & Macdonald, J. S. 2003. Stream temperatures in two shaded reaches below cutblocks and logging roads: downstream cooling linked to subsurface hydrology. Canadian Journal of Forest Research 33:1383-1396. 156 Suding, K. N., Gross, K. L., & Houseman, G. R. 2004. Alternative states and positive feedbacks in restoration ecology. Trends in Ecology and Evolution 19:49-53. Team, R. D. C. 2003. R: a language and environment for statistical computing 1.8. 1. R Foundation for statistical computing (http://www.r-project.org), Vienna, Austria, + Pastecs library (http://www.sciviews.org/pastecs). Thayer, S. A., Taylor, W. W., Hayes, D. B., & Haas, R. C. 2007. Weight of evidence for underlying dynamics of yellow perch in Saginaw Bay, Lake Huron. Ecological Modelling 206:31-40. Tong, S. T. & Chen, W. 2002. Modeling the relationship between land use and surface water quality. Journal of Environmental Management 66:377-393. Webb, B. W. & Nobilis, F. 2007. Long-term changes in river temperature and the influence of climatic and hydrological factors. Hydrological Sciences Journal 52:74-85. Webb, B. W., Hannah, D. M., Moore, R. D., Brown, L. E., & Nobilis, F. 2008. Recent advances in stream and river temperature research. Hydrological Processes: An International Journal, 22:902-918. Wehrly, K. E., Wiley, M. J., & Seelbach, P. W. 2003. Classifying regional variation in thermal regime based on stream fish community patterns. Transactions of the American Fisheries Society 132:18-38. Wehrly, K. E., Wiley, M. J., & Seelbach, P. W. 2006. Influence of landscape features on summer water temperatures in lower Michigan streams. American Fisheries Society Symposium 48:113-127. Wehrly, K. E., Wang, L., & Mitro, M. 2007. Field-based estimates of thermal tolerance limits for trout: incorporating exposure time and temperature fluctuation. Transactions of the American Fisheries Society 136:365-374. Zorn, T. G., Seelbach, P. W., & Wiley, M. J. 2002. Distributions of stream fishes and their relationship to stream size and hydrology in Michigan’s Lower Peninsula. Transactions of the American Fisheries Society 131:70-85. Zorn, T. G., Seelbach, P. W., & Wiley, M. J. 2004. Utility of species-specific, multiple linear regression models for prediction of fish assemblages in rivers of Michigan’s Lower Peninsula. Michigan Department of Natural Resources, Fisheries Research Report 2072, Ann Arbor Zorn, T. G., Seelbach, P. W., Rutherford, E. S., Wills, T. C., Cheng, S. -T., & Wiley, M. J. 2008. A Regional-scale Habitat Suitability Model to Assess the Effects of Flow Reduction on Fish Assemblages in Michigan Streams. Michigan Department of Natural Resources, Fisheries Research Report 2089, Ann Arbor 157