MODELING TEMPERATURE AND NITROGEN DYNAMICS IN MIXED LANDUSE WATERSHEDS USING A PROCESS - BASED HYDROLOGIC MODEL By Han Qiu A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of Civil Engin eering - Doctor of Philosophy 201 9 ABSTRACT MODELING TEMPERATURE AND NITROGEN DYNAMICS IN MIXED LANDUSE WATERSHEDS USING A PROCESS - BASED HYDROLOGIC MODEL By Han Qiu Hypoxia and eutrophication resulting from excessive nutrient loading are one of the most s ignificant environmental issues around the world. Although the 1972 Clean Water Act has effectively reduced point source loadings of nutrients to surface waters, controlling diffuse, nonpoint source pollution continues to be a challenge. Anthropogenic acti vities, including land use change, are considered some of the main reasons for the excessive riverine nitrogen (N) loading. Temperature, stream d ischarge, the structure of the drainage network as well as soil moisture are among the important factors influe ncing nitrogen transport and transformation in watersheds. Of particular interest is temperature, which was found to be a key factor, influencing nitrogen transformation processes; however, modeling temperature in watersheds is challenging due a large numb er of coupled processes involved. Stream thermal regimes are primarily driven by climatic conditions and influenced by a host of other factors, i ncluding topographic conditions, stream discharge, land cover near the stream and interactions with the subsurf ace. Riparian vegetation processes close to the stream banks control canopy shading, as do factors such as the spatial heterogeneity of vegetatio n density and temporal aspects of vegetation growth. Vegetation type affects stream temperature while also infl uencing the riparian microclimate including air temperature, wind speed and relative humidity. These complexities call for an integrated model th at can describe coupled hydrologic - vegetation processes. This dissertation research involves the development an d application of an integrated and fully process - oriented water - temperature - nitrogen model based on the modeling framework of PAWS (Process - based Adaptive Watershed Simulator). The integrated model was tested using data from two watersheds of different siz es and climatic conditions - the Wood Brook watershed in central England located at the Birmingham Institute of Forest Research (BIFOR) and the K alamazoo River watershed in Michigan. The phenology and surface energy modules in the coupled model were used t o quantify the impacts of vegetation processes on radiation fluxes (e.g., canopy shading and the effect of vegetation growth on optical parameter s). The integrated temperature model enabled accurate simulations of the movement and partitioning of water and thermal fluxes in stream, soil, streambed, and groundwater domains and allowed the identification of gaining and losing portions of stream reach es. Nitrogen transport and transformations on the landscape were modeled by representing multiple sources and p rocesses (fertilizer / manure application, WWTPs, atmospheric deposition, N itrogen retention and removal in wetlands and other lowland storage, t emperature - dependent transformation rates etc.) across multiple hydrologic domains (streams, groundwater, soil water). The coupled model provides a tool to examine N itrogen budgets and to quantify the impacts of human activities and agricultural practices on the riverine export of nitrogen species. iv To my beloved v ACKNOWLEDGEMENT S I would like to express my sincere gratitude to my advisor, Dr. Mantha S. Phanikumar , the completion of this dissertation would be impossible without his help. Many years later, I will still remember his valuable, insightful, and thoughtful guidance and priceless instructions throughout my Ph.D. period. His consistent enthusiasm and uncom promising attitude toward academic research are lightening my way. I appreciate that I could ha ve the opportunity to work with Dr. Mantha. I would like to thank my committee members: Dr. Kalyanmoy Deb, Dr. Shu - Guang Li , Dr. Yadu Pokhrel and Dr. Amor V. M. Ines for their valuable comments, suggestions and advice to improve this dissertation. I would also like to thank my senior lab mate, Jie Niu, for his precious advice and valuable discussions. I also thank Dr. Huasheng Liao for his sugg estions and technical instructions in simulating groundwater dynamics. I thank Drs. Stefan Krause and Philip Blaen from the University of Birmingham for sharing the data for the Millbrook watershed. Last but not the least, I would like to express my deepe st gratitude to my parents and my girlfriend, your encouragement and unconditional support endowed me with the power to persevere in pursuing unknowns and to insist on my academic career. vi TABLE OF CONTENTS LIST OF TABLES ································ ································ ················· viii LIST OF FIGURES ································ ································ ·················· x Chapter 1. Background and Motivation ································ ··························· 1 Motivation for developing a new model of nitrogen transport and transformation · 1 Motivation for developing a new watershed scale temper ature model ··············· 3 Chapter 2. Evaluating a coupled phenology surface energy bala nce model to understand stream subsurface temperature dynamics ································ ······················· 7 Introduction ································ ································ ···················· 7 Materials and Methods ································ ································ ····· 11 2.2.1. Site Description and Observational data ································ ··········· 11 2.2.2. The hydrologic model ································ ································ · 14 2.2.3. The stream temperature module ································ ····················· 19 2.2.4. The streamb ed temperature module and hyporheic exchange ·················· 32 2.2.5. The soil temperature module ································ ························· 34 2.2.6. The groundwater temperature module ································ ·············· 34 Results and Discussion ································ ································ ····· 36 2.3.1. Stream flow, groundwater head and soil moisture comparisons ················ 37 2.3.2. LAI and ET comparisons ································ ····························· 41 2.3.3. The temperature results in different hydrologic domains ························ 45 Conclusions ································ ································ ·················· 68 Chapter 3. Modeling the effects of vege tation on stream temperature dynamics in a large, mixed land cover watershed in the Great Lakes region ································ ······· 70 Introduction ································ ································ ·················· 70 Materials and Methods ································ ································ ····· 74 3.2.1. Model description ································ ································ ····· 74 3.2.2. Study sites and data sour ces ································ ·························· 78 3.2.3. Model setup ································ ································ ············· 81 Results and Discussion ································ ································ ····· 84 3.3.1. Hydrology ································ ································ ·············· 84 3.3.2. Soil moisture and soil temperature results ································ ········· 93 3.3.3. Stream temperature simulations ································ ····················· 95 3.3.3.1 Effects of resolving spatial heterogeneity in vegetation ····················· 95 3.3.3.2 Compared wit h observed stream temperature data ··························· 99 3.3.3.3 Seasonal heat budgets ································ ···························· 103 3.3.3.4 Stream temperature responses to potential deforesta tion ··················· 107 Conclusion ································ ································ ·················· 110 Chapter 4. An integrated, catchment - scale framework to model temperature - dependent nitrogen transport and transformations ································ ························· 111 vii Introduction ································ ································ ················· 111 Methods and materials ································ ································ ···· 115 4.2.1. The nitrogen observation data ································ ······················ 115 4.2.2. The model framework ································ ······························· 116 4.2.3. The Nitrogen sources ································ ································ 118 4.2. 3.1 The fertilizer/manure amount and timing: Kalamazoo River watershed · 118 4.2.3.2 The fertilizer/manure amount and timing: Wood Brook watershed ······ 119 4.2.4. Nitrogen reaction and transport in river ································ ··········· 121 4.2.5. Nitrogen Reaction and Transport in the overland flow ························· 1 23 4.2. 6. Nitrogen reaction and Transport in the vadose zone ···························· 126 4.2.6.1 The internal nitrogen cycle ································ ······················ 127 4.2.6.2 The external nitrogen cycle ································ ····················· 128 4.2.7. Nitrogen reaction and transport in the groundwater domain ··················· 131 4.2.8. Temperature dependent reaction rates ································ ············· 133 Results and discussions ································ ································ ··· 134 4.3.1. Comparisons of hydrology results between CLM4.5 and CLM4.0 ··········· 134 4.3.2. The Kalamazoo watershed nitrogen results. ································ ······ 137 4.3.2.1 Comparison of nitrate concentrations in rivers ······························· 137 4.3.2.2 Nitrate leaching ································ ································ ··· 141 4.3.2.3 Comparison of nitrate concentrations in groundwater ······················ 143 4.3.2.4 Effects of temperature on the riverine nitrogen export ····················· 145 4.3.2.5 Uncertainty analysis ································ ····························· 148 4.3. 3. The Wood B rook watershed nitrogen results ································ ····· 153 4.3.3.1 The nitrate concentration comparison in river ······························· 153 4.3.4. Nitrogen budget ································ ································ ······· 155 Conclusions ································ ································ ················· 156 Chapter 5. Conclusion ································ ································ ············· 158 APPENDIX ································ ································ ························ 160 REFERENCES ································ ································ ····················· 163 viii LIST OF TABLES Table 2.1 Model calibration parameters for hydrology ................................ .................... 18 Table 2.2 Model parameters for temperature module ................................ ....................... 18 Table 2.3 Seasonal heat flux budgets for year 2016 at the catchment outlet stream segment (segment CD in Figure 2 .1) ................................ ................................ ........ 52 Table 2.4 Performance of the streambed temperature sub - model at four sampling sites for different depths (5 cm, 10 cm, 20 cm, and 30 cm) ................................ ................... 56 Table 2.5 Thermal conductivities and leakance values parameterized for the four streambed temperature sampling sites. ................................ ................................ ..... 56 Table 2.6 Model parameters for temperature modules ................................ ..................... 62 Table 2.7 Average percentage changes ( o C x 100/ o C) of hourly stream/streambed temperatures in response to parameter perturbations ................................ ................ 63 Tab le 3.1 Performance of Stream discharge comparisons between simulated results and observations by USGS gauges ................................ ................................ .................. 86 Table 3.2 Pair - wise linear correlation coefficients for different Land use/land co ver types and soil types with GLB and MLT simulated ET outputs ................................ ........ 88 Table 3.3 Performance of the temperature comparisons between simulated results and observations using different scales of stream side land use information .................. 97 Table 3.4 Characteristics of the marked three stream segments ................................ ..... 108 Table 4.1 Fertilizer/manure nitrogen application for KRW (USDA, 2013) ................... 120 Table 4.2 Fertilizer/manure nitrogen application for Wood Brook watershed ............... 120 Table 4.3 Defaul t reaction rates used in the river t ransport module ............................... 122 Table 4.4 Monod kinetics reaction parameters ................................ ............................... 133 Table 4.5 Performance of St ream discharge comparisons between simulated results using CLM4.0, CLM 4.5 and observations by USGS gauges ................................ .......... 136 Table 4.6 Performance metrics of river nitrogen concentration time series compariso ns ................................ ................................ ................................ ................................ 140 Table 4.7 Pair - wise linear correlation coefficients for nitrate leaching values with different Land use/land cover types and soil texture ................................ .............. 140 ix Table 4.8 RMSE values of nitrate concentrations between baseline simulations and perturbated simulations (The RMSE value is calculated as the mean value of +10% and 10% for each case) ................................ ................................ .............. 149 Table 4.9 Nitrogen budget results ................................ ................................ ................... 156 x LIST OF FIGURES Figure 2.1 Map of the Wood Brook watershed with land us e and land cover. Locations of stream flow sensors, borehole sensors, streambed temperature sensors and soil moisture and temperature sensor are shown. Elevation is shown as the color gradient in the sub - map. ................................ ................................ ................................ .......... 12 Figure 2.2 Program flow diagram for the essential processes of PAWS+CLM. The modules in the blue boxes are the primary processes for the hydrologic components, the modules in the red boxes are the primary processe s of temperature simulat ions, and the modules in the black boxes are primary sub - processes. V/M/H Transfer denotes Vapor/Momentum/Heat Transfer. ................................ ............................... 19 Figure 2.3 Sketches of the radiation fluxes partitions for (a) short - w ave radiation fluxes (b) long - wave radiation fluxes (c) weighted - average heat flux computation for grid cells with mixed land use ................................ ................................ .......................... 24 Figure 2.4 Simulated and observed streamflow comparison at ba sin outlet ..................... 39 Figure 2.5 Simulated and observed transient groundwater head comparison at four different boreholes (BH), asl = above sea level ................................ ........................ 40 Figure 2.6 Comparison of simulated and observed soil moisture ................................ ..... 40 Figure 2.7 Spatial maps of 8 - day LAI: (a ) MODIS data for Sep.07 - Sep.14, 2015 (b) Simulated results for September 07 - 14, 2015 (c) MODIS data for July 05 - 12, 2016 (d) Simulated results for July 05 - 12, 2016 ................................ .................... 43 Figure 2.8 Time series comparisons between MODIS data (MOD15 product) versus simulation resul ts for 8 - day catchment averaged leaf area index ............................. 43 Figure 2.9 Spatial maps of 8 - day ET: (a) MODIS data for Sep.07 - Sep.14, 2015 (b) Simulated results for September 07 - 14, 2015 (c) MODIS data for July 05 - 12, 2016 (d) Simulated results for July 05 - 12, 2016 ................................ .................... 44 Figure 2.10 Time series comparisons between MODIS data (MOD16 product) versus simulation results for 8 - day cat chment avera ged ET ................................ ................ 44 Figure 2.11 Comparison of observed and simulated soil temperatures ............................ 46 Figure 2.12 Comparison of observed and simulated gro undwater temperatures at four borehole locations with mean annual air temperature shown. The blue bands correspond to th e uncertainty associated with ±30% changes in the hydraulic conductivity of the first groundwater layer. RMSE = root - mean - square error ........ 47 xi Figure 2.13 Stream temperature comparisons (observed versus simulated) at the basin outlet (a) Time series comparison. (b) 1:1 plot. RMSE = root - mean - square ........... 50 Figure 2.14 Photographs of the main channel showing some complex features not included in the current modeling (a) a wood debris dam (b) a bedrock outcrop withi n the channel and (c) fallen trees that provide permanent shading in some stream reaches. Photo courtesy: Dr. Mantha S. Phanikumar. ................................ .............. 51 Figure 2.15 Comparison of observed and simulated streambed temperatures at site 1 for different depths: (a) 5 cm, (b) 10 cm, (c) 20 cm, (d) 30 cm, and (e) Close - up views of the areas marked A, B and C in Figure 2.15. ................................ ....................... 54 Figure 2.16 Comparison of observed and simulated streambed temperatures at site 4 for different dep ths: (a) 5 cm, (b) 10 cm, (c) 20 cm, and (d) 30 cm. .............................. 55 Figure 2.17 Stream temperature residuals (difference between simulated and observed stream temperatures at the basin outlet) for two simula tions : baseline and constant ground - water temperature. ................................ ................................ ........................ 58 Figure 2.18 Comparison of observed and simulated streambed temperatures at Site 1 for different depths (a) 5 cm (b) 10 cm (c) 20 cm (d) 30 cm. The red line represents observed data and the blue line the baseline simulation that explicitly simulated groundwater temperature. The green line represents a simulation th at used a constant groundwater temperature based on the annual mean air tempe rature. ...................... 59 Figure 2.19 Comparison of observed and simul ated streambed temperatures at Site 4 for different depths (a) 5 cm (b) 10 cm (c) 20 cm (d) 30 cm. The red line represents observed data an d the blue line the baseline simulation that explicitly simulated groundwater temperature. The green line represents a simulation that used a constant groundwater temperature based on the annual mean air temperature. ...................... 60 Figure 2.20 Sensitivity of sim ulated stream temperature to changes in parameters listed inTable 2.6. Each parameter was changed by ±10%, ±30%, ±50% and changes in hourlystream temperature at the basin outlet relative to the baseli ne simulation results. ................................ ................................ ................................ ....................... 64 Figure 2.21 Sensitivity of simulated stream temperature to changes in parameters listed inTable 2.6. Each parameter was changed by ±10%, ±30%, ±50% and changes in hourl ystream temperature at the basin outlet relative to the baseline simulation results. ................................ ................................ ................................ ....................... 65 Figure 2.22 Daily - average, cross - sectional temperature profiles of stream, streambed and groundwater along a portion of the main stream (between points marked A and D in Figure 2.1) for different seasons. Panel (a) is a schematic of the general land use fe atures along the stream segment. The numbers 1, 2, 3, 4 in (a) denote the streambed temperature sampling l ocations. Each panel in (b) through (e) shows the temperature profile for one representative day in each of the four seasons: (b) Spring, (c) Summer, (d) Autumn, and (e) Winter. Head difference (groundwater xii head minus river stage) variations as a function of distance from point A along the main stream are plotted using the black solid lines using the second Y - axis on the right. Gaining reaches of the st ream correspond to positive head differences while losing portions are associated with negative head diffe rences. ................................ 66 Figure 3.1 Schem atic illustration of the river temperature energy components for Kalamazoo river watershed. ................................ ................................ ...................... 78 Figure 3.2 Map of the Kalamazoo River watershed. Elevation is shown as the color gradient. National Hydrography Dataset (NHD) streams, Stream temperature observation sites, U.S. Geological Survey (USGS) gauges, National Climatic Data Center (NCDC) weather stations a nd Michigan Automatic Weather Network (MAWN) stations are shown. ................................ ................................ ................... 79 Figure 3.3 Land use and land cover map of Kalamazoo River watershed. ...................... 82 Figure 3.4 Schematic illustration of using the square boxes to extract the streamside land use information. ................................ ................................ ................................ ........ 83 Figure 3.5 Stream discharge comparisons between simulated results and obser vations by USGS gauges. ................................ ................................ ................................ ........... 86 Figure 3.6 Monthly watershed - averaged ET comparisons between simulated and remotely sensed MODIS ET products. ................................ ................................ ..... 87 Figure 3.7 Spatial map of yearly averaged evapotranspiration for the Kalamazoo River watershed for the 7 - year period (2003 2009) of (a) simulated output and (b) MODIS data. ................................ ................................ ................................ ........................... 88 Figure 3.8 Plots of simulated versus observed depth to groundwater table (from Well - logic data set) for each computation grid cell. NASH is the Nash - Sutcliffe effic iency metric. ................................ ................................ ................................ ....................... 90 Figure 3.9 Spatial and temporal comparisons between MODIS product and the PAWS+CLM simulated LAI. (a) Time series comparison for four representative vegetation types. (b) Spatial map comparison for June 21 - 28, 2007 and December 19 - 26, 2007. ................................ ................................ ................................ .............. 91 Figure 3.10 10 cm Soil Moisture comparisons of simulated outputs with MAWN (Michigan Automatic Weather Network) station observations at (a) Albion and (b) MSUKBS. Sim is the simulated outputs; Obs is the MAWN station observations . . 94 Figure 3.11 10 cm Soil Temperature comparisons of si mulated outputs with MAWN (Michigan Automatic Weather Network) station observations at (a) Albion and (b) MSUKBS. Sim is the simulated ou tputs; Obs is the MAWN station observations. . 94 Figure 3.12 Performance difference of stream temperature comparisons using different scales of stream side land use information for three selected observation sites: (a) xiii USGS04107850 (b) 1:1 plot of USGS0 4107850; (c) Bear Creek, (d) 1:1 plot of Bear Creek; (e) Spring Brook; (f) 1:1 plot of Spring Brook. OBS is the observation. .... 98 Figure 3.13 Land use types and percentages for different tested land use resolutions of the 1 km stream segment for the presented three observation sites. ............................... 99 Figure 3.14 Stream temperature comparison of simulated results with observations for the eight observation sites: (a) Bear Creek, (b) Silver Creek, (c) Rice Creek, (d) Schnable Brook, (e) Spring Brook , (f) Indian Creek, (g) Marshal, (h) Allegan, USGS04107850. Sim denotes the simulated results; Obs denotes the observations. ................................ ................................ ................................ ................................ 102 Figu re 3.15 Simulated seasonal heat fluxes budgets for the eight observation sites : (a) Bear Creek, (b) Silver Creek, (c) Rice Creek, (d) Schnable Brook, (e) Sprin g Brook, (f) Indian Creek, (g) Marshal, (h) Allegan, USGS04107850. Net radiation heat flux is the sum of short wave radiation and net long wave radiation heat flux. Subsurface heat flux is the sum of streambed conduction heat flux and advective groundwater heat flux. ................................ ................................ ................................ ................. 105 Figure 3.16 Average spatially distributed (a) summer and (b) winter temperatures for the KRW stream network during the simulation period (2003 - 2010). Example reaches shown in box es are discussed in the text. ................................ ............................... 106 Figure 3.17 Temperature difference map for the KRW stream network under a model scenario of riparian deforestation. Example reaches shown in boxes are discussed in the text. ................................ ................................ ................................ .................... 109 Figure 4.1 Map of the Kalamazoo River watershed. Elevation is shown as the col or gradient. National Hydrography Dataset (NHD) streams, Storet nitrogen stations, and nitrate sam pling sites in Baas (2009), U.S. Geological Survey (USGS) gauges, National Climatic Data Center (NCDC) weather stations and Michigan Automatic Weather Network (MAWN) stations are shown. ................................ .................... 116 Figu re 4.2 Flow diagram for the essential processes of PAWS+CLM. The modules in the blue boxes are t he primary processes for the hydrologic components, the modules in the red boxes are the primary processes of temperature simulations, and the modules in the bl ack boxes are primary sub - processes. V/M/H Transfer denotes Vapor/Momentum/Heat Transfer. ................................ ................................ .......... 117 Figure 4.3 Schematic showing the major processes of the nitrogen cycle simulated. .... 118 Figure 4.4 Sketches of Boundary Layer approach. ................................ ......................... 124 Figure 4.5 Vegetation flux and pools structure of CLM 4.5 (adapted from Oleson et al., 2013). ................................ ................................ ................................ ...................... 128 Figure 4.6 stream flow comparisons between simulated results using CLM 4.0, CLM 4.5, and observations by USGS gauges. ................................ ................................ ........ 135 xiv Figure 4.7 Monthly watershed - averaged ET comparisons between simulated results using CLM 4.0, CLM4.5 and remotely sensed MODIS ET products. ............................. 136 Figure 4.8 Comparison between simulated and observed nit rate concentration time series in nine different sampling locations in Baas (2009). The simulated concentrations of nitrite, ammonia and organic nitrogen are also shown. ................................ ......... 138 Figure 4.9 Comparis on between simulated and observed nitrate, nitrite, ammonia and o rganic nitrogen concentration time series of three observation stations in Storet database. NO 3 - sim, NO 2 - sim, organicN - sim and NH 4 - sim are the simulated concentrations of nitrate - N, nitrite - N , organic - N and ammonia - N; NO 3 - obs, NO 2 - obs, organicN - obs and NH 4 - obs are the observed concentrations of nitrate - N, nitrite - N, organic - N and ammonia - N, respectively. ................................ .......................... 139 Figure 4.10 The annual a vera ge nitrate leaching map. ................................ ................... 142 Figure 4.11 Comparison of spatially distributed groundwater concentrations between observations and model results. ................................ ................................ .............. 144 Figure 4.12 Time series comparisons between observations in Baas (2009) and simulated nitrate concentrations for baseline simulation and simulation with temperature dependent reaction rates; (a) zoomed in figure of site Battle Creek at Emmett St. Dam (b) zoomed in figure of site Kalamazoo River at Augusta. NO 3 base is the baseline simulatio n, NO 3 temp is the simulation with temperature dependent reaction rates. ................................ ................................ ................................ ........................ 146 Figure 4.13 Time series comparisons between observations in Storet databse and simulated nitrate concentrations for baseline simulation and simulation with temperature dependent reaction rates. NO 3 - sim (NO 3 - sim - t), NO 2 - sim (NO 2 - sim - t), organicN - sim (organicN - sim - t) a nd NH 4 - sim (NH 4 - sim - t) are the simulated (sim ulated with temperature corrected reaction rates) concentrations of nitrate - N, nitrite - N, organic - N and ammonia - N; NO 3 - obs, NO 2 - obs, organicN - obs and NH 4 - obs are the observed concentrations of nitrate - N, nitrite - N, organic - N and ammonia - N, respectively. ................................ ................................ ................................ ....... 147 Figure 4.14 90% confidence bands of simulated nitrate concentrations resulting from uncertainties of anthropogenic nitrogen sources (fertilizer, ma nure and point sources). NO 3 - sim is the nitrate concentation of baseline simulation, NO 3 - obs is the observed nitrate concentration, uncertainty band denotes the 90% confidence band resulting from uncertainties of anthropog enic nitrogen sources. ............................ 150 Figure 4.15 90% confidence bands of simulated nitrate concentrations resulting from denitrification rate i n rivers. NO 3 - sim is the nitrate concentation of baseline simulation, NO 3 - obs is the observ ed nitrate concentration, uncertainty band denotes the 90% confidence interval resulting from denitrification rate in rivers. .............. 151 Figure 4.16 90% confidence bands of simulated nitrate concentrat ions resulting from nitrification rate in rivers. NO 3 - sim is the nitrate concentation of baseline xv sim ulation, NO 3 - obs is the observed nitrate concentration, uncertainty band denotes the 90% confidence interval resulting from nitrification rate in rivers. .................. 152 Figure 4.17 Comparison between simulated and observed stream nitrate concentration time series near the catchment outlet of WBW. The simulated concentrations of nitrite, ammonia and organic nitrogen are also shown. ................................ .......... 154 1 Chapter 1. Background and Motivation Motivation for developing a new model of nitrogen transport and transformation Hypoxia and eutrophication resulting from excessive nutrient loading are one of the most significant environmental issues around the world (Heisler et al., 2008 ; Galloway et al., 2008 ) . Harmful algal blooms (HABs ) due to eutrophication may result in oxygen depletion and mortality of aquatic species (Anderson et al., 2002) . G reenhouse gas emission and soil acidification could also be stimulated due to exc essive nitrogen emission (Pope et al., 1995 ; Galloway et al., 2003). In the Great L akes area , algal blooms frequently occur in the western part of Lake Erie and the Saginaw Bay ( Hinderer and Murray, 2011 ). T he implementation of the 1972 F ederal Clean W ater Act has effectively reduced the p oint source s loading of nutrients to surface waters in USA . However, it is still a big challenge to control the diffuse nonpoint source pollution , including agricultural fertilizer/manure use , which is considered as the major reason of eutrophication in the United States (USEPA, 1996). In China, rapid indust rialization and urbanization associated with increasing nutrient release exacerbated the situation (Le et al., 2010). Frequent occurrences of eutrophication in Lake s Tai and Chao have caused economic losses of billions of dollars (Le et al., 2010) . In Eur ope and Canada , a critical load approach has been adopted to set regional limits for acceding the potential impacts of acidic deposition re sult ing from atmospheric nitrogen loading (Burns et al., 2008) . The above issues worldwide call for well - informed wat er resources management, agricultural management and decision making. Sound management regulations and decision making require a clear understanding of the c ause - effect relationship s . However, eutrophication and nitrogen enrichment are complex set of coevo lving processes, including 2 anthropogenic , hydrologic and biogeochemical processes , etc . Anthropogenic activities are considered to be the main reason for the e xcessive r iverine nitrogen loading, including agricultural fertilizer/manure applications, land u se development and urbanization ( Carpenter et al., 1998 ; Boyer et al., 2002 ; Beman et al., 2005 ; Galloway et a l., 2008; Thomas et al., 2016 ) . In addition to anthropogenic activi ties, hydrologic processes and biological processes also influence the nitrogen reaction and transport processes at the watershed scale (Breemen et al., 2002) . Temperature, s tream discharge, the structure of drainage network as well as soil moisture are id entified as importa nt factors influencing nitrogen transformation and transport processes ( Wollheim et al., 2006 ; Helton et al., 2018 ; Schaefer and Alber, 2007 ; Miller et al., 2016 ; Porporato et al., 2003 ; 2003 ) . Of particular interest is temp era ture, which regulates the biogeochemical processes , which is found to be the key factor influencing nitrogen processes in watersheds (Schaefer and Alber, 2007 ; Miller et a l., 2016) . To elucidate the cause - effect relationship s , therefore, we need a comprehensive understanding of multiple dynamic processes. Model s are useful tools to link processes that occur at differe nt scales . Statistical models, con ceptual models/ semi - p rocess - based models, and fully distributed process - based models are the primary three types of model s developed so far to study the watershed - scale nitrogen dynamics. Statistical models , such as the SPARROW model (Schwarz et al., 2006; Robertson and Saad, 2011) , have the advantage s of efficiency and small demand of c omputational re sources, while they lac k the capability to address temporal cause - effect relationship s . Conceptual/ semi - process - based models, for example, SWAT (Arnold et al., 1998) , HSPF (Bicknell et al., 1997), I NCA (Whitehead et al., 1998; Wade et al., 2002) and LASCAM (Viney et al., 2000) , describe nitrogen processes in a holistic view while 3 simplifying some important processes, which may result in reduced accuracy . Fully process - based hydrologic models, such as MOHID (Neves, 1985; Trancoso et al., 2009) , have comprehensive descriptions of all important processes by explicitly solving the governing equations . However, they usually have a heavy computational burden and require large amount s of input data to cons truct the model . The present study aims at developi ng a process - based modeling framework that strike s a balance between model complexity and process es representation . As mentioned above, temperature plays an important role in bi ogeochemical processes that potentially impact nitr ogen transformation and transport processes. There re mains a need to develop all integrated catchment - scale temperature module to better estimate the reaction rates of nitrogen processes, which are temper ature dependent. This is the second major objective of this dissertation. Motivation for developing a new watershed scale temperature model Stream temperature is an important variable that affects ecosystem functioning and controls biogeochemical process es in aquatic systems ( Caiss ie, 2006; Allan and Castillo, 2007 ; Webb et al, 2008; Baranov et al., 2016; Folegot et al., 2017 ). Increased stream temperature can n egatively impact water quali ty and the health of aquatic ecosystems (Roth et al., 2010 ; Folegot et al., 2018 ) . Stream thermal regimes are primarily driven by climatic conditions and influenced by a host o f other factors, including t opographic conditions, stream discharge, land cover near the stream and the interactions with subsurface (Caissie, 2006 ; Hannah and Garner, 2015 ). Researchers found that riparian vegetation surrounding the river channels play a significant role in affectin g the stream temperatures (Roth et al., 2010 ; Sun et al., 2015 ; Cao et al., 2016 ; Garner et al., 2017). Additionally, spatial heterogeneities 4 and coevolution of surface and sub - surface processes make understanding the pr ocesses that drive stream te mperature dynamics a cha llenging task (Caissie et al., 2007 ; Caissie et al., 2017; Halloran et al., 2017). Primarily, two types of stream te mperature models have been developed t o date: regression models that make use of statistical linkages between meterological and geophysical conditions to predict stream temperatures ( Mohseni et a l., 1998; Benyahya et al., 2007; Jackson et al., 2017 , 2018 ), and mechanistic models based on conservation of energy that directly simulate the underlying stream temperature dynamics ( St - Hilaire et al., 2000 ; Cox and Bolte, 2007 ; Loinaz et al., 2013; see review by Dugdale et al., 2017). Mechanistic models usuall y demonstrate clear cause - effect relat ionships because of the direct descriptions of the underlying controlling processes in their governing equations ( Caissie et al., 2007 ; MacDonald et al., 2014 ; Sun et al., 2015 ; Gallice et al., 2016 ). Some mechanistic stream temperature models have incorporated the thermal advection and riverbed conduction fluxes (Haag and Luce, 2008 ; MacDonald et al., 2014 ; Gallice et al., 2016) . They considered the heat flux through the riverbed as being propotional to the di fference between river and riverbed te mpertures at a certain depth (Moore et al., 2005). However, these models used lump ed parameters to estimate the hyporheic exchange water flux which lack an explicit expression of the water - heat exchange at the interfac e of stream/GW system. Other researche rs used fully three - dimensional (3D) models (Brookfield et al., 2009) or cross se ctional two - dimensional models (Halloran et al., 2017) to simulate the integrated surface/subsurface thermal tran s port. However, due to the computational expense involved in solving fully 3D equations, these models could only be applied to relatively small portions of a river or a single reach. Therefore, there exists a 5 strong s cientific motivation to develop a water shed - scale stream temperature model that could efficiently tr ack wate r and heat fluxes coevolving in the surface and subsurface domains . Given the challenges and opportunities, t his dissertation research describes the development and application of an integrated , process - oriented coupled water - temperature - nitrogen model based on the framework of PAWS+CLM (Shen and Phanik umar, 2010; Shen et al., 2013). The mod el include s a holistic framework of nitrogen transport and transformations with multiple sources and interactions . Specifically, th e model include s the processes of nitrogen transport and transformations in the domains of streams, Groundwater , vadose zone as well as overland flow. Meanwhile, an interactive stream subsurface temperature model that takes into account the effect of ve getation processes is coupled with the nitrogen model to correct the reaction rates (e.g., denitrification in sediments) . R ar ely have previous studies applied integrated and fully process - oriented water - temperature - nitrogen model for studying the nitrogen dynamics over long period s of time at the watershed scale. The advantages of realizing this is two - fold: o n one hand, the m odel could suf ficiently take control of the integrated hydrologic, ecologic, biological and anthropogenic effects on the watershed - scale temperature and nitrogen dynamic s ; on the other hand, computational efficiency and long term predications could inform effect ive decision making. This novel model framework is evaluated by applying the model to two different watersheds with different sizes and climate conditions , i.e. the Wood Brook w atershed (WBW) located in Birmingham, UK and the Kalamazoo River watershe d (KRW) located in Michigan, USA . The other chapters of this disse r tation are organized as follows. I n 6 chapter 2, mathematical details of the development of an integrated, catchment - scale framework to model stream, soil, streambed and groundwater temperat ures under the influence of hydrologic and vegetation dynamics are presented. The he model performance was tested by applying it to the WBW with a drainage area of 3.5 km² . In C hapter 3 , the application of the developed temperature modeling framework is ex tended to the KRW with a drainage area of 5200 km² and with more than 200 tributaries. M oreover , the effects of resolving the spatial resolution of vegetation heterogeneity are evaluated to reach a better approximation of the vegetation effect with a relat ively coarse model resolution, i.e. 1000 m × 1000 m . I n chapter 4, a catchment - scale framework is developed to simulate integrated hydrologic, temperature and nitrogen processes including reactions and transport of multiple nitrogen species. Th e applicatio n s of the model framework are discussed in detail and are evaluated with multiple observations and literature values. 7 Chapter 2. Evaluating a coupled phenology surface energy balance model to understand stream subsurface temperature dynamics In this chapter, an integrated, catchment - scale framework to model stream, soil, streambed and groundwater temperatures is developed under the influence of hydrologic and vegetation dynamics in a mixed land use catchment in central England. The phenology and surf ace ener gy modules in the coupled model were used to quantify the impacts of vegetation processes on radiation fluxes (e.g., canopy shading and the effect of vegetation growth on optical parameters). The model enabled accurate simulations of the movement a nd parti tioning of water and thermal fluxes in different hydrologic domains with R 2 values of observed and simulated temperatures in the range 0.60 - 0.87. Simulated groundwater heads and stream stages allowed the identification of gaining and losing portion s of str eam reaches and the estimation of Darcy fluxes. Simulation results show significantly dampened diel streambed temperature fluctuations below 0.3 m in gaining reaches , while in losing reaches the diel fluctuations showed relatively strong fluctuatio ns below 0.3 m. The model enabled evaluation of the relative contributions of different processes to the stream thermal budget. Results indicate that net radiation was the dominant heat source while latent heat flux was the primary heat sink. The model pro vides a useful tool to explicitly simulate water and heat fluxes for analysis of temperature - dependent reaction rates in biogeochemical analyses. Introduction Previous research has shown that canopy shading and vegetation growth surrounding the stream cha nnels pl ay an important role in affecting stream temperatures ( e.g., Roth et al., 2010 ; Sun et al., 2015 ; Cao et al., 2016 ; Garner et al., 2017; Dugdale et al., 2018). However, 8 very few catchment - scale temperature models exist that explicitly consider interactions between surface and subsurface hydrologic processes and vegetation processes at scales relevant to mechanistic modeling of stream temperature (e.g., Davison et al., 2015, 2018; Sulis et a l., 2017; Loicq et al., 2018). Relatively small - scale riparian vegetation processes close to the stream banks control canopy shading, as do factors such as spatial heterogeneity of vegetation density and temporal aspects of vegetation growth (Hannah et al . , 2008). Vegetation type determines physical vegetation characteristics, such as vegetation height and leaf and branch size, which affect stream temperature while also influencing the riparian microclimate including air temperature, wind speed and relativ e humidity (Garner et al., 2017). The ability to accurately represent these processes while spanning the horizontal and vertical length scales of catchments represents a major challenge to mechanistic stream temperature modeling and the present work is an a ttempt to address these challenges. Different approaches were used to represent vegetation processes in stream temperature models in the past including the use of tunable shading coefficients (Herb & Stefan, 2011; MacDonald et al. 2014). Advances in land s urface modeling and the availability of land surface models (LSMs) in recent years opened up the possibility to accurately simulate vegetation processes by coupling hydrologic models with LSMs. For example, Li et al. (2015 ) developed a large - scale stream t emperature model within the Community Earth System Model (CESM) framework (Gent et al., 2011) by coupling the Community Land Model CLM ver 4.0 (Oleson et al., 2010) with a river routing model. The focus of their work was on studying the effects of human i m pacts such as reservoir regulation, and the model in its current form does not explicitly include interactions between streams and the subsurface domain. 9 Recently, studies of controlling processes have extended from the air / surface - water interface to t h e interface of surface - water/groundwater and surface - water/streambed ( Caissie et al., 2014 ; Caissie and Luce., 2 017; Halloran et al., 2017). In stream reaches where groundwater contribution is significant, heat fluxes from the subsurface can represent a substantial fraction of the stream energy budget ( Hannah et al., 2004). Studies have shown the effectiveness of u s ing the streambed thermal regime to quantify vertical water fluxes between the surface and sub - surface domains ( Anderson, 2005 ; Gordon et al., 2012 ; Vandersteen et al., 2015 ; Anibas et al., 2009 ). Some mechanistic stream temperature models have incorporated the thermal advection and streambed conduction fluxes ( Haag and Luce, 2008 ; MacDonald et al., 2014 ; Gallice et al., 2016 ). They considered the heat flux through the streambed as being proportional to the difference betwee n stream surface water and streambed temperatures at a given depth (Moore et al., 2005). Various approaches to estimate streambed temperatures were reported in the literature. For example, Haag and Luce (2008) introduced the concept of effective streambed t emperature as a function of two lumped parameters to estimate the streambed conduction flux. Gallice et al. (2016) assumed that the streambed temperature is the same as simulated soil temperature. MacDonald et al. (2014) used an empirical function that d e pends on the air temperature to estimate the streambed temperature. Caissie et al. (2014) used the vertical one - dimensional model to estimate the Darcy and streambed heat fluxes. Other researchers have used fully three - dimensional (3D) models ( Brookfield et al., 2009) or cross - sectional two - dimensional models (Halloran et al., 2017 ) to simulate the integrated surface/subsurface thermal t ransport . However, due to the computational expense involved in solving fully 3D eq uations these models were only be applied to relatively small portions 10 of a river or a single reach. There is a need for computationally - efficient, catchment - scale (as oppos ed to site - specific) models of coupled water and heat balance that can be used to t est hypotheses, identify parameters and evaluate the relative importance of individual processes. The objectives of this chapter are to (1) develop a catchment - scale stream subsurface temperature model that takes into account the effect of vegetation pr ocesses on radiative fluxes while explicitly simulating interactions between surface and subsurface domains (2) test the model against detailed field observations from diffe rent hydrologic domains and (3) use the model to understand the key factors that co ntrol stream thermal budgets in different seasons. In addition, questions related to scale ( the relative sizes of the catchment, stream widths and the grid sizes) as well as the ability of the model to represent thermal fluxes between hydrologic domains ac curately will be addressed. The following qu estions will be addressed: (1) i f canopy radiative fluxes are computed over grid cells whose resolution is typically larger than stream widths, is it possible to simulate stream and subsurface temperatures accura tely within an integrated modeling framework where surface and subsurface domains are coupled? (2) In a grid cell with multiple land uses, is it possible to simulate stream temperatures accurately by representing the sub - grid scale variability in surface r adiative fluxes using an area - weighted approach? (3) What is the impact of using mean air temperature as a proxy for groundwater temperature instead of explicitly simulating subsurface temperatures? (4) What are the dominant components of the stream therma l budget in different seasons? (5) Can the model correctly identify gaining and losing portions of stream reaches? 11 M aterials and Methods 2.2.1. Site Description and Observationa l data Data supporting the model building, parameterization and validation were col lected from the Wood Brook catchment at the Birmingham Institute of Forest Research (www.birmingham.ac.uk/bifor) field site in Staffordshire, UK, between March 2015 and Apri l 2017 (Figure 2. 1). The Wood Brook stream drains a 3.5 km² catchment ranging in el evation from 90 to 150 m asl and comprising of a mixture of arable farmland with juvenile and mature deciduous woodland. The catchment geology is dominated by Permotriassic sandstone, overlain by deposits of glacial till with up to approximately 10 m thick ness as well as sandy clay organic - rich top soils between 0.15 and 0.6 m thickness (Blaen et al., 2017). The Wood Brook stream was instrumented with hydrometeorological and water quality monitoring along a 1000 m study reach above the outflow of the catchm ent in an area of mature deciduous woodland (Figure 2. 1). Woodland vegetation is dominated by English oak (Quercus robur) with an understory layer of hazel (Corylus avellana ), hawthorn (Crataegus spp.), and sycamore (Acer pseudoplatanus). Adjacent to the s tream, there is additionally a presence of common alder (Alnus glutinosa), goat willow (Salix caprea) and wych elm (Ulmus glabra). Dense canopy cover results in high levels of shading along the entire reach during spring and summer months. Average stream width in the watershed is 1.2 m, and average depth is 0.3 m. The average tree height across Mill Haft is 10.6 3.8 m. The maximum (i.e. the canopy) height is 25.0 m. 12 Figure 2 . 1 Map of the Wood Brook watershed with land use and land cover. Locations of stream flow sensors, borehole sensors, streambed temperature sensors and soil moisture and tempe rature sensor are shown. Elevation is shown as the color gradient in the sub - map. 13 A stream monitoring station was established at the catchment outlet with a combined pressure transducer and thermistor (Adcon, Austria) for stage and water temperature measu rements. A stage - discharge relationship (R² = 0.89) was established from salt dilution gauging measurements (Hudson & Fraser, 2005; Blaen et al., 2017). Measurements were acquired hourly and transmitted via a telemetry system to an internet server for data storage. Vertical temperature lances (Tempcon, UK) were installed in the streambed at four locations along the study reach (Figure 2. 1). Subsurface temperatures were measured at 15 min reso lution at 5, 10, 20 and 30 cm below the streambed and stored local ly on HOBO U12 data loggers (Onset, MA, USA). Local groundwater level and water temperature was monitored hourly using Mini - Divers (Van Essen Instruments B.V., Netherlands) in four boreholes within the catchment, plus in a further borehole adjacent to the catchment maintained by the Environment Agency (England). Volumetric soil moisture content was measured at 10 cm depth every 15 min by a 5TM probe (Decagon Devices, WA, USA) situated in a yo ung woodland clearing towards the south of the catchment (elevatio n 100 m asl). Air temperature was measured at 15 min resolution using a Vaisala HMP155 probe attached to a meteorological flux tower at 1.2 m above the ground in the south of the catchment. Local precipitation was recorded using ARG100 tipping bucket rain gauges (EML, UK) in the centre of the catchment (Figure 2. 1) and approximately 200 m south of the catchment (not shown). Both rain gauges were positioned away from trees or other structures that might otherwise have impacted measurements. Additional climat e data (hourly air temperature, dew point temperature cloud cover and average wind speed) were obtained from the nearby (within 25 km) 14 Shawbury station in the Met Office MIDAS database (Met Office, 2012). For the evapotranspiration (ET) observational data , we used the level - 4, 500 m resolution, 8 - day Moderate Resolution Imaging Spectroradiometer (MODIS) Global ET product downloaded from the NASA Earth Data website (https://search.earthdata.n asa.gov/search). MODIS Leaf Area Index (LAI) data (MOD15A2H produc t based on level - 4, 500 m resolution, 8 - day LAI) were also downloaded from the same website and used to compare with simulated LAI results. Soil properties were extracted from the STATSGO so il database (Soil Survey Staff, 2017). Van Genuchten soil water re tention parameters and unsaturated hydraulic conductivities (Warrick, 2003) were thereafter estimated using pedotransfer functions as implemented in the Rosetta software (Zhang & Schaap, 201 7; Schaap et al., 2001), which uses an artificial neutral network optimized model. Groundwater hydraulic conductivities were generated based on typical hydraulic conductivities of sandstone and borehole data from pumping tests. Ordinary kriging was used fo r interpolation of point data to produce an initial hydraulic cond uctivity field. The groundwater hydraulic conductivity field was optimized during model calibration. 2.2.2. The hydrologic model The temperature module was developed within the framework of the in tegrated hydrologic model PAWS (Process - based Adaptive Watershed S imulator, Shen et al., 2016; Shen et al., 2014; Shen et al., 2013). PAWS simulates key hydrologic processes including surface and subsurface flow, channel flow, topography - induced overland flow, and soil water processes. Detailed vegetation processes and surface energy balance are solved via coupling to the Community Land Model, CLM 4. 0 (Oleson et al., 2010; Shen et al., 2011). 15 The model has been extensively tested using field observations and remotely - sensed data (Shen et al., 2014 ; Niu et al., 2014 ; Niu and Phanikumar, 2015; Safaie et al., 2017 ). The governing equations of PAWS are solved on structured grids which discretize the model into different hydrologic domains and layers. On the topmost (overland flow) layer, runoff occurs when the depth of ponding domain is in excess of the int erception depth which is governed by the diffusive wave equation. Infiltration and evaporation occur in the ponding domain; water may backfill into the ponding domain during flood conditions. A one - dimensional diffusive wave equation was used to describe t he stream flow dynamics. The stream segments are connected to the groundwater domain or the vadose zone depending on the streambed elevation and the depth of water table. Th e leakance concept (Orhan Gunduz, 2005 ) was used to compute the water flux between the groundwater domain and the stream channels. Soil water dynamics in the vadose zone are governed by the one - dimensional Richards equation with vegetation uptake representing a sin k term. Two - way dynamic interactions of different hydrologic components within PAWS+CLM (referred to as PAWS hereinafter) provide a convenient framework for studying water partitions (Niu et al., 2014) , understanding controls on hydrologic and vegetation processes ( Shen et al., 2013 ) and contaminant fate and transport ( Niu and Phanikumar, 2015 ). A grid resolution of 20m × 20m (which produced a mesh of 113×118 cells) was used to discr etize the Wood Brook catchment domain horizontally. 20 soil layers were used to dis cretize the domain between the land surface and the initial groundwater table with an initial depth of approximately 5 m. The spatial resolution was refined near the land su rface to account for the sharp gradients in fluxes appearing there. An adaptive cel l size was used for the bottom cell of the soil column, which was adjusted based on groundwater table 16 fluctuations (Shen and Phanikumar, 2010). Groundwater flow in PAWS is g overned by the quasi - ater layers were employed underneath the vadose zone to represent the aquifers. The first layer represents the unconfined aquifer with an initial thickness of 15 m, the bott om elevation of which was 20 m below the land surface. The second layer represents a deep bedrock aquifer with low hydraulic conductivities. The model was driven by hourly climate data. A nearest neighbor interpolation method was used to interpolate the sp atially distributed precipitation inputs using data from the two rain gauges. The d ifferential evolution algorithm (Price et al., 2005) was used to optimize simulated stream discharge and groundwater heads. An assumption often made in catchment models is t hat groundwater and surface water divides match perfectly (Dingman, 2015). If this assumption is violated, then groundwater can enter and leave the catchment over unknown portions of the surface water divide (which represents the catchment boundary based o n a delineation of topographic features). A challenge associated with modeling grou ndwater flow dynamics in a small catchment such as the Wood Brook is that surface water and groundwater divides may not coincide. To address this concern, we developed a reg ional - scale groundwater model of the Staffordshire area covering approximately 100 km2 containing the Wood Brook catchment and obtained steady - state groundwater head distributions in and around the catchment boundary. The regional - scale model then guided t he selection of boundary conditions (e.g. no - flow over the east boundary, constant head condition near the basin outlet etc.) for the groundwater model used in our catchment simulations. The integrated surface subsurface model used in our work couples a newly developed temperature model (described here) with an existing process - based hydrologic model that 17 incorporates complete surface and subsurface processes. F igure 2. 2 is a flow chart of the integrated model that shows how the different modules are con nected. Stream temperature dynamics were simulated using the one - dimensional advect ive heat transfer equation with multiple heat sources and sinks (described in Section 2.3). A quasi three - dimensional heat transport equation was used to simulate groundwat er temperature dynamics in the unconfined aquifer. Streambed temperature dynamics w ere simulated using a vertical one - dimensional heat transfer equation based on the assumption that lateral heat flux is of minor importance and can be ignored (Caissie et al ., 2014). The temperature module uses an explicit, two - way coupling scheme summariz ed below. While solving the stream temperature equation the latest simulated values of streambed and soil temperatures are used to approximate the conductive streambed and l ateral advective heat fluxes respectively. Meanwhile, simulated stream temperature is employed as the top boundary condition for the streambed temperature equation. The groundwater temperature module uses the bottom soil layer temperature and streambed tem perature to calculate the advective and conductive heat fluxes to groundwater. As a feedback, the groundwater temperature is used as the bottom boundary condition of soil temperature and streambed temperature equations. The initial values of stream tempera ture, streambed temperature, soil temperature and groundwater temperature are set a s the mean daily air temperature of the first day of simulation. The simulation allows sufficient time (8 months) for the model to spin up; the simulation started from Septe mber 2014 while the model calibration started from May 2015. Parameters adjusted du ring model calibration are listed in Tables 2.1 and 2.2 . Model performance was tested using detailed stream, streambed, soil and groundwater temperature measurements in a sm all mixed - use headwater catchment in central England, 18 representing spatial heteroge neity in land use and land cover patterns, as well as hydrometeorological dynamics. Table 2 . 1 Model calibration parameters for hydrology Symbo l Parameter Meaning Parameter in ( Lai and Katul, 2000) , root Efficiency function ice Scale - dependent freezing fraction parameter as in ( Niu and Yang , 2006) K 1 (m day - 1 ) First layer G roundwater Hydraulic Conductivity K 2 (m day - 1 ) Second layer G roundwater Hydraulic Conduc tivity K s (m day - 1 ) Soil Saturated Hydraulic conductivity N Van Genuchten parameter A (m - 1 ) Van Genuchten parameter l (m) Length of flow path for runoff contribution to overland flow domain h o (m) overland flow ground inception depth K r (m day - 1 ) strea mbed hydraulic conductivity Table 2 . 2 Model parameters for temperature module S ymbol Parameter Meaning k b ( W m - 1 K - 1 ) Streambed effective thermal conductivity (m) C haracteristic depth to calculate the streambed conduction heat flux C w W ind speed sheltering factor K r ( m day - 1 ) Streambed hydraulic conductivity K 1 (m day - 1 ) Hydraulic conductivity of the first groundwater layer 19 Figure 2 . 2 Program flow diagram for the essential processes of PAWS+CLM. The modules in the blue boxes are the primary processes for the hydrologic components, the modules in the red boxes are the primary processes of temperature si mulations, and the modules in the black boxes are primary sub - processes. V/M/H Transfer denotes Vapor/Momentum/Heat Transfer. 2.2.3. The stream temperature module Stream temperature is simulated using the one - dimensional heat transfer equation: ( 2. 1) where T s denotes the stream temperature, Q v (W m - 2 ) corresponds to the sum of the net heat fluxes excluding the groundwater advective heat flux, Q adv (W m - 2 ) denotes the groundwater advective heat flux, C ( J kg - 1 K - 1 ) is the specific heat, (kg m - 3 ) is the density of water, A (m 2 ), Q (m 3 s - 1 ), and W (m) are the stream segment cross - sectional area, the stream discharge, and the wetted width ( W , m), which were directly obtained from PAWS for the stream segment. Q v includes heat fluxes from di fferent components which can be summarized: ( 2. 2) 20 where Q s is the net shortwave radiation, Q l is the net longwave radiation positive toward the stream surface, Q h is the sensible heat flux, Q e is the latent hea t flux, Q f is the friction heat flux, Q b is the heat flux from streambed conduction, and Q m is the energy loss due to the melting of snow if there is snow accumulation on the stream surface. All heat flux terms in equation 2. 2 and Q adv have units of W m - 2 . Equation 2. 1 was solved using the hyb rid Lagrangian - Eulerian method based on an operator - splitting strategy (Niu and Phanikumar, 2015). Canopy shading and vegetation growth surrounding the stream channels have an obvious impact on the solar radiati on reaching and absorbed by the stream cha nnels. Radiative transfer in the presence of vegetative canopies was simulated via CLM 4.0 using the two - stream approximation of Dickinson and Sellers (Dickinson, 1983; Sellers, 1985) which provided solutions for forward and backward radiative fluxes in a n absorbing, scattering medium. The advantage of this approach is that via CLM, the model keeps track of radiative fluxes over multiple wavebands while using appropriate ecophysiological parameters to account for c hanges in optical variability between diff erent plant and tree species. Evapotranspiration (ET) is an important process that controls the riparian microclimate including water vapor, momentum and heat fluxes. The Penman - Monteith equation for ET uses the we t bulb temperature to approximate the surf ace temperature of vegetation while CLM computes the leaf temperature by solving the coupled heat transfer equations and estimates ET based on canopy and aerodynamic resistances. Canopy resistances are computed tog ether with different models for photosynth esis for C3 and C4 plant categories while aerodynamic resistances are computed using the Monin - Obukhov similarity theory (Zeng et al., 1998; Kundu et al., 2015). A potential advantage of this 21 detailed approach is t hat it provides a mechanistic description of CO 2 assimilation and hence is suitable for evaluating the impacts of land use and climate change scenarios on stream temperatures (Oleson et al., 2010; Shen et al., 2013). The two - stream approximation of Seller s (1985) and Dickinson (1983) yields the f ollowing differential equations (Bonan, 1996): (2.3) (2.4) where and denote the upward and downward diffuse radiative fluxes if assuming the incident solar flux is unity, L and S denote the exposed LAI and stem area index respectively, denotes the cosine of the zenith angle of the incident beam, and are two upscattering parameters which are applied for diffuse and direct beam radiation, is a coefficient describing the scattering effect, and K denotes the optica l depth of direct beam per unit leaf and stem area, which can be calculated as: ( 2.5 ) Here , denotes the relative projected area of leaf an d stem in the direction , which can be calculated as ( 2.6 ) where the calculations of the and are as follows: ( 2.7 ) 22 Here, denotes the departu re of leaf angles from a random distribution. For horizontal leaves, is +1, for v ertical leaves, is - 1, and is 0 for random leaves. in equations 2.3 and 2 .4 is the average inverse diffuse optical depth within a unit leaf and stem area, which can be calcula ted as: ( 2.8) where denotes the scattered flux direction. Equations 2.3 and 2.4 are used to calculate the partitions of the direct and diffuse radiation into three co mponents: (a) absorbed by the vegetation, (b) reflected by the vegetation and (c) penetrated through the vegetation for different w avebands. The optical parameters of , , in Equa tions 2.3 and 2.4 are dependent on wavelength, the calculations of which are described in Sellers (1985) and in the technical description of CLM 4.0 ( Oleson et al., 2010). Additional details including optical properties of plant functional types (PFTs) for different tree a nd shrub species as well as leaf and stem optical properties and snow properties are available in CLM documentation (Oleson et al., 2010). Details of the remaining components in equation 2.2 and the calculations of groundwater advective he at flux Q adv are described as follows. The incident shortwave (solar) radiation flux coming from the atmosphere was calculated as the sum of the direct beam and the diffuse solar fluxes: ( 2.9 ) where is the net solar radiation flux coming from the atmosphere, denotes the direct beam flux coming from the direction of the sun ( denotes the solar zenith angle 23 and denotes the wavelength) while is the diffuse solar flux (W m - 2 ). In the present work, and were calculated following the method of (Spokas and Forcella, 2006 ). The atmospheric solar radiat ion reaching the land surface was partitioned into components associated with vegetated and ground (that is, non - vegetated) surfaces. For vegetated surfaces, before solar radiation reaches the ground surface (or the water surface if trees are close to the stream), canopy will first intercept and absorb a portion of the solar fluxes as shown in Figure 2 (a). If L and S denote the exposed LAI and stem area index, K the optical depth of direct beam per unit leaf and stem area, then considering unit incident di rect and diffuse solar fluxes, the direct beam flux that is transmitted through the canopy is . Similarly, the portions of direct beam and diffuse fluxes absorbed by the canopy per unit incident flux are: ( 2.10 ) ( 2.11 ) Here the symbols and denote portions of upward and downward diffuse radiation fluxes per unit incident flux, denotes wavel ength and the superscript denotes direct beam flux coming from the direction of the sun ( is the cosine of the zenith angle of the incident beam). The fluxes most relevant for stream temperature dynamic s are the fluxes received by the ground or the stream surface below the canopy and these are and which denote the downward diffuse fluxes below the vegetation per unit incident direct beam and diffuse radiation at the top of the canopy. In equations 2.10 and 2.11 , ground albedos are denoted by the symbols and for direct and diffuse radiation, respectively. 24 Figure 2 . 3 Sketches of the radiation fluxes partitions for (a) short - wave radiation fluxes (b) long - wave radiation fluxes (c) weighted - average heat flux computation for grid cells with mixed land use The variables K , L and S are dynamically simulated in the phenology module of CLM4.0 and used in PAWS ( Shen et al., 2013 ) . Extending the above analysis for unit incident fluxes to the actual fluxes either directly measured at the site or compu ted, the solar radiation absorbed by the vegetation and the ground surface over all wavelengths was then calculated as below: ( 2.12 ) ( 2.13 ) For non - vegetated surfaces such as bare soil, open stream reaches, lakes and urban areas, all of the unit incident flux is transmitted down, therefore , and ; thereby the shortw ave radiation absorbed by the ground surface is the total incoming solar radiation absorbed by the su rfaces (Figure 2.3 (a)) : ( 2.14 ) 25 For streams the direct beam and diffuse albedos are calculated following Bonan (1996): ( 2.15 ) Multiple land use land cover (LULC) classes general ly exist within a single grid cell representing the land surface (herein referred to as the ground cell). To reduce computational effort, LULC maps wer e reclassified into plant functional types and a small number of dominant land use types were used to rep resent land use within a grid cell ( Shen et al., 2014). Considering the fine grid size used in our simulations, three dominant land uses were used in e ach grid cell in the present work. The incoming shortwave radiation absorbed by each grid cell was comput ed as the weighted sum of the radiation absorbed by the three dominant PFTs in each cell. For example, in a grid cell with trees, grass and water as th e dominant land uses as shown in Figure 2 .3 (c), the SW radiation absorbed was computed as: (2.16) where is the area fraction of the total grid cell occupied by the trees, is the SW radiation absorbed by the portion of the ground cell covered with trees and so on. This equation can be written in a more general form as below: ( 2.17 ) where Q s is the solar radia tion absorbed by the stream segment, n = 3 is the number of PFTs within the ground cell, and F ( p ) is the ar ea fraction of the pth PFT within the cell. This operation keeps control of different scenarios: a) for large stream segments, water will be the domi nant PFT within the cell; b) For small stream segments surrounded by heavy vegetation, vegetation will be t he dominant PFT; c) The impact of seasonal vegetation 26 growth cycles on the radiation fluxes absorbed by vegetation are explicitly represented using d ynamic simulations of vegetation growth cycles (via the parameters K, L and S ) which are computed in the ph enology module of CLM. Net longwave radiation Q l is equal to the downward incoming long wave radiation minus the portion emitted from the stream su rface ( Figure 2.3 (b)) : ( 2.18 ) where is the downward longwave radiation into the stream surface and is the backward longwave radiation emitted from the stream surface. The incoming longwave radiation fluxes into the stream include the longwave radiation from the atmosphe re, surrounding vegetation and topography. Several models take the air temperature as the mean temperature of the surrounding objects that are emitting the longwave radiation to the streams (e.g. Leach & Moore, 2010 ; MacDonald et al., 2014 ; Cheng & Wiley, 2016 ). Since PAWS explicitly simulates vegetation tem perature via CLM, we computed the net longwave radiation from vege tation close to the stream surfaces using the weighted sum approach for different LULC types described earlier for shortwave radiation. The computed longwave radiation flux was then applied as an input to the stream segment within the ground cell. The lo ngwave radiation from the atmosphere was estimated following the empirical equation of ( Konzelmann et al. , 1994) : ( 2.19 ) where is the downward atm ospheric longwave radiation, and are the clear sky emittance and overcast emittances respectively, n is the cloud cover fraction, is the Stefan - Boltzmann constant (5.67 ×10 - 8 W m - 2 K 4 ), p = 3 is an empirical coefficient, and T a 27 is the air temperature (K). Follow ing Konzelmann et al. (1994 ), was taken as 0.96 and was estimated using the relation ( 2.20 ) Vegetation or bare ground (bare soil, lakes or urban areas) directly take as the longwave radiation input (Figure 2 (b)). However, the longwave radiation input into the ground surface covered by vegetation, , was calculated based on the vegetation temperature T v (simulated in CLM) as below: ( 2.21 ) where is the emissivity of vegetation and was calculated as ( Oleson et al. , 2010): ( 2.22 ) where L and S are the leaf and stem area indices as before and is the average inverse optical depth for longwave radiation. Similarly, we used the longwave back radiation model from CLM based on the Stefan - Boltzmann Law , in which the upward longwave radiation from the ground, vegetation or water was calculated as: ( 2.23 ) where is a s tep function (equal to zero when the sum of exposed leaf and stem areas was less than 0.05 and one otherwise), the atmosphere is the upward longwave radiation from the vegetation/soil system when the sum of exposed leaf and stem areas was larger than 0.05, is the land use - dependent ground emissivity and and are the snow/soil surface temperatures at the current and previ ous time steps, respectively. 28 F or vegetated surfaces, the above equation becomes: ( 2.24 ) where L vg was calculated as: ( 2.25 ) For non - vegetated surfaces, equation 2.25 can be simplified as: ( 2.26 ) For streams, we omitted the first and last terms which are negligible (the last term represents changes in stream temperature within a time step) hence, t he above equation can be further simplified as: (2 .2 7 ) where T sur is the stream surface temperature (K) at the current time step and is assumed to be the same as the stream temperature T s for shallow streams, and = 0.96 is the stream water emissivity. Therefore, the net longwave r adiation flux for the ground surface was calculated as: ( 2.28 ) For the stream surface = . The weighted sum approach used for Q s (shortwave radiation) was also used to calcula te Q l combining the different land uses surrounding the stream : (2.29) In addition to the options available in CLM for the computation of latent ( ) and sensi ble 29 ( ) heat fluxes, two general expressions for latent and sensible heat fluxes are also included in the current module which are widely used in literatures ( Martin & McCutc heon , 1998 ; H erb & Stefan , 2011) : ( 2.30 ) where is the water density (m 3 kg - 1 ), L w is the latent heat of evaporation which is approximately 2.4 ×10 6 (J kg - 1 ), e w (mb) is the sa turated vapor pressure at the water surface temperature, e a (mb) is the vapor pressure of air, is the wind speed function, u w is the wind speed (m s - 1 ), a and b are empirical coefficients. Various wind functions hav e been proposed in the past and they were summarized by Edinger et al. (1974) , Theurer et al., (1984) and Martin & McCutcheon (1998) . We used the coefficients proposed by The urer et al. (1984) for coefficients a an d b (a = 2.25×10 - 9 , b = 9.40×10 - 9 ). The observed wind speed was adjusted with a wind sheltering coefficient, C w , following a similar approach adopted by Herb & Stefan (2011) : (2.31 ) where (m s - 1 ) is the observed wind speed. Commonly used formulations for calculating e w and e a suggested by the U.S. Army Corps of Engineers ( Environmental Laboratory , 1985) were used: ( 2.32 ) ( 2.33 ) where T d is the dew point temperature directly obtained from the UK Met Office MIDAS database ( Met Office , 2012) . 30 Sensible heat flux is generally considered to represent only a small component of the stream heat budget. A co mmonly used relationship between sensib le and latent heat was used to calculate convective heat flux due to sensible heat transfer: ( 2.34 ) where the Bowen ratio ( Dingman , 2015) was calculated as: ( 2.35 ) where C B is the psychrometric constant equal to 0.61 (mb K - 1 ), P a (mb) is the atmospheric pressure and P (mb) is the reference air pressure at mean sea level. We c omputed the friction heat flux using an equation from Theurer et al. (1984) and MacDo nald et al. (2014): ( 2.36 ) where S o is the stream segment slope within the grid cell, W is the stream wetted width. St reambed conduction heat flux was estimated as: ( 2.37 ) where k b (W m - 1 K - 1 ) is the effective streambed thermal conductivity and T b (K) is the temperature at a characteristic depth (m) beneath the streambed. ( Herb & Stefan , 2011) use d a characteristic depth of 1 m and k b = 1.00 W m - 1 K - 1 . Moore et al., (2005) and MacDonald et al. , (2014) , however, used 0.05 m and a k b of 2.6 (W m - 1 K - 1 ). After examining the sensitivity of model results, we used 0.5 m and k b = 2.21 (W m - 1 K - 1 ). Depending on the climatic region, snow/ice processes may control stream temperatures during winter months. If snow precipitated on the stream surface (lar ger than 0.5 kg m - 2 ) 31 and if the stream temperature was larger than the freezing temperature of water ( T f ), then the surface temperature was reset to T f . The sensible and latent heat fluxes were then recalculated using T f . Thereafter the net heat flux was recalculated using Eq. ( 2.2 ). There wi ll be energy available to melt snow if the recalculated Q net was larger than zero ( Oleson et al. , 2010): ( 2.38 ) where L f is the latent heat of fusion, equal to 3.337 10 5 (J kg - 1 ); W sno denotes the snow accumulation (kg m - 2 ). Heat flux from groundwater advection , Q adv , is given by the following equations ( Caissie et al., 2014) : ( 2.39 ) where q gw (m 2 s - 1 ) is the groundwater seepage volume per segment length (simulated by PAWS hydrology modules) and calculated using the leakance concept ( Gunduz & Aral , 2005) , T soil is the ave rage soil temperature (in the top 6 - 8 layers). The 6 - 8 soil layers covered the ranges of stream depths for most of the channel segments in this study, as the water levels during our simulation period fluctuated between 0.8 - 2.8 m below the ba nk elevatio n . The continuity equation can be written as (ignoring surface runoff): ( 2.40 ) Expanding the partial derivatives on the left - hand side of Eq. ( 2. 1) using the product rule, 32 ass uming smooth variation of A , Q and T s along the stream segment and combining Eq. ( 2.39 ) and ( 2.40 ), the following equation can be obtained ( Moore et al., [2005]; Kurylyk et al., [2016]; Gallice et al., [2016]): ( 2. 41) The temperature difference term in the right - hand side of Eq. ( 2.41 ) drops out for losing reaches, such that the advective groundwater heat flux will not directly cause the stream temperature change. 2.2.4. The streambed temperature modul e and hyporheic ex change A one - dimensional (vertical) advection conduction heat transport equation was used (Anderson, 2005; Caissie et al., 2014) to model the dynamics of streambed temperatures: ( 2.42 ) where z denotes the vertical coordin ate (positive downward), k b is the effective thermal conductivity (W m - 1 K - 1 ) of the saturated water - sediment matrix, T b is the streambed temperature (K), is the vertical Darcy flux (m s - 1 ) where a p ositive value represents down welling, c m (J kg - 1 K - 1 ) is the specific heat of soil water matrix and m (kg m - 3 ) is the density of solid - water matrix. In this study, the stream depths were relatively shallow and vertical water flux was the dominant compone nt, therefore v z was estimate d as , where is the groundwater flux entering or leaving the channel (computed in PAWS) , W is the wetted width of the channel. The volumetric heat capacity m c m for the solid - water matrix was estimated as (Caissie et 33 al., 2014): ( 2.43 ) where n is the porosity of the streambed material , and represent the density (kg m - 3 ) and specific heat (J kg - 1 K - 1 ) of the streambed solid material, respectively . and for the streambed material (mainly sandstone) were estimated as 2650 (kg m - 3 ) and 920 (J kg - 1 K - 1 ), respectively. Similar ly, k b was estimated using the relation: (2.44 ) where k w and k g ( W m - 1 K - 1 ) are the thermal conductivities of water and the aquifer material, respectively. We used a value of 0.25 for the porosity and 0.59 and 2.75 (W m - 1 K - 1 ) for k w and k g , respectively, the latter being a typical thermal con ductivity value for sandstone (Incropera and DeWitt, 1996), the main substrate in the research area. With the above properties, k b was calculated as 2.21 W m - 1 K - 1 . Extending to depths of 5 - 6 m below the surface, the temperature fluctuations in the st reambed are usually small relative to sur face soil or stream temperature fluctuations (Caissie & Luce, 2017) . Therefore, after examining the sensitivity of model results to the streambed depth we used a depth of 6 m and discretized this regio n into 12 layers. The streambed layers were finer near the stream bottom and became coarser while approaching the groundwater layers to account for the steeper temperature fluctuations near the stream. An implicit upwind scheme (Phanikumar & McGuire, 2010 ) was used to solve the streambed conduction equation with temperatures in the stream and the unconfined aquifer serving as the top and bottom boundary conditions respectively. 34 2.2.5. The soil temperature module PAWS us es the soil temperature module implemente d in CLM 4.0. Briefly, CLM uses the unsteady heat conduction equation for vertical soil heat transfer: ( 2.45 ) where c is the volumetric snow/soil heat capacity (J m - 3 K - 1 ), is the soil thermal conductivity (W m - 1 K - 1 ) which is heterogeneous based on the soil constitutes and saturation extent (Oleson et al., 2010), z is the depth in the vertical direction (m) and T is the soil temperature (K). We used 20 soil layers to discr etize the soil column and modified the no heat flux bottom boundary condition used in CLM to include heat flux from the groundwater domain by using groundwater temperature as the bottom boundary condition. The top boundary condition was a heat flux conditi on from the overlying atmosphere into the surface snow/soil layer. The heat flux into the snow/soil surface also included the net radiation and sensible and latent heat fluxes following descriptions in earlier sections. Th e calculations of net radiation h eat fluxes coming into the soil surface followed the same steps described in sections 2. 2.2 . The latent and sensible heat fluxes on the soil surface were calculated using an aerodynamic resistance method based on the Monin - Obukhov similarity theory ( Kundu et al. , 2015) in CLM. This equation was numerically solved using the Crank - Nicholson method and additional details are available in the CLM 4.0 documentation ( Oleson et al. , 2010). 2.2.6. The groundwater temperature module Groundw ater temperature is often approxi mated as the mean annual air temperature and used in stream temperature modeling (e.g., MacDonald et al. , 2014). Since PAWS uses 35 process - based descriptions of flow and transport, a general 2 - D advection - dispersion equation ( Domenico & Schwartz , 1998) was used to model groundwater temperature: ( 2.46 ) where T gw t (s) is time, (kg m - 3 ) and c b (J kg - 1 K - 1 ) are the density and spe cific heat of aquifer - water matrix , (kg m - 3 ) and c w (J kg - 1 K - 1 ) are the density and specific heat of the water, q (m s - 1 ) is the Darcy flux vector directly obtained from the PAWS groundwater flow module, Q soil is the soil hea t flux exchange between soil and groundwater and (W m - 3 ) includes the heat transfer due to advection and conduction through the streambed , k e (W m - 1 K - 1 ) is the effective thermal conductivity : ( 2.47 ) where k o (W m - 1 K - 1 ) is the bulk thermal conductivity, (m) is the thermal dispersivity term and (m s - 1 ) is the absolute value of the Darcy flux vector. Following earlier research that ind icated negligible impacts of thermal dispersion in modeling groundwater surface water interac tions ( Ingebritsen & Sanford , 1998; Hopmans et al. , 2002; Vandersteen et al ., 2015) we assumed = 0 in this study. After computing the temperature of the groundwater domain, heat exchange between the soil and groundwater domains w as estimated as: ( 2.48 ) where T soil - is the soil temperature of the second bottom layer and A ( m 2 ) is the grid cell area. For grid cells that contain stream segments, we also considered heat exchange between streambed and groundwater. 36 The calculation of can be expressed as (Brookfield et al., 2009; Therrien et al., 2010) : ( 2.49 ) where temperature if the groundwater domain is gaining water from the stream or the groundwater temperature if the groundwater domai n is losing water to the stream; (K) is the bottom layer streambed temperature. For the grid cells in the groundwater domain without any interaction with stream segments, = 0. Based on the vertical exte nt of the soil column and the first groundwater layer depth used in the present study, the groundwater temperature simulated here accounted for an average temperature corresponding to approximately 5 - 20 m aquifer beneath the ground surface elevation. T gw was then used as the bottom boundary condition for the stre ambed temperature simulations. The model was run for a period of approximately two years (May 2015 April 2017) using a uniform step size of 20 m in the horizontal x - and y - directions with typica l computer run times of 1.5 days on a workstation with Inte l i7 core processor. Results and Discussion The performance of the PAWS model was assessed by comparing model results with different types of observed data including point measurements at the fie ld sites as well as remotely sensed data for the whole catc hment. Multiple model performance metrics were used to assess model performance including the Nash - Sutcliffe efficiency coefficient (NASH) ( Nash & Sutcliffe , 1970) , the coefficient of determination (R 2 ), and the root mean squared error (RMSE). The order of presentation is as follows. T he ability of the model is first tes ted to reproduce observed water fluxes in the catchment by showing comparisons 37 between observed and simulated streamflows, transient groundwater heads and soil moisture. Comparisons are then show n vegetation p roce sses since vegetation growth and shading directly control stream temperatures while ET fluxes control both water and energy fluxes within the catchment. Therefore, comparisons of simulated LAI and ET with remotely - sensed data (MODIS) are presented next . Su bsequently, the performance of the temperature model is shown in different hydrologic domains including soil and groundwater temperatures, stream temperatures and streambed temperatures. 2.3.1. Stream flow , groundwater head and soil moisture comparisons A com pari son of simulated and observed streamflows at the outlet of the catchment shows that, the model successfully captured the diel fluctuations of the stream discharge (Figure 2.4 ). The NASH value for the comparison of hourly stream discharge is 0.71 and t he R MSE value is 0.01. The stream discharge was slightly underestimated during the period May 2015 to August 2015, whereas overestimation was observed during March 2016 to July 2016. Due to the small size of the stream channels, usually less than 3 m wide , in accurate representations of the channel geometry may have introduced uncertainties and overestimated peaks in the simulated stream discharge. Some unmeasured spatially - heterogeneous parameters, such as streambed hydraulic conductivity ( K r , m day - 1 ) and vegetation interception depths may have introduced uncertainties as well. Considering these uncertainties and the relatively small values of stream discharge, the simulated results provided an acceptable ( Krause et al ., 2005; Legates et al ., 1999) descrip t ion of the observed data. 38 Figure 2.5 shows the simulated and observed groundwater head dynamics at four boreholes whose locations are shown in Figure 1. The R 2 values for the four groundwater sites are 0.87, 0.78, 0.64 and 0.32, respectively. Overall, the model simulations were able to reproduce the amplitude and the general trend of the observed groundwater heads. At BH12, although the simulated groundwater heads were generally in the range of observed values, the decreasing trend during the simulation pe r iod was not adequately captured. This is attributed to various factors including limited borehole data and uncertainties associated with boundary conditions for the groundwater flow model. The simulated soil moisture results (Figure 2.6 ) generally capture d the trend and fluctuations measured by the soil moisture sensor. One obvious mismatch occurred during a dry period in July 2016. The observed soil moisture approached the value 0.07, while the lowest simulated soil moisture value remained at 0.14. This i s likely the result of a scale mismatch between simulations (an area of a grid cell is 400 m 2 while the sampling volume of a moisture sensor is of the order of a few cubic centimeters). Considering the uncertainties and scaling issues associated with model parameterization, we conclude that the simulated water fluxes can be used to serve as the basis for our temperature simulations in the catchment. 39 Figure 2 . 4 Simulated and observed streamflow compar ison at basin outlet 40 Figure 2 . 5 Simulated and observed transient groundwater head comparison at four different boreholes (BH), asl = above sea level Figure 2 . 6 Comparison of simulated and observed soil moisture 41 2.3.2. LAI and ET comparisons Figure 2.7 shows spatial comparisons of 8 - day - average LAI at two selected 8 - day periods: September 07 September 14, 2015 and July 05 July 12, 2016 (individual pixels fr om the original MODIS data can be seen in this figure due to the relatively small size of the catchment). For both periods, the two pixels close to the southern - most portion near the catchment outlet and th e pixel in the northwest area show relatively hig h LAI for both simulated results and the MOD15A2H product, and these areas correspond to land use types with heavy portions of deciduous trees. However, during September 07 September 14, 2015, several pixe ls with predominantly agricultural land use porti ons are slightly overestimated by simulated results relative to MOD15A2H data. This may be due to discrepancies between phenological parameters utilized in C3 crop of CLM and uncertainties associated with MO DIS product. Figure 2. 8 shows a time series compa rison of 8 - day catchment - averaged LAI with simulation results for the same period. The general trends of simulated LAI and observed MODIS data are in good agreement: higher LAI values in growing season and l ower LAI in other seasons. The simulated catchmen t average LAI values remain around 2 - 2.6 during growing seasons, while observed catchment averaged LAI (MOD15A2H) values represent a slightly wider range (1.6 - 2.7) during the growing seasons. It should also be noted that some parts of the MODIS pixels are outside the simulation domain, and this mismatch in domains may have introduced uncertainties into the comparison. Spatial ET distributions over two 8 - day periods ( September 07 September 14, 2015 and July 05 July 12, 2016 ) based on MODIS data and simulated results are in good agreement, see Figure 2. 9. The ET values are relatively high in the northwest region and the outlet 42 areas of the catchment for both MODIS data and simulated results during the two 8 - day periods. This is primarily because of t he heavy portions of deciduous trees in these areas. As discussed in previous work, PAWS usually represents a better spatial heterogeneity of ET due to its more detailed sub - cell land use information than re motely sensed MODIS product ( Shen et al., 2013; N iu et al., 2014). To quantitatively evaluate the ET values, we compare the catchment - averaged 8 - day ET time series between MODIS data and simulated results , see Figure 2. 10 . Due to cloud contamination and other possible data quality issues, MODIS observati ons are no t continuously available during the simulation period. Annual ET cycles are matched well between MODIS data and simulated results with an R 2 value of 0.89. However, the winter and early spring ET values were overestimated by the simulated results . Probable explanation for this mismatch could include a combination of several factors: a) different ET algorithms used between MODIS (Penman - Monteith) and PAWS (resistance approach based on the two - big leaf model ( Dai et al. , 2004) with explicit calculat ion of lea f temperatures and aerodynamic and canopy resistances); b) inexact parameterization of vegetation phenology parameters for the land use categories, for example, we used the broadleaf deciduous tree category to represent the deciduous woodland as described before; c) well - known uncertainties in the MODIS ET product ( Mu et al. , 2011). Overall, the ET magnitude and spatial variations in ET agree well between simulated results and MODIS observations. 43 Figure 2 . 7 Spatial maps of 8 - day LAI: (a) MODIS data for Sep.07 - Sep.14, 2015 (b) Simulated results for September 07 - 14, 2015 (c) MODIS data for July 05 - 12, 2016 (d) Simulated results for July 05 - 12, 2016 Figure 2 . 8 Time series comparisons between MODIS data (MOD15 product) versus simulation results for 8 - day catchment averaged leaf area index 44 Figure 2 . 9 Spatial maps of 8 - day ET: (a) MODIS data for Sep.07 - Sep.14, 2015 (b) Simulated results for September 07 - 14, 2015 (c) MODIS data for July 05 - 12, 2016 (d) Simulated results for July 05 - 12, 2016 Figure 2 . 10 Time series comparisons between M ODIS data (MOD16 product) versus simulation results for 8 - day catchment averaged ET 45 2.3.3. The temperature results in different hydrologic domains The comparison between observed and simulated soil temperature is shown in Figure 2. 11 for a de pth of 12 cm. Soil te mperature and moisture were measured at the same location (Figure 1). Downward heat fluxes into the soil surface were calculated under conditions of typical deciduous forest cover, the dominant land use in this part of the study area. An R 2 value of 0.84 and an RMSE value of 1.15 indicate that a good agreement between simulated and observed temperature time series was obtained (Figure 2. 11 ). The upper envelope of diurnal fluctuations was well captured by the model, while the lower bounds were overestimated. Possible reasons for thi s discrepancy include uncertainties associated with parameterization of soil thermal conductivities and /or heat capacities and instrument errors as well as scale mismatch. Simulated groundwater temperatu res at four sites (Figure 2.12 , with locations of bo reholes marked in Figure 2.1 ) generally captured the observed trend and amplitudes of measurements. The blue bands are from the results of a sensitivity analysis and correspond to an uncertainty of ±30% i n the hydraulic conductivity of the first groundwate r layer which represents the unconfined aquifer. Additional sensitivity analysis results are presented in the Supporting Information. The R 2 between daily simulated groundwater temperature results and obs ervations for the four BHs are 0.78, 0.83, 0.74 and 0.59; RMSE values are 0.23, 0.26, 0.28 and 0.35, respectively. It should be noted that the simulated groundwater temperature results are from a single groundwater layer with an average aquifer depth betwe en 5 - 20 meters, while temperature sensors were locat ed at depths of 10.18 m, 8.45 m, 6.30 m and 8.11 m for the four BHs respectively. The differences in depths may have caused the shift of temperature phases as noted by other researchers ( Vandersteen et al . , 46 2015 ). This could be the reason for the slight d isagreement in the timing of groundwater temperature minima between simulated and observed data for the four BHs. Although much smaller compared with stream or air temperatures, groundwater temperatures s till showed monthly variation. Our results indicate that specifying the bottom boundary condition based on simulated groundwater temperatures improves overall model performance compared to using the annual mean air temperature or the no flux boundary condi tion. For comparison, the constant value of annual m ean air temperature is also plotted using dashed lines in Figure 2.12. Figure 2 . 11 Comparison of observed and simulated soil temperatures 47 Figure 2 . 12 Comparison of observed and simulated groundwater temperatures at four borehole locations with mean annual air temperature shown. The blue bands correspond to the uncertainty associated with ±30% changes in the hydraulic conductivity of the first groundwater layer. RMSE = root - mean - square error Comparisons of observed and simulated stream temperature results in the form of time series (Figure 2.13 a) and a 1:1 plot (Figure 2.13 b) show that diurnal fluctuatio ns were adequately reproduced by the model. The temperature residuals between simulated a nd observed temperature are also presented in Figure S5 (Supporting Information). The R 2 value for the hourly comparison at the basin outlet is 0.87, and the RMSE valu e is 1.32. However, the model underestimated the peak stream temperatures in April by app roximately 1.5 - 2 o C compared with observations. Several factors (or a combination of factors) could have contributed to this mismatch: a) underestimated incoming sol ar radiations during April; b) overestimated vegetation shading effect during April due t o inaccurate phenological parameters. In the absence of direct measurements of solar radiation reaching the stream surface and runoff contributions, it is difficult to pinpoint the exact reason for the April peak temperature mismatch. It is relatively diff icult to simulate 48 the temperature dynamics in a small stream due to shallow stream depths and the relatively small heat capacity. MacDonald et al. (2014 ) indicate that stream temperature simulations were sensitive to the stream wetted depth information used as input to the model and that simulated temperatures were more sensitive to decreased rather than increased wetted width. Inaccurate descriptions of channel geometry may have introduced errors int o the simulated wetted width and depth, which probably introduced uncertainties into the estimation of net heat flux. The deciduous tree category used to represent rip arian tropic CLM 4.0 which is the closest approximation we could find for the riparian vegetation in the Wood Brook catchment which is mostly composed of English Oak and oth er trees (see methods section). It is possible that some of the mismatch between observed and simulated stream temperatures is due to mismatches in phenology of these tree species from the generic deciduous tree in CLM. In addition, it is note d that sever al complexities at the field sites are not incorporated into the current version of the m odel. These include features such as wood debris dams that produced ponding, bedrock outcrops that created drop - offs producing small local waterfalls within the channe l, as well as fallen trees that provided permanent shading in some reaches (Figure 2.14 ). Seasonal heat flux budgets for the year 2016 were calculated at the catchment outlet stream segment (the reach marked CD in Figure 2.1 ) for different components, as s hown in Table 2.3 . We found that net radiation was the dominant heat source and latent h eat flux was the primary heat sink, while the annual net sensible heat flux was close to zero. Runoff contribution was found to be negligible during the simulation per iod and hence not shown in Table 2.3 . Although relatively small compared to the net radia tion heat fluxes, streambed 49 c onduction heat fluxes play a role in damping the diurnal and seasonal temperature amplitudes driven by atmospheric effects ( MacDonald et al., 2014; Cox & Bolte , 2007 ) . The streambed heat conduction fluxes represented heat sources fo r sustaining warm stream temperatures during winter and served as heat sinks for cooling the streams during summer (Table 2.3 ; also Hannah et al ., 2004). Friction heat fluxes, which represe nt the highest potential friction flux, were the smallest portion (< 1 W m - 2 ) of the budget due to the small stream size and low flow rates. Similarly, groundwater advective heat fluxes were of small magnitude year round, and this finding is consistent wit h the results of ( Caissie & Luce , 2017) who noted that lateral groundwater advective heat fluxes are generated only from rapid flow events when water entering the stream channel did not have adequate time to reach equilibrium temperature before mixing . 50 Figure 2 . 13 Stream temperature comparisons (observed versus simulated) at the basin outlet (a) Time series comparison. (b) 1:1 plot. RMSE = root - mean - square 51 Figure 2 . 14 Photographs of the main channel showing some complex features not included in the current modeling (a) a wood debris dam (b) a bedrock outcrop within the channel and (c) fallen trees that provide permanent shading in some stream reaches. Photo courtesy: Dr. Mantha S. Phanikumar. 52 Table 2 . 3 Seasonal heat flux budgets for year 2016 at the catchment outlet stream segment (segment CD in Figure 2.1) Unit (W m - 2 ) Spring Summer Fall Winter Net Radiation 19.00 35.84 7.32 3.32 Latent heat flux - 17.46 - 25.85 - 14.02 - 9.21 Sensible heat flux - 0.79 4.51 - 0.74 - 3.61 Friction heat flux 0.96 0.43 0.77 0.78 Streambed conduction 2.65 - 9.10 - 2.97 6.74 Groundwater advective heat flux 1.37 - 2.23 - 1.03 1.44 Figu res 2.15 and 2.16 show comparisons of observed and simulated streambed temperatures at sites 1 and 4 (locations marked in Figure 1) for different depths (5cm, 10cm, 20cm, 30cm); a close - up view of the regions marked A, B and C in Figure 11 is shown in Figu re S7 (Supporting Information). The R 2 and RMSE values, shown in Table 2 marked The thermal conductivities and heat capacities at the four measurement sites were slightly adjusted to acknowledge spatial dependence of these parameters (Table 2.5 ). The comparison is shown for 50 days for which observed data are available . Caissie et al. (2014) found that the diel variations of streambed temperature were no longer vis ible at depths greater than 70 cm for the two streams they sampled in New Brunswick, Canada. They estimated the vertical (upwelling) Darcy flux with values of 2.5 and 5.1 mm/ hour for the two streams, respectively. Among the four sampling sites, at sites 1 and 2 which represented gaining stream reaches, the diurnal fluctuations and amplitude at 20 cm and 30 cm were strongly attenuated due to the vertical movement of groundwater . In contrast, at 53 sites 3 and 4, the diurnal fluctuations were not visibly dampene d even at a depth of 30 cm, which was evidence of losing reaches. The average v z values from the PAWS simulations during the 50 - day period for the four sites were - 3.63, - 2.1 3, 1.63 and 0.46 mm/hour respectively. Positive v z values indicate losing stream r eaches while gaining stream reaches are characterized by negative values of velocity. The simulated v z values are consistent with the above observation that sites 1 and 2 are in gaining reaches while sites 3 and 4 represent losing reaches. There is field e vidence of a clay layer of unidentified spatial extent near sites 3 and 4 and the presence of this layer was expected to inhibit groundwater surface water interactions. How ever, detailed subsurface characterization is needed to understand the spatial ext ent of the clay layer and its thickness to rule out interactions between domains and such data are not available at this time. It is possible that the clay layer is not conti nuous, or the layer may have little to zero thickness in places allowing water to move to the overlying region as shown by the negative v z v alue at site 1. Our results demonstrate the ability of the one - dimensional heat transfer model to simulate streamb ed temperature fluctuations accurately at various depths. However, the simulated s treambed temperatures at 20 cm and 30 cm at site 1 slightly overestimated the diurnal fluctuations relative to the observed data, which may be due to underestimated v z values and/or inaccurate values of streambed thermal properties. 54 Figure 2 . 15 Comparison of observed and simulated streambed temperatures at site 1 for different depths: (a) 5 cm, (b) 10 cm, (c) 20 cm, (d) 30 cm, and (e) Close - up views of the areas marked A, B and C in Figure 2.15 . 55 F igure 2 . 16 Comparison of observed and simulated streambed temperatures at site 4 for different depths: (a) 5 cm, (b) 10 cm, (c) 20 cm, and (d) 30 cm. 56 Table 2 . 4 Performance of the streambed temperature sub - model at four sampling sites for different depths ( 5 cm, 10 cm, 20 cm, and 30 cm ) 5 cm 10 cm 20 cm 30 cm Locations Simulation R 2 RMSE R 2 RMSE R 2 RMSE R 2 RMSE Site1 Baseline 0.75 0.84 0.79 0.53 0.86 0.24 0.87 0.24 Constant GW T 0.73 0.91 0.75 0.58 0.81 0.32 0.76 0.34 Site2 Baseline 0.82 0.85 0.83 0.54 0.84 0.26 0.80 0.23 Constant GW T 0.79 0.88 0.81 0.59 0.79 0.23 0.73 0.34 Site3 Baseline 0.74 0.92 0.77 0.58 0.78 0.46 0.81 0.23 Constant GW T 0.74 0.93 0.76 0.58 0.75 0.50 0.78 0.27 Site4 Baseline 0.73 1.04 0.74 0.80 0.81 0.54 0.82 0.42 Constant GW T 0.72 1.04 0.75 0.80 0.78 0.57 0.79 0.46 Table 2 . 5 Thermal conductivities and leakance values paramet erized for the four streambed temperature sampling sites. Sites Streambed conductivity, K r (m day - 1 ) v z (mm hour - 1 ) Streambed thermal conductivity, k b (W m - 1 o C - 1 ) Site 1 0.105 - 3.63 2.55 Site 2 0.098 - 2.13 2.21 Site 3 0.131 1.63 2.21 Site 4 0.112 0. 46 2.45 57 We examine the effects of using mean annual air temperature, a constant value during the simulation period, as a proxy of groundwater temperature in the simulation relative to the baseline (i.e., explicit simulation of groundwater temperature). Stream temp erature results at the basin outlet are presented as the temperature residuals (simulated minus observed) for the two cases (Figure 2.17 ). The model performance at the basin outlet using these two methods are not significantly different. One maj or reason is due to the dominant control of stream surface heat fluxes (radiation heat fluxes plus sensible and latent heat fluxes). In addition, the basin outlet is located within a losing portion of the stream whi ch the effect of groundwater is expected to be even smaller. However, for streambed temperature results, the performance differences are expected to be larger especially at deeper depths due to the proximity to groundwater. Comparisons of streambed tempe ratures at two representative sites si te 1 (a gaining reach) and site 4 (a losing reach) are shown in Figures 2.18 and 2.19 for the two cases (baseline simulation and constant groundwater temperature simulation). 58 Figure 2 . 17 Strea m temperature residuals (difference between simulated and observed stream temperatures at the basin outlet) for two simulations: baseline and constant ground - water temperature. 59 Figure 2 . 18 Comparison o f observed and simulated streambed temperatures at Site 1 for different depths (a) 5 cm (b) 10 cm (c) 20 cm (d) 30 cm. The red line represents observed data and the blue line the baseline sim ulation that explicitly simulated groundwater temperature. The gr een line represents a simulation that used a constant groundwater temperature based on the annual mean air temperature. 60 Figure 2 . 19 Compar ison of observed and simulated streambed temperatures at Site 4 f or different depths (a) 5 cm (b) 10 cm (c) 20 cm (d) 30 cm. The red line represents observed data and the blue line the baseline simulation that explicitly simulated groundwater temperature. The green line represents a simulation that used a constant groun dwater temperature based on the annual mean air temperature. 61 The R 2 and RMSE values for the streambed temperature comparisons are tabulated in Table 2 .4 mean annual air t emperature as groundwater temperature. The baseline simulation gener ally produced a better match with observed streambed temperatures. The differences are more obvious for gaining reaches (Site 1 and Site 2) and at greater depths (depths of 20 cm and 30 c m). Explicitly simulated groundwater and streambed temperatures provi de a more realistic method than using the mean annual air temperature for calculating streambed conduction fluxes. However, the two - point gradient method used to approximate the streambe d conduction heat flux may introduce some errors especially when the vertical water flux is larger than 5 mm h - 1 ( Caissie & Luce , 2017) . A sensitivity analysis was conducted to understand how simulated stream and streambed temperatures respond to changes in model paramet ers. Although a detailed examination of sensitivity and equifinality of parameters is beyond the scope of the present work, we examined the effect of changing five model parameters on sim ulated temperatures. The parameters include streambed hydraulic condu ctivity ( ), the wind speed sheltering factor ( ), streambed thermal conductivity ( ), the characteristic dist the two - point gradient method and the hydraulic conductivity of the first groundwater layer which represents the unconfined aquifer ( ) , as tabulated in Table 2.6 . Parameters (with the e by relative to the baseline simulation while keeping all other parameters as the calibrated values . For the Stream temperatur e results (Figure 2.20 , Tab le 2.7 ) indicate that of the five parameters examined, simulated stream temperatures are most sensitive to the streambed hydraulic conductivity 62 and that this parameter affects simulated temperatures at the hourly time scale. The other parameters in decrea sing order of importance after are: C w k b , and finally K 1 (Figure 2.20 ). Sensitivity analysis results for streambed temperatures followed a trend that is similar to stream temperatures although the magnitudes of changes are different (Figure 2.21 , Table 2.7 ). The diffe rences in temperatur es at streambed site s 1 and 4 in ter ms of average hourly temperature percent change for the selected parameters are similar . The simulated hourly temperature results are most sensitive to K r (streambed hydraulic conductivity) among the selected parameters. Perturbations of K r influence excha nge of both water and heat fluxes between str eam and the groundwater domain thus affect ing the stream discharge via v z , and the streambed heat conduction simultaneously. The hourly temperature results are moderately sens itive to C w , k b . The parameter C w primarily influences the latent heat flux while k b and directly affect the approximation of stream bed heat conduction flux. The temperature results have the smallest sensitivity in response to perturbations of gro undwater hydraulic conducti vities. Table 2 . 6 Model parameters for temperature modules S ymbol Parameter Meaning k b ( W m - 1 K - 1 ) Streambed effective thermal conductivity (m) C haracteristic depth to calcu late the streambed conduction heat flux C w W ind speed sheltering factor K r ( m day - 1 ) Streambed hydraulic conductivity K 1 (m day - 1 ) Hydraulic conductivity of the first groundwater layer 63 Table 2 . 7 Average percentage changes ( o C x 100/ o C) of hourly stream/streambed temperatures in response to parameter perturbations average percent change of hourly temperature (%) parameter change by ±10% C w k b K r K 1 S tream 1.478 3.152 4.770 3.345 0.136 Streambed : site 1 5cm 1.276 3.155 4.092 3.358 0.130 Streambed: site 1 10cm 1.191 3.039 3.653 3.230 0.128 Streambed: site 1 20cm 1.129 2.890 3.278 3.065 0.125 Streambed: site 1 30cm 1.080 2.767 3.052 2. 920 0.125 Streambed: site 4 5cm 1.381 3.447 4.558 3.720 0.127 Streambed: site 4 10cm 1.320 3.256 4.021 3.539 0.126 Streambed: site 4 20cm 1.271 3.044 3.542 3.334 0.130 Streambed: site 4 30cm 1.219 2.879 3.252 3.151 0.136 parameter change by ±30% C w k b K r (=0.25m) K 1 S tream 5.293 3.189 7.388 3.634 0.198 Streambed: site 1 5cm 4.719 3.374 6.217 3.627 0.164 Streambed: site 1 10cm 4.424 3.441 5.545 3.476 0.154 Streambed: site 1 20cm 4.158 3.457 5.171 3.291 0.141 Streambed: site 1 30cm 3.957 3.435 4.999 3 .122 0.133 Streambed: site 4 5cm 5.195 3.489 6.217 3.943 0.152 Streambed: site 4 10cm 4.804 3.323 5.545 3.774 0.146 Streambed: site 4 20cm 4.448 3.135 5.171 3.572 0.167 Streambed: site 4 30cm 4.187 2.999 4.999 3.385 0.199 parameter change by ±50% C w k b K r K 1 S tream 7.886 4.196 16.468 5.633 0.510 Streambed: site 1 5cm 6.851 5.151 14.203 5.556 0.423 Streambed: site 1 10cm 6.470 5.906 12.879 5.371 0.390 Streambed: site 1 20cm 6.214 6.427 12.105 5.170 0.369 Streambed: site 1 30cm 6.011 6.637 11. 662 4.947 0.354 Streambed: site 4 5cm 7.222 4.712 15.148 5.719 0.407 Streambed: site 4 10cm 6.838 4.593 13.700 5.662 0.359 Streambed: site 4 20cm 6.559 4.422 12.834 5.571 0.378 Streambed: site 4 30cm 6.324 4.337 12.353 5.412 0.429 64 Figure 2 . 20 Sensitivity of simulated stream te m perature to changes in parameters listed in Table 2.6. Each parameter was changed by ±10%, ±30%, ±50% and changes in hourlystream temperature at the basin outlet relative to th e baseline simulation results. 65 Figure 2 . 21 Sensitivity of simulated stream te m perature to changes in parameters listed in Table 2.6 . Each parameter was changed by ±10%, ±30%, ±50% and changes in hourlystr eam temperature at the basin outlet relative to the baseline simulation results. 66 Figure 2 . 22 Daily - average, cross - sectional temperature profiles of stream, streambed and groundwater along a portion of the main stream (between points marked A and D in Figure 2.1) for different seasons. Panel (a) is a schematic of the general land use features along the stream segment. The numbers 1, 2, 3, 4 in (a) deno te the streambed temperature sampling locations. Each pa nel in (b) through (e) shows the temperature profile for one representative day in each of the four seasons : (b) Spring, (c) Summer, (d) Autumn, and (e) Winter . Head difference (groundwater head minus river stage) variations as a function of distance from point A along the main stream are plotted using the black solid lines using the second Y - axis on the right. Gaining reaches of the stream correspond to positive head differences while losing portions are associated with negative head differences. 67 As show n in Figure 2.22 , the cross - sectional daily - average temperature profiles of stream, streambed and groundwater along a portion of the main stream (between points A and D in Figure 2. 1) are presented fo r one day in each of the four seasons. The stream temper atures are shown in the topmost layer while the bottom layers represent the streambed temperatures with the lowermost layer representing groundwater temperatures. The daily average head differences b etween groundwater heads and stream stages are superimpo sed on the figures to illustrate the gaining or losing portions of the stream. Stream bank elevations and land use scenarios along the main streams are also shown in Figure 13 (a). Gaining reaches of the stream correspond to positive head differences (seco nd Y - axis on the right) while losing portions of the stream are associated with negative head differences. It could be seen from Figures 13 (c) and 13 (d) that on July 15, 2016 and October 15, 2016, t he daily average stream temperatures with crop land use nearby are slightly higher than portions surrounded by deciduous trees. However similar phenomenon is not obvious on April 15, 2016 and January 15, 2017. This could be attributed to the seasonal shift of phenological properties of the two vegetation classe s, such that the differences of shading effects between deciduous trees and crops have a greater effect in summer and fall than in spring and winter. Overall, the losing and gaining portions of the st ream are not altered for the four selected days for the year considered in this analysis and appear to be stable features over the simulation period. This aspect will be examined in detail in future work. There is a clear shift of the temperature divide ( shown in yellow color in the figures) within the profile s as the stream reach switches between gaining and losing reaches for all four days. Near the catchment outlet where sites 3 and site 4 are located, the temperature divide approaches the groundwater t emperature, 68 indicating a losing stream, while at the loc ations of site 1 and site 2, the temperature divide approaches the stream temperature, indicating a gaining stream. These phenomena are consistent with previous discussions. The analysis of streambed temperature stratifications and time series have become a standard approach for quantifying the vertical exchange fluxes between groundwater - surface water domains ( Caissie et al ., 2007; Vandersteen et al ., 2015; Krause et al ., 2011 ). Simulated streambe d temperature results from process - based models could be used to gain further insights into key processes and parameters (for example, to optimize reach scale streambed hydraulic conductivity values). Therefore, integrated process - based models are useful t ools to gain insights into multiple processes that affec t thermal dynamics in catchments. Conclusions This chapter provided a catchment - scale framework to simulate coupled hydrologic and stream subsurface thermal transport processes. When simulating soil and streambed temperatures, significant computational expense could be saved by reducing the three - dimensional he at transport equations into a combination of two - (groundwater) and one - dimensional (streambed) equations with limited loss of physics in the catchment examined in our work. Process - based simulati ons make it possible to manipulate individual processes, ev aluate the impacts of different heat fluxes, and to understand the sources of errors. As with other process - based hydrologic models, PAWS requires a large amount of hydrological and geophysical dat a. The weighted - average method of flux computation for gri d cells with mixed land use proved to be an efficient way to estimate the (below canopy) net radiation fluxes reaching the streams while reducing the necessary parameters and observations. This ap proach can be used to evaluate the potential impacts of cli mate and 69 land use changes on stream temperature dynamics and to identify best strategies from the point of protecting steam habitat ( Mohseni et al., 2003; Jackson et al ., 2018). Future work will te st the ability of this method to simulate stream temperatur e dynamics in large river basins. It is possible that the estimated radiation fluxes may need to be recomputed at smaller scales and aggregated to the coarse grid scale. Future work will also evalu ate model sensitivity to parameters and vegetation scenario s as well as snow / ice dynamics during winter months. Scaling effects of heat transport processes in different landscape units will also be investigated by applying the model to large watersheds with more complex stream networks and land use in different parts of the world. 70 Chapter 3. Modeling the effect s of vegetation on stream temperature dynamics in a large , mixed land cover watershed in the Great Lakes region This study aim s at quantifying th e effects of resolving the spatial variability of vegetation on the temperature simulations using a process based hydrologic model, i.e. PAWS+CLM, in a 5200 km 2 mixed - landuse watershed in Michigan, USA , i.e. the Kalamazoo River watershed ( KRW , Figure 3.1 ). The model explicitly solves the equations that govern the stream discharge and temperature processes. By incorporating interactions of thermal and water fluxes at the interfaces of air/surface and surface/subsurface, the model enable s comprehensive repres entation of the stream thermal dynamics. The streamside landuse effects on the stream temperature are explicitely represented by coupling with the CLM phenological module. The streamside land use information is rescaled using di fferent nested resolutions t o represent the land use effects on stream temperature; we find that a resolution of 90 m land use best represent the land use effects while simulating stream temperature . The simulation cover s a 7 - year period from 2003 to 2009 , during which period the hyd rologic and temperature results are compared with multiple observations. Stream thermal budgets and r esponses of stream temperature to vegetation scenarios including potential deforestation effects are reported. Introduction St ream temperature is a key ec ological variable (Caissie, 2006 ; Allan and Castillo, 2007 ) , mediating aquatic meta bolism and nutrient retention and transformation (Schaefer and Alber, 2007; Starry et al., 2005) , as well as influenc ing the abundance and distribution of aquatic species in cluding fishes (Ebersole et al., 2001; Wehrly et al., 2003 ) . Understanding 71 controls on stream temperature is necessary to predict how streams will respond to the warming climate (Schindler, 2001; Leibow itz et al., 2014). While primarily driven by climatic conditions, stream thermal regimes are also shaped by multiple factor s including near - stream (riparian) land cover, stream - groundwater interactions, discharge, and channel geomorphology (Caissie, 2006) . Land use development and deforestation have been fou nd to affect stream temperature regimes, generally increasing diel variation (Johnson and Jones, 2000) , and maximum summer temperatures are of particular interest when they could negatively impact coldwa ter fish habitat (Buisson et al., 2008; Hrachowitz et al., 2010; Kelleher et al., 2012). Moreover, streams vary in their thermal responses as a function of depth, which determines their heat capacities (Sun et al., 2015), as well as the relative influence of groundwater inputs and hyporheic exchange (Gordon et al., 2012; Caissie and Luce, 2017 ; Qiu et al., 2019 ) . M odeling approaches have been widely used to predict stream temperature responses to climate and land - use change (Jackson et al., 2017; Gallice e t al., 2016 ; Sun et al., 2015; MacDonald et al., 2014 ). Generally, s tream temperature models can be categorized in to two groups: statistical models ( Jackson et al., 2017; Benyahya et al., 2007 ) and process - based hydrologic models (PBHMs) ( Qiu et al., 2019; Gallice et al., 2016; Sun et al., 2015; MacDonald et al., 2014 ). Statistical models have the advantages of simplicity and efficiency due to the relatively small data demands and computational expense; however, they lack the ability to offer insights into fundamental and the often coevolving roles of multipl e processes, or to extrapolate beyond the range of conditions under which they are calibrated . PBHMs explicitly simulate the interactiv e physical processes that drive stream temperature dynamics by solvi ng t he governing equations of mass, momentum and ener gy 72 conservation . Time - marching simulations enable PBHMs to evaluate temporal variations, track individual processes , and understand the controlling factors (Beven, 2002). Various PBHMs have been develop ed to study stream temperature dynamics for different research purposes. For example, MacDonald et al. (2014) developed a watershed - scale river temperature model for studyi ng the temperature dynamics in mountain ous region s by solving the one - dimensional heat advection equation; Sun et al. (2015) built a rive r reach - scale , process - oriented stream temperature model by focusing on the riparian land - use effects . Recently, those models have inco rporated more detailed process r ep resentations to incorporate groundwater lateral advection and streambed conduction heat fluxes. For example, Gallice et al. (2016) assumed proportional relationship s between stream bed conduction heat flux and temperature difference between the stream and t he streambed at a certain depth. Haag and Luce (2008) estimate d the streambed conductio n flux by using the concept of two lumped parameters. However, PBHMs are usually data - demanding t o construct models that ensure process fidelity (Beven, 2002). Particula rly at the watershed scale, spatial heterogeneit y within the watershed could bring unce rtainties in estimating the stream discharge as well as water temperature (Beven, 1993 ; Beven, 2001; Beven, 2002). T o date, some PBHMs either simplified the simulation o f stream discharge by using empirical equations, e.g. , MacDonald et al. (2014), or coul d only be applied to small watersheds due to the fully three - dimensional nature of the models and the consequent high computational expense and data requirements. Ripar ian vegetation plays an important role in controlling stream temperature via shading (R oth et al., 2010; Moore et al., 2005) and by changing the riparian microclimate. Land use / land cover maps (represented using plant functional types or PFTs) generally contain a 73 much finer level of detail than one can include within model grid cells. To m ake computations tractable, these maps are often reclassified to include only a small number of dominant land cover types within a model grid cell, replacing land cover types that occupy only a small area within the grid cell with the dominant land cover i n that cell. Since the outcome of this reclassification depends on the size of the grid cell, the impacts of representing land cover in a hydrologic model depend on the grid resolution, and stream temperature results are particularly sensitive to these det ails. In large watersheds, grid cell sizes tend to be relatively large as well to make computations practical. One way to provide an adequate representation of riparian vegetation in large watersheds is to reclassify riparian land cover at a smaller nested resolution relative to the original grid size to upscale the vegetation effects to the size of the larger grid cell (i.e., by using several nested grid cells within a l arger cell; for example, ten 100 x 100 m nested cells can be created within a single 1 x 1 km grid cell). However, a systematic assessment of such a reclassification method and the effects of the nested grid cell resolutions on simulated stream temperature s is not reported in the literature. In this chapter, we apply the PAWS model ( P rocess - based A daptive W atershed S imulator, Shen and Phanikumar, 2010; Shen et al., 2013; Shen et al., 2014; Qiu et al., 2019a) , a PBHM of intermediate complexity, to predict streamwater temperatures over space and time across a humid temperate watershed with seas onally present, deciduous vegetation canopies in the riparian zones . The newly - developed stream temperature model (Qiu et al., 2019) uses the radiation transfer modules in CLM 4.0 (Dickinson, 1983; Sellers, 1985; Oleson et al., 2010) to partition radiation fluxes among riparian vegetation canopies, back radiation and fluxes received by stream surfaces . The dynamic seasonalit y of vegetation phenol ogical properties such as leaf area index (LAI) 74 and stem area index (SAI ) and their effects on radiation fluxes w ere simulated for different plant functional types in the vegetation module of CLM . Our study watershed is considerably l arger than a small experimental watershed where the model has been previously validated, and thus our study developed the model to run on data with a coarser spatial resolution. Research questions that will be addressed in this chapter include: 1. How can the riparian vegetation shading effect be approximated using relatively coarse grid resolutions ( e.g. 1 km 2 ) used to model large waters heds? 2. Is there an optimal grid resolution for representing the vegetation effects in large watersheds? 3. What are the dominant controls on stream temperature regimes over the seasons? and 4. What is the effect of riparian deforestation on temperature i n large vs. small streams? Materials and Methods 3.2.1. Model description PAWS is a process - based, distributed hydrological m odel that uses the finite volume method to solve the governing partial differential equations for different hydrologic domains based on structured grids (Shen and Phanikumar, 2010; Shen et al., 2013; Shen et al., 2014). The PAWS model structure allows effi cient coupling of surface and subsurface processes (Shen and Phanikumar, 2010) by simplifying the fully three - dimensional (3 - D), variab ly - saturated subsurface model using a combination of the one - equation (for the vertical soil column ) and quasi - 3 - D saturated groundwater equation for the aquifers, assuming that lateral soil water exchanges with the stream channels ar e negligible (Shen et al., 2013). The coupling of PAWS (Shen and Phanikumar, 2010) with the land surface model CLM (Community Land Model, version 4; Oleson et al., 2010) allows for a comprehensive representation of hydrologic and vegetation processes 75 inclu ding land surface processes, subsurface processes, and interactions among different hydr ologic domains (Shen et al., 2013) . Applications and validation of the PAWS+CLM model have been widely reported for various watershed s around the world (Shen and Phanik umar, 2010; Shen et al., 2013; Shen et al., 2014; Niu et al., 2015; Safaie et al., 2017; Niu et al., 2017) , but the model has only been applied to simulate stream temperatures by Qiu et al. (2019a ). Details of the coupled PAWS+CLM stream temperature model are described in Qiu et al. (2019a ) and therefore only a brief summary is presented her e. S imulation domains are discretized into different lateral grid cells and vertical layers. T he t opmost layer represents the land surface layer or overland flow layer in which surface runoff is computed based on the diffusive wave equation. Surface energy balances, evaporation, infiltration and snow processes are computed in the overland layer. Beneath the land surface layer are the vadose zone layers which are governed b y the Richards equation. Two groundwater layers (unconfined and confined) are connect ed to the bottom most layer of the vadose zone and are governed by the quasi - 3D groundwater segments crisscross the overl and flow layer and are connected to the first groundwater layer which represents the unc onfined aquifer . St ream discharge is governed by the one - dimensional diffusive wave equation exchanging fluxes with the overland flow and groundwater layer s . Several ea rlier studies demonstrated that one - dimensional models adequately describe solute and thermal transport in rivers ( Gallice et al., 2016; MacDonald et al., 2014; Anderson and Phanikumar, 2011; Phanikumar et al., 2007; Shen and Phanikumar, 2009). Therefore, the stream temperature module solves the one - dimensional he at transport equation (Equa tion 3.1, Figure 3.1) that includes multiple heat sources/sinks 76 such as short wave radiation, long wave radiation, back radiation, latent and sensible heat s , and hea t exc hange from subsurface, as shown below: (3.1) where T s (K) is the stream temperature; Q v (W m - 2 ) is the sum of the net heat fluxes excluding the groundwater advective heat flux ; Q adv ( W m - 2 ) is the advective heat flux due to groundwater flow; A (m 2 ) and W (m) are the cross - sectional area and wetted width of the stream segments respectively, Q (m 3 s - 1 ) denotes stream discharge and C (J kg - 1 K - 1 ) and (kg m - 3 ) are the specific heat and density of water , respectively . Q v includes multiple heat sources/sinks and can be expressed as: (3.2) where Q s (W m - 2 ) is the net short wave radiation reaching the s tream surface; Q l (W m - 2 ) is the net long wave radiation; a nd Q h (W m - 2 ) and Q e (W m - 2 ) are the sensible and latent heat flux es ; Q f (W m - 2 ) is the friction heat flux, Q b (W m - 2 ) is the stream bed heat conduction flux; Q m ( W m - 2 ) is the heat flux used to melt the snow when there is snow accumulation in the channel ; and d (m) i s the wetted depth of the river. Detailed calculations of these heat sources/sinks in e quation s 3. 1 and 3.2 are described in ( Qiu et al., 2019a ). The hybrid Lagrangian - Eulerian meth od based on an operator - splitting strategy was employed to solve e quation ( 3. 1) ( Niu and Phanikumar , 2015) . Following MacDonald et al., (2014) and Leach and Moore (2010) , the temperatures of stream junction segments are updated at the end of each time step based on the temperatures of the upstream and downstream segments using the mixing model concept as shown in e quation ( 3 .3 ): 77 (3 .3 ) where (K) and (K) are the stream temperatures at the upstream and downstream junction segments at the end of time step n ; (m 3 s - 1 ) and (m 3 s - 1 ) are the upstream and downstream discharge rate s ; (K) is the updated junction segment temperature at the end of time step n . The stream bed temperature module solves the one - dimensional vertical heat transfer equation based on the assumption that the lateral heat transfer within the streambed is negligi ble (Caissie et al., 2014) . The stream temperature serves as the top boundary condition of the streambed equation and the average unconfined groundwater temperature is employed as the bottom boundary condition. PAWS solves for temperature in the groundwate r domain; however, since detailed groundwater temperature data are not available for the KRW, we used the mean annual air temperature during 2003 - 2009 (10.07 o C ) as the unconfined groundwater temperature based on previous successful model applications base d on this assumption ( Cox and Bolte, 2007 ; Leach and Moore, 2010). Additionally, PAWS uses the soil/snow temperature simulation in CLM which is based on the solution of the one - dimensional he at conduction equation (Oleson et al., 2010) . Snow melt is simula ted if the snow layer temperature is above the freezing point. 78 Figure 3 . 1 Schematic illustration of the river temperature energy components for Kalamazoo river watershed . 3.2.2. Study sites and data source s Located in the southwest L ower P eni nsula of Michigan, USA, the Kalamazoo River Watershed (KRW) drains 5,200 km 2 of glacial deposits that deliver water to the river primarily (>70%) via groundwater flow paths (Rheaume 1990; Figure 3.2). The KRW has an av erage annual precipitation of approximately 960 mm with about half falling as snow, and increasing annual snowfall is observed from the head waters to the outlet due to the effect of Lake Michigan (Wesley, 2005). Mean daily air temperature s for the recent decade range from approximately - 24 to 3 8 . The land surface elevation ranges from 175 to 380 m.a.s.l. Agricult ure (dominated by corn and soybeans) is the primary land use which 79 occupies approximately 47% of the watershed , followed by forest (21%), open land (9 %), and urban (7%), see Figu re 3.3 . Figure 3 . 2 Map of the Kalamazoo River watershed. Elevation is shown as the color gradient. National Hydrography Dataset (NHD) streams, Stream temperature observation sites, U.S. Geological Survey ( USGS) gauges, National Climatic Data C enter (NCDC) weather stations and Michigan Automatic Weather Network (MAWN) stations are shown. The PAWS model integrat ed multiple sources of Geographical Information System (GIS) data. Shen et al. (2014) elaborat ed the data integration algorithms for constructing the PAWS model, and here we briefly describe the primary data sources. The overland layer was built with the 30 m resolution National Elevation Dataset (NED ) , from the United States Geological Survey (US GS , https://nationalmap.gov/elevation. html ) for calculations of topographic processes (e.g. surface runoff and lowland storage etc.). The National 80 Hydrography Dataset (NHD, https://nationalmap.gov/hydr o.html) was overlaid on the landscape layer to extract profiles of river reaches. The river n etwork was organized by sequentially ranking the river reaches from upstream to downstream. For the land use and land cover (LULC) information, we used the Integrated Forest Monitoring Assessment and Prescription (IFMA P) 30 - m resolution raster data set pro vided by the Michigan Department of Natural Resources (MDNR, 2010). Lying underneath the overland layer are the soil layers from the Soil Survey Geographic Database (SSURGO, Soil Survey Staff , 2019 ) of the U.S. Depart ment of Agriculture, Natural Resources Conservation Services (NRCS). V an Genuchten soil parameters wer e thereafter processed using pedotransfer functions provided in Rosetta (Schaap et al., 2001). The simulation was driven by hourly climate data (e.g., pr ecipitation, solar radiation, air temperature and wind speed etc.) from the National Climatic Data Center (NCDC) of the National Oceanic and Atmospheric Administration (NOAA) and th e Michigan Automated Weather Station Network (MAWN , now called EnviroWeathe r, http:// www.enviroweather.msu.edu ) ; l ocations of the mete orological stations are shown in Figure 1 . Spatial data on depth to groundwater, based on static water levels in water sup ply wells throughout the watershed, were obtained from the State of Michiga logic data base ( https://secure1.state.mi.us/wellogic/ ) . Soil temperature data from two MAWN stations ( the Albion station and the MSU Kellogg Biological Station, MSUKBS ) are used to test the model . The nearest neighbor me thod was used to interpolate stat ions data to model grids. Onset Stowaway XTI Model 2 sensors were used for stream temperature data collection (accuracy of better than 0.6 o C). 81 3.2.3. Model setup We used a grid reso lution of 1000 ×1000 m for horizontal discretization , which produced a mesh of 101×150 grids for the entire watershed. 20 vertical layers were used to simulate the vadose zone dynamics and 2 layers for the groundwater domain (unconfined and confined aquife rs). The streams were discretized as 1 km segments for solving the governi ng stream discharge and water temperature equations. The model incorporated up to level - 5 rivers to represent a more realistic channel network and to reduce model uncertainties resu lt ing from reduced channel density (Wang and Wu, 2013) . To better represen t the vegetation effect of the riparian land cover , we resampled the fraction s of stre am segment PFTs by extracting finer - scale riparian PFTs along each 1 km stre am segment. As shown in Figure 3 , we extracted the land cover information in fine nested grid cells and recalculated the land cover categories and p ortions by adding the fraction s of different PFTs in the square boxes. Three domin ant PFTs were selected including the stream portion as a PFT of open water to calculate the weighted sum radiation fl uxe s (Qiu et al., 2019a) . Since all grid cells adjacent to stream channels will contain open water as a land cover, during land cover recla ssification a minimum portion of water should be maintained to ensure that the effects of vegetation on water are adequa tely represented. Based on sensitivity analysis, a minimum portion of 20% water wa s assigned to underscore the stream PFT representing o pen water . That is, if water is incorporated into the three dominant PFTs , we will choose max{default wat er portion with in the three PFTs , 20%} as the water portion; otherwise, the water portion is assigned as 20%, and the other three dominant land covers occupy the re maining 80%. We evaluated the impact of representing riparian vegetation by comparing the simulated results using different nested grid 82 resolutions within a larger 1 x 1 km grid cell and comparing the simulation results with observed data. The differential evolution algorithm (Price et al., 2005) was employed to tune the model parameters for improving the model performance in predicting stream discharge, groundwater heads and evapotranspiration (ET) rates . Parameters and their physical meanings were summarized in Shen et al., (2013) . Figure 3 . 3 Land use and l and cover map of Kalamazoo River watershed . 83 Figure 3 . 4 Schematic illustration of using the square boxes to extract the streamside land use information . 84 Results and Discussion In this section, we as sess model performance by comparing simulated results with observations. The primary model performance metrics employed are Nash - Sutcliffe efficiency coefficient (NASH, Nash and Sutcliffe, 1970), the coefficient of determination (R 2 ), and the root mean squ ared error (RMSE). 3.3.1. Hydrology Model performance evaluated using observed datasets for stream discharge, ground water heads and evapotranspiration during a 7 - year simulation period (2003 - 2009) has been previously reported for the KRW b y Qiu et al. (2019b) a lthough the earlier results were based on a slightly different set of calibration methods and parameters. In the earlier paper multi - site calibration was used with a different set of parameters in each sub - watershed to address the res earch questions in tha t paper while a single set of parameters are used for the entire watershed in this paper. Here we first briefly describe model results for hydrology in this section. Additionally, we present the comparisons of leaf area index (LAI) be tween simulated result s and MODIS data as vegetation growth cycles are important for simulating stream temperature. The simulated stream discharge results during the 7 - year simulation period are compared with USGS obser vations for four different gauging st ations in the watershe d, see Figure 3.5 . The NASH and RMSE values are tabulated in Table 3. 1 . The NASH and RMSE values tabulated in Table S 1 show satisfactory model performance. The simulated discharge near the watershed outlet, i.e. , the location of USGS gauge 04108660 (Kalama zoo River at New Richmond, Michigan) , matched the amplitude and fluctuations of observed streamflows fairly well. However, some peak st r eam discharge values are overestimated by the model. 85 It should be noted th at rather than purely ca librating to the strea m discharge, the model performance was constrained by calibrating to stream discharge, ground water heads and monthly ET simultaneously. The stream discharge performance may be offset because of compensation errors among different hydrologic components ( se e, for example, Anderton et al., 2002; Beven and Freer, 2001). Additionally, inaccurate represen tat ion of field heterogeneities and uncertainties within both input data and observations could also lea d to the imperfect comparisons for stream discharge ( Bev en, 1993; Beven, 2001; Anderton et al., 2002; Beven and Freer, 2001). Overall, the simulated stream discharge results are considered acceptable considering all sources of uncertainty . As shown in Figure 3.6 , the monthly simulated ET results (watershed aver age) are compared with MODIS satellite observations obtained from NASA earth data search engine (ht tps://search.earthdata.nasa.gov). The mon thly ET fluctuations match very well between simulated results and MOD16 product. In some summer months, simulated r esults overestimate the ET, while during winter months the model outputs underestimate ET relative to MOD16 data. This phenomenon could be attributed to the different algorithms used in MOD16 product and PAWS+CLM. MOD16 product is derived based on the Penm an - Monteith formulation (Mu et al., 2011), while a resistance approach based on the two - big leaf mo del (Dai et al., 2004) is used in PAWS+CLM to calculate ET. 86 Figure 3 . 5 Stream discharge comparison s between simulated results and observations by USGS gauges . Table 3 . 1 Performance of Stream discharge comparisons between simulated results and observations by USGS gauges USGS station N o. NASH RMSE 0410301 0 0.61 2.04 04106000 0.63 10.18 04108600 0.64 1.79 04108660 0.76 20.35 87 Figure 3 . 6 Monthly watershed - averaged ET c omparison s between simulated and remotely sensed MODIS ET produc t s . 88 Figure 3 . 7 Spatial map of yearly averaged evapotranspiration for the Kalamazoo River watershed for the 7 - year period (2003 2009) of (a) simulated output and (b) MODIS data. Table 3 . 2 Pair - wise linear correlation coefficients for different Land use/land cover types and soil types with GLB and MLT simulated ET outputs Land Use/Land Cover Soils Forest Grass Crops Urban Wetland Sand Clay Organic Matter 0.28 - 0.05 - 1.99 - 0.24 0.27 0.53 - 0.50 0.08 p 0 0.00 0.00 0.00 0.00 0.00 0.00 0.00 is the Correlation Coefficient. p < 0.05 indicates statistical significance. 89 Figures 3.7 (a) and 3.7 (b) show the annual - average spatial maps of ET based on the simula tions and MODIS16 data respectively. Annual average ET of the 7 - year simulation period is 583.43 mm yr - 1 , which is comparable to the MODIS value, 559.89 mm yr - 1 . The spatial maps of ET from simulations and MODIS data generally follow a similar pattern. A l inear correlation analysis was performed on the spatially distributed ET simulated values against the LU LC types and soil types and the results are summarized in Table 3.2 . In this table, is the pair - wise linear correlation coefficient, and p is the probability for testing the null hypothesis of no correlation against the alternative hypothesis that there is a significant correlation. If p is small (i.e., < 0.05), then the correlation i s significantly - - 0.24) are negatively correlated to ET. The grass portion within KRW tends to show median ET values among all the LULC types and shows the least cor - 0.04) with the spatially distributed ET values and oscillates slightly between negative and positive with the two different calibration methods. High positive correlations of ET values correlation with - 0.50) indicate that the ET is also closely related to soil infiltration capacity. The PAWS model outputs resolved the ET heterogeneity better than did the remotely sensed MODIS data. The ma jor land cover in northwestern KRW is forest and there a re many lakes and reservoirs located in the middle of the watershed. Therefore, we expect high ET values within this area as shown in simulated ET maps. The south - central areas of Kalamazoo (where the MODIS data are blank) are urban areas, which correspond to the low ET values in PAWS output (colored blue). The comparison between simulated depth to groundwater table and well - logic data is shown in Figure 3 .8 . A NASH 90 value of 0.90 shows that a good agre ement was obtained between si mulated land surface dept h to water table and observed data. There is a slight trend for the simulated depth to water table to move below the 1:1 line which indicates, on average, higher simulated groundwater table depth compa red with well - logic data. This discrepancy could be at tributed to several reasons, incl uding potentially overestimated recharge, underestimated groundwater contribution to rivers and uncertainties in the hydraulic conductivities. Figure 3 . 8 Plots of simulated versus observed depth to groundwater table (from Well - logic data set) for each computation grid cell. NASH is the Nash - Sutcliffe efficiency metric. 91 Figure 3 . 9 Spatial and temporal comparisons between MODIS product and the PAWS+CLM simulated LAI. (a) Time series comparison for four representative vegetation type s. (b) Spatial map comparison for June 21 - 28, 2007 and December 19 - 26, 2007. 92 Additionally, we compa re the spatial and temporal LAI values between simulated results and level - 4, 8 - day, MODIS data (MOD15, https://search.earthdata.nasa.gov ), see Figure 3.9 . Figure 3.9 (a) shows the temporal LAI compari sons between MOD15 and simulated results for four representative vegetation cover during 2005 - 2007, which the representative vegetation occupy more than 85% of the simulated grid cell. Model results show abrupt LAI increa se in growing season and LAI decrea se in wilting season for all the four representative vegetations, while MODIS data show r elatively smooth change. This dif ference is most obvious for C3 grass. Moreover, the model overestimates the duration of high LAI va lues of C3 grass compared with MODI S data. This discrepancy may be due to the phenological parameters utilized in C3 crop of CLM and uncertainties associated with MODIS product. Several outliers are observed for Corn during the growing season of 2005. In this period, the model results sho w LAI in the range of 3 ~ 4 while several observed LAI values are above 5. These outliers may be attributed to the uncertainties of the remotely sensed MODIS data, including cloud contamination or other data quality issu es. Figure 3.9 (b) shows the compar ison of spatially distributed LAI for two 8 - day period in June and December, respectively. During June, the spatial pattern of LAI is mostly matched between simulated results and observations, including the high LAI value s in the northwest and some parts w ithin the middle of the watershed. Forest is the major land cover in northwestern KRW which tends to show high LAI values during growing season. Meanwhile, the lake effect may increase the length of growing season from ~ 150 days at the eastern end of KRW to ~ 180 days near the watershed mouth and Lake Michigan, which could also contribute to high LAI. Two blank areas in the south - central part of KRW of MODIS data are urban areas, where the model generates very low LAI va lues. Although 93 the dominant land use is urban for these two areas, some minor portion of land use contributions may still produce positive LAI, for example, the grass and trees in parks. This information is reflected in the data integration algorithm of PAWS+CLM, which assimilate data at fine resolutions and upscale to grid resolution by preserving the major information (Shen et al., 2014). During December, for the most part both the simulated and MODIS data show nearly zero LAI values. The northwest par t in December show s relative ly high er LAI values which are primarily coming from evergreen forest. 3.3.2. S oil moisture and soil temperature results Figures 3.10 show s the 10 cm soil moisture comparisons at two MAWN stations. It should be noted that the observed data represent a point measurement as data were collected using a Campbell Scientific CS616 water content reflectometer (WCR) whereas our simulated results represent an average of a grid cell domain with area of 400×400 m 2 . At station Albion (Figure 3.10 (a)), simulated soil moistures show almost the same trend comparing with observations but generally lower in winter and higher in summer. For example, around February 2004, the simulated soil moisture values are below 0.1 while the observed soil moist ure values are between 0.2 and 0.25. Arou nd July 2005, the simulated soil moistures are above 0.1 while the observations are slightly lower. At Michigan State University Kellogg Biological station (MSUKBS) (Figure 3.10 (b)), the relatively higher soil mois tures could not be simulated accurately i n February 2009, which may due to the underestimated rainfall intensity during this period. 94 Figure 3 . 10 10 cm Soil Moisture comparisons of simulated outputs wit h MAWN (Michigan Automatic Weather Networ k) station observations at (a) Albion and (b) MSUKBS. Sim is the simulated outputs; Obs is the MAWN station observations. Figure 3 . 11 10 cm Soil Temperature comparisons of simulated outputs with MAWN (Mich igan Automatic Weather Network) station observations at (a) Albion and (b) MSUKBS. Sim is the simulated outputs; Obs is the MAWN station observations. 95 The simulated soil temperatures at a depth of 10 cm are compared with observations at two MAWN stations in Figure 3.11 . The comparison at station MSUKBS is only shown from December 2007 to July 2009 based on data availability. Except in the coldest period s, the daily fluctuations of soil temperature are captured quite well by the model, with NASH values of station (MSUKBS), respectively. The amplitude of simulated summer soil tempe ratures matched well with observations at Albion station as shown in Figure 3.11 (a) although at statio n MSUKBS, the amplitudes of peak summer temperatures are slightly underestimated, which may be due to overestimated vegetation shading effect at this sta tion , see F igure 3.11 ( b). In the coldest winter periods, the observed temperatures show nearly constan t values around 0 o C, whereas the simulated results show much lower temperatures. These discrepancies are due to the known limitations in describing the movement of water through frozen soil in CLM 4.0 and are associated with soil hydraulic and thermal pro perty schemes used while simulating freeze - thaw cycles (Swenson et al., 2012; Chen et al., 2018; Fang et al., 2016). Fortunately, accurate simulation of the coldest soil temperatures is not important to the goals of this study because as long as the surfac e soils are frozen there is neither infiltration nor overland runoff. 3.3.3. S tream temperature simulation s 3.3.3.1 Effects of resolving spatial heterogeneity in veget ation We evaluate the effect s of resolving the spatial variability in vegetation on stream temperature by testing the model performance using different nested land use resolutions within model grid cells near the streams. The sizes of the nested grids with in the 1 km x 1 km model grid cell considered are: 6 0 m × 6 0 m, 9 0 m × 9 0 m, 2 5 0 m × 2 5 0 m, and 1000 96 m × 1000 m, respectively. The stream temperature comparisons are performed at eight observation al sites when observed data are available during the simula tion period. Locations of the observati on sites are shown in Figure 3. 1 . A ll the comparisons are shown at hourly scale, except for the Kalamazoo River station near Allegan ( USGS # 04107850 ) , for which the comparison of daily average temperatures is reported due to data availability . The performance metrics for the tests are tabulated in Table 3.3 . Results with land use resolutions of 60 m × 60 m, 90 m × 9 0 m show similar performance, and are generally better than model performance ba sed on resolutions of 250 m × 250 m, and 1000 m × 1000 m . At site Marshal, 90 m × 9 0 m showed the best performance among the f our test cases. Three observation sites are selected to present the performance difference, as shown in Figure 3.12 . The land use PFTs and percentages for the stream segments where the three observational sites are located are shown in Figure 3.13 . As the resolution of the nested grids used to represent land use within a model grid cell changes, the dominant land use fractions within the model grid cell change producing changes in stream temperature due to changes in riparian climate. For example, at the Al legan site (USGS # 04107850), deciduous broadleaf is the most dominant streamside land use for all the four resolutions tested whi le the 250 m resolution dataset produced the largest portion of deciduous broadleaf (41.8%). Open water proportions for 60 m a nd 90 m resolution streamside land use datasets exceeded the minimum threshold of 20%, while open water portions with resolutions of 250 m and 1000 m are estimated as 20%. As a result , the summer temperature at this site was underestimated using 250 m and 1000 m streamside land use information. In contrast, at site Spring Brook, the summer temperature was overestimated using land us e information from 250 m and 1000 m nested grids. This could be accounted from 97 relatively small percentage of deciduous broadl eaf, thus relatively small summer shading effect, using 250 m and 1000 m streamside land use. At site Bear Creek, dominant land u se types and percentages for all the four different resolutions are almost the same. Accordingly, the simulated stream tempera ture results did not show significant differences. All subsequent results of stream temperature presented in the paper are based o n the optimum resolution of 90 m x 90 m for resampling the land use information within a model grid cell of size 1 km x 1 km. Table 3 . 3 Performance of the temperature comparisons between simulated results and observations using different scales of stream side land use information Site Bear Creek Schnable Creek Silver Creek Indian Cr eek Scale NASH RMSE NASH RMSE NASH RMSE NASH RMSE 60m 0.85 1.88 0.83 2.38 0.78 2.58 0.85 1.94 90m 0.85 1.86 0.84 2.36 0.78 2.61 0.85 1.94 500m 0.85 1.89 0.81 2.53 0.76 2.74 0.83 2.12 1000m 0.85 1.89 0.83 2.44 0.77 2.68 0.83 2.14 Site Marshall Rice C reek Spring Creek USGS04107850 Scale NASH RMSE NASH RMSE NASH RMSE NASH RMSE 60m 0.85 1.81 0.82 1.56 0.85 2.05 0.89 1.67 90m 0.86 1.70 0.83 1.43 0.86 2.04 0.89 1.63 500m 0.85 1.79 0.83 1.48 0.79 3.04 0.83 2.73 1000m 0.84 2.09 0.80 1.68 0.79 3.06 0.84 2.76 98 Figure 3 . 12 Performance difference of stream temperature comparisons using different scales of stream side land use information for three selected observation sites: (a) USGS04107850 (b) 1:1 plot of USGS04107850; (c) Bear Creek, (d) 1:1 plot of Bear Creek; (e) Spring Brook; (f) 1:1 plot of Spring Bro ok. OBS is the observation. 99 Figure 3 . 13 Land use types and percentages for different tested land us e resolution s of the 1 km stream segment for the presented three observation sites. 3.3.3.2 Compared with observ ed stream temperature data Comparisons of simulated and observed s tream temperature for the eight observation sites are shown in Figure 3.14 . For the o bservation sites at Indian Creek (Figure 3.14 f) and Marshal l (Figure 3.14 g), simulated stream temperature results show negative temperature values during some days in winter while the observed temperature remain s constant around 0 o C , which is similar to th e discrepancy in soil temperature comparisons presented earlier. This is probably due to overestimation of subzero soil temperatures in winter as noted earlier which also influence stream temperatures via exchange fluxes. In addition, the high turbulence o f flowing water tends to prevent freezing (Mohseni and Stefan, 1999). Alt hough friction heat flux due to fluid friction at the riverbed and the banks is included in the modeling, detailed effects of turbulence were not explicitly modeled. Previous researc h suggested that neglecting the energy associated with s nowmelt could imp act stream temperature simulations (Caissie et al., 2007). Although the PAWS model (via CLM ) takes 100 the energy of snowmelt into consideration as described in (Qiu et al., 2019a) , freez ing and thaw ing of stream ice are not explicitly simulated. Future effort s will focus on improving the representations of stream ice and phase change processes. At the Silver Creek site, the winter temperature is underestimated by the model from December 2008 to February 2009 ( Figure 3.14b ). Underestimated warm heat fluxes fro m the subsurface during winter period s may be responsible for this mismatch. At the Schnable Brook site, the amplitude of diurnal temperature fluctuations is overestimated by the mode l which may be due to inaccurate representation of stream geometry and/or inaccurate stream discharge outputs ( Figure 3.14d ). Diurnal stream temperature fluctuations are found to be sensitive to stream geometries (width and depth), e.g. MacDonald et al., ( 2014) . The stream width values were estimated from report ed values (Wesley, 2005) and Google E arth, which may have introduced uncertainties in to the stream width representation . In addition, stream discharge influences the stream heat capacity and the adve ctive heat fluxes. Although the stream discharges were calibrated to several USGS gauges, it is possible for mismatches to exist between simulated discharge and actual discharge in uncalibrated tributaries due to a combination of factors. These factors may include inaccu rate representation of the heterogeneities in strea m properties (such as stream depths and widths) as well as processes not simulated (such as tile drains) . At the Marshal l site, the summer temperature of 2008 is slightly underestimated by the model (Figure 3.14g) . In contrast, at the Indian Creek site, t he summer temperature of 2009 is slightly overestimated (Figure 3.14f). Possible explanations for the se discrepancies in summer temperatures could include the inaccurate estimation of vegeta tion shading effect s and/or the heat fluxes from the subsurface. M oreover, the assumption of a minimum 20% of water fraction within the grid cell may 101 not exactly describe the reality, for example, forest streams may be fully shaded by broadleaf trees in su mmer. It should also be noted that the observed data are point mea surements while the model results represent average values of temperature in a 1 - km stream seg ment. Discrepancies could occur if the vegetation covers around the sampled points are significa nt ly different from the average cover along the 1 - km stream segmen t. Despite these general observations and caveats, we conclude that a good overall agreement is noted between observations and model results (Figure 3.14). 102 Figure 3 . 14 Stream temperature comparison of simulated results with observations for the eight observation sites: (a) Bear Creek, (b) Silver Creek, (c) Rice Creek, (d) Schnable Brook, (e) Spring Brook, (f) Indian Creek, (g ) Marshal, (h ) Allegan, USGS04107850. Sim denotes the simulated results; Obs denotes the observations. 103 3.3.3.3 Seasonal heat budgets As shown in Figure 3.15 , average seasonal heat flux budgets during the simulation period are calculated at the temperature sampli ng sites . For all the eight sampling sites, the dominant heat source is net radiation and the primary heat sink is latent heat , which is consistent w ith discussions in the literature (e.g. MacDonald et al., 2014; Leach and Moore, 2010) . Sensible heat flux is a heat sink in s pring, autumn and w inter, while it shifts to a heat source in s ummer. Heat flux from the subsurface is an important component and switch es from a source in spring and w inter to a sink in summer and f all during the si mul ation period; ther efore s ubsurface heat fluxes serve as a heat buffer to cool down the high summer temperatures and warm up the chilly winter temperatures. This brings out the importance of investigating the connections between surface and ground water t emperatures and thei r co - evolving responses in a warming climate (Kurylyk et al., 2013). Friction heat flux and snow melt heat flux occupy a small portion of the heat budget and show relatively less impact on the stream thermal regimes in this study. Give n the similar latitud e s of the eight sampling sites, incoming solar radiation and longwave radiation flux es can be expected to be similar for all sites . However, distinct seasonal net radiation heat fluxes are shown for the eight sampling sites. This bring s out the significant role the vegetation shading effect plays i n altering the net radiation fluxes received by the stream surface in different reaches . Figure 3.16 shows average spatially distributed stream temperatures for the KRW during the simulation period for summer (Fi gure 3.16(a)) and winter (Figure 3.16(b)) , respectively. Five distinct locations in the watershed are marked using rectangular boxes for further examination of the spatial patterns . At location (1), the two highlighted stream segments 104 are connected to seve ral lakes which are less affected by the vegetation shading effect, thus the summer temperature are relatively high. Similarly, at location (2), the streamside land use is dominated by urban and open land, high summer temperature is al so shown as the veget ation shading effect becomes less important . In contrast, l ocation (3) shows relatively cool summer temperature s and low winter temperature s , which are primarily due to streamside heavy f orest proportion (Broad leaf deciduous forest tr ee). The monitoring site Spring Brook is located at location (4) which shows relatively low summer temperature and a relatively warm winter temperature. The higher subsurface heat flux contributions at the site are the primary reason to stabilize the temperatures in tributari es of location (4), see also Figure 3.15 . In contrast, at location (5), simulated results show relatively high temperature s in summer and average frozen temperatures in winter, which can be attributed to the relatively small subsur face heat flux contributi ons. 105 Figure 3 . 15 Simulated seasonal heat fluxes budgets for the eight observation sites : (a) Bear Creek, (b) Silver Creek, (c) Rice Creek, (d) Schnable Brook, (e) Spring Brook, (f) Indian Creek, (g) Marshal, (h) Allegan, USGS04107850. Net radiation heat flux is the sum of short wave radiation and net long wave radiation heat flux. Subsurface heat flux is the sum of streambed conduction heat flux and advective groundwater h eat flux. 106 Figure 3 . 16 Average spatially distributed (a) summer and (b) winter temperatures for the KRW stream network during the simulation period (2003 - 2010). Example reaches shown in boxes are discussed in the text. 107 3.3.3.4 Stream temper ature responses to potential deforestation To evaluate the potential effects of deforestation on high summer temperatures, a modeling scenario was conducted by replacing 10% of the dec iduous broadleaf forest with C3 grass for river segments that contained a riparian land cover of more than 10% deciduous broadleaf forest. This change is only incorporated into calculations of stream incoming radiation heat fluxes, while preserving calcula tions of other hydrologic variables such as river discharge and ET. A ma p of the June July average daily stream temperature difference (simulated temperature from this scenario minus the baseline temperature) is shown in Figure 3. 1 7 . As expected, small a nd lower - order (usually headwater) streams are more sensitive to potenti al deforestation effects due to their relatively small heat capacities in comparison with large and higher - order rivers. The extreme temperature increase could reach 6 7 o C . Three example stream segments that initially have more than 10% riparian deciduo us broadleaf forest are highlighted to understand the temperature changes under the deforestation scenario, marked as boxes 1 - 3 in Figure 3.17 . Results for three relevant hydro logic properties the average stream width, average stream depth, and the baseflo w index (BFI), which is the ratio of groundwater contribution to total stream discharge volume during June July are shown in Table 2. Segment 1 shows the highest temperature increase among the three segments. It has moderate stream width and depth value s, but the lowest BFI values. Lacking the buffering effect from inputs of cooler groundwater, this stream segment is most vulnerable to potential deforestation. Segment 2 is o ne of the headwater streams of the Kalamazoo River and has a BFI of more than ha lf during June July. However, the small stream size increased its sensitivity to potential 108 deforestation, and as a result the temperature increased by 4 - 4.5 o C. Segment 3 shows the least temperature increase; a relatively high BFI value and stream siz e render this stream least sensitive to potential deforestation effects. In these scenarios, we used the mean annual air temperature as the groundwater temperature, as explaine d earlier. However, the shallow groundwater temperatures could become responsive to short - term air temperature and land surface temperature changes (Kurylyk et al., 2013), which may in turn influence the surface water temperature. And ultimately the ground water temperature in the aquifers that supply much of the stream discharge will respond to the warming climate, though likely over decadal time scales in this watershed. It will be interesting to further explore the coevolving effects of land cover change and climate change on stream and groundwater temperature regimes. Table 3 . 4 Characteristics of the marked three stream segments Stream ID Average Width (m) Average Depth (m) Baseflow index 1 4.3 1.68 18% 2 2.6 1.12 46% 3 4.8 1.86 45% 109 Figure 3 . 17 Temperature difference map for the KRW stream network under a model scenario of riparian deforestation. Example reaches shown in boxes are discussed in the text. 110 Conclusion This study extended the application of the PAWS model and its newly developed temperature module to a 5400 km 2 mixe d land use watershed. We incorporated publicly available spatial data sets, tested the optimal spatial resolution for representing riparian vegetation effects, and validated our process - based simulations with field measurements of soil and stream water temperature regimes. We applied the model to quantitively evaluate the effects of riparian vegetation on the stream temperatures, accounting for land use and cover and seasonal vegetation phenolo gy . The model results suggested overall satisfactory performance considering several sources of uncertainties. Results of this work indicate the possi bility to extend the application of the developed model to evaluate potential impacts of climate and land cover changes on the stream temperature regimes, which will be increasingly important to identify effective strategies for fisheries and aquatic ecosy stem management. . 111 Chapter 4. An integrated, catchment - scale framework to model temperature - dependent ni trogen transport and transformation s This chapter describes the development of a process - based, integrated hydrology, temperature and nitrogen transpor t and transformation s model. This integrated modeling approach enabled real - time (or near real - time) fore casts of nitrogen concentrations in different hydrologic domains by seamlessly integrat ing with a mechanistic, grid - based hydrologic model. T he model performance was evaluated by applying the model to two watersheds in different parts of the world, i.e. , t he Wood Brook watershed in the U.K. and the Kalamazoo River watershed in the U.S.A. Introduction Nitrogen loading and transport in river basin s are closely related to several environmental issues such as hypoxia and eutrophication (Heisler et al., 2008) . Excessive nitrogen in waterbodies stimulates the growth of harmful algal blooms (HABs) which may result in oxygen depletio n and mortalit y of aquatic species (Anderson et al., 2002) . In addition to deterior ated water quality, extensive nitrogen enrichment may stimulate greenhouse gas emission s , lead to soil acidifica tion, and induce production of tropospheric ozone and aerosols which may cause respiratory illness (Pope et al., 1995 ; Galloway et al., 2003). Given the rapid growth of industrial and e conomic development, negative impacts result ing from excessive nitrogen emission have been considered as one of the biggest world - wide pollution problems (Heisler et al., 2008; Galloway et al., 2008). Increasing riverine nitrogen loading is usually related to anthropogenic activities, e.g. agricultural practices including the application of fertilizers and manures (Carpenter et al., 1998 ; Boyer et al., 2002 ; Beman et al., 2005 ;Ga lloway et a l., 2008), land use change and 112 urbanization (Thomas et al., 2016) . Interacting with these activities are the regional climate conditions, hydrological p rocesses and biogeochemical processes (Breemen et al., 2002) . Particularly, temperature was found to play an important role in affecting the proportion of riverine nitrogen export to the total catchment nitrogen input (Schaefer and Alber, 2007 ; Miller et al., 2016) . Positive relationship between temperature a n d denitrification was commonly found by researches, e.g. Holmes et al., (1996), Strauss et al., (2002), Starry et al., (2005), Schaefer and Alber (2007). Additionally, the structure of drainage network, stream discharge and the stream size were found to b e important factors influencing stream nitrogen transformation processes (Wollheim et al., 2006 ; Helton et al., 2018) . S oil moisture was identified as a nother important hydrologic factor in controlling plant nitrogen uptake, as well as nitrogen mineralization and denitrification processes in the vadose zone ( Porporato et al., 2003 ; . Regarding those import ant hydrologic and topographic controls on the nitrogen fate and transport processes, understanding the r iverine nitrogen export mechanisms requires a joint analysis of the regional climate, catchment characteristics, anthropogenic processes such as land u se development and fertilizer applications, as well as hydrologic and biogeochemical processes (Breemen et al., 2002) . Therefore, there is a ne ed to develop a tool that could be effectively and efficiently utilized to simulat e nitrogen reaction s an d transport in different hydrologic domains, predict riverine nitrogen expo rt based on multiple nitrogen sources and inform sustainable water resources , agricultural management and decision making. Meanwhile, the tool would be utilized to improve the abili ty to understand the joint functioning of the water - temperature - nitrogen system under different climate conditions, catchment characteristics and influ enc es of anthropogenic activities. 113 Previously, t hree primary types of catchment - scale models have been de veloped that link hydrological and nitrogen reaction and transport processes: statistical models, conceptual or semi - distributed models and fully distr ibuted process - based (or mechanistic) models. Statistical models, such as the SPARROW model (Schwarz et al., 2006; Ro bertson and Saad, 2011) , estimate the riverine nitrogen export by using statistically significant factors including climate and topographic conditions , land use types, soil pr operties and drainage densities. While statistical models are usually efficient for quantitative analysis due to their light computational expense and low data demand, they lack the capability to identify temporal cause - effect rela tionship s and are unable to address the underlying mechanics driv ing the interactions among different pro cesses. Conceptual or semi - process based models, e.g. SWAT (Arnold et al., 1998) , HSPF (Bicknell et al., 1997), INCA (Whitehead et al., 1998; Wade et al., 2002 ) and LASCAM (Viney et al., 2000), trace the water and nitrogen species in different hydrologic domains while conceptualiz ing some of the driv ing processes using empirical equations or lumped parameters. For example, the SWAT model simplifies the processe s in the subsurface by using empirical equations, which lack the capability to track nitrogen species in the subsurface aquifers. While the INCA model includes nitrogen reactions in different hydrologic domains, it lacks the ability to track the nitrogen v ariations in space due to simplifying assumptions that solve only time - dependent ordinary differential e quations . Some f ully process - based mechanistic models, e.g. MOHID (Neves, 1985; Trancoso et al., 2009), use time - marching simulations by solving the go verning mechanistic equations. Explicitly simulated water and nitrogen fluxes provide the ability to eval uate the impacts of both point and non - point sources, track individual processes, and to understand the interactions among different hydrologic 114 domains . However, MOHID solves the thre e - dimensional Richards equation and demand s high computational expense s , especially when applied to large - scale watershed s. There remain opportunities for integration of hydrologic and biogeochemical processes into model frameworks that could efficiently couple surface and subsurface processes and be adaptively applied to diffe rent catchments, for developing flexible model structures that could strike a balance between model complexity and computational expense. The coupled PAWS+CLM model (Shen and Phanikumar, 2010; Shen et al., 2013) has readily incorporated the land surface pr ocesses, subsurface processes, biogeochamical processes and the interactions among different processes. With little loss of physics, PAWS+CLM significantly reduces the computational demand by simplifying the fully three - dimensional (3 - D) subsurface model t o one - - 3 - D saturated gr oundwater domain (Shen et al., 2013). The integration of PAWS+CLM with transport model has been effectively applied for simulating the bacteria fate and transport in Red Cedar River watershed, Michigan (Niu and Phanikumar, 2015). PAWS is recently coupled w ith century - based nitrogen modules in CLM4.5 to update the nitrogen and carbon dynamics with user - defined fertilizer nitrogen flux rather than using default fertilizer flux in CLM4.0 (Oleson et al., 2010 ; Oleson et al., 2013). Given the challenges and op portunities, there remain strong scientific and technological motivations to build a watershed - scale model with integrated hydrologic, thermodynamic and biogeochemical processes based on curren t framework of PAWS+CLM. This innovati ve integrated modeling a pproach should be able to make real - time (or near real - time) forecasts of nitrogen concentrations at different hydrologic domains by tracking the nitrogen reaction s and fate at different hydrol ogic components. 115 The main objectives of this chapter are to: 1 ) report the development of an integrated, catchment - scale framework based on the original PAWS+CLM framework to model nitrogen reaction s and transport processes under the influence of hydrolog ic, biogeochem ical and anthropogenic impacts ; 2 ) test the newly coupled biogeochemistry modules of century based nitrogen model in CLM4.5 by investigating the nitrogen leaching results 3) evaluate the influence of temperature - dependent reaction rates for n itrification and denitrification within streams and streambeds. 4) Test t he model performance by applying the model to two different watersheds with different sizes and climate conditions and by investigating the interactions of multiple hydrologic and bio geochemical processes. The findings of this research are expect ed to aid the management of agricultural and aquatic ecosystem s . Methods and materials The detailed hydrologic conditions of the two tested watersheds, i.e. the Kalamazoo river watershed (KRW) and the Wood brook watershed (WBW) have been extensively discu ssed in former chapters. In this section, after briefly introduc ing the nitrogen observation data, we describe the model framework of nitrogen processes and the nitrogen sources for the two wat ersheds. 4.2.1. The nitrogen observation data The nitrogen observation data are primarily obtained from Baas (2009) and the National Water Quality Monitoring Council ( http s://www.waterqualitydata.us/portal/ ), i.e. the Storet database. The locations of the samplin g points and monitoring stations are shown in Figure 4.1 . For WBW, the nitrate concentrations are measured at the catchment outlet using an OPUS UV spectral sensor (Blaen et al., 2017). 116 Figure 4 . 1 Map of t he Kalamazoo River watershed. Elevation is shown as the color gradient. National Hydrography Dataset (NHD) streams, St oret nitrogen stations , and nitrate sampling sites in Baas ( 2009), U.S. Geological Survey (USGS) gauges, National Climatic Data Center (NC DC) weather stations and Michigan Automatic Weather Network (MAWN) stations are shown. 4.2.2. The model framework The nitrogen reaction and transport module was developed within the framework of the integrated hydrologic model PAWS (Process - based Adaptive Watersh ed Simulator, Shen et al., 2016; Shen et al., 2014; Shen et al., 2013) and its bacteria transport model (Niu and Phanikumar, 2015) and the newly developed temperature model (Qiu et al. , 2019a) as described in Chapter 2. Figure 4. 2 shows the flow chart of how the different modules are coupled. Primarily, s tream nitrogen dynamics were simulated using the one - dimensional reaction and transport equation s considering multiple species, i.e. Nitrate ( NO 3 - N ), nitrite 117 ( NO 2 - N ), ammonia ( NH 4 - N ) and organic nitrogen. A quasi three - dimensional reaction and transport equation was used to simulate the groundwater nitrogen dyn amics in the unconfined aquifer . A two - dimensional reaction and transport equation was employed to simulate the nitrogen dynamics in overland flow. The nitrogen processes in the soil column are simulated by coupling PAWS with the Community Land Model version 4.5 ( Ole son et al., 2013). I mportant processes that are simulated in this framewor k are illustrated in Figure 4. 3 . In t his se ction, we will des cribe the nitrogen sources, the model framework and governing equations in detail for different hydrologic domains. Figure 4 . 2 F low diagram for the essential processes of PAWS+CLM. The modules in the blue boxes are the primary processes for the hydrologic components, the modules in the red boxes are the primary processes of temperature simulat ions, and the modules in the black boxes are primary sub - processes. V/M/H Transfer denotes Vapor/Momentum/Heat Trans fer. 118 Figure 4 . 3 Schematic showing the major processes of the nitrogen cycle simulated . 4.2.3. The Nitrogen sources Major n itrogen sources include atmospheric deposition, biological fixation, fertilizer and manur e application, septic tanks and point sources (Breemen et al., 2002; Chapra, 2008 ). The PAWs model is coupled with CLM4. 5 to a ccount for atmospher ic deposition and biological fixation of nitrogen which will be described in a later section . Here, we describ e how the amount and spatial distribution of nitrogen sources from fertilizer/manure are treated in this study. 4.2.3.1 T h e fertilizer/manure amount and timing: Kalamazoo River watershed For determining the timing and spatial distribution of intra - annu al fertilize r applications , we used data from United States Department of Agriculture (USDA, 2013) w hich has fertilizer applica tion records for certain important crops including corns , soybeans and 119 wheats . Crop rotation between corn and soybean, although not shown in Table 4.1, is considered by referring to the USDA CropScape ( https://nassgeodata.gmu.edu/CropScape/ ) map to reduce the nitrogen expense and to increase corn yields. We used the manure amount in Luscz et al. (2015 ) which estimated spatially distribute d annual net nitrogen input of manure applica t ions for several watersheds in Michigan including the Kalamazoo River watershed. The monthly application amount s are thereafter distributed into daily amount s by assumin g fertilizer /manure are applied duri ng continuous three days without rainfall. The amount of nitrogen fertilizer/manure used is averaged over the computation al grid cells based on the types of crop land. The applications of fertilizers/manures for Kalamazoo watershed are performed mainly in April, May, June and October. Table 4.1 lists the average nitrogen amount applied to the corn, soybean an d wheat farmland during April, May, June and October for the Kalamazoo River waters hed from 2003 - 2010. 4.2.3.2 The fertilizer/manure amount and timing: Wood Brook watershed For WBW, the date and amount of each fertilizer/manure application are recorded and meas ured for each crop field during the year 2015~ 2017, which cover the whole simulati on period in this study. These records are directly used in the study. Table 4.2 lists the dates and total nitrogen amounts from fertilizer/manure applied to WBW. 120 Ta ble 4 . 1 Fertilizer/manure nitrogen application for KRW (USDA, 2013) Nitrogen applied from Fertilizer/Manure (kg/ha) corn soybean wheat April 79 . 06 3 . 89 0.0 0 May 1543.3 0 0.0 0 10 .20 June 31 . 51 1 . 46 0.0 0 O ctober 17 . 00 0. 9 3 121 . 94 Table 4 . 2 Fertilizer/manure nitrogen application for Wood Brook watershed Date Nitrogen applied from Fertilizer/Manure (k g/ ha ) 03/18/2016 3 5 .48 04/04/2016 4 0 .86 06/ 1 5 /2016 48.20 07 /10/2017 33 . 60 03/15/ 2017 34 . 00 121 4.2.4. Nitrogen reaction and transport in river The one - dimensional channel transport and reacti ons involving multiple nitrogen species (i.e. NO 3 - N, NO 2 - N, NH 4 - N and organic - N ) follow the advection - dispersion - reaction equat ions below ( Gunduz, 2004 ; Chapra, 2008; Alam and Dutta, 2016) : (4.1) (4.2) (4.3) (4.4) where D L is the longit udinal dispersion coefficient (L 2 T - 1 ); is the mean velocity of the river segment (LT - 1 ); A is the cross - sectional area of the stream segment; K r1 is the organic N accumulation rate due to algal growth ; K r2 is the decay rate of or ganic N ( T - 1 ) , K r3 is the loss rate of ammonia ( T - 1 ) , K r4 is the nitrification rate ( T - 1 ) , K r5 is the denitrification rate (T - 1 ) and q 1 , q 2 , q 3 , q 4 are source or sink terms (L 2 T - 1 ) for organic - N, NH 4 - N, NO 2 - N and NO 3 - N, respectively. D L is est imated using empirical equation from Fischer et al., (1979): (4.5) w here W is the width of the stream segment (L) ; g is the gravity acceleration (LT - 2 ) ; d is the depth of the stream segment (L) an d S ( - ) is the slope of the channel . Mulholland et al., (2008) studied the biotic uptake and denitrification using nitrogen stable isotope tracer experiments across 72 streams and 8 regions including 9 rivers in Michigan. They obtained one regression equa tion depicting the relationship between denitrification 122 rate, river depth a nd NO 3 - concentration. In order to incorporate the effects of stream depth and nitrate concentration on the denitrification rate, we use the regression equation generated from a pow er law relationship by Mulholland et al., (2008): (4.6) where V f (L T - 1 ) is the N uptake velocity , - 0.493 and - 2.975 are coeffi cients estimated from regression analysis. The denitrification rate (T - 1 ) can be expressed as: (4. 7 ) where h is the river depth (L), K r5 is the denitrification rate (T - 1 ) as in equation 4.4. However, this experimental study is taken under the stream temperature conditions within the range of 12~ 25 for ge nerating this regression equation (Mulholland et al., 2008). A temperature dependent modifie r will be used and evaluated to adjust the denitrification rate. We use typical literature values (Chapra, 2008) as the default values of all other reaction rates, see Table 4.3 , and the rates will be modified and evaluated with a temperature dependent fac tor in a later section. The denitrification rate will also be evaluated in the uncertainty analysis section. Table 4 . 3 Default reaction rates used in the river transport module Reaction Rates Values K r1 (mg l - 1 day - 1 ) 0.10 K r2 (day - 1 ) 0.05 K r3 (day - 1 ) 0.02 K r4 (day - 1 ) 0.01 K r5 (day - 1 ) Mulholland et al., 2008 123 The initial nitrate concentrations were set as 1.0 0 m g/ L for all t he river reaches in KRW; while the initial concentrations of NO 2 - N , NH 4 - N and organic - N are set as 0.1 mg/L uniformly in the river network . For the WBW, the initial concentrations of n itrate are set as 5 mg/L which are the average observed nitra te output concentrations at the catchment outlet. For the boundary condition, the N eumann boundary ( zero concentration gradient ) condition is specified at the external ends of each river reach for both watersheds. Additionally, internal boundary condition need to be specified to guarantee the concentration continuity at the river junctio n segments. The mixing model concept is adopted to calculate the concentrations at river junction segments: (4. 8 ) where (M L - 3 ) and ( M L - 3 ) are the stream nitrogen concentrations at the upstream and downstream junction segments at the end of time step n ; (m 3 s - 1 ) and ( L 3 T - 1 ) are the up stream and downstream flow rate s at the junction segments; ( M L - 3 ) is the updated junction segment concentration at the end of time step n which will be used to update and at the junction segment. 4.2.5. Nitrogen Rea ction and Transport in the overland flow Since ammonia is volatile and NO 2 - N is unstable i n surface runoff, NO 3 - N and organic nitrogen are the only nitrogen specie s simulated in surface runoff. When surface ponding accumulates, solute transfer will happen between the top layer soil solution and the runoff water. Current theories describing the transfer process between soil solution and the surface runoff can be generalized into two categories: 1) the boundary layer approach and 2) the mixing layer approach (Shi et al., 2011 ; Rumynin, 2015) . The boundary layer approach 124 assumes the transfer of solut es between the topsoil layer and runoff water is diffusion driven and dependent on the concentration gradient. The mixing layer approach assumes that an active thin layer exists at the interface of runoff water and topsoil water, in which the solute mixing and transport processes happen. The depth of the thin layer, is a parameter dependent on the soil properties and rainfall intensity which may need experimental data to calibrate (Tong et al., 2010) . Figure 4 . 4 Sketches of Boundary Layer approach . In this study, the boundary layer approach (Figure.4. 4 ) is proposed to be applied for the purpose of reducing the number of parameters. Based on the boundary layer approach, the mass transfer fl ux J d can be calcu lated as: ( 4. 9 ) where, C s (M L - 3 ) is the top layer soil nitrate concentration; C o (M L - 3 ) is the cross - sectional averaged surface runoff solute concentration, k e (L s - 1 ) is th e mass transfer coefficient between top layer soil solution and surface runoff. Thereafter, the overland solute transport equation can be expressed as two dimensional ad vection - dispersion equation (Deng et al., 2005 ) : 125 (4.10) where k is the de cay rate [T - 1 ]; i i s the infiltration rate [L T - 1 ]; D x and D y are di spersion coefficients in the x - and y - direction s [L 2 T - 1 ]. The transport of n itrogen in overland flow is usually dominated by advection under aerobic condition s , the nitrate nitrogen and organic nitrogen co uld be c onsidered as conservative when they are transported overland . Therefore, the nitrate nitrogen and organic nitrogen are the only two species tracked in overland flow, the dispersion coefficient and the decay coefficient k are assumed as zero when it is applied in this re search ( Creed and Beall, 2009 ; Rumynin, 2015 ) . However, w etlan ds, ponds and lakes in the landscape make a great contribution to the nitrogen transport and transformation processes (Cheng and Basu, 2017). They serve as important nitrog en sinks that regulate the nitrogen retention processes over its transport across th e landscape. The transient storage concept is adopted in this study to account for nitrogen retention and denitrification processes in wetlands and ponds, which can be expr essed as ( Cheng and Basu, 2017) : (4.1 1 ) (4.1 2 ) w here Q ( L 3 T - 1 ) is steady state flow rate en te ring or leaving the water column, Q gw ( L 3 T - 1 ) is the groundwater inflow or outflow rate, C o ( M L - 3 ) is the Nitrogen concentration in the inflow of surface runoff , C w (M L - 3 ) is the n itrogen concentration in the lowland storage (wetlands and ponds) and the outflow of surface runoff , K den (T - 1 ) is th e den itrification rate in the sediment, C r ( M L - 3 ) is the n itrogen concentration in the sediment water, A is the 126 wetted contact area (L 2 ), d r i s the effective reactive depth ( L ), and ( T - 1 ) is the mass exchange rate coefficient. The values of Q , Q gw , V w and A are calculated in PAWS+CLM and its lowland storage module (Shen et al., 2013). d r is used as 30 mm based on the estimation of active denitrification zone depth ( Cheng and Basu, 2017 ; H arvey et al., 2013 ). The mass exchange rate coefficient is estimate d as : (4.1 3 ) where D is the dispersion coefficient at the streambed active layer which can be estimated as 5×10 - 6 cm 2 /s. T he denitrification rate for nitrate is estimated using empirical relationship generated by Cheng and Basu (2017): (4.1 4 ) ater residence time which can be calculated as V w / Q . Similarly, the denitrification rate for the total nitrogen is estimated as: (4.1 5 ) In this study, the total nitrogen i s assumed as the sum of nitrate and organic nitrogen. Equation 4.11 and Equation 4.12 are solved using an implicit finite difference method (Appendix A) . At each time step, the nitrate concentration is first solved, followed by the total nitrogen concentra tion. The organic nitrogen concentration is obtained by subtracting the nitrate concentration from the total nitrogen concentration. 4.2.6. Nitrogen reaction and Transport in the vadose zone The Century model (Parton et al.,1988) incorporated in CLM 4.5 enables user defined fertilizer applications in top soil layers rather than using default fertilizer fluxes in CLM 4.0 ( Oleson et al., 2010; Oleson et al.,201 3), which provide no fl exibility to simulate 127 anthropogenic activities. In this se ct ion, we briefly introdu ced the nitrogen processes in the vadose zon e of the Century - based model in CLM 4.5 which is now coupled with PAWS. 4.2.6.1 The internal nitrogen cycle Here, we briefly describe the internal nitrogen dynamics of plant - litter - soil system of Century based model empl oyed in this study. CLM is fully prognostic with respect to all carbon and nitrogen state variables in the vegetation, litter, and soil organic matter domains (Thornton and Zimmermann, 2005) . The processes of vegetation growth and litterfall are al so prognostic, which are simulated in the phenology module of CLM for each plant fu nctional type (PFT). Leaf, live stem, dead stem, live coarse root, dead coarse root, and fine root pools are tracked with different pools, see Figure 4. 5 . Two additional s torage pools, one short term storage pool and one long term storage pool, are assoc iated with these pools. The newly assimilated carbon from photosynthesis and nitrogen from available mineral nitrogen are firstly allocated to meet the maintenance respirat ion costs. The vegetation tissue mass, nitrogen concentrations and temperature are three major factors that determine t he total carbon and nitrogen demand for maintenance respiration. The extra carbon flux coming from photosynthesis after supplying the mai ntenance respiration can be applied to support new vegetation growth. Potential car bon supply to new growth is determined based on the C : N ratios for each tissue type and related parameters. CLM uses a synthesis of microcosm decomposition studies based o n radio - labeled substrates to parameterize the pools sizes and reaction rates. More over, CLM uses rates modifiers to account for the soil moisture and temperature limitations on the reaction rates. 128 Figure 4 . 5 Vegetation flux and pools structure of CLM 4.5 (adapted from Oleson et al., 2 013) . 4.2.6.2 The external nitrogen cycle CLM directly couples the in terna l cycling of nitrogen within the plant litter soil organic matter system with the external nitrogen processes and important sources and sinks, i.e. the atmospheric deposition, denitrific ation, nitrification, plant uptake etc., as shown in Figure 4. 3 . Here, we primarily describe the external nitrogen processes used in the Century ni trogen model incorporated in CLM, which includes divided NH 4 - N and NO 3 - N pools. A single variable is used in CLM to represent the total atmospheric nitrogen deposition which may include a combination of wet and dry deposition of different nitrogen species 129 (Oleson et al., 2013). In the Century - based model, it is assumed that all the d eposited nitrogen enters int o the mineral nitrogen pool. Meanwhile, the biological nitrogen fixation is explicitly simulated in CLM by assuming it is a function of net primary production. Similar with atmospheric nitrogen deposition, the biological fixed nitrogen also enters into th e mineral nitrogen pool. The Century based nitrification and denitrification processes (Parton et al., 1996, 2001; del Grosso et al., 2000) as incorporated in CLM 4.5 are used in this study. Based on this approach, the nitrific ation rate of NH 4 + to NO 3 - i s modified by several environmental conditions, including soil temperature, soil moisture, and pH: (4.1 6 ) where R nitra,p is the potential nitrification rate for the p th PFT, k nitra denotes the max imum nitrification rate which is assumed as 0.1 day - 1 following Parton et al., (2001), f (H 2 O), f (T) and f (pH) are nitrification rate modifiers relative to soil moisture, soil temper ature and soil pH value, respectively. The soil moisture modifier is calcul ated as: (4.1 7 ) where is the soil water potential as level j , is a lower limit for soil water potential which is set as - 10 MPa , is the saturated soil water potential which is calculated using the multivariate regression mo del in CLM. The temperature modifier in CLM is calculated using a Q 10 approach: ( 4.1 8 ) w here T soil,j is the soil temperature at level j , T ref is the reference temperature which is 25 130 in CLM, Q 10 equals 1.5 which is a constant number. Since pH is not solved in CLM, a constant pH val ue of 6.5 is used in CLM following Parton et al., (1996). The potential denitrification rate is assumed to happen only in anoxic conditions, and is limited by the nitrate concentration and carbon assumption rates: (4.1 9 ) where R denitr,p denotes the potential denitrification rate; f ([NO3]) and f ( decomp ) are the nitrate and carbon related rate modifiers; frac anox is the anoxic fraction of the soil which is calculated using the anoxic microsite formulation from Arah and Vinten (1995) : (4.20) where a , , , , and are constant coefficients, R denotes the typical pore space radius of soil moisture , V denotes the oxygen consumption rate, C denote s the oxygen concentration, is the pore space filled with water, is the pore space filled with air and is the ratio of oxygen diffusivity in water to oxygen diffusivity of in air. All these parameters are calculated in CLM to describe the fractions of anoxic pore space in the soil column. The fertilizer /manure nitrogen is applied in the top soil layers . I t was assumed the fertilizer was fully mixed in the active mixing zone (i.e. the top 10 mm soil layer) which is consistent with the CLM - century bas ed model setting. Nitrogen leaches into groundwater when excessive soil mineral nitrogen remains after plant up take, immobilization, and denitrification. The amount of leached nitrogen is dependent on the concentration of dissolved inorganic (mineral) nitr ogen in the soil solution. PAWS+CLM follows the leaching scheme of CLM in which the leached N species include NO 3 - N and Organic N; NH 4 - N is assumed to be adsorbed onto mineral surfaces and unaffected by leaching 131 (Oleson et al., 2013). The leached nitrogen enters unconfined groundwater layer via a groundwater nitrogen source term. 4.2.7. Nitrogen reaction and t ransport in the groundwater domain The t raced nitrogen species in groundwater are consistent with species in river domain which contain organic - N, NO 3 - N, NO 2 - N and NH 4 - N. General two - dimensional solute reactive - transport equations involving multiple n itr ogen species in the fully saturated groundwater domain are proposed to be adopted as follows ( Zheng and Bennett, 2002 ; Le e, et.al. 2009): (4.21) (4.22) (4.23) (4. 2 4 ) where R 1 , R 2 , R 3 , R 4 are retardation factors (dimensionless); q x , q y are Darcy velocities (LT - 1 ); D xx , D xy , D yy , D yx are dispersion coefficient tensors (L 2 T - 1 ); is the 132 porosity(dimensionless); K g1 , K g2 , K g3 , K g4 are reaction rates (T - 1 ); q s1 , q s2 , q s3 , q s4 are source or sink terms . For the degradation rate of organic nitrogen, i.e. K g4 , we use the same constant number of 0.05 (day - 1 ) as we used in river. The Mo nod reaction kinetics (Molz et al., 1986; Lu et al., 2009 ) are used to estimate reaction rates following the method of Lee et al., ( 2006) involved in equation: ( 4. 2 5 ) ( 4 .2 6 ) ( 4.2 7 ) where (T - 1 ) is the maximum substrate utilization rate for nitrification of ammonium - nitrogen to nitrite nitrogen, (T - 1 ) i s the maximum substrate utilization rate for nitrification of nitrite - nitrogen to nitrate - nitrogen, X 1 , X 2 and X 3 are the concentrations of autotrophic ammonia - oxidizing, nitrite - oxidizing biomass and heterotrophic biomass , respectively; , , , , are the half saturation constants for NH 4, NO 2 , O 2 , CH 2 O and NO 3 , respectively; k b1 , k b 2 , k b3 and are the inhibition constants of the ammonia - oxidizing biomass, nitrite - oxidizing biomass, heterotrophic biomass and oxygen, respectively. All those parameters values are obtained from the Vasse Research Station simulation by Lee, et.al. ( 2009 ), see Table 4.4. How e ver, since the concentrations of oxygen ([ O 2 ]), BOC ( [CH 2 O] ) and biomass ( X 1 , X 2 , X 3 ) are not solved in this study, we 133 use constant values for these concentrations which are average steady state values from Lee, et.al. ( 2009 ). Table 4 . 4 Monod kinetics reaction parameters Parameter Value Parameter Value (day - 1 ) 1 k b1 (mg/l) 1 (day - 1 ) 3.5 k b2 (mg/l) 1 (day - 1 ) 10 k b3 (mg/l) 0.5 (mg/l) 0.1 X 1 (mg/l) 0.1 (mg/l) 0.3 X 2 (mg/l) 0.002 (mg/l) 0.5 X 3 (mg/l) 5 (mg/l) 0.77 [O 2 ] (mg/l) 1 (mg/l) 6 [CH 2 O] (m g/l) 10 (mg/l) 0.01 For KRW, we use the average nitrate concentration during the simulation period (2003 ~ 2010) obtained from wellogic data archived by DEQ ( State of Michigan, 2016 ) as the initial concentration for nitrat e. For multiple wells within the same simulation grid, the average value of these records is used. For WBW, h owever , a constant value of 5 mg/L is used as the initial nitrate concentration for groundw ater due to data limitation s . Typical literature conce ntrations for other N species from (Lee, et.al. 2009) are used as initial conditions for other N species for both watersheds. 4.2.8. Temperature dependent reaction rates The Arrhenius equation ( Whitehead et al., 1998; Chap ra, 2008; Wade et al., 2002) is adopted in this study to adjust the reaction rates in different hydrologic domains other than the vadose zone which is already considered in CLM: (4. 2 8 ) 134 where R 0 is the reaction ra tes at a default temperature of 20 , is a simplified temperature dependent modifier of the reaction rates. In this study, we use a value of 1.047 following the INCA model (Whitehead et al., 1998; Wade et al., 2002) to adjust the temperature dependent reaction rates. Inwood et al., (2007 ) and Li et al., (2017) studied one tributary of KRW, i.e. the Little Rabbit Creek, and found that t he den itrification rate exponentially decreas ed with streambed depth . Therefore, we use the temperature of the first layer streambed which cover s a depth of 0 ~ 10 cm to adjust the denitrification rate. O ther reaction rates in the river domain are adjusted directly using the simulated river temperature . For the reactions in groundwater, the simulated groundwater temperature is used to adjust the reaction rat es. Results and discussions The model performance was tested using simulated results a gainst observed data. Several different metrics are employed to quantitatively evaluate the model performance, including the correlation coefficient (R 2 ) , Nash - Sutcliffe efficiency coefficient (NASH) ( Nash & Sutcliffe , 1970) and root mean square error (RMSE). 4.3.1. Comparisons of h ydrology results b etween CLM4.5 and CLM4.0 Figure 4. 6 shows the stream discharge time series comparisons between simulated results using CLM 4.0 , CLM 4.5 , and observations by USGS gauges of KRW . The NASH and RMSE values of the stream discharge comparisons are tabulated in T able 4.5. The model results almost show the same performance between CLM 4.0 and CLM 4.5 . The major reason is that the biogeochemistry module is the only updated module by us ing CLM 4.5 century - based model. Other dominant hydrologic processes, such as th e stream flow, evapotranspiration and groundwater processes etc., are kept the same as previous version 135 of PAWS+CLM. Therefore, the monthly watershed average ET also preserve t he similar magnitude and trend after coupling with CLM 4.5 century - based model , see Figure 4.7 . These model performance, are acceptable to serve as the basis for our nitrogen simulations in the catchment considering all the uncertainties discussed in prev ious chapters. Figure 4 . 6 stream flow comparisons between simulated results using CLM 4.0, CLM 4.5, and observations by USGS gauges . 136 Figure 4 . 7 Monthly watershed - averaged ET comparisons between simulated results using CLM 4.0 , CLM4.5 and remotely sensed MODIS ET products . Table 4 . 5 Performance of Stream discharge compa risons between simulated results using CLM4.0, CLM 4.5 and observations by USGS gauges CLM 4.0 CLM4.5 USGS stat ion No. NASH RMSE NASH RMSE 04103010 0.61 2.04 0.6 0 2 .22 04106000 0.63 10.18 0.61 10.39 04108600 0.64 1.79 0.65 1.37 04108660 0.76 20.35 0.72 21.47 137 4.3.2. The Kalamazoo watershed nitrogen results. 4.3.2.1 C omparison of nitrate concentration s in river s Figure 4. 8 a nd Figure 4.9 show the river nit rogen concentration comparisons between simulated results and the observations. The R 2 values and RMSE values for the nine sampling locations are tabulated in Table 4. 6 , which show comparable performance with other watershed - scale nitrogen models, e.g. Lu et al ., (2017), Yang et al., (2018). The performance metrics for the site USGS 04102810 from the Storet database are not shown since only several observation points are within our simulation period. At two sampling locations of Battle Creek, the model could cap ture several nitrate peak values during the year 2006. In the year 2005, the range of spring observed nitrate concentrations are reproduced by the model whereas the summer nitrate concentrations are overestimated. At t he other seven sampling locations of the Kalamazoo River , the model could reproduce the nitrate fluctuations reasonably well, especially during the year 2006. The timings of several peaks are well captured by the model which showed the capability of the mo del to accurately predic t the timing of nitrogen flush - out. However , some of the peak values are either overestimated or underestimated by the model which could be accounted from multiple uncertaint y sources including the model input and parameter uncertai nties. Accurate knowledge of agricult ural practices is lacking, uncertainties of the timing and amount of fertilizer applications as well as their spatial distributions are a major source of model uncertainties (Yang et al., 2018). The unce rtainties asso ciated with model parameters and the model results of other hydrologic variables, e.g. , stream discharge, soil moisture could also affect model performance ( Porporato et al., 2003 ). The simulated concentrations of organic nitrogen, nitrite and ammonia all follow the range of simul ation s, 138 see Figure 4.9. The concentrations of ammonia and nitrite occupy a very small portion in the total nitrogen concentrations, with an average portion of 0.98% and 2.73%, respectively, for the three observation sites. Figure 4 . 8 Comparison between simulated and observed nitrate concentration time series in nine different sampling locations in Baas (2009) . The simulated concentrations of nitrite, ammonia and organic nitrogen are also shown. 139 - Figur e 4 . 9 Comparison between simulated and observed nitrate , nitrite, ammonia and organic nitrogen concentration time series of three observation stations in Storet database . NO 3 - sim, NO 2 - sim, organicN - sim and NH 4 - sim are the simulated concentrations of nitrate - N, nitrite - N, organic - N and ammonia - N; NO 3 - obs, NO 2 - obs, organicN - obs and NH 4 - obs are the observed concentrations of nitrate - N, nitrite - N, organic - N and ammonia - N, respectively. 140 Table 4 . 6 Performance metrics of river nitrogen concentration time series comparisons Baseline simulation After temperature correction Sampling Location R 2 RMSE R 2 RMSE Battle Creek at Emmett St. Dam 0. 70 0.32 0. 7 4 0.27 Battle Creek 0. 69 0.31 0. 72 0.26 Kalamazoo river at Battle Creek 0.67 0.26 0.74 0.24 Kalamazoo river at Comstock 0. 70 0.21 0. 71 0.21 Kalamazoo river at Caulkins Dam 0.6 8 0.25 0.6 6 0.27 Kalamazoo river at east of BC 0.62 0.34 0.6 2 0.3 5 Kalamazoo river at Augusta 0 .8 4 0.22 0.8 4 0.23 Kalamazoo river at Kalamazoo 0.62 0.32 0.64 0.29 Kalamazoo river at Galesburg 0.78 0.19 0.7 8 0. 20 USGS 04105707 (Nitrate) 0.76 0.24 0.77 0.23 USGS 04105707 (Nitrite) 0.68 0.01 0.69 0.01 USGS 04105707 (Ammonia) 0.79 0.01 0.77 0.01 U SGS 04105707 (Organic N) 0. 59 0.1 9 0. 60 0.1 8 USGS 04108660 (Nitrate) 0.81 0.24 0.81 0.23 USGS 04108660 (Nitrite) 0.78 0.01 0.80 0.01 USGS 04108660 (Ammonia) 0.73 0.01 0.72 0.01 USGS 04108660 (Organic N) 0.66 0.18 0.65 0.19 Table 4 . 7 Pair - wise linear correlation coefficients for nitrate leaching values with different Land use/land cover types and soil texture Land Use/Land Cover Soils Forest Grass Crops Urban Wetland Sand Clay Organic Matter - 0. 64 0. 2 5 0.84 - 0. 82 - 0. 5 7 0.5 3 - 0. 34 0. 11 p 0 .00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 is the Correlation Coefficient. p < 0.05 indicates statistical significance. 141 4.3.2.2 N itrate leaching Nitrate leaching happens when excessive nitrogen percolates into t he unsaturated groundwater layer after plant uptake and denitrification. As shown in Figure 4. 10 , the spatially distributed annual average nitrate leaching rate was comput ed during the simulation period. Several low nitrate leaching z one s are located in t he northwestern part and t wo south - central areas of KRW. The northwestern part of KRW is primarily forest occupied and many open lakes and reservoirs are in this area. Low fertilizer usage is applied in this region such that the nitrate leakage is expected to be low. Two south - cent ral areas of KR are urban areas, which are also expected to have low fertilization. We performed a linear correlation analysis on the spatially distributed nitrate leaching results against the LULC types and the soil texture , and the results are summarized in Table 4. 7 . In this table, the pair - wise linear correlation coefficient is expressed with , and p denotes the probability for testing the null hypothesis of no correlation against the alternative hypothesis that there is a significant correlation. Smal l p (i.e., < 0.05) indicates the correlation is important . As expected, crops are highly po sitively correlated with nitrate leaching values 84 ) due to the heavy agriculture fertilizer use. U - 0. 82 ) and forest areas ( - 0. 64) are negatively correlated to nitrate leaching values, which is consistent with previous discussio n . In addition to low fertilizer use, the pavements in urban areas may prevent the water from percolating into groundwater which reduces the nitrate leaching. However, point sources and surface runoff produced in the urban areas may accelerate the transpor t process of nitrogen into the river (Thomas et al., 2016). The sand soil texture tends to show a h igh positive correlation with nitrate leaching values 53 ) whereas clay soil tends to show relatively high negative correlatio n - 0. 34 ) with nitrate leaching values, wh ich 142 indicate s that the nitrate leaching is also closely related to soil infiltration capacity . The soil texture and infiltration cap acity influence the soil water movement which are also highly correlated with ET and groundwater recharge (Qiu et al., 2019 b). The sandier the soil, the more potential for the nitrate water percolates into groundwater (Lin et al., 2001; Wick et al., 2012) . Figure 4 . 10 The annual average nitrate leaching map . 143 4.3.2.3 C omparison of nitrate concentration s in groundwater The comparison between annual average spatial map s of simulated groundwater nitrate concentrat ion and observed nitrate concentration from the DEQ wellogic data is shown in Figure 4. 11 . The simulated annual average groundwater nitrate concentration map preserves the major pattern of the observed groundwater nitrate concentration map. Several high n itrate concentration areas are located in the upper ( south east) portions of the watershed which are dominated by agriculture land use. The blank pixels in the observed map are due to a lack of observations. In the middle of KRW, the simulated groundwater concentrations show relatively high concentrations with yellow and red color while only several pixels of obser vations are available . Low nitrate concentrations are shown in the northwestern part of KRW for both observations and simulations, which is consi stent with the nitrate leaching map. For the two urban areas, the observations show some high concentrations (r ed spots) while the simulations do not. This discrepancy might result from anthropogenic nitrogen point sources which are not explicitly included in the model. The aquifer material of Kalamazoo is mainly glacial deposits which are composed of outwash sand and gravel with a relatively high average hydraulic conductivity of approximately 18 m/day (Wesley, 2015, Qiu et al., 2019 b). The nitrogen tran sport processes are therefore relatively high due to active groundwater movement. This can be reflected from si mulated results: despite the groundwater of forest and urban areas receive extremely small amount of nitrate leaching, their nitrate concentratio ns could remain within the range of 0.5~4 mg/l. 144 Figure 4 . 11 Comparison of spatially distributed groundwater concentrations between observations and model results . 145 4.3.2.4 Effects of temperature on the riverin e nitrogen export The time series comparisons between observations and simulated n itr ogen concentrations for baseline simulation and the simulation with temperature - dependent reaction rates are shown in Fig ure 4. 12 and Figure 4.13 . The performance metrics of the comparison after applying the temperature - dependent reaction rates are als o tabulated in Table 4. 6 . It can be seen from T able 4. 6 that after applying the temperature dependent reaction rates, the performance of nitrate concentrations at the two sam pling sites of Battle Creek was improved with R 2 values increased and RMSE values decreased . However, in the main river, performance improvement s are not consistently shown for all the sampling sites. The obvious change we can tell is during the summers, as shown in the zoomed - in Figures 4. 1 2 (a) and 4 .1 2 (b), that the simulate d nitrate concentrations decreased after applying the temperature dependent reaction rates. In contrast, the nitrate concentrations i ncreased during winters (zoomed - in figures not sh own). It is consistent with previous literature and discussions that warmer temperature s will enhance the nitrogen removal processes in the catchment and riverine network (Schaefer and Alber, 2007 ; Miller et al., 2016) . The model performance for the organic nitrogen, nitrite and ammonia are not obviously changed after applying the temperature dependent reactions rates at the Storet stations. One major reason for this could be the large time gap of the observations . 146 Figure 4 . 12 T ime series comparisons between obser vations in Baas (2009) and simulated nitrate concentrations for baseline simulation and simulation with temperature dependent reaction rates ; ( a ) zoomed i n figure of site Battle Creek at Emmett St. Dam (b) zoomed in figure of site Kalamazoo River at Augusta . NO 3 base is the baseline simulation, NO 3 temp is the simulation with temperature dependent reaction rates. 147 Figure 4 . 13 T ime series comparisons between obser vations in S toret databse and simulated nitrate concentrations for baseline simulation and simulation with temperature dependent reaction rates . NO 3 - sim (NO 3 - sim - t), NO 2 - sim (NO 2 - sim - t), or ganicN - si m (organicN - sim - t) and NH 4 - sim (NH 4 - sim - t) are the simulated (simulated with temperature corrected reaction rates) concentrations of nitrate - N, nitrite - N, organic - N and ammonia - N; NO 3 - obs, NO 2 - obs, organicN - obs and NH 4 - obs are the observed concentrations o f nitrate - N, nitrite - N, organic - N and ammonia - N, respectively. 148 4.3.2.5 U ncertaint y analysis T o quantitatively estimate the model sensitivity result ing from the parameter uncertaint y and nitrogen sources , two important parameters, i.e. the denitrifi cation and nitri fication rate s in rivers , as well as the anthropogenic nitrogen sources were perturbed by relative to the baseline simulation . Figure 4. 1 4 , Figure 4.15 and Figure 4.16 show the results of the uncertainty analysis , in which the up per and lower limits of shaded area s correspond to the maximum and minimum concentration results due to change of parameter or anthropogenic nitrogen source s . The observed nitrate concentrations are mostly covered by th e uncertainty band s of the simulated nitrate concentrations for all the sampling sites. Although the band s cover a range of around 0. 1 - 0.5 mg/l, the fluctuations of the band s preserve the similar trend of the baseline simulation, which indicates the consis tency of the model behavior. We calcu late the RMSE values between baseline simulations and the simulations after applying the perturbations of parameters or anthropogenic sources. Larger RMSE values means larger deviation s between baseline simulations and perturb ed simulations, which indicat es larger sensitivity. As shown in Table 4.8, t he model is most sensitive to uncertainties resulting from nitrogen sources, followed by the denitrification rate and the nitrification rate in rivers. In addition to uncer tainty in the magnitude of sources, t here are also uncertainties associated with the location and timing of manure and fertilizer application which are not known precisely. These additional uncertainties are probably responsible for some observed peaks in data that were not captured by the mo del. 149 Table 4 . 8 RMSE val u es of nitrate concentrations between baseline simulations and perturbated simulations (The RMSE value is calculated as the mean value of +10% and 10% for each case) Anthropogenic Sources Denitrification rate Nitrification rate Sampling Location RMSE RMSE RMSE Battle Creek at Emmett St. Dam 0. 49 0. 4 2 0. 28 Battle Creek 0. 52 0. 46 0.3 3 Kalamazoo river at Battle Creek 0. 3 6 0. 33 0.2 8 Kalamazoo river at Comstock 0. 18 0. 15 0. 14 Kalamazoo river at Caulkins Dam 0.2 1 0.2 1 0.2 0 Kalamazoo river at east of BC 0.3 0 0. 27 0. 24 Kalamazoo river at Augusta 0. 39 0.2 6 0.2 1 Kalamazoo river at Kalamazoo 0. 37 0. 27 0. 27 Kalamazoo river at Galesburg 0. 26 0.1 7 0.1 6 USGS 04105707 (Nitrate) 0.25 0.23 0.22 USGS 04108660 (Nitrate) 0.22 0.20 0.20 150 Figure 4 . 1 4 90 % confidence bands of simulated nitrate concentration s result ing from uncertainties of anthropogenic nitrogen source s (fertilizer , manure and point sources ) . N O 3 - sim is the nitrate concentation of baseline simulation, NO 3 - obs is the observed nitrate concentration, uncertainty band den otes the 90% confidence band result ing from uncertainties of anthropogenic nitrogen sources. 151 Figure 4 . 15 90 % confidence bands of simulated nitrate concentrations result ing from denitrification rate in r ivers. N O 3 - sim is the nitrate concentation of baseline simulation, NO 3 - obs is the observed nitrate concentration, uncertainty band denotes the 90% confidence interval result ing from denitrification rate in rivers . 152 Figure 4 . 16 90 % confidence bands of simulated nitrate co ncentrations result ing from nitrification rate in rivers. N O 3 - sim is the nitrate concentation of baseline simulation, NO 3 - obs is the observed nitrate concentration, uncertainty band d enotes the 90% confidence interval result ing from nitrification rate in rivers . 153 4.3.3. The Wood Brook watershed nitrogen results 4.3.3.1 The nitrate concentration comparison in river Figure 4.1 7 shows the comparison between simulated and observed stream nitrate concen trations at the catchment outlet of WBW. T he simul ated concentrations of nitrite, ammonia and organic nitrogen are also shown although no observed data are available . The R 2 and RMSE of the nitrate concentration comparison are 0.69, and 1. 09 , respectively , which shows acceptable performance. While t he ti ming of nitrate peaks could be captured by the model before September 2016, the nitrate peaks in fall and winter are not well reproduced by the model . The simulated peaks are consistent with the timing of f ertilizer /manure applications, the fertilizer application dates are during spring and summer of 2016 in WBW (Table 4.2). However, the observations still show nitrate peaks during fall and winter . Unlike the nitrate sampling in the Kalamazoo watershed, whi ch was bas ed on wet chemistry, the data from the WBW were obtained using a nitrate sensor. Some of t he nitrite concentrations obtained from a similar sensor were negative and therefore not used to test the model. Large uncertainties are associated with the sensors. In addition, some of the mismatch might be primarily due to the nitrogen retention a nd storage processes in rivers and/ or land surface which could strongly influence the timing of nitrogen flush - out ( Alexander et al., 2009; Covino et al., 2010; Y e et al., 2012) . Even though the current model considers the denitrification processes in wetlands and lowland, the rainfall related nitrogen activation processes are not represented. It is reported in the same catchment that the nitrogen flush - out could be strongl y affected by the intensity of the storm events and time - d ependent nitrogen source zone activations ( Blaen et al., 2017 ) . 154 Figure 4 . 17 Comparison between simulated and observed stream nitrate c oncentration time series near the catchment outlet of WBW . The simulated concentrations of nitrite, ammonia and organic nitrogen are also shown. 155 4.3.4. Nitrogen budget Table 4. 9 lists the simulated nitrog en budget results for the two tested watersheds. The tota l nitrogen input from anthropogenic sources is counted as 100%. The groundwater nitrogen contribution to rive rs for the KRW ( 22.91 %) is almost the same as that of the WBW ( 22.74 %). T he riverine nitrogen output ratio of WBW ( 39.45 %) is much higher than that of KRW ( 27.26 %). This could also be seen from the river nitrate concentrations that the average outlet nitrate concentration of Wood Brook (~ 5 mg/L) is much higher than the nitrate co ncentrations in the sampled two rivers of KRW (0.5 mg/L - 2.5mg/L). The net removal nitrogen amount in river network of KRW occupies 6.37 % of the total anthropogenic nitro gen input of the watershed which is almost double the removal ratio of WBW ( 3.24 %). T he mean annual air temperature of WB W during the simulation period is a round 10 .5 o C which is 0.43 o C higher than the mean annual temperature of KR W ( ~ 10.07 o C ). Based on previous discussions, high er watershed temperature s could enhance the nitrogen rem oval processes in the catchment, where as the warmer WBW has a lower nit rogen removal ratio in the river network . The primary reason could be t he more complex river network and longer water - nitrogen retention time of KRW , which could increase the duration of denitrification (Wollheim et al., 2006 ; Helton et al., 2018) . 156 Table 4 . 9 Nitrogen bud get results Kalamazoo River watershed Wood Brook watershed N amount (kg/year) Percentage N amount (kg/year) Percentage Total anthropogenic N input 2.37×10 7 100.00% 29188.47 100.00% GW nitrogen contribution to river 5.43×10 6 22.91% 6639.78 22.74% O VN nitrogen contribution to river 2.54×10 6 10.71% 6421.48 22.00% Riverine export N 6.46×10 6 27.26% 11514.02 39.45% Net removal in river network 1.51×10 6 6.37% 947.24 3.24% Net removal in wetlands 1.37×10 6 5.78% 687.14 2.35% Conclusions A catchment - s cale framework was developed to simulate integrated hydrologic, temperature and nitrogen processes. This innovative modeling framework enables real - time (or near real - time) forecasts of nitrogen concentrations at different hydrologic domains. The model pe rformance was evaluated by applying the model to two different watersheds with different size and climate conditions. Extensive comparisons between modeling results and observations indicated the capability of this developed modeling framework to be effect ively utilized in simulating watershed - scale temperature and nitrogen processes in different hydrologic domains. Strongly positive correlation was identified between groundwater nitrate leaching amount with crop land use in KRW, whereas the fo rest and unba n land use showed negative correlation. The detailed spatial information of ground water nitrate concentrations could be reproduced by PAWS+CLM which are within reasonable range of observations . The use of t emperature - dependent reaction rates confirmed the positive r elationship between nitrogen removal rate in the 157 catchment and temperature (Schaefer and Alber, 2007 ; Miller et al., 2016) . Nitrogen budget results indicated that more co m plex river network and longer water retention time in KRW increased the nitrogen removal ratio in the rivers compared with WBW (Wollheim et al., 2006 ; H e lton et al., 2018) . Moreover, results of model uncertainty analysis showed the sensitivity of the model performance in response to the parameter and nitrogen source uncertainties, which confirmed the robustness of the model. 158 Chapter 5. Conclusion This dissertation research developed a process - based, integrated hydrology - temperature - nitrogen model that can be used to make real - time forecasts of river temperature and nitrogen concentrations at the watershed scale. The developed model frame work was ap plied to two watersheds with different sizes and climate conditions and to explore the interactions of multiple hydrological and biogeochemical processes. Extensive comparisons between simulated results and observations indicated good model perf ormance. De tailed representations of important hydrologic, thermal and nitrogen processes enabled the model to track and quantify fluxes of water, heat and nitrogen species across hydrologic domains at different scales. The model framework improved the und erstanding of individual processes and enabled the identification of controlling processes using budgets of water, heat and nitrogen. Detailed descriptions of the interactions between surface - subsurface and at stream bed layers could be used to gain furthe r insights of the functions of hydrological and biogeochemical processes at the hyporheic zone. The streambed temperature stratifications and fluctuations are dependent on the vertical exchange water fluxes between groundwater - surface water domains, whic h can be us ed to investigate streambed processes and parameters. Potential deforestation in the KRW underscored the importance of protecting the cold fish habitats, especially at source water zone with relatively small stream size. The simulated results of riverine n itrogen output in KRW are found to be more sensitive to anthropogenic nitrogen sources than the denitrification and nitrification rates in the river network. The nitrate leaching amount was strongly positively correlated with crop land use in KR W, whereas negative correlation was found for the forest and urban land use. In addition, warmer temperatures 159 facilitated the natural nitrogen removal processes in the stream network. These findings underscore the importance of protecting aquatic ecosystem s and they could inform management of agricultural practices. It would also be interesting to project the simulations to future scenarios with potential climate change and land use change. 160 A PPENDIX 161 A PPENDIX . Solution of the transient storage equat ions for wetlands The transient storage equations in Chapter 4 (equations 4.11 and 4.12) can be solved using an implicit finite difference method as explained in following procedures . Equation 4.12 can be discretized using an implicit Euler method to deriv e an expression for : (A.1) Collecting term for , we get: (A.2) Such that: (A.3) Simplifying equation A.3 by substituting the constants with P 1 and P 2 , can be expressed as: (A.4) where , . Equation A.1 can be discretized as: (A.5) Substituting in equation A.5 with equation A.4, we can obtain: 162 (A.6) Re ar ranging the items in equation A.6 : (A.7) Equation A.7 can be solved implicitly by arranging it in the tri - diagonal form: (A.8) where A i = 0, , C i = 0, . 163 REFERENCES 164 REFERENCES Alam, M.J., Dutta, D., 2016. A Sub - Catchment Based Approach for Modelling Nutrient Dynamics and Transport at a River Basin Scale. Water Resour Manage 30, 5455 5478. https://doi.org/10.1007/s11269 - 01 6 - 1500 - x Allan, J.D., Castillo, M.M., 2007. Stream Ecology: Structure and Function of Running Waters, second ed. Springer, Dordrecht, The Netherlands. Alexan der, R.B., Smith, R.A., Schwarz, G.E., 2000. Effect of stream channel size on the delivery of nit rogen to the Gulf of Mexico. Nature 403, 758. https://doi.org/10.1038/35001562 Alexander, R. B., Böhlke, J. K ., Boyer, E. W., David, M. B., Harvey, J. W., Mulholland, P. J., Seitzinger, S. P., Tobias, C. R., Tonitto, C. and Wollheim, W. M . , 2009. Dynamic modeling of nitrogen losses in river networks unravels the coupled effects of hydrological and biogeochemical processes, Biogeochemistry, 93(1), 91 116, doi: 10.1007/s10533 - 008 - 9274 - 8 . Anderson, D.M., Glibert, P.M., Burkholder, J.M., 2002. Harmful algal blooms and eutrophication: Nutrient sources, compositi on, and consequences. Estuaries 25, 704 726. https://doi.org/10.1007/BF02804901 Anderson, M. P., 2005. Heat as a Ground Water Tracer, Ground Water, 43(6), 951 968. doi:10.1111/j.1745 - 6584.2005.00052.x An derson, E. J. and Phanikumar, M. S., 2011. Surface storage dynamics in large rivers: Comparing three - dimensional particle transport, one - dimensional fractional derivative, and multirate transient storage models, Water Resour. Res., 47(9), W09511, doi: 10.1029/2010WR010228 Anderton, S., Latron, J. and Gallart, F.: Sensitivity analysis and multi - response, multi - criteria evaluation of a physically based distributed model, 2002. Hydrol. Process., 16(2), 333 3 53, doi: 10.1002/hyp.336 . Anibas, C. , Fleckenstein, J.H., Volze, N., Buis, K., Verhoeven, R., Meire, P., Batelaan, O., 2009. Transient or steady - state? Using vertical temperature profiles to quantify groundwa ter surface water exchange. Hydrol. Process. 23, 2165 2177. doi:10.1002/hyp.728 Ani bas, C., Schneidewind, U., Vandersteen, G., Joris, I., Seuntjens, P., Batelaan, O., 2016. From streambed temperature measurements to spatial - temporal flux quantification: u sing 165 the LPML method to study groundwater surface water interaction. Hydrol. Process . 30, 203 216. doi:10.1002/hyp.10588 Arah, J. and Vinten, A., 1995. Simplified models of anoxia and denitrification in aggregated and simple - structured soils. European Jour nal of Soil Science 46(4):507 - 517. Arnold, J.G., Srinivasan, R., Muttiah, R.S., Wil liams, J.R., 1998. Large Area Hydrologic Modeling and Assessment Part I: Model Development1. JAWRA Journal of the American Water Resources Association 34, 73 89. https://doi.org/10.1111/j.1752 - 1688.1998.tb05961.x Baranov V., Lewandowski, J. & Krause, S. , 2016 . Bioturbation enhances CO2 emission from temperate lakes in a warming world, Biol. Lett, 201 6 12 201 60448; doi: 10.1098/rsbl.2016.0448. Baas, Dean G., 2009. Inferrring dissolved phosphorus cycling in a TMDL watershed using biogeochemistry and mixed linear models. Ph.D. dissertation, Michigan State University. Benyahya, L., Caissie, D., St - Hilaire, A., Ouarda, T.B.M.J., Bobée, B., 2007. A Review of Statistical Water Temperature Models. Canadian Water Resources Journal / Revue canadienne des ressources hydriques 32, 179 192. doi:10.4296/cwrj3203179 Beman, J.M., Arrigo, K.R., Matson, P.A., 2005. Agricultu ral runoff fuels large phytoplankton blooms in vulnerable areas of the ocean. Nature 434, 211. https://doi.org/10.1038/nature03370 Beven, K., 1993. Prophecy, reality and uncertainty in di stributed hydrol ogical modelling, Advances in Water Resources, 16(1), 41 51, doi: 10.1016/0309 - 1708(93)90028 - E Beven, K.J., 2001. Uniqueness of place and process representations in hydrological m odelling. Hydrol . Earth Syst. Sci. 4, 203 213. https://doi.org/10.5194/hess - 4 - 203 - 2000 Beven, K. and Freer, J., 200 1 . Equifinality, data assimilation, and uncertainty estimation in mechanistic modell ing of complex e nvironmental systems using the GLUE methodology, Journal of Hydrology, 249(1 4), 11 29, doi: 10.1016/S0022 - 1694(01)00421 - 8 , 2001. Beven, K., 2002. Towards an alternative blueprint for a physical ly based digitally simulated hydrologic response modelling system, Hydrol. Process., 16(2), 189 206, doi: 10.1002/hyp.343 . Bicknell, B.R., Imhoff, J.C., Kittle, J.L., Jr., Donigian, A.S., Jr., and Johanson, R .C., 1997, Hydrological Simulation Program -- Fortran, User's manual for version 11: U.S. Environmental Protection Agency, National Exposure Research Laboratory, Athens, Ga., EPA/600/R - 97/080, 755 p. 166 Bothe, H.; Ferguson, S.; Newton, W. E. Biology of the Nit rogen Cycle; Elsevier Science, 2006. Blaen, P.J., Khamis, K., Lloyd, C., Comer - Warner, S., Ciocca, F., Thomas, R.M., MacKenzie, A.R., Krause, S., 2017. High - frequency monitoring of catchment nutrient exports reveals hig hly variable storm event responses a nd dynamic source zone activation. J. Geophys. Res. Biogeosci. 122, 2017JG003904. https://doi.org/10.1002/2017JG003904 Bonan, G. B., 1996. A land surface model (LSM v ersion 1.0) for ecological, hydrolog ical, NCAR/TN - 417+STR, National Center for Atmospheric Research, Boulder, Colo. Boyer, E.W., Goodale, C.L., Jaworski, N.A., Howarth, R.W., 2002. Anthropogenic nitrogen source s and relationships to riverine nitrogen export in the northeastern U.S.A. Biogeochemistry 57 58, 137 169. https://doi.org/10.1023/A:1015709302073 Breemen, N. van, Boyer, E.W., Goodale, C.L., Jaworski, N.A., Paustian, K ., Seitzinger, S.P., Lajtha, K., May er, B., Dam, D. van, Howarth, R.W., Nadelhoffer, K.J., Eve, M., Billen, G., 2002. Where did all the nitrogen go? Fate of nitrogen inputs to large watersheds in the northeastern U.S.A. Biogeochemistry 57 58, 267 293. https://doi.org/10.1023/A:1015775225913 Birgand, F., Skaggs, R.W., Chescheir, G.M., Gilliam, J.W., 2007. Nitrogen Removal in Streams of Agricultural Catchments A Literature Review. Critical Reviews in Environmental Science and Technology 37, 381 487. https://doi.org/10.1080/10643380600966426 Brookfield, A.E., Sudicky, E.A., Park, Y. - J., Conant, B., 2009. Thermal transport modelling in a fully integrated surface/subsurface framewor k. Hydrol. Process. 23, 2150 2164. doi:10.1002/hyp.7282 Bu isson, L., Thuiller, W., Lek, S., Lim, P. and Grenouillet, G., 2008. Climate change hastens the turnover of stream fish assemblages, Global Change Biology, 14(10), 2232 2248, doi: 10.1111/j.13 65 - 2486.2008.01657.x Bur ns, D. A., Blet t, T., Haeuber, R., and Pardo, L.H., 2008. Critical loads as a policy tool for protecting ecosystems from the effects of air pollutants, Front. Ecol. Environ , 6, 156 159, doi:10.1890/070040 Caissie, D., 2006. The thermal regime of rivers: a review. Freshwater Biology 51, 1389 1406. https://doi.org/10.1111/j.1365 - 2427.2006.01597.x Caissie, 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. doi:10.1016/j.jhydrol.2007.01.008 167 Caissie, D., Kurylyk, B.L., St - Hilaire, A., El - Jabi, N., MacQuarrie, K.T.B., 2014. Streambed temperature dynamics and corresponding heat fluxes in small streams e xperiencing seasonal ice cover. Journal of Hydrology 519, 1441 1452. doi:10.1016/j.jhydrol.2014.09.034 Caissie, D., Luce , C.H., 2017. Quantifying streambed advection and conduction heat fluxes. Water Resour. Res. 53, 1595 1624. doi:10.1002/2016WR019813 C alow, P., Petts E., 1992. The Rivers Handbook, Volume 1, First edition, Wiley - Blackwell, Oxford, United Kingdom. Cao, Q. , Sun, N., Yearsley, J., Nijssen, B., Lettenmaier, D.P., 2016. Climate and land cover effects on the temperature of Puget Sound streams . Hydrol. Process. 30, 2286 2304. doi:10.1002/hyp.10784 Carpenter, S.R., Caraco, N.F., Correll, D.L., Howarth, R.W., Sha rpley, A.N., Smith, V.H., 1998. Nonpoint Pollution of Surface Waters with Phosphorus and Nitrogen. Ecological Applications 8, 559 568. https://doi.org/10.1890/10510761(1998)0 08[0559:NPOSWW]2.0.CO;2 Chapra, S.C., 2008. Surface Water - Quality Modeling. Waveland Pr Inc, Long Grove, Ill. Conley, D.J., Paerl, H. W., Howarth, R.W., Boesch, D.F., Seitzinger, S.P., Havens, K.E., Lancelot, C., Likens, G.E., 2009. Controlling Eutrophica tion: Nitrogen and Phosphorus. Science 323, 1014 1015. https://doi.org/10.1126/ science.1167755 Chen, J., Wen, J., Tian, H., Zhang, T., Yang, X., Jia, D., & Lai, X., 2018. A Study of Soil Thermal and Hydraulic Properties and Parameterizations for CLM in the SRYR, Journal of Geophysical Research: Atmospheres, 123, 8487 8499. DOI: 10.1 029/2017JD028034. Cheng, F. Y. and Basu, N. B., 2017. Biogeochemical hotspots: Role of small water bodies in landscape n utrient processing, Water Resources Research, 53(6), 5038 5056, doi: 10.1002/2016WR 020102 Covino, T., McGlynn, B. and Baker, M., 2010. Separating physical and biological nutrient retention and quantifyin g uptake kinetics from ambient to saturation in successive mountain stream reaches, Journal of Geophysical Research: Biogeosciences, 115(G4), doi: 10.1029/2009JG001263 . Cowling, E.B., Erisman, J.W., Smeu lders, S.M., Holman, S.C., Nicholson, B.M., 1998. Optimizing air quality management in Europe and North America: Justification for integrated management of both oxidized and reduced forms of nitrogen. Environmental Polluti on, Nitrogen, the Confer - N - s First International Nitrogen Conference 1998 102, 599 608. https://doi.org/10.1016/S0269 - 7491(98)80088 - 2 168 Cox, M.M., Bolte, J.P., 2007. A spatially explicit network - based model for estimating stream temperature distribution. Environmental Modelling & Software, Special section: Environmental Risk and Emergency Management 22, 502 514. doi:10.1016/j.envsoft.2006.02.011 Creed, I. F. and Beall, F. D.: Distributed topograp hic indicators for predicting nitr ogen export from headwater catchments, Water Resources Research, 45(10), doi: 10.1029/2008WR007285 , 2009. Dai, Y., Dickinson, R.E. & Wang, Y.P., 2004. A two - big - leaf mod el for canopy temperature, photosy nthesis, and stomatal conductance, J. Climate, 17, 2281 2299, doi: 10.1175/1520 - 0442(2004)017<2281:ATMFCT>2.0.CO;2 . Davison, J.H., Hwang, H. - T., Sudicky, E.A. & Lin, J.C., 2015. Coupled atmospheric, land surface, and subsurface modeling: Exploring water and energy feedbacks in three - dimensions, Advances in Water Resources, 86, 73 85. del Grosso, S.J., et al. 2000. General model for N 2 O and N 2 gas emissions from soils due to dentrification. Global Biogeochem. Cycles 14(4):1045 - 1060. Dickinson, R. E. (1983), Land surface processes and climate Surface albedos and energy balance, Adv. Geophys., 25, 305 353. Dingman, S. L., 2015. Physical Hydrology, Third Edition, Waveland Pres s, Long Grove, IL - Iturbe, I., 2003. Hydrologic controls on soil carbon and nitrogen cycles. II. A case study. Advances in Water Resour ces 26, 59 70. https://doi.org/10.1016/S0309 - 1708(02)00095 - 7 Dugdale, S.J., Bergeron, N.E., St - Hilaire, A., 2013. Temporal variability of thermal refuges and water temperature patterns in an At lantic salmon river. Remote Sensing of Environment 136, 358 373. doi:1 0.1016/j.rse.2013.05.018 Dugdale S.J, Hannah D.M. & Malcolm I.A., 2017. River temperature modelling: a review of deterministic approaches and future directions, Earth Science Reviews, 1 75, 97 - 113 doi: 10.1016/j.earscirev.2017.10.00 Dugdale S.J., Malcolm I.A., Kantola K. & Hannah D.M., 2018. Stream temperatures under contrasting riparian forest cover: an examination of thermal dynamics, microclimates and heat exchanges, Science of the T otal Environment, 610, 1375 - 1389, doi: 10.1016/j.scitotenv.2017.08.198 Ebersole, J. L., Liss, W. J. and Frissell, C. A., 2001. Relationship between stream temperature, thermal refugia and rainbow trout Oncorhynchus mykiss abundance in arid - 169 land streams in the northwestern United States, Ecology of Freshwater Fish, 10(1), 1 10, doi: 10.1034/j.1600 - 0633.2001.100101.x Environmental Laboratory, 1985. CE - QUAL - R1: A Numerical One - Dimensional Model of - 82 - 1, US Army Engineer Waterways Experiment Station, Vicksburg. Epelde, A.M., Antiguedad, I., Brito, D., Jauch, E., Neves, R., Garneau, C., Sauvage, S., Sánchez - Pérez, J.M., 2016. Different modelling approaches to evaluate nitrogen transport and turnover at t he watershed scale. Journal of Hydrology 539, 478 494. https://doi.org/10.1016/j.jhydrol.2016.05.066 Fang, X., Luo, S., Lyu, S., Chen, B., Zhang, Y., Ma, D. and Chang, Y., 2016. A Simulation and Validation of CLM during Freeze - Thaw on the Tibetan Plateau, Adv. Meteorology, Article ID: 9476098, DOI: 10 .1155/2016/9476098. Ferrari, M. R., Miller, J. R. and Russell, G. L., 2007. Modeling changes in summer temperature of the Fraser River during the next century, Journal of Hydrology, 342(3), 336 346, doi: 10.1016/j.jhydrol.2007.06.002 . Fischer, H. B ., List, E.J., Koh, R.C.Y., Imberger, J., and Brooks, N.H., 1979. Mixing in Inland and Coastal Waters, pp. 483, Academic Press, San Diego, CA. Folegot, S., Hannah, D. M., Dugdale, S. J., Kurz, M. J., Drummond, J. D., Klaar, M. J. & Krause , S., 2017 . Low f low controls on stream thermal dynamics, Limnologica. doi: 10.1016/j.limno.2017.08.003 Folegot, S., Krause, S., Mons, R., Hannah, D.M. & Datry, T., 2018. Mesocosm experiments reveal the direction of groundwater surface water exchange alters the hyporheic refuge capacity under warming scenarios, Freshwater Biol., 63:165 177. https://doi.org/10.1111/fwb.13049 Garner, G., Malcolm, I.A., Sadler, J.P., Hannah, D.M., 2017. The role of riparian vegetation density , channel orientation and water velocity in determining river temperature dynamics. Journal of Hydrology, 553, 471 - 485, doi:10.1016/j.jhydrol.2017.03.024 Gallice, A., Bavay, M., Brauchli, T., Como la, F., Lehning, M., Huwald, H., 2016. StreamFlow 1.0: an e xtension to the spatially distributed snow model Alpine3D for hydrological modelling and deterministic stream temperature prediction. Geosci. Model Dev. 9, 4491 4519. doi:10.5194/gmd - 9 - 4491 - 2016 G alloway, J.N., Aber, J.D., Erisman, J.W., Seitzinger, S.P., Howarth, R.W., Cowling, E.B., Cosby, B.J., 2003. The Nitrogen Cascade. BioScience 53, 341 356. https://doi.org/10.1 641/0006 - 3568(2003)053[0341:TNC]2.0.CO;2 170 Galloway, J.N., T ownsend, A.R., Erisman, J.W., Bekunda, M., Cai, Z., Freney, J.R., Martinelli, L.A., Seitzinger, S.P., Sutton, M.A., 2008. Transformation of the Nitrogen Cycle: Recent Trends, Questions, and Potenti al Solutions. Science 320, 889 892. https://doi.org/10.1126 /science.1136674 Gent, P. R. et al. (2011). The community climate system model version 4, J. Clim., 24, 4973 - 4991 Glibert, P., Anderson, D., Gentien, P., Granéli, E., Sellner, K., 2005. The Glob al, Complex Phenomena of Harmful Algal Blooms. Oceanography 18, 136 147. https://doi.org/10.5670/oceanog.2005.49 Gordon, R.P., Lautz, L.K., Briggs, M.A., McKenzie, J.M., 2012. Automated calculation o f vertical pore - water flux from field temperature time seri es using the VFLUX method and computer program. Journal of Hydrology 420, 142 158. doi:10.1016/j.jhydrol.2011.11.053 Gunduz, O., 2004. Coupled flow and contaminant transport modeling in large watersheds. Dissertation Thesis, Georgia Institute of Technolog y, Georgia Institute of Technology. Haag, I., Luce, A., 2008. The integrated water balance and water temperature model LARSIM - WT. Hydrol. Process. 22, 1046 1056. doi:10.1002/hyp.6983 Haag, I. & Luce, A., 2008. The integrated water balance and water temper ature model LARSIM - WT, Hydrol. Processes, 22, 1046 1056. doi:10.1002/hyp.6983 Halloran, L.J.S., Andersen, M.S., Rau, G.C., 2017. Investigation of the thermal reg ime and subsurface properties of a tidally affected, variably saturated streambed. Hydrol. Pro cess. 31, 2541 2555. doi:10.1002/hyp.11197 Hannah, D.M., Malcolm, I.A., Soulsby, C., & Youngson, A.F., 2004. Heat exchanges and temperatures within a salmon spaw ning stream in the Cairngorms, Scotland: seasonal and sub - seasonal dynamics. River Res. Appl. 20 (6), 635 652. http://dx.doi.org/10.1002/rra.771 . Hannah, D.M., Malcolm, I.A., Soulsby, C., & Youngson, A.F., 2008. A comparison of forest and moorland stream microclimate, heat exchanges and thermal dyn amics. Hydrol. Process. 22 (7), 919 940. http://dx.doi.org/10.1002/hyp.7003 Hannah D.M., Malcolm I.A. & Bradle y, C., 2009. Seasonal hyporheic temperature dynamics over riffle bedforms, Hydrological Proces ses, 23, 2178 - 2194 DOI: 10.1002/hyp.7256 Hannah D.M. and Garner G., 2015. River water temperature in the United Kingdom: changes over the 20th century and possib le changes over the 21st century, Progress in Physical Geography, 39, 68 92, doi: 10.1177/0309 133314550669 171 Harvey, J. W., Böhlke J.K., Voytek M.A., Scott D., and Tobias C. R., 2013. Hyporheic zone denitrification: Controls on effective reaction depth and contribution to whole - stream mass balance: Scaling hyporheic flow controls on stream denitrif ication, Water Resour. Res.,49,6298 6316, doi:10.1002/wrcr.20492. Heisler, J., Glibert, P., Burkholder, J., Anderson, D., Cochlan, W., Dennison, W., Gobler, C., Dortch, Q., Heil, C., Humphries, E., Lewitus, A., Magnien, R., Marshall, H., Sellner, K., Stoc kwell, D., Stoecker, D., Suddleson, M., 2008. Eutrophication and Harmful Algal Blooms: A Scientific Consensus. Harmful Algae 8, 3 13. https://doi.org/10.1016/j.hal.2008.08.006 Helton, A.M., Poole, G.C., Meyer, J.L., Wollheim, W.M., Peterson, B.J., Mulholland, P.J., Bernhardt, E.S., Stanford, J.A., Arango, C., Ashkenas, L.R., Cooper, L.W., Dodds, W.K., Gregory, S.V., Hall, R.O., Hamilton, S.K., Johnson, S.L., McDowell, W.H., Pott er, J.D., Tank, J.L., Thomas, S.M., Valett, H.M., Webster, J.R., Zeglin, L., 2011. Thinking outside the channel: modeling nitrogen cycling in networked river ecosystems. Frontiers in Ecology and the Environment 9, 229 238. https:// doi.org/10.1890/080211 Helton, A.M., Hall, R.O., Bertuzzo, E., 2018. How network structure can affect nitrogen removal by streams. Freshw Biol 63, 128 140. https://doi.org/10.1111/fwb .12990 Herb, W.R. & Stefan, H.G., 2011. Modified equilibrium temperature models for cold - water streams, Water Resour. Res., 47, W06519. doi:10.1029/2010WR009586 Hinderer J.M. and M.W. Murray , 2011. Feast and famine in the Great Lakes how nutrients and in vasive species intera ct to overwhelm the coasts and starve offshore waters Technical Report , National Wildlife Federation. Holmes, R.M., Jones, J.B., Fisher, S.G., Grimm, N.B., 1996. Denitrification in a nitrogen - limited stream ecosystem. Biogeochemistry 33, 125 146. https://doi.org/10.1007/BF02181035 Homer, C.G., Dewitz, J.A ., Yang, L., Jin, S., Danielson, P., Xian, G., Coulston, J., Herold, N.D., Wickham, J.D., and Megown, K., 2015, Completion of the 2011 National Land Cover Database for the conterminous United States - Represe nting a decade of land cover change information . Photogrammetric Engineering and Remote Sensing, v. 81, no. 5, p. 345 - 354 Bristow, K.L., 2002. Indirect estimation of soil thermal properties and water flux using heat pulse pro be measurements: Geometry and dispersion effects, Water Resour. Res., 38, 7 1. doi:10.1029/2000WR000071 Hrachowitz, M., Soulsby, C., Imholt, C., Malcol m, I. A. and Tetzlaff, D., 2010. Thermal regimes in a large upland salmon river: a simple model to ident ify the influence of 172 landscape controls and climate change on maximum temperatures, Hydrological Processes, 24(23), 3374 3391, doi: 10.1002/hyp.7756 Hudson, R. & Fraser, J., 2005. The mass balance (or dry in jection) method, Streamline Watershed Manage., Bull., 9, 6 - 12. Incropera, F.P. & D. P. DeWitt, D.P., 1996. Introduction to Heat Transfer, Third ed. Joh n Wiley & Sons, 801p. Ingebritsen, S.E. & Sanford, W.E. (1999). Groundwater in Geologic Processes, Cambridge University Press, Cambridge, UK. Inwood, S. E., Tank, J. L. and Bernot, M. J. , 2007. Factors Controlling Sediment Denitrification in Midwestern S treams of Varying Land Use, Microb Ecol, 53(2), 247 258, doi: 10.1007/s00248 - 006 - 9104 - 2 . Jackson, F. L., Hannah, D.M., Fryer, R. J., Millar, C. P., Malcolm, I. A., 2017. Development of spatial regre ssion models for predicting summer river temperatures from landscape characteristics: Imp lications for land and fisheries management. Hydrol. Process. 31, 1225 1238. doi:10.1002/hyp.11087 Jackson, F.L., R.J. Fryer, D.M. Hannah, C.P. Miller and I.A. Malcol m, 2018. A spatio - temporal statistical model of maximum daily river temperatures to infor m the management of Scotland's Atlantic salmon rivers under climate change. Sci. Tot. Environment, 612, 1543 - 1558, https://doi.org/10.1016/j.scitotenv.2017.09.010 Johnson, S. L. and Jones, J. A.: Stream temperature responses to forest harvest and debris flows in western Cascades, Oregon, 2000. Can. J. Fish. Aquat. Sci., 57(S2), 30 39, doi: 10.1139/f00 - 109 . Kelleher, C., Wagener, T., Goos eff, M., McGlynn, B., McGuire, K. and Marshall, L.. 2012. Investigating controls on the thermal sensitivity of Pennsylvania streams, Hydrological Processes, 26(5), 771 785, doi : 10.1002/hyp.8186 . Krause, P., hydrological model assessment, Advances in Geosciences , 5, 89 - 97 Krause S., Hannah D.M. & Blume T., 2011. Interstitial pore - water temperature dynamics across a pool - riffle - pool sequence, Ecohydrology, 4, 549 563 DOI: 10.1002/eco.199 Krause S., Hannah, D.M. & Blume T., 2011. Heat transport patterns at pool - riffle sequences of a UK lowland stream, Ecohydrology, 4(4), 549 - 563, doi: 10.1002/eco.199 Kundu, P.K., Cohen, I.M. & Dowli ng, D.R., 2015. Fluid Mechanics, 6th edition, pp. 928, Aca demic Press, New York. 173 Kurylyk, B. L.; Moore, R. D.; MacQuarrie, K. T. B., 2016. Scientific briefing: quantifying streambed heat advection associated with groundwater surface water interactions. Hyd rological Processes, 30, 987 992, doi: 10.1002/hyp.10709 . Leach, J.A., Moore, R.D., 2010. Above - stream microclimate and stream surface energy exchanges in a wildfire - disturbed riparian zone. Hydrol. Process . 24, 2369 2381. doi:10.1002/hyp.7639 Lee, Mee - Sun, Kang - Kun Lee, Yunjung Hyun, T. Prabhakar Clement, and David Hamilton. Modelling 192, no. 1 2 (February 15, 2006): 143 59. https://doi.org/10.1016/j.ecolmodel.2005.07.013 . Legates, D. R. & McCabe Jr., G. J. , 1999 - of - measures in hydrologic and hydroclimatic model validati on, Water Resour. Res., 35, 1, 233 - 241 Leibowitz, S. G., Comeleo, R. L., Wigington Jr., P. J., Weaver, C. P., Morefield, P. E., Sproles, E. A. and Ebersole, J. L., 2014. Hydrologic landscape classification evaluates streamflow vulnerability to climate ch ange in Oregon, USA, Hydrol. Earth Syst. Sci., 18(9), 3367 3392, doi: 10.5194/hess - 18 - 3367 - 2014 Li, G ., Jackson, C. R. & K. A. Kraseski., 2012. Modeled riparian stream shading: Agreement with field measurements and sensitivity to riparian conditions. J. Hydrol., 428 429, 142 151. doi:10.1016/j.jhydrol.2012.01.032 Li, H. - Y. , Leung, L. R., Tesfa, T., Vois in, N., Hejazi, M., Liu, L., Liu, Y., Rice, J., Wu, H. & Yang, X., 2015. Modeling stream temperat ure in the Anthropocene: An earth system modeling approach, Journal of Advances in Modeling Earth Systems, 7, 1661 - 1679, doi: 10.1002/2015MS000471. Lin, B. - L., Sakoda, A., Shibasaki, R., Suzuki, M., 2001. A Modelling Approach to Global Nitrate Leaching Caused by Anthropogenic Fertilisation. Water Research 35, 1961 1968. https://doi.org/10.1016/S004 3 - 1354(00)00484 - X Loicq P., Moatar, F., Jullian Y., Dugdale S.J. & Hannah D.M., 2018. Improving representation of riparian vegetation shading in a regional stream temperature model using LiDAR data, Science of the Total Environment, 624, 480 - 490 DOI: 10.1 016/j.scitotenv.2017.12.129 Loinaz, M. C., Davidsen, H. K., Butts, M. & Bauer - Gottwein, P., 2013 . Integrated flow and temperature modeling at the catchment scale, J. Hydrol., 495, 238 251. https://doi.org/10.1016/j.jhydrol.2013.04.039 174 Lu, M. - C., Chang, C . - T., Lin, T. - C., Wang, L. - J., Wang, C. - P., Hsu, T. - C., Huang, J. - C., 2017. Modeling the terrestr ial N processes in a small mountain catchment through INCA - N: A case study in Taiwan. Science of The Total Environment 593 594, 319 329. https://doi.org/10.1016/j.scitotenv.2017.03.178 Luscz, E.C., Kendall, A.D., Hyndman, D.W., 2015. High resolution spatially explicit nutrient source models for the Lower Peninsula of Michigan. Journal of Great Lake s Research 41, 618 629. https://doi.org /10.1016/j.jglr.2015.02.004 MacDonald, R.J., Boon, S., Byrne, J.M., 2014. A process - based stream temperature modelling approach for mountain regions. Journal of Hydrology 511, 920 931. doi:10.1016/j.jhydrol.2014.02.009 Martin, J.L., McCutcheon, S.C., 19 98. Hydrodynamics and Transport for Water Quality Modeling. CRC Press. Meisner, J. D., Rosenfeld, J. S. and Regier, H. A., 1988. The Role of Groundwater in the Impact of Climate Warming on Stream Salmonines, Fisheries, 13(3), 2 8, doi: 10.1577/1548 - 8446(1988)013<0002:TROGIT>2.0.CO;2 . Met Office, UK, 2012, Met Office Integrated Data Archive System (MIDAS) Land and Marine Surface Stations Data (1853 - current), NCAS British Atm ospheric Data Centre, Access date: 2017/04/03. http://catalogue. ceda.ac.uk/uuid Michigan Department of Natural Resources (MDNR), 2001. IFMAP/GAP Lower peninsula land cover, Available at <.http://www.mcgi. state.mi.us/mgdl/?rel=thext&action=thmname&cid5&cat=Land+Cover +2001> (retrieved 28.11.16). Miller, M.P., Tesorier o, A.J., Capel, P.D., Pellerin, B.A., Hyer, K.E., Burns, D.A., 2016. Quantifying watershed - scale groundwater loading and in - stream fate of nit rate using high - frequency water quality data. Water Resour. Res. 52, 330 347. https://doi.org/10.1002/2015WR017753 Mohseni, O., Stefan, H.G., Erickson, T.R., 1998. A nonlinear regression model for weekl y stream temperatures. Water Resour. Res. 34, 2685 2692. doi:10.1029/98WR01877 Mohseni, O. and Stefan, H.G., 1999 . Stream temperature/air temperature relationship: a physical interpretation, J. Hydrology, 218, 128 141. Mohseni O., Stefan H.G., Eaton, J.G ., 2003. Global warming and potential changes in fish habitat in US streams. Climatic Change 59(3): 389 409. 175 Moor e, R. D., Sutherland, P., Gomi, T. & Dhakal, A., 2005. Thermal regime of a headwater stream within a clear - cut, coastal British Columbia, Can ada, Hydrol. Processes, 19, 2591 2608. doi:10.1002/hyp.5733 Mu, Q., Zhao, M. & Running, S.W., 2011. Improvements to a MODIS global terrestrial evapotranspiration algorithm, Remote Sens. Environ., 115, 1781 - 1800. doi:10.1016/j.rse.2011.02.019. Mulholland , Patrick J., Ashley M. Helton, Geoffrey C. Poole, Robert O. Hall, Stephen K. Hamilton, Bruce J. Peterson, Jennife 20 08): 202 5. https://doi.org/10.1038/nature06686 . Nash, J. E. & Sutcliffe, J.V., 1970. River flow forecasting through conceptual models part I A discussion of principles, J. Hydrol., 10, 282 290. doi:10 .1016/0022 - 1694(70)90255 - 6 National Climatic Data Center (NCDC) , Available at < http://www.ncdc.noaa.gov/oa/climate/climatedata.html#daily > (accessed 06.01.16). Neitsch, S.L., Arnold, J.G., Kiniry, J.R., Williams, J.R., 2011. Soil and Water Assessment Tool Theoretical Documentation Version 2009. Grassland, Soil and Water Research Laboratory, Agricultural Research Service and Blackland Research Center. Texas Agricultural Experim ent Station, College Station. Texas, USA Neves, R.J.J., 198 5. Etude experimentale et modelisation des circulations trasitoire et residuelle dans l´Estuaire du Sado PhD. Thesis. University of Liege. Niu, J., Shen, C., Li, S. - G. & Phanikumar, M.S., 2014. Q uantifying storage changes in regional Great Lakes watershed s using a coupled subsurface - land surface process model and GRACE, MODIS products, Water Resour. Res., 50, 7359 7377. doi:10.1002/2014WR015589 Niu, J., Phanikumar, M.S., 2015. Modeling watershed - scale solute transport using an integrated, process - based hy drologic model with applications to bacterial fate and transport. Journal of Hydrology 529, Part 1, 35 48. doi:10.1016/j.jhydrol.2015.07.013 Niu, J., Shen, C., Chambers, J. Q., Melack, J. M. and Riley, W. J., 2017: Interannual Variation in Hydrologic Budg ets in an Amazonian Watershed with a Coupled Subsurface Land Surface Process Model, J. Hydrometeor., 18(9), 2597 2617, doi: 10.1175/JHM - D - 17 - 0 108.1 Oleson, K.W., Lawrence, D.M., Gordon, B., Flanner, M .G., Kluzek, E., Peter, J., Levis, S., Swenson, S.C., Thornton, E., Feddema, J.,2010. Technical Description of version 4.0 of the Community Land Model (CLM). NCAR Technical Note, Rep. NCAR/TN - 478þ STR, National Center for Atmospheric Research, Boulder, Colo . 176 Oleson, K., and Coauthors, 2013. Technical description of version 4.5 of the Community Land Model (CLM). NCAR Technical Note, Rep. NCAR/TN - 503+STR, 420 pp, National Center for Atmospheric Resear ch, Boulder, Colo. Parton, W., Stewart, J. and Cole, C., 1 988. Dynamics of C, N, P And S in Grassland Soils - A Model. Biogeochemistry 5:109 - 131. Parton, W.J. et al. 1996. Generalized model for N2 and N2O production from nitrification and denitrificatio n. Global Biogeochemical Cycles 10(3):401 - 412. Parton, W.J. et al. 2001. Generalized model for NOx and N2O emissions from soils. J. Geophys. Res. 1 06(D15):17403 - 17419. Phanikumar, M. S., Aslam, I., Shen, C., Long, D. T. and Voice, T. C.2007. Separating surface storage from hyporheic retention in natural streams using wavelet decomposition of acoustic Doppler current profiles, Water Resour. Res., 43( 5), W05406, doi: 10.1029/2006WR005104 . Phanikumar, M.S . and McGuire J.T. , 2010. A Multi - species reactive transport model to estimate biogechemical rates based on single - well push - pull test data, Computer s & Geoscie nces, 36(8), pp. 997 - 1004 Pope, C.A., Thun, M.J., Namboodiri, M.M., Dockery, D.W., Evans, J.S., Speizer, F.E., Heath, C.W., 1995. Particulate air pollution as a predictor of mortality in a prospective study of U.S. adults. Am. J. Respir. Crit. Care Med. 151, 669 674. https://doi.org/10.116 4/ajrccm/151.3_Pt_1.669 - Iturbe, I., 2003. Hydrologic controls on soil carbon and nitrogen cycles. I. Modeling scheme. Advances in Water Resources 26, 45 58. https://doi.org/10.1016/S0309 - 1708(02)00094 - 5 Price, K. V., Storn, R. M. & Lampinen, J.A. , 2005. Differential Evolution: A Practical Approach to Global Optimization, Berlin: Springer. Qiu, H., Blaen, P., Comer Warner, S., Hannah, D. M., Krause, S. and Phanikumar, M. S, 2019a. Evaluating a coupled phenology surface energy balance model to understand stream subsurface temperature dynamics in a mixed - use farmland catchment, Water Resources Research, 0(ja), do i: 10.1029/2018W R023644 . Qiu, H., Niu, J. and Phanikumar, M. S.: Quantifying the space time variability of water balance components in an agricultural basin using a process - based hydrologic model and t he Budyko framework, 2019b. Science of The Total Environment, 676, 176 189, doi: 10.1016/j.scitotenv.2019.04.147 . 177 Robertson, D.M., Saad, D.A., 2011. Nutrient Inputs to the Laurentian Great La kes by Source and Watershed Estimated Using SPARROW Watershed Models 1. JAWRA Journal of the American Water Resources Association 47, 1011 1033. https://doi.org/10.1111/j.1752 - 1688.2011.00574.x Roth, T.R., Westhoff, M.C., Huwald, H., Huff, J.A., Rubin, J.F., Barrenetxea, G., Vetterli, M., Parriaux, A., Selker, J.S., Parlange, M.B., 2010. Stream Temperature Response to Three Riparian Vegetation Scenarios by Use of a Distributed Temperature Vali dated Model. Environ. Sci. Technol. 44, 2072 2078. doi:10.1021/es902 654f Royer, T.V., Tank, J.L., David, M.B., 2004. Transport and fate of nitrate in headwater agricultural streams in Illinois. J. Environ. Qual. 33, 1296 1304. Rumynin, V.G., 2015. Overlan d Flow Dynamics and Solute Transport. Springer. Rumynin, V.G., 2015 . Overland Flow Dynamics and Solute Transport, Theory and Applications of Transport in Porous Media. Springer International Publishing. Runkel, R.L., 2007. Toward a transport - based analys is of nutrient spiraling and uptake in streams. Limnology and Oceano graphy: Methods 5, 13. Safaie, A., Wendzel, A., Ge, Z., Nevers, M.B., Whitman, R.L., Corsi, S.R. & Phanikumar, M.S., 2016. Comparative evaluation of statistical and mechanistic models of Escherichia coli at beaches in southern Lake Michigan, Environmental Science & Technology, 50, Safaie, A., Dang, C., Qiu, H., Radha, H. & Phanikumar, M.S., 2017. Manifold methods for assimilating geophysical and meteorological data in Earth system models and their components, J. Hy drol., 544, 383 396. doi:10.1016/j.jhydrol.2016.11.009 Schaefer, S.C., Alber, M., 2007. Temperature controls a latitudinal gradient in the proportion of watershed nitrogen exported to coastal ecosystems. Biogeochemistry 85, 333 346. https://doi.org/10.1007/s10533 - 007 - 9144 - 9 Schindler, D. W.2001. The cumulative effects of climate warming and other human stresses on Canadian freshwaters in the new millennium, Can. J. Fish. Aquat. Sci., 58(1), 18 29, doi: 10.1139/f00 - 179 . Schlesinger, W.H., Reckhow, K.H., Bernhardt, E.S., 2006. Global change: The nitrogen cycle and rivers. Water Resour. Res. 42, W03S06. https://doi.org/10.1029/2005WR004300 Schwarz, G. E., A. B. Hoos, R. B. Alexander, and R. A. Smith, 2006. SPARROW: concepts and user documentation for the SPARROW surface water - quality model. Techniques and Methods Report, Book 6, Chapter B3. U.S. Geologi cal Survey, Reston, Virginia, USA. 178 Sellers, P. J. (1985), Canopy reflectance, photosynthesis and transpiration, Int. J. Remote Sens. 6, 1335 1372. , https://doi.org/10.1080/01431168508948283 Schaap , M. G., Leij, F. J. & van Genuchten, M. T., 2001. Rosetta: a computer program for estimating soil hydraulic parameters with hierarchical pedotransfer functions. J. Hydrol., 251, 163 176. doi:10.1016/S0022 - 1694(01)00466 - 8 Shen, C., Phanikumar, M.S., 2010. A process - based, distributed hydrologic model based on a large - scale method for surface subsurface coupling. Advances in Water Resources 33, 1524 1541 . doi:10.1016/j.advwatres.2010.09.002 Shen, C., Niu, J., Phanikumar, M.S., 2013. Evaluating controls on coupled hydrologic and vegetation dynamics in a humid continental climate watershed using a subsurface - land surface processes model. Water Resour. Res. 49, 2552 2572. doi:10.1002/wrcr.20189 Shen, C., J. Niu, and K. Fang, 2014. Quantifying the effects of d ata integration algorithms on the outcomes of a subsurface land surface processes model, Environ. Model. Software, 59, 146 161. doi:10.1016/j.envsoft.2 014.05.006 Shen, C., Riley, W.J., Smithgall, K.M., Melack, J.M. & Fang, K., 2016. The fan of influence o f streams and channel feedbacks to simulated land surface water and carbon dynamics, Water Resources Research, 52, (2), pp. 880 902 Shi, X., Wu, L., C hen, W., Wang, Q., 2011. Solute Transfer from the Soil Surface to Overland Flow: A Review. Soil Science S ociety of America Journal 75, 1214 1225. Smith, R.A., Schwarz, G.E., Alexander, R.B., 1997. Regional interpretation of water - quality monitoring data. Water Resour. Res. 33, 2781 2798. https://doi.org/10.102 9/97WR02171 Soil Survey Staff, Natural Resources Conservation Service, United States Department of Agriculture, U.S. General Soil Map (STATSGO2). Ava ilable online at: https://sdmdataaccess.sc.egov.usga.gov. Accessed 03/01/2017. Starry, O.S., Valett, H.M ., Schreiber, M.E., 2005. Nitrification rates in a headwater stream: influences of seasonal variation in C and N supply. Journal of the North American Benthological Society 24, 753 768. https://doi.org/10.189 9/05 - 015.1 State of Michigan, 2016. Wellogic digital water well dataset, Michigan GIS open data portal, State of Michigan. URL http://gis.michigan.opendata.arcgis.com/datasets? q=wellogic&sort by=relevance. Accessed on June 21, 2016. St - Hilaire, A., Mor in, G., El - Jabi, N., Caissie, D., 2000. Water temperature modelling in a small forested stream: implication of forest canopy and soil temperature. Can. J. Civ. Eng. 27, 1095 1108. doi:10.1139/l00 - 021 179 Sulis, M., Williams, J.L., Shreshta, P., Diederich, M., Simmer, C. & Maxwell, R.M., 2017. Coupling groundwater, vegetation, and atmospheric processes: A comparison of two integrated mod els, J. Hydrometeorology, vol. 18, pp. 1489 - 1511, doi: 10.1175/JHM - D - 16 - 0159.1 Subehi, L., Fukushima, T., Onda, Y., Mizugaki, S., Gomi, T., Kosugi, K., Hiramatsu, S., Kitahara, H., Kuraji, K. and Terajima, T. 2010. Analysis of stream water temperature cha nges during rainfall events in forested watersheds, Limnology, 11(2), 115 124, doi: 10.1007/s10201 - 009 - 0296 - 2 . Sun, N., Yearsley, J., Voisin, N., Lettenmaier, D.P., 2015. A spatially distributed model for the assessment of land use impacts on stream temperature in small urban watersheds. Hydrol. Process. 29, 2331 2345. doi: 10.1002/hyp.10363 Swenson, S.C., Lawrence, D.M., and Lee, H., 2012. Improved simulation of the terrestrial hydrol ogical cycle in permafrost regions by the Community Land Model, Journal of Advances in Modeling Earth Systems, vol. 4, M08002. doi:10.1029/201 2MS000165. Testa, J.M., Li, Y., Lee, Y.J., Li, M., Brady, D.C., Di Toro, D.M., Kemp, W.M., Fitzpatrick, J.J., 201 4. Quantifying the effects of nutrient loading on dissolved O2 cycling and hypoxia in Chesapeake Bay using a coupled hydrodynamic biogeochemic al model. Journal of Marine Systems 139, 139 158. ht tps://doi.org/10.1016/j.jmarsys.2014.05.018 Therrien R, McLaren R.G., Sudicky E.A., Panday S. M., 2010. HydroSphere: A Three - Dimensional Nume rical Model Describing Fully - Integrated Subsurface and Surface Flow and Solute Transport. Groundwater Simulations Group: Waterloo. Theurer, F. B., K. A. Voos, K.A., and W. J. Miller, 1984. Instream water temperature model. Instream flow information paper 16, U.S. Fish. Wildl. Serv., FWS/OBS - 85/15. Thomas, Z., Abbott, B.W., Troccaz, O., Baudry, J., Pinay, G., 2016. Proximate and ultimate controls on carbon and nutrient dynamics of small agricultural catchments. Biogeosciences 13, 1863 1875. https://doi.org/https://doi.org/10.5194/bg - 13 - 1863 - 2016 Thornton, P.E., Rosenbloom, N.A., 2005. Ecosystem model spin - up: Estimating steady state conditions in a coupled terrestrial carbon and nitro gen cycle model. Ecological Modelling 189, 25 48. https://doi.org/10.1016/j.ecolmodel.2005.04.008 Tong, J. - X., Yang, J. - Z., Hu, B.X., Bao, R., 2010. Experimental study and mathematical modell ing of soluble chemical transfer from unsaturated/saturated soil to surface runoff. Hydrol. Process. 24, 3065 3073 . https://doi.org/10.1002/hyp.7722 Trancoso, A.R., Braunschweig, F., Chambel Leitao, P., Obe rmann, M., Neves, R., 2009. An advanced modeling tool for simulating complex river systems. Sci. Total Environ. 407, 3004 3016. http://dx.doi.org/10.1016/j.scitotenv.2009.01.015. 180 USDA, ARMS Farm Financial and Crop Production Practices. Available online at: https://www.ers.usda.gov/data - products/arms - farm - financial - and - crop - production - practices . Accessed 12/ 01/2013. USEPA, 1996. Environmental indicators of water quality in the United States. EPA 841 - R - 96 - 002, USEPA, Office of Water (4503F), U.S. Government Printing Office, Washington, D.C. Vandersteen, G., Schneidewind, U., Anibas, C., Schmidt, C., Seuntjen s, P. & Batelaan, O., 2015. Determining groundwater - surface water exchange from temperature - time series: Combining a local polynomial method with a maxim um likelihood estimator, Water Resour. Res., 51, 922 939. doi:10.1002/2014WR015994 Van Vliet, M. T. H. , Yearsley, J. R., Franssen, W. H. P., Ludwig, F., Haddeland, I., Lettenmaier, D. P. & Kabat, P., 2012. Coupled daily streamflow and water temperature m odelling in large river basins, Hydrol. Earth Syst. Sci., 16, 4303 4321. Viney, N.R., Sivapalan, M., D eeley, D., 2000. A conceptual model of nutrient mobilisation and transport applicable at large catchment scales. Journal of Hydrology 240, 23 44. https://doi.org/10.1016/S0022 - 1694(00)00320 - 6 Wang, D. and Wu, L.2013. Similarity of climate control on base flow and perennial stream density in the Budyko framework, Hydrol. Earth Syst. Sci., 17(1), 315 324, doi: 10.5194/hess - 17 - 315 - 2013 . Wade , A.J., Durand, P., Beaujouan, V., Wessel, W.W., Raat, K.J., Whitehead, P.G., Butterfield, D., Rankinen, K., Lepisto, A., 2002. A nitrogen model for Europ ean catchments: INCA, new model structure and equations. Hydrol. Earth Syst. Sci. 6, 559 582. https://doi.org/10.5194/hess - 6 - 559 - 2002 Wallach, R., Grigorin, G., Rivlin (Byk), J., 2001. A comprehensiv e mathematical model for transport of soil - dissolved chemicals by overland flow. Journal of Hydrology 247, 85 99. https://doi.org/10.1016/S0022 - 1694(01)00365 - 1 Warrick, A.W., 2003. Soil Water D ynamics, Oxford University Press, UK, pp. 416. Webb, B.W., Hannah, D.M., Moore, R.D., Brown, L.E. & Nobilis, F., 2008. Recent advances in stream and river temperature research, Hydrol. Process. 22, 902 918, DOI: 10.1002/hyp.6994. Wehrly, K. E., Wiley, M. J. and Seelbach, P. W.2003. Classifying Regional Variation in Thermal Regime Based on Stream Fish Co mmunity Patterns, Transactions of the American Fisheries Society, 132(1), 18 38, doi: 10.1577/1548 - 8659(2003)132<0018:CRVITR>2.0.CO;2 , 2003. 181 Wesley, J.K., 2005. Kalamazoo Ri ver Assessment. Fisheries Division Special Report. Michigan Department of Natural Resources, Michigan, USA. Available online at http://www.michigan.gov/dnr/0,4570,7 - 153 - 10364_52259_19056 - 46270 -- ,00.html Whitehead, P.G., Wilson, E.J., Butterfield, D., 1998. A semi - distributed Integrated Nitrogen model for multiple sour ce assessment in Catchments (INCA): Part I model structure and process equations. Science of The To tal Environment 210 211, 547 558. https://doi.org/10.1016/S0048 - 9697(98)00037 - 0 Wick, K., He umesser, C., Schmid, E., 2012. Groundwater nitrate contamination: Factors and indicators. J Environ M anage 111, 178 186. https://doi.org/10.1016/j.jenvman.2012.06.030 Wollheim, W.M., Vörösmarty , C.J., Peterson, B.J., Seitzinger, S.P., Hopkinson, C.S., 2006. Relationship between river size and nutrient removal. Geophys. Res. Lett. 33, L06410. https://doi.org/10.1029/2006GL025845 Wool, T.A., Am brose, R. B., Martin, J. L. & Comer, E.A., 2001. Water Quality Analysis Simulation Program (WASP) Ve Yang, X., Jomaa, S., Zink, M., Fleckenstein, J.H., Borchardt, D., Rode, M., 2018. A New Fully Distributed Model of Nitrate Transport and Removal at Catchment Scale. Water Resources Research 54, 5856 5877. https://doi.org/10.1029/2017WR022380 Ye, S., Covino, T. P., Sivapalan, M., Basu, N. B., Li, H. - Y. and Wang, S. - W.: Dissolved nutrient retention dynamics in river networks: A modeling investigation of tra nsient flows and scale effects, Water Resources Research, 48(6), doi: 10.1029/2011WR010508 , 2012. Zeng, X., Zhao, M. & Dickinson, R.E. , 1998 . Intercomparison of bulk aerodynamic algorithms for the comput ation of sea surface fluxes using TOGA COARE and TAO Data, J. Clim., 11, 2628 2644. Zhang , Y. & Schaap, M.G. , 2017 . Weighted recalibration of the Rosetta pedotransfer model with improved estimates of hydraulic parameter distributions and summary statistics (Rosetta3), J. Hydrol., 547, pp. 39 - 53, doi: 10.1016/j.jhydrol.2017.01.004. Zheng, C., Bennett, G.D., 2002. Applied Contaminant Transport Modeling. Wiley Interscience, 621 pp.