Modeling Landscape Change Effects on Stream Temperature Using the Soil and Water Assessment Tool

Stream temperature is one of the most important factors for regulating fish behavior and habitat. Therefore, models that seek to characterize stream temperatures, and predict their changes due to landscape and climatic changes, are extremely important. In this study, we extend a mechanistic stream temperature model within the Soil and Water Assessment Tool (SWAT) by explicitly incorporating radiative flux components to more realistically account for radiative heat exchange. The extended stream temperature model is particularly useful for simulating the impacts of landscape and land use change on stream temperatures using SWAT. The extended model is tested for the Marys River, a western tributary of the Willamette River in Oregon. The results are compared with observed stream temperatures, as well as previous model estimates (without radiative components), for different spatial locations within the Marys River watershed. The results show that the radiative stream temperature model is able to simulate increased stream temperatures in agricultural sub-basins compared with forested sub-basins, reflecting observed data. However, the effect is overestimated, and more noise is generated in the radiative model due to the inclusion of highly variable radiative forcing components. The model works at a daily time step, and further research should investigate modeling at hourly timesteps to further improve the temporal resolution of the model. In addition, other watersheds should be tested to improve and validate the model in different climates, landscapes, and land use regimes.


Introduction
Stream temperature is an important water quality parameter that affects physical and chemical processes in streams [1]. Higher stream temperatures in river systems represent a growing concern worldwide and can affect the habitat and life spans of fish [2,3]. According to Eaton and Scheller [4], some fish species will disappear from the water body, if stream temperature transcends an upper limit. Using historical data ranging from 30 to 100 years, Kaushal et al. [5] reported that stream temperatures have been increasing throughout the United States at a rate of 0.009-0.077 _C/year, with a significant increase in the western United States. In particular, stream temperatures in the Pacific Northwest have reached historical records-at times, they have exceeded the lethal limit of 21.1 _C for some aquatic species such as salmon. For example, in the summer of 2015, the river temperature in the Columbia River in the State of Washington reached the level of 24.5 °C and led to the death of 235,000 sockeye salmon out of the total 507,000 that passed through the Bonneville Dam [6].
Similarly, Marys River (Hydrologic Unit Code 17090003) is a tributary of the Willamette River in Northwest Oregon and has experienced increasing temperatures over the last few years. The OregonDepartment of Fish and Wildlife (ODFW) conducted a stream survey on flow conditions between the years of 1991 and 1993 and observed that the maximum stream temperature reached the maximum limit of 17.8 °C (ODFW). However, recently, increasing levels of human activities have resulted in even higher water temperatures. For example, high water temperatures between 21.1 °C and 26.7 °C were observed from June to August for some tributaries of Marys River (Marys River Watershed Council), which has resulted in several United States Environmental Protection Agency (USEPA) 303(d) listings for temperature exceedances. These trends may be related to changing climatic drivers as well as land use practices (e.g., harvest of timber, increasing barren lands and clear-cutting areas throughout the watershed) and landscape changes (e.g., urban and agricultural development).
Amidst large-scale landscape and land use changes, preservation of riparian buffers can increase stream shading, thereby helping regulate water temperature along stream reaches [7]. Stream shading intercepts and absorbs a large portion of solar radiation before it reaches the water surface, resulting in less thermal energy that reaches and is stored in streams, which indirectly helps to cool stream temperatures. Brown et al. [8] conducted a study in the Alsea watershed in Oregon along the coast range to study the impact of shading on stream temperatures before and after clear cuts in the watershed. They found that clear-cutting resulted in stream temperature increases of 7.8 • C one year after the cuts. Bond et al. [9] investigated the impact of riparian reforestation on summer stream temperatures in the Salmon River in northern California, and they found that partial reforestation lowered stream temperatures by 0.11-0.12 °C/km and by 0.26-0.27 °C/km for full reforestation.
Since increasing water temperatures have remained a major concern in many watersheds, many models have been proposed to simulate stream temperatures at time scales varying from minutes to months (see Ficklin et al. [10] for a brief review). For this paper, we specifically focus on a semi-distributed mechanistic watershed model called the Soil and Water Assessment Tool (SWAT) [11], which has been extensively used to evaluate the effects of landscape and land use changes on different hydrologic components. Ficklin et al. [10] improved the original stream temperature model within SWAT, which was a linear regression model by Stefan and Preud'homme [12] that correlates thermal energy exchange of air temperature to water temperature. Ficklin et al. [10] developed a daily-scale model for stream temperature prediction by integrating multiple climate and hydrological components, including snowmelt, surface runoff, lateral flow, groundwater flow, and finally air temperature. Several studies have found the Ficklin et al. model [10] produces more realistic simulation results compared to the linear regression proposed by Stefan and Preud'homme [12] (Barnhart et al., 2014;Ficklin et al., 2012;Ficklin et al., 2014 [10,13,14]). However, the model developed by Ficklin et al. [10] does not explicitly account for the different types of radiation that affect thermal energy of water systems, and only recent work has attempted to improve the model by incorporating select radiative components [15]. In general, thermal energy added and removed from any water system consists of incoming radiation that adds thermal energy to the water and results in increasing stream temperatures. This incoming radiation mainly consists of solar radiation coming from the sun, atmospheric longwave radiation, landcover longwave radiation, convection, and evaporation. In contrast, backscattered radiation removes thermal energy and helps to cool temperatures. This radiation consists of emitted longwave radiation from the water surface as well as convection and evaporation. Convection and evaporation radiative components can either add or remove energy from any water body depending on stream temperatures and climate conditions, specifically air temperature, humidity, and wind speed. This paper highlights model development to explicitly incorporate multiple thermal radiation components into the Ficklin et al. [10] stream temperature model within SWAT. These thermal radiation equations are used within the widely used HEATSOURCE model [16], but until now, these equations have yet to be incorporated into SWAT. HEATSOURCE differs from SWAT because it is a reach-based stream temperature model, whereas SWAT is a watershed model. This means that SWAT simulates upland processes in addition to in-stream processes. HEATSOURCE can model stream temperatures at hourly timesteps and requires site-specific data (e.g., shading, canopy structure, stream morphology) that is oftentimes not available over the entire spatial extent of watersheds. Conversely, SWAT simulates hydrologic components and stream temperatures throughout a watershed using a daily time step and utilizes generally obtainable input data, such as spatially distributed precipitation, temperature, elevation, land use, and soil type. Users may prefer to use SWAT instead of HEATSOURCE when site-specific data is unavailable or when the study goal is to determine the effect of alternative land management scenarios on stream temperatures throughout large, heterogeneous watersheds.
The paper is organized as follows. First, the study area and the SWAT model setup are described. Then, three different SWAT stream temperature models are analyzed, including Stefan and Preud'homm [12] air temperature regression, a mechanistic model by Ficklin et al. [10], and an extension of the Ficklin et al. [10] model in which we specifically incorporate radiative components. We calibrate SWAT for hydrologic discharge in the Marys River watershed, and we compare the three stream temperature models to examine their relative performance for multiple sub-basins with different land use/land cover. We demonstrate the utility of our results by comparing simulations for sub-basins within primarily forested and agricultural landscapes.

Study Area
The Marys River watershed, shown in Figure 1, is located in the Pacific Northwest of the United States (Hydrologic Unit Code (HUC) 17090003) and is part of the Willamette River basin (HUC 170900) in Oregon. It is one of five major river systems located on the western side of the Willamette River. The area of the watershed is of 782 square kilometers. The highest point of the watershed is at Marys Peak at an elevation of 1280 meters above sea level, and the lowest point is in Corvallis, Oregon, where Marys River drains into the Willamette River at an elevation of 76 meters above sea level. The climate of the watershed in the winter season is mild and wet, with an average winter temperature of 5 °C and rainfall during the winter ranges from 1000 mm downstream of the watershed to more than 2500 mm at the highest elevation upstream of the watershed. In general, the watershed tends to be dry, sunny, and warm throughout the summer (Marys River Watershed Council). It has an average summer temperature of 17.5 °C. High rainfall intensity results in high stream discharge during winter and spring and mean annual flows of 12-13 m 3 /s. However, during the summer, flows are generally very low, and discharge sometimes drops below one cubic meter per second. Base flow is a major contributor to the flow of the river, where 61-70% of the total stream flow comes from groundwater contributions [17].
The watershed is divided into three different land use categories: Forest, agricultural, and urban. Most of the watershed (65% of the total area) consists of forest, which is largely located along the western portion of the watershed. In these mixtures of deciduous and evergreen forests, small streams flow over beds of gravel and cobbles with high velocities due to steep slope gradients. Flow leaving the forested region then enters agricultural land in the Willamette Valley that consists mainly of cultivated crops, hay, pasture, wheat, and grass seed production. The streams in this region flow on sand and silt with mild slope gradients, resulting in decreased flow velocities. Furthermore, urban areas are situated further downstream within the Willamette Valley (e.g., the cities of Philomath and Corvallis), and the stream flows over mostly flat to nearly flat gradients. Stream velocities decrease significantly as Marys River enters Philomath, Oregon, and then continues eastwards into Corvallis, Oregon, until it meets the Willamette River at the lowest point located in the watershed.

Soil and Water Assessment Tool (SWAT)
The Soil and Water Assessment Tool (SWAT) is a semi-distributed watershed model that is designed to predict the impact of management on water, sediment, and agricultural chemical yields in gauged and ungauged watersheds [11]. In this study, SWAT was used to simulate the hydrologic and stream temperature dynamics within the Marys River watershed and to evaluate a model extension to the stream temperature model developed by Ficklin et al. [10]. The Marys River watershed was divided into smaller sub-watersheds using pre-defined drainage boundaries and a 10-meter digital elevation model, and then the sub-watersheds in the SWAT model were further divided into smaller units called hydrologic response units (HRUs) using ArcSWAT, a toolbox within ArcGIS for SWAT, with a HRU percentage threshold of 5%. Each HRU is a unique combination of land-use, soil type, and topographic slope and represents the basic unit for conducting mass balances and hydrologic flow in SWAT. The area of the 46 predefined sub-basins varies from 35 square kilometers for the largest sub-basin to 3.0 square kilometers for the smallest. The average sub-basin area is 17 square kilometers. The watershed slope was divided into two categories: (1) A steep gradient area located mostly within the forested regions in the western portion of the watershed, and (2) the nearly flat region located within the Willamette valley, east of the watershed where the cities Philomath and Corvallis are located.
SWAT's input data types include spatial GIS input files such as a Digital Elevation Model (DEM), a land use land cover layer, and a soil layer [18]. Input data needed to delineate the watershed including the DEM, sub-basins, and stream layers in addition to necessary land use and soil SSURGO (Soil Survey Geographic Database) layers to build the HRUs were acquired from United States Department of Agriculture [19]. Three weather stations were used as climate forcings: Corvallis Water Bureau (CWB) COOP ID of (351877), Hyslop weather station, which is also known as Oregon State University weather station (OSU) COOP ID of (351862), and finally Corvallis municipal airport (KCVO) weather station. Weather data of the Corvallis Water Bureau (CWB) and Hyslop weather station were obtained from National Oceanic and Atmospheric Administration (NOAA) for 2005 to 2014. Weather data for Corvallis Water Bureau station included only precipitation and minimum and maximum air temperature. The Hyslop weather station and Corvallis municipal airport stations included data for precipitation, minimum and maximum temperature, wind speed, and humidity for the period of 2005 to 2014. SWAT was used to simulate flow and stream temperature throughout the Marys River watershed for the period 2005-2014. This includes the period 2010-2014 when observations for stream temperature were available.

Model 1: Linear Regression-
The default SWAT stream temperature model uses a linear relation between air temperature and stream temperature developed by Stefan and Preud'homme [12] to calculate stream temperature in the Mississippi River basin, as shown in Equation (1): T water is the average daily water temperature (°C), and T water is the average daily water temperature (°C). Stream temperatures predicted from the above equation will always be higher than air temperature, which is generally a fair assumption for small streams with shallow water depths where stream temperature is primarily controlled by air temperature. However, this may not be necessarily true for streams influenced by snowmelt, surface runoff, and groundwater contributions [10].

Model 2: A Mechanistic Approach Involving Air Temperature and
Hydrological Flows-Ficklin et al. [10] developed a mechanistic stream temperature model within SWAT by combining air temperature (heat exchange) and hydrological inputs (flow mixing) including different hydrological parameters, surface runoff, lateral flow, snowmelt, and groundwater contributions. The Ficklin et al. [10] stream temperature model discretizes stream temperature determination into three components: (1) Within the subbasin, (2) contribution of upstream sub-basins to the targeted sub-basin, and (3) finally heat exchange between air temperature and the stream.
The first part of the stream temperature calculation within the sub-basin (Equation (2)) calculates the local temperature based on a mixing of surface runoff, lateral flow, groundwater, and snowmelt temperatures within the sub-basin flowing to the main stream: T w,local = α 0.1 Sub snow + β T gw Sub gw + λ T air, lag Sub surq + Sub latq Sub wyld (2) T air,lag is average daily air temperature with a lag (°C), and α, β, and λ are calibration coefficients that relate the relative contribution of the hydrologic components to local water temperature (dimensionless). Sub snow is the snowmelt contribution in sub-basin (m 3 /day), Sub gw is the groundwater contribution in sub-basin (m 3 /day), Sub surq is the surface runoff in the sub-basin (m 3 /day), Sub latq is the lateral soil flow in sub-basin (m 3 /day), and Sub wyld is the water yield in the sub-basin combining all of the above hydrological inputs (m 3 /day).
The second part of the Ficklin et al. [10] calculates the effect of upstream sub-basin flow on stream temperature, as shown in Equation (3): T waterinitial is the stream temperature adding the effects of flow within the sub-basin, T w,local was calculated previously, T w,upstream is the water temperature of streams entering the subbasin (°C), and Q outlet is the stream flow discharge at the outlet of sub-basin (m 3 /day).
The final step is to calculate the stream temperature by including the effect of air temperature: T water = T waterinitial + K T air − T waterinitial (TT) if T air > 0 T water = T waterinitial + K T air + ε − T waterinitial (TT) if T air < 0 (4) T water is the final stream temperature of water (°C) for a given sub-basin, T air is the average daily air temperature (°C), K is the bulk coefficient of heat transfer (1/h), TT is travel time of water through the sub-basin (hour), and finally ε is air temperature addition coefficient (for when air temperature drops below zero).
This mechanistic stream temperature model requires calibration coefficients α, β, γ, k, lag as well as annual groundwater temperatures as inputs. All of the other inputs needed to run the model are provided by SWAT. Groundwater temperature can be estimated from weather data provided as the annual average air temperature, and it is often taken 1-2 °C higher than the average annual air temperature [20].

Model 3: A Mechanistic
Approach Involving Air Temperature, Hydrological Flows, and Radiative Components-Changes in stream temperature are affected by heat and mass transfers [16] that are dependent on channel morphology, hydrology, and stream vegetation, which provides shading near streams. Vegetation especially helps in cooling temperatures by intercepting and absorbing incoming solar radiation. The mechanistic model introduced in the last section accounts for flow transfer and mixing of various hydrologic components, including sub-basin surface runoff, lateral flow, snowmelt, and ground water, which oftentimes help to reduce temperatures, depending on the season. However, the model utilizes a bulk heat coefficient to account for radiative heat exchange between the air-water interface and does not account for land cover or vegetation near streams. The dependency of the model on the air-water correlation of heat exchange can thereby lead to over-prediction of stream temperatures. Water temperature change related to heat transfer is a function of several sources of radiative heat exchange, as shown in Equation (5): where, Φ total is the net radiation exchange and is equal to the direct and diffuse solar radiation Φ SR as well as the longwave-atmosphere, longwave-landcover, convection, and evaporation components.
Direct and diffusive solar radiation represent the largest sources of incoming thermal energy into streams. Longwave radiation from different sources also plays a role in increasing and decreasing temperatures: Atmospheric and land-cover longwave radiation add energy to water volumes and increase water temperature, while longwave radiation from the water surface emits radiation from the water surface to the atmosphere to cool streams. Energy loss due to evaporation is considered a larger contributor of decreasing stream temperatures when the required energy is met to change the water phase from liquid to gas. Overall, convection or the air-water interface is considered a very small portion of the total energy budget. Groundwater flux also helps to decrease stream temperatures when added as a thermal input.
These radiative components have been included in the HEATSOURCE model but have not been incorporated explicitly within spatially distributed watershed models such as SWAT.
Each of the components are defined as follows.
Φ SR = H day − 0.5 × H day 1 − e −k × LAI Φ SR is the amount of solar radiation reaching the water surface, H day is the incident total solar radiation per day (MJ/m 2 ·day), k is the light extinction coefficient, and LAI is the leaf area index. Atmospheric solar radiation is calculated as follows: Φ longwave-atmosphere = 0.96 × ε atm × σ × (T air + 273.15) 4 Φ longwave-atmosphere is the longwave radiation emitted from the atmosphere, ε atm is the emissivity of the atmosphere (unitless), σ is the Stefan-Boltzmann constant (MJ·m −2 ·day −1 ·K −4 ), and T air is the average air temperature per day ( • C). The atmospheric solar radiation depends on the emissivity of the atmosphere (ε atm ), in which the percent of cloudiness and type of landuse can increase or decrease atmospheric emissivity. The longwave radiation emitted from landcover Φ longwave-landcover is dependent on the view to the sky θ VTS (unitless) and can be calculated as follows: Φ longwave-landcover = 0.96 × 1 − θ VTS × 0.96 × σ × T air + 273.15 4 (8) The last component of longwave solar radiation is the water surface longwave: Φ longwave-water surface = ε W × σ × T s + 273.15 4 (9) Φ longwave-water surface surface is the longwave radiation emitted from the water surface, ε w is the water emissivity taken as 0.97, and T s is the stream temperature (°C). The evaporation from the water surface is the most effective component in decreasing the thermal energy stored in water: Φ evaporation = ρ × L e × E (10) ρ is the density of water (kg/m 3 ), L e is the latent heat of vaporization (MJ/kg), which is calculated as L e = 2.501 − 2.361 × 10 −3 × T air . E is the evaporation rate of the water surface (m/day) and is calculated using a mass transfer method: f(w) × (e sw − e aw ), where f(w) is the wind function a + b × w that depends on coefficients a and b (mbar −1 ) and the wind speed w (m/s) measured 2 m above the water surface [16]. Finally, e sw is the saturation vapor pressure of water (mbar), and e aw is the vapor pressure of water (mbar).
The convection radiation component is calculated using the previously calculated evaporative flux Φ evaporation and Bowen's ratio BR: The change of stream temperature due to thermal energy flux is calculated as follows: Here T w−TD°C day is the temperature change generated from the thermal components, Φ total is Here, T w−TD the net driver (MJ/m 2 ·day), C w (MJ/kg·C) is the specific heat capacity of water, and d w (m) is the depth of water in the channel, which is estimated by the SWAT model.
The final stream temperature is calculated by replacing the second term in Equation (4) by the new generated stream temperature, as follows: T water =T waterinitial + T w−TD (13) Note that, as mentioned above, the majority of the components are calculated within the SWAT model automatically.
Overall, this model is useful because it utilizes all of the distributed information (mechanistically simulated via SWAT) regarding stream temperature and its relation to hydrologic components within a networked watershed, yet it also explicitly incorporates radiative energy exchange at the surface of the stream.

Model Calibration/Validation Methodology
SWAT was used to simulate daily hydrologic discharge at each of the sub-basins within the Marys River watershed from 2010-2014 in order to match observed discharge and stream temperature data. The SWAT model was manually calibrated for stream flow between 2010-2014 using the United StatesGeological Survey (USGS) (14171000) Philomath flow gauge, which is located 6.7 km southwest of where the Marys River meets the Willamette River, covering a 394 km 2 .
The model was only calibrated without validation due to the availability of observations for flow since the available observations only included a period of less than 10 years. Based on the data availability, the model was calibrated for the period of January 2010 to December 2014.
The Nash Sutcliffe efficiency (NSE; Nash and Sutcliffe (1970)) criterion Equation (14) and Pearson's product moment correlation coefficient (1999; Equation (15)) were used to evaluate hydrologic model efficiency: O is the observed value, S is the model prediction, O avg is the overall observed mean, and S avg is the overall simulated mean. The NSE values range from −∞ to one; a NSE value of less than 0.5 designates an 'unsatisfactory' model, while a NSE value above 0.75 is considered a 'very good' model [21]. R 2 values range from zero to one, with zero indicating a nonlinear relationship between the observed and predicted value and one indicating a perfect fit and a linear relationship between the observed and simulated variables.
The stream temperature models were manually calibrated using root mean square error (RMSE) values as well as percent bias (PBIAS). According to Chai et al. [22], RMSE is widely used as a statistical metric tool to assess performance of models. RMSE can be calculated as follows: Where O is observed value, S is model predicted value, and n is the total number of the points. A RMSE value of zero indicates a perfect fit. According to Singh et al. [23], values of RMSE less than half of the standard deviation of the measured data can be taken as acceptable for the model evaluation.
In addition to RMSE, PBIAS was also used to evaluate the models. PBIAS measures the average tendency of the simulated data to be larger or smaller than their observed counterparts [24]. Zero is the optimal value of PBIAS, and a low absolute value implies an accurate model. According to Gupta et al. [24], positive PBIAS values indicate a model underestimation bias, and negative PBIAS values indicate a model overestimation bias. PBIAS can be calculated as follows: A set of seven calibration parameters were selected and manually modified to calibrate SWAT for hydrology (Table 1), and five parameters were manually modified to calibrate SWAT for stream temperature ( To compare the differences in model performance, kernel density estimates were calculated using R software. This nonparametric technique is similar to using histograms to highlight the differences between model simulations and observed data for the three tested models.

Hydrology Calibration in SWAT
As mentioned previously, SWAT was used to simulate daily hydrologic discharge at each of the sub-basins within the Marys River watershed for 2005-2014, which included 2010-2014, when observations for stream temperature were also available. The NSE and R 2 values of the default model's simulated flow (no calibration) were −0.37 and 0.50, respectively, indicating an unsatisfactory model. From Figure 2, it is clear that the uncalibrated model overpredicts the peaks during storm events. Also, the model's responses to each rain event are very rapid, and the water loss rates are excessive, which results in zero flow for late summer periods.
Manual calibration resulted in model improvement of NSE values from −0.37 to 0.72. This designates a 'good' model according to Moriasi et al. [21], since it is >0.65. The R 2 value of the model increased to 0.80. Figure 3 shows that the calibrated model matches both the base flow and peaks well, whereas all of the peaks were mainly over predicted by the default model. The default model could not capture the very low stream flows for some summer days and resulted in zero flow (Figure 2), but the calibrated model ( Figure 3) fixed the lowflow problem and improved the results. Figure 2 shows that the maximum simulated peak for the uncalibrated model was around 500 m 3 /s, whereas the calibrated model ( Figure 3) reduced this value to match the peak observed around 300 m 3 /s. The rapid response of the main channel to any storm event led to the over-prediction of peaks even for small rain storms in the uncalibrated model. Also, the lateral flow travel time parameter helped to slow the response to storm events and smoothed the hydrograph. The other major problem was associated with the excessive loss of water in a short period of time after a rapid response to any storm; water in the uncalibrated model was lost instantaneously and therefore resulted in zero flow for late summer days. The effective hydraulic conductivity in the main channel (CH_K2) was used to prevent the excessive loss from the main channel and managed to eliminate the zero flow days. Also, this parameter helped to smooth and eliminate the transient fluctuations in the hydrograph. The SCS curve number for moisture conditions II (CN2) was used to reduce or increase the simulated peaks to match the observed hydrograph. The canopy interception parameter for specified land cover types (CANMX) as well as the minimum and maximum snowmelt factors (SMFMN, SMFMX) all helped to increase the evaporation rate. The high surface runoff surge was decreased using the soil evaporation compensation factor (ESCO), the SCS curve number for moisture conditions II (CN2), and the threshold depth of water in the shallow aquifer (GWQMN).

Stream Temperature Calibration in SWAT
After a satisfactory hydrologic calibration was performed, manual calibration was performed for two of the three stream temperature models that will be tested in the Marys River using SWAT. The first model is the linear regression model from Stefan and Preud'homme [12] and was not calibrated. The second model is the Ficklin et al. [10] model, and the third model is our extension to the Ficklin et al. [10] model. As shown in Table 2, we manually calibrated five parameters to best match the second and third models to the stream temperature observations between 2010 and 2014. Note that only summer stream temperature observations were available and that only a single set of calibration parameters were used for the entire watershed.

Stream Temperature Model Comparison
After satisfactory hydrologic and stream temperature calibrations were performed, three models were tested to simulate stream temperature in the Marys River using SWAT. The first was the linear regression model from Stefan and Preud'homme [12], the second the Ficklin et al. [10] model, and the third is our extension to the Ficklin et al. [10] model that incorporates HEATSOURCE radiative forcing components to the Ficklin et al. [10] model. For the remainder of this paper, these models will be referred to as Model 1, Model 2, and Model 3, respectively. Table 3 shows RMSE and PBIAS results for these three models using daily data. Kernel density estimates for the differences between simulated and observed stream temperatures for the four sub-basins are plotted in Figure 4.
In general, Model 1 performed the worst among all three models (see Table 3 and Figure 4), and Model 2 outperformed Models 1 and 3 for all sub-basins. Model 1 consistently overestimated stream temperatures for all of the sub-basins. This is apparent in Figure 4, where the distributions of the differences between the simulated and observed temperatures are shifted from zero for Model 1. The coefficients of Model 1 guarantee that the simulated stream temperature will be above average daily air temperature values when the average air temperature is less than 20 • C. Therefore, inclusion of cold water from groundwater or upstream sources is not captured in this model, thus resulting in overestimations. Model 2, which is the calibrated Ficklin et al. [10] model, showed improvements compared to the linear model (Model 1) in both Table 3 and Figure 4, which agrees with previous studies (Barnhart et al., 2014;Ficklin et al., 2012;Ficklin et al., 2014) [10,13,14]. This is presumably because Model 2 incorporates hydrologic components, including groundwater upstream temperatures, in addition to an air-heat exchange transfer coefficient. Model 3, which replaced the simple air-heat exchange transfer coefficient from Model 2 with explicitly calculated radiative components, shows similar distributions with Model 2 ( Figure  4), yet the performance values of Model 2 are generally better. This is likely due to the high variability of the radiative components included in Model 3, which will be discussed further in the next section.

Land Cover Effects on Stream Temperature
We now compare Models 2 and 3 to demonstrate that the incorporation of radiative components is able to simulate the influence of land use and land cover on stream temperatures (e.g., forested vs. agricultural regions). Stream temperature simulations for two sub-basins-a forested area (sub-basin 8) and an area dominated by agriculture with low vegetation (sub-basin 36)-are shown in Figure 5.
Model 2 simulates nearly identical stream temperatures for both forested and agricultural sub-basins. Conversely, by including the various radiative components into the model, Model 3 simulates consistently increased stream temperatures associated with the agricultural sub-basin. This is especially apparent during the early summer months. To examine differences in the models further, Figure 6 compares the net radiative components of Model 3 Equation (5) with the K-component second term in Equation (4) of Model 2. Figure 6 shows that the net radiation as calculated using Model 3 (orange lines) changes according to the primary land use cover for a given sub-basin in SWAT. For example, agricultural sub-basins have larger incoming (positive) radiation contributions that help to increase stream temperatures, whereas forested areas have small or negative radiative effects, depending on the season, due to increased LAI, reduced solar radiation reaching the water surface, and evaporative fluxes. Model 2 uses a bulk coefficient of heat transfer in the second term of Equation (4). This reflects a convection component of the net radiative balance, but it does account for the other radiative energy terms, including solar radiation, atmospheric longwave, land surface longwave, water surface longwave, and evaporation. Therefore, it is not able to capture cover-related differences in net radiation and therefore changes in stream temperature due to landscape and land use changes. Figure 6 also shows that the net radiative driver as calculated in Model 3 has much higher variability than the K-component used in Model 2. We found that the high variability (i.e., noise) is mainly due to SWAT's estimation of solar radiation as well as the evaporation calculations, which are not explicitly accounted for in Model 2. This likely led to the reduced model performance exhibitedwhen comparing model simulations to observed data. Figure 7 compares these data further by plotting the observed stream temperature as well as simulations using Models 2 and 3 for both the forested (sub-basin 8) and agricultural (subbasin 36) sub-basins.
The mean (standard deviation) differences between the stream temperature simulations for the forested sub-basin (8) (Table 3). However, Model 2 is symmetric about the 1:1 line shown in Figure 7 and does not capture the positive bias of agricultural stream temperatures that is shown by the observed values as well as in the simulated values of Model 3. Model 3 captures the increases in stream temperature associated with the lack of forest cover in the agricultural sub-basin, yet Model 3 overpredicts this effect, and the variation in the simulated values are much greater than Model 2 or the observed data.
Overall, Model 3 may be useful for simulating the watershed-scale impacts of land use conversion from forest to agriculture on stream temperatures; however, model improvement -potentially through improved calibration-is needed to better match observed data. In addition, while Model 3 is able to simulate impacts of land use on stream temperature, this advantage is produced with the trade-off that radiative components exhibit much higher variability, and this noisy fluctuation (which can be seen directly from Figures 5-7) further decreases RMSE values (Table 3).

Conclusions
This study sought to explicitly incorporate radiative forcing components into an existing mechanistic, semi-distributed stream temperature model Ficklin et al. [10] using SWAT. Ultimately, stream temperature is controlled by different climate components including humidity, wind speed, evaporation, and solar radiation besides air temperature and is also heavily dependent on hydrologic components including surface flow, groundwater, and snow melt processes. Our new model leverages the Ficklin et al. [10] stream temperature model, which accounts for hydrological and climatological components and discretizes the prediction of stream temperature into three parts: (1) Streamflow within a sub-basin, (2) contributing upstream sub-basin hydrologic components, and (3) accounting for the heat exchange between air and water surface, which can be described as a convection term. The extended model replaces the convection K-component term within the Ficklin et al. [10] model with a more comprehensive characterization of radiative energy terms, including solar radiation, atmospheric longwave, land surface longwave, water surface longwave, evaporation, and finally convection drivers. The extended model was used along with the Ficklin et al. [10] model and the linear regression model to simulate stream temperatures within agricultural, forested, and mixed sub-basins within the Marys River watershed. Results showed that all models performed reasonably well, and the Ficklin et al. [10] model outperformed the others. However, the extended model was capable of simulating differences between stream temperatures associated with agricultural and forested watersheds that reflected observed data, although the differences were overestimated. The reduced performance of the extended model that included radiative components might be able to be improved by further calibration; yet, the high variability of the radiative terms is also limiting. For example, the model relies on incoming solar radiation as well as wind velocity, which are difficult to represent over large spatial scales and feature high variability. Alternative formulations of the radiative components should be considered in future work. Bogan et al. [1] suggests using a shading factor instead of a leaf area index, which may lead to more accurate results, assuming shading information is available for the watershed of interest. In addition, alternative estimations for solar radiation, or perhaps any of the radiative energy terms, could improve the model and should be pursued. Overall, incorporating radiative components into the Ficklin et al. [10] stream temperature provides a new mechanism for simulating the effects of alternative land uses on stream temperature within SWAT. This will be especially useful for land managers and decision makers when considering alternative land management scenarios and conservation strategies using SWAT.   Kernel density estimates for the differences between observed and simulated stream temperatures as calculated using the three models for four sub-basins (a-d). The vertical line at zero indicates simulations that exactly match the observed data. Comparison of the thermal energy impacts on stream temperatures for Models 2 and 3 in both primarily forested and agricultural sub-basins. Comparison of observed and simulated stream temperatures for agricultural and forested sub-basins. The observed data (black) show a consistent increase in stream temperatures for the agricultural sub-basin (sub-basin 36). Model 2 does not capture consistent changes in stream temperature associated with land cover type. Model 3 simulates increased temperatures for the agricultural sub-basin, but it overestimates the effect compared to the observed data.