A River Temperature Model to Assist Managers in Identifying Thermal Pollution Causes and Solutions

: Thermal pollution of rivers degrades water quality and ecosystem health, and cities can protect rivers by decreasing warmer impervious surface stormwater inﬂows and increasing cooler subsurface inﬂows and shading from riparian vegetation. This study develops the mechanistic i-Tree Cool River Model and tests if it can be used to identify likely causes and mitigation of thermal pollution. The model represents the impacts of external loads including solar radiation in the absence of riparian shade, multiple lateral storm sewer inﬂows, tributaries draining reservoirs, groundwater ﬂow, and hyporheic exchange ﬂow in dry weather steady ﬂows and wet weather unsteady ﬂows. The i-Tree Cool River Model estimates the shading e ﬀ ects of the riparian vegetation and other features as a function of heights and distances as well as solar geometry. The model was tested along 1500 m of a New York mountain river with a riparian forest and urban areas during 30 h with two summer storm events in 2007. The simulations were sensitive to the inﬂows of storm sewers, subsurface inﬂows, as well as riparian shading, and upstream boundary temperature inﬂows for steady and unsteady conditions. The model simulated hourly river temperature with an R 2 of 0.98; when shading was removed from the simulation the R 2 decreased 0.88, indicating the importance of riparian shading in river thermal modeling. When stormwater inﬂows were removed from the simulation, the R 2 decreased from 0.98 to 0.92, and when subsurface inﬂows were removed, the R 2 decreased to 0.94. The simulation of thermal loading is important to manage against pollution of rivers.


Introduction
Excessive river temperatures are detrimental to water quality and ecosystem health [1]. River warming can result from increased inflows of warm point and non-point source discharges, decreased inflows of cool sub-surface waters, removal of riparian shade, increased air temperatures, and changes in channel substrate and depth that increase absorption, conduction, and convection in heat transfer [2][3][4]. River thermal pollution is often associated with discharges of coolant water used by industry, but it is also associated with land-use change, including urbanization, river impoundment, channel management, and regulation [5,6]. River temperature is a critical water quality parameter for riverine systems, that affects the saturation of dissolved oxygen [7,8], kinetic reactions and resulting pollutant concentrations [9,10], and fish distribution, metabolism, growth, reproduction, and mortality [11,12]. Urbanization can elevate river temperatures through changes in riparian land cover which affects shade on the water surface, through river morphology which affects water depth, surface area, and velocity, and through flow connectivity with groundwater, stormwater, and other point and non-point source inflows [13][14][15][16][17][18]. When precipitation strikes hot impervious surfaces of urban areas, this generates warmer stormwater relative to the temperature of river water [19][20][21]. designed to allow for flexible shading factor algorithms, steady and unsteady flow, as well as other heat and mass transfer processes. The i-Tree Cool River Model is an open-source tool written in C++, and its package contains the routines and an executable file for running the code, which can be downloaded from http://www.itreetools.org/research_suite/coolriver. The model executable is called at the command line along with a configuration extensible markup language (XML) file, which includes the required initial information. The i-Tree Cool River Model C++ algorithms can be edited and recompiled with Visual Studio 2017 Community Edition or later, which is freeware. The simulation output includes the simulated river temperature and the heat fluxes.
The objectives of this paper are to present the theory of the i-Tree Cool River Model, to apply the model in a case study with unsteady stormwater inflows, and to evaluate the importance of the heat and mass transfer processes. To that end, following the model development, the manuscript provides a model testing to address the application of the model. The science questions are: When analyzing sources of thermal pollution, and possible mitigation scenarios, what is the relative contribution of (a) storm sewer inflows, (b) subsurface inflows of groundwater and hyporheic exchange, (c) riparian shading and weather, on the accuracy of simulated river temperature?

Heat Flux Formulation
The i-Tree Cool River Model simulates an advection-dispersion equation with inflows and heat fluxes following Martin et al. [32] ∂T w ∂t where T w is the cross-sectional averaged river temperature ( • C), t is time (s), U is the reach average flow velocity (m/s), x is river distance (m), D L is the dispersion coefficient (m 2 /s), R h is the heat flux reaction term, also known as heat transfer [3,27], and R i is the reaction term of the external inflows. When R i is combined with the advection and dispersion terms in Equation (1), they are collectively referred to as mass transfer [3,27]. The R h and R i are defined as R i = Q W T W + Q GW T GW + Q Hyp T Hyp + Q SS T SS Q i + Q GW + Q Hyp + Q SS − T W,t−1 (3) where Φ net is the net exchange of thermal energy (W/m 2 ), ρ is the water density (kg/m 3 ), C p is the specific heat capacity of water (J/kg • C), y is the average water column depth (m), Q is discharge (m 3 /s), T is water temperature ( • C), and subscripts W is the river flow, GW is groundwater flow, Hyp is hyporheic exchange, and SS is stormwater inflow, t − 1 is prior time step. River velocities, dispersion, and inflows are calculated using standard methods, described in Supplementary Materials Section S1. The subsurface inflows distinguish between hyporheic and groundwater inflows due to their different environmental processes and use a separate mathematical formulation for each term. For surface inflows, users can include tributaries in place of storm sewers, and assign an unlimited number of surface inflows for each cross section. The net exchange of thermal energy is defined as in Boyd et al. [27] as where the Φ is the heat flux (W/m 2 ), and subscripts net is the net heat flux at the water surface, longwave is the longwave radiation flux at the water surface, shortwave is the shortwave radiation at the water surface, latent is the latent heat flux from evaporation, sensible is the sensible heat flux representing the convective thermal flux from the water surface, and sediment is the bed sediment heat flux representing conduction forcing at the water column interface. The longwave radiation flux in Equation (4) is composed of positive downward fluxes from the atmosphere and land cover over the water surface, and a negative upward flux from the waterbody to the air, following the approach of Boyd et al. [27] where Φ atmospheric longwave is the atmospheric flux (W/m 2 ), Φ landcover longwave is the land cover flux (W/m 2 ), and Φ back longwave is the back-to-air flux (W/m 2 ). Atmospheric longwave radiation is a function of air temperature and exposure from the river surface to the atmosphere, called the view-to-sky factor (f ), calculated using Boyd et al.
where T air is air temperature ( • C), the ε atm is the emissivity of the atmosphere (0 to 1), σ is the Stefan-Boltzmann constant (5.6696 × 10 −8 , W/m 2 K 4 ), and min(f 1 , f 2 , f 3 ) is the minimum of the three view-to-sky factors (0 to 1), where f 1 represents building effects, f 2 represents vegetation effects, and f 3 represents topographic effects ( Figure 1). The emissivity of the atmosphere ε atm is calculated using [33] ε atm = 1.72( 0.1e a T air + 273.2 ) where e a is the actual vapor pressure (mbar), and C L is the cloudiness, which ranges from 0 for a clear sky to 1 for full cloud cover [34].
convective thermal flux from the water surface, and sediment is the bed sediment heat flux representing conduction forcing at the water column interface. The longwave radiation flux in Equation (4) is composed of positive downward fluxes from the atmosphere and land cover over the water surface, and a negative upward flux from the waterbody to the air, following the approach of Boyd et al. [27] where Tair is air temperature (°C), the atm ε is the emissivity of the atmosphere (0 to 1), σ is the Stefan-Boltzmann constant (5.6696 × 10 −8 , W/m 2 K 4 ), and min(f1, f2, f3) is the minimum of the three view-to-sky factors (0 to 1), where f1 represents building effects, f2 represents vegetation effects, and f3 represents topographic effects ( Figure 1). The emissivity of the atmosphere atm ε is calculated using [33] where ea is the actual vapor pressure (mbar), and CL is the cloudiness, which ranges from 0 for a clear sky to 1 for full cloud cover [34].
The view-to-sky factors value of 1 indicates a full unobstructed sky view [27,35,36]. The general sky-view-factor formula for fi is computed for each cross-section based on Chen et al. [14] 2 1 Figure 1. Shading of the river surface. A cross-sectional view, in which BSA is the building shading angle, VSA is the vegetation shading angle, and TSA is the topographic shading angle. h building , h tree , and h bank are building, vegetation, and bank heights respectively. D building , D canopy , and D bank are building to the bank, canopy to the bank, and bank.
The view-to-sky factors value of 1 indicates a full unobstructed sky view [27,35,36]. The general sky-view-factor formula for f i is computed for each cross-section based on Chen et al. [14] where i indicates the object at that cross-section, where 1 = building, 2 = vegetation, or 3 = topography; and SA is the shade angle (radians), computed as and h c is the combined height of the objects above the water (e.g., if a tree is set on a hill, h c = h tree + h bank ), and max (D i ) is the maximum distance from all objects at that cross-section to the edge of the water.
The land cover longwave radiation in Equation (5) also uses the view-to-sky factors. The land cover radiation represents the land cover, e.g., vegetation such as trees' influence on water temperature, and the model by default sets land cover temperature equal to atmospheric temperature, following the approach of Boyd et al. [27] Φ landcover The waterbody to air radiation term in Equation (5) is a function of water temperature, representing heat flux emitted from the water surface, following the approach of Boyd et al. [27] Φ back longwave = −0.96σ(T w + 273.2) 4 (11) where the T W is the river temperature ( • C). See the Supplementary Materials, Section S2 for the methods used to find the remaining right-hand side terms in Equation (4), which are short wave radiation, latent heat flux, sensible heat flux, and bed sediment heat flux. Table S1 lists the 10 input files required by i-Tree Cool River, and names and describes the parameters in each of the files.

Study Area and Model Inputs
The i-Tree Cool River Model's accuracy in representing thermal loading was tested in unsteady state (i.e., wet weather) using unpublished data from 11 to 12 June 2007 for a 1500 m reach of Sawmill Creek, in Tannersville, New York (42.1955 N, 74.1339 W, WGS84). Sawmill Creek is a second-order mountainous river with varying watershed land use, starting in forests and transitioning to urban land. At the end of the Sawmill Creek study reach, the time of travel was approximately 30 min and the upstream watershed area is 8.16 km 2 , which includes a nested urban watershed of 1.8 km 2 draining to the river in storm sewers (Figure 2a,b). Sawmill Creek flow at the upstream boundary was estimated using stage-discharge relations, monitoring stage with pressure transducers (manufactured by Global Water Instruments) at the upstream and downstream stations (0 m and 1500 m respectively) and in storm sewers. Stage was converted to discharge using the Manning equation, with stage converted to channel area and hydraulic radius using geometry relations, and the Manning roughness coefficients estimated from pebble counts Wolman [37] at each cross section by Crispell et al. [29]. We installed compound weir plates in the sewers and used the weir plate manufacturer equations to convert stage to discharge for the storm sewer inflows. Observed storm flows were corroborated with simulated flows by the i-Tree Hydro model (Yang,et al. [38]), calibrated to match the estimated baseflow.
Rainfall occurred twice during 12 June 2007 the first time with a 2 h duration totaling 3.3 mm and the second time with a 3 h duration totaling 8.4 mm. The storm sewer inflows were active during dry and wet weather, in dry weather due to illicit connections draining buildings, and in wet weather due to storm runoff. Crispell et al. [29] monitored river temperatures and storm sewer drainage temperatures at 12 river stations in the Sawmill Creek reach and the two inflow locations at 10-minute intervals using ibutton temperature data loggers, which were used to set boundary conditions. The temperature monitoring ibuttons in Sawmill Creek reach were strategically placed and considered representative of the reach temperature, capturing the influence of stormwater inflows after they had distance to mix with the channel water [29]. In the upstream section of the reach, between the cross section at station 0 m and a station at 600 m, the observed average river temperature increased by 0.03 • C per 100 m (3.3%), and between the cross section at station 900 m and a station 1500 m, temperature increased by 0.008 • C per 100 m (1%). In the middle section of the reach, between the cross section at station 600 m and station 900 m, the temperature increased by 0.1 • C per 100 m, three times the rate of the upstream reach, an increase attributed to the warming effect of the Tannersville's storm sewer inflow (Table S2).
The i-Tree Cool River was also tested in steady state mode, e.g., no rainfall events, for a 475 m reach of Meadowbrook Creek (43.0306 N, 76.0680 W, WGS84), a first order and urbanized river in the city of Syracuse, New York (Figure 3a,b). Flow at the upstream boundary, cross section survey data, and river temperature at 30 monitoring locations at 5-minute intervals were provided by Glose et al. [31], who used these data to develop HFlux. We simulated the 5-day period of 13-19 June 2012. times the rate of the upstream reach, an increase attributed to the warming effect of the Tannersville's storm sewer inflow (Table S2).
The i-Tree Cool River was also tested in steady state mode, e.g., no rainfall events, for a 475 m reach of Meadowbrook Creek (43.0306 N, 76.0680 W, WGS84), a first order and urbanized river in the city of Syracuse, New York (Figure 3a,b). Flow at the upstream boundary, cross section survey data, and river temperature at 30 monitoring locations at 5-minute intervals were provided by Glose et al. [31], who used these data to develop HFlux. We simulated the 5-day period of 13-19 June 2012.   The i-Tree Cool River Model simulated Sawmill Creek using input data from multiple sources. Specification of hourly weather data, including air temperature T air , dew point temperature T dew , shortwave radiation S in , fraction cloudiness C, and wind speed U wind were obtained from the National Solar Radiation Data Base (NSRDB) station ID#1227776, located 23 km from the study site. The NSRDB provides satellite estimated surface radiation at 30 min intervals. The shading factor, SF, in Equation (S2) and view-to-sky coefficients, f in Equation (S4) were estimated at 1 m intervals along Sawmill Creek using the TTools algorithm from observations of riparian vegetation and aerial images of the study area [29]. The river base width and bank slope were obtained from field surveys at each cross section [29], which defined the irregular pattern of river widening and narrowing. The simulated Sawmill Creek reach was delineated into 15 segments considering the locations where the temperature was observed, with segment lengths no greater than 100 m, which resulted in a simulation timestep of 0.5 seconds to satisfy the i-Tree Cool River Model stability criteria. The simulations represented the observed conditions, as well as alternative scenarios to determine model sensitivity to shading, subsurface inflow, and the calculated upstream boundary condition, which are sometimes difficult to obtain. Our calculated upstream boundary condition was derived with Mohseni et al. [30] Equation (12), a non-linear regression between air temperature and river temperature.
In the equation, the coefficient α is the estimated maximum stream temperature, γ is a measure of the steepest slope of the function, and β represents the air temperature at the inflection point.   Our observed upstream boundary condition was obtained from ibutton thermistor measurements. We analyzed the simulated and observed river temperatures for each of the cross sections averaged with respect to time to obtain a 30-hour average at each of the 12 cross sections. Simulations were written hourly for each cross section, which can be written at any timesteps for each meter of the river. We ran this simulation using Equations (5) to (11) for longwave radiation, Equation (S2) for shortwave radiation, Equation (S11) for latent heat, Equation (S12) for sensible heat, and Equation (S15) for the sediment heat.  The i-Tree Cool River model simulated hourly water temperatures were not significantly different than the observed, for reach averaged data, based on the p-values calculated using a pairedsamples t-test and the α = 0.05 (See Table S3 for more details). The 30-hour average simulated and observed river temperature, along the entire 1500 m Sawmill Creek reach, increased by 0.4 °C at a slope of approximately 0.02 °C per 100 m, but with longitudinal variation in that slope. For the 1500 m reach, the model had a root-mean-square error (RMSE) of 0.03 °C and a coefficient of determination (R 2 ) of 0.98 with a p-value of 0.87 which was greater than the α of 0.05. The model simulated the relatively rapid increase in water temperature recorded by the sensors, between cross section 600 m and 900 m, corresponding to the reach with storm sewer inflow. This relatively rapid increase in temperature leveled at station 900 m, which is the first station downstream of the last storm sewer outfall. Relatively warm water in the Tannersville storm sewer entering Sawmill Creek between cross sections at station 600 m and at station 900 m was a major driver of the i-Tree Cool River Model forecasting a rise in river temperature during both wet and dry weather conditions ( Figure S1a). During the wet weather, a total of 7 hours, the rate that simulated river temperature increased from station 600 m to station 900 m at a rate of 0.32 °C per 100 m ( Figure S1b), much steeper than during dry weather. During dry weather, the simulated river temperature from station 600 m to station 900 m increased at a rate of 0.04 °C per 100 m ( Figure S1c), approximately 12% of the wet weather slope. The i-Tree Cool River algorithms for shading, net groundwater discharge, hyporheic exchange, and upstream boundary condition temperature influenced the simulation of longitudinal river temperature and the model goodness of fit. In all of the scenarios, the calculated paired t-test p-values were smaller than the α = 0.05, rejecting the null hypothesis, H0 of a significant difference between the means of the simulated and observed reach averaged river temperatures (See Table S4 for more details). When the shading algorithm was disabled, i.e., no shade was simulated, the i-Tree Cool River Model overestimated the river temperature for all cross sections by 0.34 °C, at a rate of 0.02 °C per 100 m, for the 30-hour period, 11 to 12 June 2007 ( Figure S2), and the model RMSE increased to 0.36 °C and the R 2 decreased to 0.88.

Model Evaluation
Diurnal sinusoidal patterns of simulated and observed river warming and cooling were driven by the heat balance but disrupted by abrupt pulses of inflow due to warm runoff during the two storm events on 11   The i-Tree Cool River model simulated hourly water temperatures were not significantly different than the observed, for reach averaged data, based on the p-values calculated using a paired-samples t-test and the α = 0.05 (See Table S3 for more details). The 30-hour average simulated and observed river temperature, along the entire 1500 m Sawmill Creek reach, increased by 0.4 • C at a slope of approximately 0.02 • C per 100 m, but with longitudinal variation in that slope. For the 1500 m reach, the model had a root-mean-square error (RMSE) of 0.03 • C and a coefficient of determination (R 2 ) of 0.98 with a p-value of 0.87 which was greater than the α of 0.05. The model simulated the relatively rapid increase in water temperature recorded by the sensors, between cross section 600 m and 900 m, corresponding to the reach with storm sewer inflow. This relatively rapid increase in temperature leveled at station 900 m, which is the first station downstream of the last storm sewer outfall. Relatively warm water in the Tannersville storm sewer entering Sawmill Creek between cross sections at station 600 m and at station 900 m was a major driver of the i-Tree Cool River Model forecasting a rise in river temperature during both wet and dry weather conditions ( Figure S1a). During the wet weather, a total of 7 hours, the rate that simulated river temperature increased from station 600 m to station 900 m at a rate of 0.32 • C per 100 m ( Figure S1b), much steeper than during dry weather. During dry weather, the simulated river temperature from station 600 m to station 900 m increased at a rate of 0.04 • C per 100 m ( Figure S1c), approximately 12% of the wet weather slope. The i-Tree Cool River algorithms for shading, net groundwater discharge, hyporheic exchange, and upstream boundary condition temperature influenced the simulation of longitudinal river temperature and the model goodness of fit. In all of the scenarios, the calculated paired t-test p-values were smaller than the α = 0.05, rejecting the null hypothesis, H 0 of a significant difference between the means of the simulated and observed reach averaged river temperatures (See Table S4 for more details). When the shading algorithm was disabled, i.e., no shade was simulated, the i-Tree Cool River Model overestimated the river temperature for all cross sections by 0.34 • C, at a rate of 0.02 • C per 100 m, for the 30-hour period, 11 to 12 June 2007 ( Figure S2), and the model RMSE increased to 0.36 • C and the R 2 decreased to 0.88.
Diurnal sinusoidal patterns of simulated and observed river warming and cooling were driven by the heat balance but disrupted by abrupt pulses of inflow due to warm runoff during the two storm events on 11 and 12 June 2007. The Sawmill Creek mean temperature, the average of measurements at the 12 cross sections, diurnally peaked at 15.8 • C by 15:00 h June 11, 2007 (Figure 5a), two hours after the peak in shortwave radiation (Figure 5b). By 20:00 hours, shortwave radiation has declined to 0, net radiation became negative, and river temperature has decreased from the peak of 15.8 • C to 15.1 • C. A storm event at 21:00 h 11 June 2007, and then again at 01:00 h of 12 June 2007, generates inflow of warmer water from the upstream and storm sewer, creating a temporary increase in temperature, which disrupts the sinusoidal pattern in cooling toward the diurnal minimum temperature at 05:00 h of 12 June 2007. Dry weather extends through the remainder of the simulation, and at 06:00 h of 12 June 2007, the increasing shortwave and thus net radiation reestablish heat flux as the main driver of the increasing river temperature, which peaks at 16.4 • C at 14:00 h. There was no significant difference between the simulated and observed time averaged river temperature datasets based on the paired-samples t-test and α = 0.05 (See Table S5 for more details). The model simulations of the abrupt pulses in river temperature during the wet weather, extending from hour 20 of day 1 to hour 3 of day 2, had a Nash-Sutcliffe Efficiency (NSE) coefficient of 0.9 and a p-value of 0.80 (Figure 5a). The magnitude of the simulated river temperature changes due to the inflow of stormwater from the Tannersville storm sewer system was 0.3 • C during the first storm and 0.4 • C during the second storm. The model simulations had their poorest fit with observed river temperatures during a 6 h period on 12 June 2007, between 03:00 and 09:00 h, centered at sunrise, when it overestimated the river temperature by an average of 0.13 • C.
has declined to 0, net radiation became negative, and river temperature has decreased from the peak of 15.8 °C to 15.1 °C. A storm event at 21:00 h 11 June 2007, and then again at 01:00 h of 12 June 2007, generates inflow of warmer water from the upstream and storm sewer, creating a temporary increase in temperature, which disrupts the sinusoidal pattern in cooling toward the diurnal minimum temperature at 05:00 h of 12 June 2007. Dry weather extends through the remainder of the simulation, and at 06:00 h of 12 June 2007, the increasing shortwave and thus net radiation reestablish heat flux as the main driver of the increasing river temperature, which peaks at 16.4 °C at 14:00 h. There was no significant difference between the simulated and observed time averaged river temperature datasets based on the paired-samples t-test and α = 0.05 (See Table S5 for more details). The model simulations of the abrupt pulses in river temperature during the wet weather, extending from hour 20 of day 1 to hour 3 of day 2, had a Nash-Sutcliffe Efficiency (NSE) coefficient of 0.9 and a p-value of 0.80 (Figure 5a). The magnitude of the simulated river temperature changes due to the inflow of stormwater from the Tannersville storm sewer system was 0.3 °C during the first storm and 0.4 °C during the second storm. The model simulations had their poorest fit with observed river temperatures during a 6 h period on 12 June 2007, between 03:00 and 09:00 h, centered at sunrise, when it overestimated the river temperature by an average of 0.13 °C. The river temperature simulated by the i-Tree Cool River Model was a function of spatially and temporally varying contributions of groundwater, hyporheic exchange, and storm sewer inflow. Analysis of these components to thermal loading can assist in developing pollution mitigation or river restoration scenarios. In cross sections without storm sewer inflow and in the absence of rain events, the river temperature was predominantly determined by river flow from the upstream reach, and groundwater only contributed approximately 1%, while hyporheic exchange contributed approximately 10% (Figure 6a; 350 m, and 1100 m). In cross sections with storm sewers, even in the absence of rain events, when the storm sewers discharged flow from illicit connections, they The river temperature simulated by the i-Tree Cool River Model was a function of spatially and temporally varying contributions of groundwater, hyporheic exchange, and storm sewer inflow. Analysis of these components to thermal loading can assist in developing pollution mitigation or river restoration scenarios. In cross sections without storm sewer inflow and in the absence of rain events, the river temperature was predominantly determined by river flow from the upstream reach, and groundwater only contributed approximately 1%, while hyporheic exchange contributed approximately 10% (Figure 6a; 350 m, and 1100 m). In cross sections with storm sewers, even in the absence of rain events, when the storm sewers discharged flow from illicit connections, they contributed approximately 25% of the flow influencing the cross section river temperature (Figure 6a; 800 m), which warmed the river water (see Table S2 reflecting warmer average temperature of the storm sewer in the dry weather). During wet weather, there was inflow from the storm sewer due to the flows from impervious areas, and this inflow contributed approximately 50% of the flow thereby influencing the river temperature (Figure 6b, 800 m) and provided the thermal load to the river system. Integrated along the 1500 m of Sawmill Creek reach, the contribution of groundwater summed to 15% of the total river volume, while the hyporheic exchange, which flows in and out within each sub-reach associated with a cross section, averaged approximately 10% of inflow in each sub-reach.
contributed approximately 25% of the flow influencing the cross section river temperature ( Figure  6a; 800 m), which warmed the river water (see Table S2 reflecting warmer average temperature of the storm sewer in the dry weather). During wet weather, there was inflow from the storm sewer due to the flows from impervious areas, and this inflow contributed approximately 50% of the flow thereby influencing the river temperature (Figure 6b, 800 m) and provided the thermal load to the river system. Integrated along the 1500 m of Sawmill Creek reach, the contribution of groundwater summed to 15% of the total river volume, while the hyporheic exchange, which flows in and out within each sub-reach associated with a cross section, averaged approximately 10% of inflow in each sub-reach. In addition to the analysis of unsteady simulations in Sawmill Creek, the i-Tree Cool River Model performance was analyzed for the steady state condition in Meadowbrook Creek for [13][14][15][16][17][18][19] June 2012. The i-Tree Cool River Model simulated the time averaged river temperatures at 30 cross sections with an RMSE of 0.2 °C. We combined these 30 cross sections into reach averaged river temperature data to examine the diurnal pattern driven by the heat balance (Figure 7a). There was no significant difference between the simulated and observed reach averaged river temperatures based on a paired-samples t-test and α = 0.05 (See Table S6   In addition to the analysis of unsteady simulations in Sawmill Creek, the i-Tree Cool River Model performance was analyzed for the steady state condition in Meadowbrook Creek for 13-19 June 2012. The i-Tree Cool River Model simulated the time averaged river temperatures at 30 cross sections with an RMSE of 0.2 • C. We combined these 30 cross sections into reach averaged river temperature data to examine the diurnal pattern driven by the heat balance (Figure 7a). There was no significant difference between the simulated and observed reach averaged river temperatures based on a paired-samples t-test and α = 0.05 (See Table S6

Sensitivity Analysis
A sensitivity analysis of the i-Tree Cool River Model examined the fluctuation in simulated temperature with changes in input data to identify the most sensitive parameters, which is useful when considering impacts of environmental change. The sensitivity analysis was performed for steady and unsteady simulations. We used global sensitivity analysis to identify the most important parameters and coordinated this analysis with that for the Meadowbrook Creek reach in summertime, by Glose et al. [31], noting both models are based on Heat Source [27]. Glose et al. [31] used an observed boundary conditions and our equations 1, 2, 4, 5, 6, 12, 21, and 22, and varied discharge by ±10%, groundwater temperature ±15%, varied shading factor and view-to-sky factor by ±20%. They identified groundwater as the most sensitive parameter, with a ±0.2 °C change on average stream temperature. We replicated this sensitivity analysis to the 1500 m Sawmill Creek reach, confirming these sensitivities. We then extended the analysis in the 1500 m Sawmill Creek reach to consider varying parameters of storm sewer temperature, sediment temperature, and recorded boundary conditions temperature by ±15% ( Figure S3) and varying parameters of substrate hydraulic conductivity (SHC), cloudiness factor (Cl), and groundwater discharge (GW) by ±20% ( Figure S4). When storm sewer temperature was varied by ±15%, the reach-averaged temperature changed by 1.65% (0.27 °C). When sediment temperature was varied by ±15%, the reach-averaged temperature changed by 0.3%. When upstream boundary conditions temperature was varied by ±15%, the reachaveraged temperature changed by 9.5%. When substrate hydraulic conductivity was varied by ±20%, the reach-averaged temperature changed by 0.15%. When cloudiness factor was varied by ±20%, the reach-averaged temperature changes by 0.02%. When groundwater discharge factor was varied by ±20%, the reach-averaged temperature changes by 0.03%. Based on this analysis, the most sensitive

Sensitivity Analysis
A sensitivity analysis of the i-Tree Cool River Model examined the fluctuation in simulated temperature with changes in input data to identify the most sensitive parameters, which is useful when considering impacts of environmental change. The sensitivity analysis was performed for steady and unsteady simulations. We used global sensitivity analysis to identify the most important parameters and coordinated this analysis with that for the Meadowbrook Creek reach in summertime, by Glose et al. [31], noting both models are based on Heat Source [27]. Glose et al. [31] used an observed boundary conditions and our Equations (1), (2), (4)-(6), (10), (11), (S2), (S11) and (S12), and varied discharge by ±10%, groundwater temperature ±15%, varied shading factor and view-to-sky factor by ±20%. They identified groundwater as the most sensitive parameter, with a ±0.2 • C change on average stream temperature. We replicated this sensitivity analysis to the 1500 m Sawmill Creek reach, confirming these sensitivities. We then extended the analysis in the 1500 m Sawmill Creek reach to consider varying parameters of storm sewer temperature, sediment temperature, and recorded boundary conditions temperature by ±15% ( Figure S3) and varying parameters of substrate hydraulic conductivity (SHC), cloudiness factor (Cl), and groundwater discharge (GW) by ±20% ( Figure S4). When storm sewer temperature was varied by ±15%, the reach-averaged temperature changed by 1.65% (0.27 • C). When sediment temperature was varied by ±15%, the reach-averaged temperature changed by 0.3%. When upstream boundary conditions temperature was varied by ±15%, the reach-averaged temperature changed by 9.5%. When substrate hydraulic conductivity was varied by ±20%, the reach-averaged temperature changed by 0.15%. When cloudiness factor was varied by ±20%, the reach-averaged temperature changes by 0.02%. When groundwater discharge factor was varied by ±20%, the reach-averaged temperature changes by 0.03%. Based on this analysis, the most sensitive model parameters, ranked in order of importance, are upstream boundary conditions, storm sewer temperature, sediment temperature, substrate hydraulic conductivity, groundwater discharge, and cloudiness (additional sensitivity analysis is presented in supplementary materials Figures S3-S6).

Discussion
The i-Tree Cool River Model simulated the warming effects of the many potential sources of thermal pollution, including radiation fluxes and urban runoff, both dry weather illicit connections and wet weather stormwater. The model also shows how groundwater and hyporheic exchange inflows can provide a cooling effect, providing a comprehensive approach to assessing and perhaps mitigating thermal loading. To determine which factors are most effective in such management, this discussion provides some perspective on the effects of each warming and cooling effect.
The impact of urban runoff on the average temperature was rapid, within 1 hour of the onset of precipitation, and caused a temperature increase of 0.3 • C for the first storm, and 0.4 • C for the second storm of 12 June 2007. The rapid and large change in river temperature can be attributed to the short duration event, which Herb et al. [1] suggest when rainfalls only last 2 to 3 h will make the largest impact on raising stormwater and water temperature. The urban storm sewer area was approximately 21% of the watershed drainage area and had 35% impervious cover, which contributed to a relatively large volume of flashy, warm, stormwater response. Relative differences between air and water temperatures contribute to the warming, as noted by Herb et al. [1]; for the Sawmill Creek reach the average air temperature was 21.2 • C and average dew point temperature was 18.0 • C, both warmer than the average river temperature which was 15.2 • C. Even though the two rainfall events occurred at night during 12 June 2007, when solar radiation was not present, the prior day averaged 20% cloud coverage, allowing 80% of summer shortwave radiation to reach the small albedo impervious surface.
Simulation of the effects of the nighttime stormwater thermal load in riverine receiving waters on Sawmill Creek contributes additional data and tools to the investigation and management of the urban heat island. A common signature of the urban heat island is elevated nighttime air temperatures in urban areas relative to rural areas, due to physical differences affecting solar heating, such as albedo and thermal capacity, and anthropogenic heat sources [39]. Daytime insulation is a common driver of thermal loading of receiving waters [40], but for 12 June in Tannersville, the daytime solar heating of impervious area did not cause thermal loading of the river until the nighttime wet weather event. The nighttime precipitation landed on warm impervious surfaces, retaining much of their daytime elevated surface temperatures due to high capacitance, and this surface warmth was conducted into the stormwater entering the relatively cool, rural origin, receiving water. The effect of urban heat islands on rivers was studied by Somers et al. [41], who noted a 1.6 • C higher warm season temperature in urban rivers than forested rivers, and 8 • C greater spatial variation in urban rivers than in rural river temperatures along a 1 km transect. During a daytime storm event affecting all rivers, the temperatures in urban rivers rose as much as 4 • C, compared with a negligible rise in temperature in the forested rivers [41]. Nelson et al. [26] forecasted the thermal impact of individual storm events and found storm-induced river temperature surged by 3.5 • C for the warm season in urban watersheds near Washington, DC, USA; with drainage areas averaging 8 km 2 .
Proper simulation with the i-Tree Cool River Model of unsteady flows and their thermal pollution of receiving waters requires consideration of model goals and limitations. Typical model goals are either model inter-comparison for contrasting scenarios, such as varying impervious or tree cover, or model simulation for hindcasting or forecasting. In cases of model inter-comparison, the model has fewer limitations and the model physics will allow users to consider changes in river temperature for changes in study site conditions; model simulation has accuracy constrained by the accuracy of inputs as well as a model epistemic error [31,42]. This project attempted to improve accuracy of model simulations in Sawmill Creek by obtaining accurate data of the storm volumes and temperatures entering at the upstream and storm sewer locations along the boundary, using ibutton sensors, which are widely used for river temperature monitoring [31,43,44]. In cases where point or diffuse sources enter at multiple, unspecified locations along the river channel, such as with groundwater seeps, ibuttons may be inefficient and a better monitoring approach may involve using distributed temperature sensing system [8,35] for high-frequency time series, or from forward-looking infrared radar [45] temporally coarser data. An alternative to monitoring storm sewer inflow temperatures is estimating those temperatures using models of impervious runoff [25,29], and upstream boundary conditions can be estimated using air temperature records with the Mohseni et al. [30] non-linear regression.
Groundwater and hyporheic exchange were significant factors of temperature regulation during wet and dry weather. The section of Sawmill Creek simulated by the i-Tree Cool River Model had groundwater flow rates of approximately 0.0024 m 3 /s per 100 m, or 1% of flow, and riverbed longitudinal slopes that generated 10% contributions of hyporheic exchange ( Figure 6). This combined subsurface flow, through the model inflow routines, contributed a cooling effect for the simulated summer period in 11-12 June 2007. Removing these inflows from the simulation caused the model to achieve RMSE of 0.18 • C and R 2 of 0.94, less important than the cooling through shading and the heat flux routines, when removed generated a RMSE of 0.36 • C and R 2 of 0.88. In winter, when river temperature is typically below subsurface water temperatures, this inflow would likely contribute a warming effect, as observed in other rivers by Risley et al. [46] and Kurylyk et al. [47]. From a survey of other studies, the relative contribution of groundwater and hyporheic exchange inflow with river water varies by site conditions and time. Poole et al. [48] working in mountain rivers, with bed slopes above 2%, also found the surface water received a larger volumetric inflow from hyporheic exchange than from groundwater, while Glose et al. [31] working in valleys with approximately 1% slopes did not identify significant hyporheic exchange and set groundwater as the only subsurface source of inflow.
Riparian shading from tree canopy, hillslope, and buildings provided the only land-based reduction in shortwave radiation and the view to the sky for the river, which influenced the longwave radiation. We used model inter-comparison simulations to contrast a scenario with and without shading and determined shading cooled river temperatures by an average of 0.34 • C during the 30 hours period. The landscape contribution to shading varied longitudinally along the reach, and at cross-section 9 was primarily from building shade, while upstream at cross sections 1 to 5 was primarily forest; hillslope topography provided minimal shading at this site. Shading is a concern in river thermal loading, and shallow and slow moving water is more vulnerable to such warming, and others have modeled this effect. Sun et al. [18] simulated 6 years along 6 separate reaches ranging in length from 85 to 1185 m of Mercer Creek in Washington State, and determined that tree and hillslope shading reduced the annual maximum temperatures by 4 • C. Roth et al. [49] simulated three cloud-free summer days in August 2007 along a 1260 m section of the Boiron de Morges River in southwest Switzerland, and determined riparian shading, by decreasing shortwave radiation, decreased daily average water temperatures by 0.7 • C. Guoyuan et al. [50] demonstrated predictions of shade from riparian vegetation (e.g., the Chen et al. [14] method used in i-Tree Cool River) were sensitive to the interaction of river azimuth and latitude, with E-W rivers in low latitudes benefiting least from riparian shade. Lee et al. [51] recommended for effective reduction in shortwave radiation, riparian areas utilize shading angles of 70 • (1.22 radian) and view-to-sky factors smaller than or equal to 0.22. In the 1500 m Sawmill Creek reach, less than 10% of the view-to-sky factors were smaller or equal to 0.22 and shading angles averaged 50 • , and therefore additional thermal management opportunities are present.
The i-Tree Cool River Model was designed to assist river managers assess mechanistic causes of thermal pollution using free, open source, relatively simple algorithms in order to negotiate the balance between complexity and accuracy. While the model requires several input files, many of these can be obtained from publicly available data, site surveys, or estimation approaches; for model inter-comparison studies the accuracy of input data become less critical than in forecast simulations. The number of input files required by the model is comparable to other mechanistic models simulating river temperature, such as HFlux, HSPF, and QUAL2K, which require approximately 25 to 40 parameters, spatial data of river geometry and riparian features, and time series data describing the weather and discharge. Obtaining these inputs is a potential limitation of the i-Tree Cool River Model, and methods to obtain or estimate these input files are discussed above.

Conclusions
In this study, we developed the one-dimensional mechanistic i-Tree Cool River Model to simulate river temperature considering a combination of advection, dispersion, heat flux, and inflow processes. The i-Tree Cool River Model has the ability to analyze the impacts of external loads including multiple lateral storm sewer inflows, groundwater flow, and hyporheic exchange flow in steady and unsteady flows. The i-Tree Cool River Model estimates the shading effects of the riparian vegetation and other features as a function of heights and distances as well as solar geometry. The model performance was tested in steady and unsteady modes for the Meadowbrook reach in Syracuse, New York and Sawmill Creek in Tannersville, New York, respectively. The i-Tree Cool River Model performed satisfactorily in both simulations. The model can be used to conduct thermal pollution analysis of urban areas and investigate land cover and hydrology-based mitigation methods. The simulated river temperature of the i-Tree Cool River Model can be used for other environmental models, such as urban development models, atmospheric models, climate change models, and hydrology models.
Supplementary Materials: The following are available online at http://www.mdpi.com/2073-4441/11/5/1060/s1, Table S1: List of the input files required for the simulation process of the i-Tree Cool River Model, Table S2: Observed water temperatures for three reaches of Sawmill Creek and for the Tannersville storm sewer, during 11 and 12 June 2007, as the average for all time steps during the dry or wet weather conditions, Table S3: Statistical analysis (paired t-test) of the reach averaged observed and simulated river temperature in Sawmill Creek for the (a) original condition including both wet and dry weather (b) wet weather, and (c) dry weather, Table S4: Statistical analysis (paired t-test) of the reach averaged observed and simulated river temperatures using the scenarios for the (a) no shading effect, (b) no groundwater and hyporheic exchange inflows, and (c) calculated boundary condition, Table S5: Statistical analysis (paired t-test) of the observed and simulated river temperatures in Sawmill Creek, between 12:00 h of 11 Jun 2007 to 17:00 h of 12 June 2007, Table S6: Statistical analysis (paired t-test) of the observed and simulated river temperatures in Meadowbrook Creek, for 13-19 June 2012, Figure S1: Time averaged observed and simulated river temperature in Sawmill Creek for the (a) original condition including both wet and dry weather (b) wet weather, and (c) dry weather, Figure S2: Time averaged observed and simulated river temperatures using the scenarios for the (a) no shading effect, (b) no groundwater and hyporheic exchange inflows, and (c) calculated boundary condition, Figure S3: Simulated time averaged river temperatures along the 1500 m Sawmill Creek reach for the original condition (Base) and for conditions with ±15% changes in (a) storm sewer temperature (T SS ), (b) sediment temperature, and (c) boundary conditions temperature, Figure S4: Simulated time averaged river temperature along the 1500 m Sawmill Creek reach for the original condition (Base) and for conditions with ±20% changes in (a) substrate hydraulic conductivity (SHC), (b) cloudiness factor (Cl), and (c) groundwater discharge (GW), Figure S5: Fluctuations of shading factors and daily average shortwave radiation along the 1500 m Sawmill Creek reach. The shading factors denoted by a triangle are measured at each of the 12 monitoring stations, and the minimum and maximum shading factors were selected from the 5 m interval set of shading factors measured between each station, Figure S6: Temperature differences between the observed and simulated river temperature when using Mohsni et al. [30], ∆T calc versus recorded, ∆T recorded boundary conditions, for (a) nighttime and (b) daytime.