MODELING THE MOVEMENT OF WATER，BACTERIA AND NUTRIENTS ACROSS HETEROGENEOUS LANDSCAPES IN THE GREAT LAKES REGION USING A PROCESS-BASED HYDROLOGIC MODEL By Jie Niu A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of Environmental Engineering - Doctor of Philosophy 2013 ABSTRACT MODELING THE MOVEMENT OF WATER， BACTERIA AND NUTRIENTS ACROSS HETEROGENEOUS LANDSCAPES IN THE GREAT LAKES REGION USING A PROCESS-BASED HYDROLOGIC MODEL By Jie Niu The development and application of process-based hydrologic models (PBHMs) continues to be a topic of significant interest to the hydrologic community. Although numerous studies have applied PBHMs at small scales ranging from plot and field scales to smallwatershed scales, the application of PBHMs to understand large-scale hydrology remains a topic that is relatively unexplored. Understanding controls on large-scale hydrology is key to climate change assessments and effective water resources management; therefore, to quantify the nature and magnitude of fluxes in regional Great Lakes watersheds, we use a new distributed hydrologic model (PAWS+CLM). Here we describe the application of the model to several large watersheds in the State of Michigan including the Grand River, Saginaw Bay, Kalamazoo and Red Cedar River watersheds and evaluate model performance by comparing model results with different types of data including point measurements of streamflows, groundwater heads, soil moisture, soil temperature as well as remotely-sensed datasets for evapotranspiration (ET) and land water thickness equivalent (GRACE). We then report a budget analysis of major hydrologic fluxes and compute annual-average fluxes due to infiltration, ET, surface runoff, sublimation, recharge, and groundwater contributions to streams etc. as percentage of precipitation and use this information to understand the inter-annual variability of these fluxes and to quantify storage in these large watersheds. After testing the model for its ability to describe hydrologic fluxes and states, we describe the development of solute transport models at the watershed scale by using a mechanistic, reactive transport modeling framework in which the advection, dispersion and reaction steps are solved using an operator-splitting strategy. The solute transport models are tested extensively using available analytical solutions for different hydrologic domains and then applied to describe transport with surface – subsurface interactions and to describe the fate and transport of fecal indicator bacteria such as Escherichia coli in the Red Cedar River watershed in Michigan. Following the successful application of the bacterial fate and transport model, we describe detailed reactive transport modules for predicting the levels of nutrients (N and P). The models are tested using available field observations for the Kalamazoo River watershed in Michigan. The watershed-scale fate and transport modules are expected to aid management by quantifying the impacts of upstream watershed influences on water quality in downstream receiving water bodies such as lakes and oceans. Together with the flow modules they represent a comprehensive suite of process-based models to describe the terrestrial hydrologic cycle coupled with vegetation/land surface and biogeochemical processes. To my parents, my wife, and my lovely daughter iv ACKNOWLEDGEMENTS I would like to express my sincere appreciation to my advisor, Dr. Mantha S. Phanikumar, for his support and diligent, insightful and incredibly careful guidance throughout the time of my dissertation research. This dissertation would have been impossible without his help. His persistent enthusiasm, relentless professionalism and uncompromising attitude toward academic research are always excellent examples for me. I am grateful that I had the opportunity to work with him during my Ph.D. studies. I would like to thank my committee members: Dr. Dave Long, Dr. Irene Xagoraraki, Dr. Jiaguo Qi and Dr. Shu-Guang Li for their invaluable comments, suggestions and advice to improve this thesis. I would like to thank previous lab mates Chaopeng Shen for his advice solving problems and the model development, and the current lab mates Guoting Kang and Han Qiu for their assistance with data processing. I also thank Dr. Huasheng Liao in Dr. Shu-Guang Li’s group for his advice with groundwater data processing. Finally I want to thank Miss Ye Yao, my wife for her great effort for my thesis formatting, and dedicated support for my Ph.D. study. v TABLE OF CONTENTS LIST OF TABLES··········································viii ·········································· LIST OF FIGURES ·········································· ········································· ix Chapter 1 Introduction ······································· 1 ······································· 1.1. Importance of quantifying water budgets and storage··················· ·················· 1 1.2. Motivation for a new transport model ···························· ··························· 4 Chapter 2 Quantifying Storage in Regional Great Lakes Watersheds Using GRACE, MODIS Products and Coupled Subsurface – Land Surface Process Models ········· 8 ········· 2.1. Introduction ··········································· 9 ··········································· 2.2. Methodology ·········································· ·········································10 2.2.1. The PAWS model ·················································· 10 ·················································· 2.2.2. Description of Sites················································· 12 ················································· 2.2.3. Data Sources ·······················································16 ······················································ 2.2.4. Simulation details···················································19 ·················································· 2.2.5. Storage and Water Budget Calculations ································19 ······························· 2.3. Results and Discussion····································21 ···································· 2.3.1. Stream flow and groundwater head comparisons ························ 21 ························ 2.3.2. Soil moisture and soil temperature comparisons ·························22 ························ 2.3.3. Evapotranspiration Comparisons with MODIS data ······················27 ····················· 2.3.4. Budget analysis ···················································· 30 ···················································· 2.3.4.1 Linear correlation analysis of annual average results ············· ············30 2.3.4.2 Storage and partition ································ ·······························33 2.3.4.3 Budyko curves ···································37 ··································· 2.4. Conclusions ··········································· ··········································42 Chapter 3 Modeling Solute Transport at the Watershed Scale Using a Process-Based, Distributed Hydrologic Model with Application to Bacterial Fate and Transport Modeling44 3.1. Introduction ··········································· ··········································45 3.2. Methods ············································· ············································45 3.2.1. 2D Overland Flow Transport Equation ·································45 ································ 3.2.2. 1D Channel Flow Equation ·········································· 47 ·········································· 3.2.3. Interactions between overland flow and channel flow···················· 53 ···················· 3.2.4. Transport in the Vadose zone ·········································56 ········································ 3.2.5. Two-dimensional groundwater transport equation ······················· 59 ······················· 3.2.6. Particle tracking scheme and mass balance ····························· 61 ····························· 3.2.7. Microorganisms release from manure and loss rate ······················ 67 ······················ 3.3. Results and Discussion····································72 ···································· 3.3.1. One-dimensional channel transport ····································72 ··································· 3.3.2. Two-dimensional overland flow transport······························ 75 ······························ 3.3.3. One-dimensional solute transport in a soil column ·······················80 ······················ vi 3.3.4. Groundwater transit time simulation ···································87 ·································· 3.3.5. Plot scale simulations of manure-borne fecal coliforms ···················94 ·················· 3.3.6. Modeling Escherichia coli fate and transport in the Red Cedar River watershed ·································································97 ································································ Chapter 4 Modeling Nutrient Transport in Regional Great Lakes Watersheds····· 102 ····· 4.1. Introduction ·········································· ········································· 102 4.2. Carbon model ········································ 104 ········································ 4.2.1. Carbon in litter pool ················································105 ··············································· 4.2.2. Carbon in humus pool ··············································106 ············································· 4.2.3. Carbon in microbial biomass pool ····································107 ··································· 4.3. Nitrogen model ········································ ······································· 108 4.3.1. Nitrogen in litter pool ·············································· 111 ·············································· 4.3.2. Nitrogen in humus pool ·············································112 ············································ 4.3.3. Nitrogen in microbial biomass pool ·································· 112 ·································· 4.3.4. Nitrogen mineralization and immobilization rates ······················ 113 ······················ 4.3.5. Ammonium and Nitrate·············································116 ············································ 4.4. Phosphorus model ······································ ····································· 118 4.4.1. Organic phosphorus ················································121 ··············································· 4.4.1.1 Organic phosphorus stored in the living biomass pool ( P )········ V ······· 121 4.4.1.2 Phosphorus stored in plant litter ( PL ) ····················· ···················· 124 4.4.1.3 Phosphorus in microbial biomass pool ( PM ) ················ 125 ················ 4.4.1.4 Phosphorus in humus pool ( PH )························ ······················· 126 4.4.2. Inorganic phosphorus ·············································· 127 ·············································· 4.5. Model results ········································· ········································ 128 4.5.1. Nitrogen test case ··················································129 ················································· 4.5.2. Phosphorus test case ··············································· 134 ··············································· 4.5.3. Nitrogen and Phosphorus model application to Kalamazoo watershed ·····140 ···· 4.5.3.1 Study area ····································· ···································· 140 4.5.3.2 Nitrogen results ·································· ································· 143 4.5.3.3 Dissolved phosphorus results ·························· ························· 146 4.6. Conclusions ·········································· ········································· 149 Chapter 5 Conclusion ······································ 151 ······································ REFERENCES ············································ ··········································· 153 vii LIST OF TABLES Table 2.1 Processes and governing equations of the hydrological model PAWS ................11 Table 2.2 NCDC weather stations ........................................................................................17 Table 2.3 LULC coverage percentages in the Grand River and Saginaw Bay watersheds ..18 Table 2.4 USGS gauging stations .........................................................................................18 Table 2.5 Spearman correlation coefficients ........................................................................32 Table 2.6 Annual budgets (unit: mm/yr)...............................................................................34 Table 2.7 Subbasin characteristics of the Grand River watershed in relation to Budyko pairs .......................................................................................................................................40 Table 4.1 Parameters representative of the climate, soil, vegetation and nitrogen pool.....129 Table 4.2 Parameters values for soil, climate, vegetation, carbon and phosphorus cycles that are representative of a Brazilian Cerrado ecosystem. ..................................................136 viii LIST OF FIGURES Figure 2.1 Map of GR and SB watersheds. Elevation is shown as the color gradient. NHD streams, USGS gauges, groundwater stations and NCDC weather stations are shown. ( For interpretation of the references to color in this and all other figures, the reader is referred to the electronic version of this dissertation.) .........................................................13 Figure 2.2 Land Use and Land Cover map for GR and SB watersheds. ..............................15 Figure 2.3 Simulated and observed streamflow comparison for GR watersheds at 2 different USGS gauges. USGS gauge station ID and R values are shown in the each subplot. ..................................................................................................................................23 Figure 2.4 Simulated and observed streamflow comparison for SB watersheds at different 2 USGS gauges. USGS gauge station ID and R values are shown in the each subplot. ........24 Figure 2.5 Simulated and observed transient groundwater head comparison at Holt, Michigan (Grand River watershed). .....................................................................................25 Figure 2.6 Observed versus simulated depth to water table. The observed data are obtained by kriging data in the WELLOGIC database. ........................................................26 Figure 2.7 Soil moisture (a) and soil temperature (b) comparison at Belding. .....................27 Figure 2.8 Monthly simulated and MODIS ET data comparison. R 2 value is shown in the figure. ..............................................................................................................................29 Figure 2.9 Five-year averaged spatial map of simulated and MODIS ET............................30 Figure 2.10 Cumulative storage changes for GR (a) and SB (b) watersheds. ......................36 Figure 2.11 Comparison of water thickness anomaly between gravity-based satellite measurements (GRACE data) and PAWS model descriptions.............................................36 Figure 2.12 Budyko curve for the GR watershed. ................................................................38 Figure 2.13 The effect of groundwater flow in shifting the Budyko pairs. ..........................42 Figure 3.1 Method of characteristics at junctions for advection step ...................................51 Figure 3.2 Treatment of junctions .........................................................................................52 Figure 3.3 The process calculating river/land exchange. ......................................................54 ix Figure 3.4 1D channel transport simulation with constant velocity but three different values of dispersion...............................................................................................................74 Figure 3.5 Contour plots of 2D overland transport showing analytical and numerical solution at the same time step. ..............................................................................................77 Figure 3.6 2D overland transport results showing the breakthrough curves at different downstream locations (a) and the longitudinal concentration distribution at different times (b) from the center of the patch source. ................................................................................79 Figure 3.7 1D soil transport showing relative concentrations for step increase as a function of time (a), step increase as a function of pore volumes (b), set decrease as a function of pore volumes (c), and pulse as a function of pore volumes (d). ........................83 Figure 3.8 1D soil transport results showing effect of apparent diffusion, velocity, and displacement length on relative concentration plotted as a function of pore volumes. ........86 Figure 3.9 Conceptual model sketch for flow in a uniformly recharged, unconfined, horizontal aquifer used for developing the general analytical solution. The left-hand boundary is impermeable, and the flow discharges through the right-hand fixed-head boundary. The travel time is calculated between the two arbitrary points: xi (the departure point at the water table) and x. ..............................................................................88 Figure 3.10 Groundwater travel time results ( K  104 m/s) showing model comparison of elapsed transit time (a), the corresponding water table profile and streamlines (b), and the velocity vector field (c) as function of distance. .............................................................91 Figure 3.11 Groundwater travel time results ( K  103 m/s) showing model comparison of elapsed transit time (a), and the corresponding water table profile and streamlines as function of distance. ..............................................................................................................92 Figure 3.12 Groundwater travel time results ( K  105 m/s) showing model comparison of elapsed transit time (a), and the corresponding water table profile and streamlines as function of distance. ..............................................................................................................93 Figure 3.13 Plot scale simulation results showing the accumulative runoff (a) and the bacteria concentration (b) as function of time. .....................................................................96 Figure 3.14 Map of RCR watershed. Elevation is shown as color gradient. NHD streams, USGS gauges, WWTP, CAFO and NCDC weather stations are shown. .............................98 Figure 3.15 Comparison of observed and simulated E. coli concentration at three monitoring sites on Red Cedar river. ..................................................................................100 x Figure 4.1 Schematic representation of the soil carbon budget. Carbon pools are represented by boxes with bold text while fluxes between pools are represented by arrows with italicized text. Carbon stored in living biomass (as well as the flux due to FIRE) is not considered explicitly in this model. ..............................................................................108 Figure 4.2 Systematic sketch of nitrogen cycles in the soil ................................................110 Figure 4.3 Schematic representation of the components of the nitrogen cycles. ................111 Figure 4.4 Schematic representation of the soil and biomass phosphorus budget. Organic and Inorganic phosphorus pools are separated in two groups by boxes with bold text while fluxes between pools are represented by arrows with italicized text. .......................121 Figure 4.5 Temporal dynamics of soil moisture (a), carbon litter (b), carbon humus (c), carbon biomass (d), ammonium (e), and nitrate (f) simulated for the case of the broadleafed savanna at Nylsvley. The broken lines represent the average values or the range of values observed at Nylsvley. The average values are reported whenever observations are available. .............................................................................................................................133 Figure 4.6 Simulated rates of litter decomposition (a), net mineralization (b), nitrate uptake (c) and nitrate leaching (d) in the broad-leafed savanna at Nylsvley. The broken lines represent the average values observed at Nylsvley. The nitrate leaching (d) is reported to be negligible at Nylsvley. .................................................................................134 Figure 4.7 Simulated long term dynamics of net mineralization (a), phosphorus uptake (b), litter decomposition (c), and leaching of PI beneath the rooting zone (d). The dashed lines represent the mean observed fluxes for this system. ..................................................138 Figure 4.8 Long term temporal dynamics for simulated state variables in a Brazilian Cerrado ecosystem: (a) soil moisture; (b) carbon in the litter pool; (c) carbon in the microbial biomass; and, (d) available inorganic phosphorus in soil. ..................................139 Figure 4.9 Map of Kalamazoo watershed showing the location of USGS gauging, NCDC stations and nitrogen sampling sites et. al. ..........................................................................141 Figure 4.10 Land use and land cover map for Kalamazoo watershed. ...............................142 Figure 4.11 Streamflow comparison for USGS 04108660 in Kalamazoo watershed ........143 Figure 4.12 Nitrate concentration at four sites in Kalamazoo River. .................................145 Figure 4.13 Dissolved phosphorus concentration at four sites in Kalamazoo River. .........148 xi Chapter 1 Introduction 1.1. Importance of quantifying water budgets and storage The Great Lakes account for 18% of the World’s and 90% of the United States’ freshwater supply [Sea Grant College Program, 1985]. The Great Lakes are an important source of fresh water, recreational resource and transportation corridor for the Midwestern United States and Canada [Cherkauer and Sinha, 2010]. The Great Lakes and their connecting waterways support large commercial shipping in the United States [Adam et al., 1995]; so the declining water levels in the Great Lakes are a source of major concern. Estimates of runoff (related to watershed storage) and other components of the hydrologic cycle are important to understand and quantify lake level changes. Studies with a focus on hydrologic water budgets provide the scientific basis for managing the water usage effectively, especially in the presence of climate change. According to previously reported research, winters in the Great Lakes region are getting shorter, annual average and high temperatures in summer are on the rise, and the frequency of heavy rain events has been increasing [Christensen et al., 2007; Kling et al., 2003; Trapp et al., 2007; Wuebbles and Hayhoe, 2004]. It is important to understand how key components of the hydrologic cycle are changing in response to climate change, relative to current conditions. However, except for a few major components of the water cycle, considerable variability exists in estimates of current water budgets due to differences in model conceptualizations of hydrologic processes (e.g., lumped parameter versus process-based models). Regional groundwater flow is a major component of the water budget globally and also in the hydrologic cycle in the Great Lakes region. However not all watershed models explicitly describe groundwater 1 processes; therefore, there is a need to quantify water budgets using models that include subsurface processes. While a number of climate change related studies are underway, there is an urgent need to quantify current water budgets, comparing estimates from different types of hydrologic models in order to assess and compare inter-model variability in water budget estimates to projected future changes due to climate change. The hydrologic water budget represents the balance between precipitation, evapotranspiration, surface run-off, and total water storage change, and includes all water pathways at or near the land surface, including channel flow, infiltration and recharge. There are several ways to estimate hydrologic budgets, e.g., using field observations and data for different hydrologic components, or via real-time satellite products [Gao et al., 2010]. Among these water budget components, channel flow can be observed by stream flow gauges. With dense weather station network and satellite-based estimates, precipitation can be quantified for this climatic region. The challenge remains for estimating other spatially varying components of the hydrologic budget at large scales. For example, information on how the largest storage components within the watershed are changing has important implications for water resources management. In addition, due to the spatial heterogeneity and sampling errors, the existing observational data alone may not be sufficient to provide accurate water budget information for research interests [Pan and Wood, 2006; Roads et al., 2003]. Spatially-distributed hydrologic fluxes from model reanalysis have become more popular (e.g. NLDAS [Mitchell et al., 2004]). Well-tested, process-based hydrologic models, which are based on conservation principles of mass, momentum and energy, may be better suited for providing detailed estimates of water 2 budgets and provide information on how storage components are changing within the watershed. One of the objectives of this thesis is to fill this important knowledge gap for two important Great Lakes watersheds. Recently, process-based models became an important tool for studying the hydrologic water budgets [Berbery et al., 2003; Luo et al., 2005; Roads et al., 2003; Roads et al., 1999; Shen et al., 2013; Shen and Phanikumar, 2010; Su and Lettenmaier, 2009]. Storing water is one of the important hydrologic functions of a watershed [Black, 1997; Sivapalan, 2006]. Watershed storage is a key variable used in watershed classification [Wagener et al., 2007] and is intimately tied to surface and subsurface stormflows and runoff ratios in a watershed [Sawicz et al., 2011; Sayama et al., 2011]. Storage is also a key variable that controls the coevolution of hydrologic and biogeochemical states of a watershed. Recent research highlights the link between watershed storage and nutrient levels [Seilheimer et al., 2013] and among nutrient (specifically, nitrogen) levels, net primary production and hydrologic fluxes such as ET [Dickinson et al., 2002; Shen et al., 2013]. Understanding changes in watershed storage is important for climate change studies as storage is directly related to the resilience of a watershed [Tague et al., 2008]; however, quantifying storage is difficult since many watershed models do not explicitly simulate subsurface process and thus do not adequately describe the past storage trajectory and the storage-discharge relations due to complex subsurface heterogeneity [Beven, 2006; Sayama et al., 2011]. The application of physics-based watershed models to large river basins is a topic of increasing interest. Combining gravity-based satellite measurements [Rodell et al., 2007], MODIS products for evapotranspiration and leaf area index and ground-based 3 measurements for other key hydrologic variables (e.g., soil moisture, streamflows, groundwater heads, temperature) opens up the possibility to test physics-based watershed models for their ability to describe changes in watershed storage in regional watersheds, however, such studies are limited to the best of our knowledge. We quantify storage in two regional Great Lakes watersheds using a physics-based watershed model and independent (satellite and ground-based) measurements of key hydrologic fluxes. Specifically, we seek to answer questions such as: How is storage changing in regional Great Lakes watersheds? What are the major components of the hydrologic cycle as described by well-tested process-based models? What is the long-term relation between different components, as expressed, for example, by the Budyko curve? How do storage and groundwater flow influence the long-term relation between different hydrologic components inn regional Great Lakes watersheds? 1.2. Motivation for a new transport model In the Great Lakes region the amounts of nutrients and organic material entering the lakes have increased with intensified urbanization and agriculture. Increased beach closures due to microbiological pollution, harmful algal blooms and eutrophication and drinking water related issues continue to be causes for concern highlighting the need for accurate transport models. According to NRDC [Dorfman and Rosselot, 2011] the number of beach closings and advisories in 2010 were the second highest in 21 years indicating that beachgoers continue to be at risk of exposure to contaminated water. The assessment and management 4 of non-point sources of microbial pollution, in particular, is an issue of great interest. Pathogens that have the potential to infect humans can be divided into the categories of bacteria, protozoans and viruses [Jamieson et al., 2004]. The United States Environmental Protection Agency (USEPA) recommends that E. coil be used as the principal indicator organism in freshwaters, instead of fecal coliforms. Research has shown that E. coli densities are strongly correlated with swimming-associated gastroenteritis [USEPA, 2001]. Models that have been developed to simulate exclusively landscape microbial pollution processes include MWASTE [Moore et al., 1989], COLI [Walker et al., 1990] and SEDMOD [Fraser et al., 1998]. These models successfully simulate the manure-borne bacteria releasing process and the transport of bacteria through runoff, however, the die-off due to solar radiation and other environmental factors was not considered, and there was no interaction among different hydrologic components in these models. Other models have also been developed to simulate survival and transport of fecal indicator bacteria in receiving waters such as lakes [Canale et al., 1993] and rivers [Wilkinson et al., 1995]. These models do not simulate the physical hydrologic processes including flow through channel networks which are key factors influencing transport. Their method is not suitable for large scale watershed transport where major areas of the watershed lack long-term observations of flow and bacterial concentrations. Models that incorporate both landscape and in-stream microbial processes include the Soil and Water Assessment Tool (SWAT) [Sadeghi and Arnold, 2002], a watershed model developed by [Tian et al., 2002], a modified SWAT model [Cho et al., 2012], which considered solar radiation associated bacterial die-off, and a coastal watershed fecal coliform fate and transport model based on 5 Hydrologic Simulation Program Fortran (HSPF) [Rolle et al., 2012]. These models tend to use conceptual representations of the groundwater and vadose zone compartments, ignoring the complexity of the flow in the subsurface domain. The subsurface flow is an integral component of the hydrologic cycle and in shallow water table conditions groundwater controls soil moisture and provides sources of water for evapotranspiration (ET) [Shen and Phanikumar, 2010]. Therefore, the subsurface flow domain should be fully considered and integrated into the transport problem. Surface water eutrophication resulting from nonpoint source phosphorus and nitrogen inputs is the most common water quality problem in the United States (USEPA, 1996). Point sources of pollution, including phosphorus, from industrial discharges and municipal waste water treatment plants are easily identified and quantified. Point source phosphorus loading to surface waters has been greatly reduced since the implementation of the 1972 Federal Clean Water Act. However, control of diffuse nonpoint source pollution has been less successful and is considered the main cause of eutrophication in lakes, streams and coastal areas in the United States. A variety of models have been developed to simulate hydrological processes, nutrient transport through surface runoff, soil infiltration, and groundwater flow, as well as in-stream nutrient processes at different scales. Examples include SWAT [Arnold et al., 1998], HSPF [Donigian et al., 1995], INCA [Whitehead et al., 1998; Wade et al., 2002] and LASCAM [Viney et al., 2000]. Some of these models provide long-term, daily simulation of nutrient load in large catchments [Arnold and Fohrer, 2005; Viney et al., 2000], whereas others are event-based and are clearly unsuitable for long-term continual predictions. Although the long-term continuous nutrient models are capable of 6 providing accurate results, a large number of parameters cannot be obtained from field measurements and must instead be determined through model calibration. An additional constraint to model development and verification is that water quality and hydrometeorological data are rarely simultaneously collected in a satisfactory resolution. A process-based, watershed-scale reactive transport modeling framework was developed in this work to describe fluxes of water, conservative solutes and bacteria exported from watersheds to downstream receiving water bodies such as lakes and oceans. Major motivation for developing the new transport model includes the following: (1) the development of mechanistic, grid-based models is expected to facilitate seamless integration of watershed models with similar (i.e., grid-based) lake and ocean models to make real-time (or near real-time) forecasts of water quality. (2) By describing coupled hydrologic and nutrient (C,N,P) processes in one integrated model that involves multiple hydrologic domains, we expect to address interdisciplinary research questions that involve ecosystem functions and services (e.g., the impact of land use and climate change on net primary productivity, evapotranspiration, water yield etc). With these objectives in mind, two models to describe carbon-nitrogen cycles and carbon-phosphorus cycles in the soil, were adapted and coupled to the PAWS+CLM hydrologic model, to describe the nutrients concentrations in the soil and use the transport routing schemes to simulate nitrate and dissolved phosphorus concentration observed in the streams in Kalamazoo watershed in Michigan. 7 Chapter 2 Quantifying Storage in Regional Great Lakes Watersheds Using GRACE, MODIS Products and Coupled Subsurface – Land Surface Process Models As a direct measure of watershed resilience, watershed storage is important for understanding climate change impacts on water resources. In this paper we estimate storage and water budgets for two of the largest watersheds in the State of Michigan, USA (the Grand River and the Saginaw Bay watersheds) using remotely-sensed data and a processbased hydrologic model that includes detailed representations of subsurface and land surface processes. Model performance is evaluated using ground-based observations (streamflows, groundwater heads, soil moisture, soil temperature) as well as satellite-based estimates of evapotranspiration and watershed storage. We also applied the model on a series of benchmark problems and compared with the results generated by other hydrologic models. We use the model to compute the annual-average fluxes due to evapotranspiration, surface runoff, recharge and groundwater contribution to streams. We also analyze the impacts of land use and land cover (LULC) and soils types on annual hydrologic budgets using correlation analysis. We compare model results of storage with GRACE product equivalent water thickness and analyze water budgets using the Budyko framework for subbasins based on topography, land use land cover and subsurface storage capacity. This work provides new estimates of water budgets and storage in Great Lakes watersheds and the results are expected to aid in the analysis and interpretation of the current trend of declining lake levels, in understanding projected future impacts of climate change as well as in identifying appropriate climate adaptation strategies. 8 2.1. Introduction Recently there has been a surge of interest in the application of physics-based watershed models that explicitly describe surface (e.g., St. Venant equations for channel flow) and subsurface (e.g., Richards and/or Darcy equation-based) processes. The PAWS (Processbased Adaptive Watershed Simulator) model used in this thesis uses a coupling method to link surface and subsurface processes efficiently in order to facilitate the application of process-based models to large watersheds. The model has been tested extensively using available analytical solutions, laboratory data, streamflow, site-measured and satellitebased observations for medium and large watersheds [Shen et al., 2013; Shen and Phanikumar, 2010]. This chapter is organized as follows. After briefly describing the PAWS model and data sources, two of the largest watersheds in the State of Michigan, the Grand River and the Saginaw Bay watersheds, are used to examine water budgets and storage components in this paper over a five year period (2003-2007). Stream flow and groundwater heads are compared to the observed data from USGS gauges. Simulated soil moisture and soil temperature fields are also compared to the available monitoring data. In addition, evapotranspiration (ET) estimates based on MODIS data are compared with model-based ET estimates and comparisons are shown as both time-series of spatial averages and 2D spatial fields. PAWS model descriptions of storage changes are then compared with gravity-based satellite measurements of storage change (GRACE data reported as water thickness anomaly). A statistical linear correlation analysis is also performed to examine 9 the relation of the hydrologic budget components to the physical characteristics of the watersheds including land use and land cover, and soils type. An analysis of water budgets based on the Budyko curve framework is then reported to understand the long-term relationship between key hydrologic variables. 2.2. Methodology 2.2.1. The PAWS model The major hydrologic processes simulated by PAWS [Shen and Phanikumar, 2010] and the solution schemes are given in Table 1. The discretized differential equations for overland, soil water and groundwater flows are solved using structured grids. Each model river is associated with a real-world, named stream. The river network, which is extracted from National Hydrography Dataset (NHD), is explicitly modeled. A flexible strategy that computes the exchange fluxes between land grid and river cells allows river and land domains to be discretized independently. Two sub-domains - the ponding water layer and the flow domain are used to separate the water on the ground. Water can only move between cells (including river cells) in the flow domain, while the infiltration and evaporation happen in the ponding layer. Runoff from ponding layer to the flow domain is estimated based on Manning's equation. A layer of snow, which is quantified by the snow cover fraction and Snow Water Equivalent (SWE), may also exist on the ground depending on the climatic condition. A mosaic approach is used to model the sub-grid heterogeneity within each land cell. Several generic plant functional types (PFT) [Thornton and 10 Zimmermann, 2007] represent model classes and are used to re-classify the land use data. Rainfall and snowmelt are joined and add to the ponding layer after canopy interception. In PAWS, for each land cell there is one soil column modeled. Water and energy can be exchanged between the soil column and the unconfined aquifer. Through a layer of aquitard, water can percolate further down to deeper aquifers and contribute to or gain from river water via riverbed materials. Water can also be directly extracted from the aquifer layers by pumping activities. Table 2.1 Processes and governing equations of the hydrological model PAWS Process Governing Equations Snowfall Accumulation and melting SNICAR model (Toon et al., 1989) in CLM 4.0 Canopy interception Bucket model in CLM, storage capacity related to total Leaf Area Index (tLAI) Depression Storage Specific depth Runoff Manning's formula + Kinematic wave formulation + Coupled to Richards equation Infiltration / Exfiltration Processes coupled to Richards equation Overland flow 2D Diffusive wave equation Overland / Channel Exchange Weir formulation Channel network Dynamic wave or diffusive wave Evapotranspiration Penman Monteith + Root extraction Soil Moisture Richards equation Lateral Groundwater Flow Quasi-3D Recharge / Discharge (Vadose zone / Groundwater interaction) Non-iterative coupling inside Richards equation Stream / Groundwater interaction Conductance / Leakance Vegetation Processes CLM ver. 4.0 11 A simple growth cycle formulation for vegetation was applied in the original version of PAWS. The coupling of CLM (version 4) to PAWS [Shen et al., 2013] provides more detailed, process-based vegetation dynamics [Thornton and Zimmermann, 2007]. The soil hydrology and river routing routines in PAWS replace the corresponding procedures in CLM. Except due to snowmelt and freeze-thaw phase change, soil moisture states are unchanged in CLM. The resulting ET fluxes, soil water / ice content are passed back to PAWS for the calculations of flow processes after completing CLM calculations. More details about this coupling are available in [Shen et al., 2013]. 2.2.2. Description of Sites The Grand River (GR) watershed is the second largest watersheds in Michigan (Figure 2.1). It is located in the middle of Michigan's Lower Peninsula. The GR is the longest river in Michigan, which runs 420 km across several counties before emptying into Lake Michigan. 2 The watershed drains an area of 14,431 km including 18 counties and 158 townships. The average flow in the GR during the 5-year simulation period was 108 m3/s. The basin mainly contains flat terrain and many swamps and lakes. The Saginaw Bay (SB) watershed is the largest drainage basin in Michigan (Figure 2.1), draining approximately 15% of the total land area in the state. The watershed includes a part of 22 counties and contains the largest contiguous freshwater coastal wetland system in the United States. It is located in the east central portion of Michigan's Lower Peninsula. 12 Figure 2.1 Map of GR and SB watersheds. Elevation is shown as the color gradient. NHD streams, USGS gauges, groundwater stations and NCDC weather stations are shown. ( For interpretation of the references to color in this and all other figures, the reader is referred to the electronic version of this dissertation.) 13 The Saginaw River is a 35 km river formed by the confluence of tributaries at southeast of Saginaw. It flows northward into the Saginaw Bay of Lake Huron just northeast of Bay 2 City. The area of the SB watershed is 22,556 km and over half of the land use in the region is agricultural. The land use in both watersheds is shown in Figure 2.2. 2 For both watersheds, we use a 1 km cell size for spatial discretization which results in a grid size of 195 170 for GR and 201 210 for SB. 20 vertical layers are used to discretize the vadose zone and to solve the Richards equation. Two layers are used for the fullysaturated groundwater domain – the top layer represents the unconfined aquifer while the bottom layer represents the confined aquifer. The discreization resolution for the river cells is 1000m also. More details about the discretization and additional comparisons with data are reported in [Shen et al., 2013]. For temporal discretization, the base time step applied for the whole watershed is 1 hour, 10 minutes for the over-land flow and channel flow, 24 hours for the saturated groundwater flow, and adaptive time steps ranging from 10 minutes to 1 hour are applied for the unsaturated zone 14 Figure 2.2 Land Use and Land Cover map for GR and SB watersheds. 15 2.2.3. Data Sources Various data sources are needed for process-based models to simulate different hydrologic components. 30m resolution National Elevation Dataset (NED) from USGS (http://ned.usgs.gov) was used for the Digital Elevation Model (DEM) input for both watersheds to calculate the slope and for surface runoff routing computations. The daily based meteorological elements, pan evaporation, precipitation, snowfall, maximum temperature, minimum temperature, wind speed, and relative humidity acquired from National Climatic Data Center (NCDC, stations listed in Table 2.2) of the National Oceanic and Atmospheric Administration (NOAA) were used as weather data inputs for the model. In our study regions there are 14 weather stations available. The nearest neighbor interpolation scheme was used for spatial interpolation of data. For land use and land cover (LULC) data we used the Integrated Forest Monitoring Assessment and Prescription (IFMAP) dataset that has a 30m resolution. Using the land use classification percentile values reported in the literature e.g., [Jia et al., 2001; Lu and Weng, 2006], we transformed the land use classification provided by the Michigan Department of Natural Resources [MDNR, 2010] into several different model classes. The eight main MDNR land use class coverage percentiles for both watersheds are listed in Table 2.3. Soil Survey Geographic (SSURGO) database from U.S. Department of Agriculture, Natural Resources Conservation Service was used as model input. For the channel flow network, we used National Hydrography Dataset (NHD) from U.S. Geological Survey (USGS). The stream flow and groundwater head data for comparison were observed by a number of USGS gauge stations (listed in Table 2.4) within the two watersheds. There are a few 16 stations with soil moisture and soil temperature data available from the Enviro-Weather Network (data available at http://www.enviroweather.msu.edu/homeMap.php) and these datasets were also used for model comparison. Table 2.2 NCDC weather stations Stations Number Stations Name Latitude Longitude 200146 Alma 43.38 -84.67 201818 Corunna 2NE 43.00 -84.07 202395 East Lansing 4 S 42.67 -84.48 202437 Eaton Rapids 42.52 -84.65 203333 Grand Rapids 42.88 -85.52 203429 Greenville 2 NNE 43.20 -85.25 203664 Hastings 42.65 -85.28 204078 Ionia 2 SSW 42.95 -85.08 204150 Jackson Reynolds FLD 42.27 -84.47 204320 Kent City 2 SW 43.20 -85.77 204641 Lansing Capital City AP 42.78 -84.58 204944 Lowell 42.93 -85.33 205712 Muskegon CO AP 43.17 -86.23 209006 Williamston 3NE 42.72 -84.25 17 Table 2.3 LULC coverage percentages in the Grand River and Saginaw Bay watersheds GR (%) SB (%) Needle leaf evergreen temperate tree 2.44 3.71 Broadleaf deciduous temperate tree 19.64 25.48 Broadleaf deciduous temperate shrub 0.63 1.91 C3 non-arctic grass 5.98 8.39 Crops 57.91 47.84 Impervious 5.30 3.73 Water 1.29 1.25 Wetland 6.81 7.69 Table 2.4 USGS gauging stations Gauge Number Gauge Name Latitude Longitude 04114000 Grand River at Portland 42.8564 -84.9122 04116000 Grand River at Ionia 42.9720 -85.0692 04117500 Thornapple River near Hastings 42.6159 -85.2364 04119000 Grand River at Grand Rapids 42.9645 -85.6764 04149000 Flint River near Fosters 43.3084 -83.9536 04155000 Pine River at Alma 43.3795 -84.6556 04155500 Pine River near Midland 43.5645 -84.3692 04157000 Saginaw River at Saginaw 43.4128 -83.9630 GR SB 18 2.2.4. Simulation details There are considerable uncertainties associated with parameters such as riverbed conductance, van Genuchten parameters for the soils, etc. There are more than 10 calibrated parameters which can be referred to Shen et al.[2013] similarly. We apply global multipliers (multiply the same value to all spatially varied variables in each grid cells) to adjust some parameters within a certain range to improve match with observations while maintaining the spatial structure. Typical parameter values used in the model are described T in the st S (t )   t 0   1 N i 1 Pi (t )  ETi (t )  Qout i (t ) udy of Shen and Phanikumar N [2010] for the Red Cedar River watershed, a sub-watershed of the GR watershed. For the SB watershed, we did not adjust any parameters. Instead all model parameter multipliers for the GR watershed were used for the SB watershed to assess model performance and to test whether parameters from a geographically similar watershed will produce acceptable 2 results. The model performance is evaluated by coefficient of determination or R , Nash_sutcliff coefficient, absolute percent bias and root mean squared error. 2.2.5. Storage and Water Budget Calculations To estimate storage for the two watersheds based on water balance components, we use two different equations – the first equation is based on the water balance equation for large watersheds while the second equation is written to facilitate comparisons with GRACE data. The following hydrologic mass-balance equation, appropriate for large watersheds and for long periods of time [Gupta, 2008] was used: 19 T S (t )  t 0  1 N i 1 Pi (t )  ETi (t )  Qout i (t ) N  (2.1) where S is an integrated measure of the change of watershed storage, P is precipitation, Qout is river outflow, and ET is evapotranspiration. The equation is used to compute cumulative storage over the entire watershed (i=1 to N denote the model grid cells) over the period of integration (T). Equation 1 emphasizes the water balance near the land surface although also contains the groundwater contribution to streams. The other equation for storage emphasizes the subsurface components while including water on the land surface as well. Equation 2 is used to compute storage (water thickness anomaly) for comparison with GRACE data: S (t )  GW (t )  SW (t )  (GWU  GWC  VDZ )  (OVN  CS  SWE  rh) (2.2) where GW represents the water stored in the subsurface that includes the unconfined (GWU) and confined (GWC) aquifers and water in the vadose zone (VDZ); OVN represents all surface water components including the ponding domain, low land and depression storage and overland flow components; CS denotes canopy storage; SWE is snow water equivalent and rh represents channel flow. 20 2.3. Results and Discussion 2.3.1. Stream flow and groundwater head comparisons Data from eight USGS gaging stations were used for comparing with the 5-year model simulation results for the GR watershed. The locations of the USGS gauges are shown in 2 Figure 2.1 as cyan circles. The R values for all eight comparisons, shown in Figure 2.3, ranged from 0.57 to 0.72. The same parameter set estimated for the GR watershed simulation was applied for the 5-year SB watershed simulation. Stream flow data observed 2 at other four USGS gauge stations located in SB were used for comparison. The R values were still above 0.5 for several USGS gauges shown in the right panels of Figure 2.4 that show the comparison for the SB watershed. We note that, in all cases, the base flow is generally described very well. Possible reasons that may have contributed to the mismatch in streamflows include the underestimation of streamflow peaks during snowmelt periods due to measurement uncertainty of snowfall depths at stations, and some deficiencies of the snow energy balance routine of current CLM model as it does not account for the effect of snow cover fraction in heat exchange, which would result in errors in snow melt timing and we have found this to be a major suspect for some under-estimation of spring peaks. For USGS 4116000, 4119000 and 4114000, the peaks in spring 2007 and generally after that are not well captured. Examining the weather station records indicate that the precipitation records of two important stations, Williamston 3NE (209006) and Eaton Rapids (202437) suffered significant loss of data coverage after 2005. Eaton Rapids station has only 9 and 4 months of valid precipitation records for 2006 and 2007 respectively, while Williamston 21 has 9, 6 and 2 months for 2005, 2006 and 2007 respectively. The missing records are filled by borrowing records from neighboring stations, which still cannot prevent the decrease in model efficiency. Besides lack of weather stations with valid/reliable data records for the 2 simulation period also caused the reduction of R of stream flow comparison for the SB watershed. The availability of transient groundwater head observations is limited in our study area. Thus we only have one groundwater head comparison in the GR watershed 2 (Holt, Michigan) shown in Figure 2.5. The comparison (R is 0.8) represents a convincing match considering the fact that transient head comparisons are sensitive to local agricultural, municipal and industrial pumping which are difficult to quantify accurately. The simulated depths to water table (average from 2003 to 2007) are compared with the measured values from the WELLOGIC database in Figure 2.6. As can be seen from the spatial fields and the 2 R value 0.56, an overall agreement is noted between simulated and observed water table depths. 2.3.2. Soil moisture and soil temperature comparisons Soil temperature and soil moisture data are obtained from several Enviro-Weather stations within the two watersheds. We compared soil moisture and soil temperature at 10 cm depth for Enviro-Weather station located at Belding (Kropf farms, latitude 43.1125, longitude 85.3122, elevation 257m) in the GR watershed. The comparisons were presented in Figures 2.7 (a) and (b) respectively. It should be noted that our simulation results represent an average over the 1 km2 grid cell, whereas the observed data essentially represent a point measurement as data were collected using a soil moisture TDR (time domain reflectometry) 22 probe. In view of these differences in scales between the measurement and the model, the comparisons are considered acceptable. Figure 2.3 Simulated and observed streamflow comparison for GR watersheds at different 2 USGS gauges. USGS gauge station ID and R values are shown in the each subplot. 23 Figure 2.4 Simulated and observed streamflow comparison for SB watersheds at different 2 USGS gauges. USGS gauge station ID and R values are shown in the each subplot. 24 Figure 2.5 Simulated and observed transient groundwater head comparison at Holt, Michigan (Grand River watershed). 25 Figure 2.6 Observed versus simulated depth to water table. The observed data are obtained by kriging data in the WELLOGIC database. 26 Figure 2.7 Soil moisture (a) and soil temperature (b) comparison at Belding. 2.3.3. Evapotranspiration Comparisons with MODIS data One of the most important but difficult to measure hydrologic component is ET. The ET comparison data were downloaded from the MODIS Global Evaportranspiration Project (MOD16), which is a part of NASA/EOS project to estimate global terrestrial 27 evapotranspiration from earth land surface by using satellite remote sensing data (http://www.ntsg.umt.edu/project/mod16). The MOD16 global ET datasets are regular 12 2 km land surface ET datasets for the 109.03 million km global vegetated land areas at 8day, monthly and annual intervals. The dataset covers the time period 2000 – 2010. Model ET result was both temporally and spatially averaged through the whole computational domain and each month during two years. MODIS data were extracted for the same watershed area and time period. The comparison is presented in Figure 2.8 in the form of time series. The seasonal trend is clear: ET reaches its maximum value during 2 summer and a minimum value during winter as expected. The R value of the two-year comparison shown in Figure 2.8 is 0.85, which represent a good match with a similar trend between the observed and the simulated, except for the deviation at the higher values in July and August. The possible reason for the deviation is attributed to the different ET algorithms and water stress formulations used between CLM and MODIS estimates. PAWS+CLM uses a resistance approach to describe ET based on the two-big-leaf model [Dai et al., 2004], while the MODIS product is based on the Penman-Monteith formulation [Mu et al., 2011; Mu et al., 2007]. 28 Figure 2.8 Monthly simulated and MODIS ET data comparison. R 2 value is shown in the figure. A comparison (observed versus simulated) of the spatial maps of five-year-average ET for the GR and SB watersheds is shown in Figure 2.9. The mean annual average ET value is about 600 mm/year, which is comparable to the MODIS evapotranspiration data. The highest ET values about 1000 mm/year occur due to evaporation from open water bodies, while some of the lowest ET values (below 300 mm/year) are found for urban areas (there is no MODIS data available for urban areas, thus these areas are blanked out in the MODIS ET map) and the places with low sand percentage in soil. The ET heterogeneity due to LULC is discussed in the next section. 29 Figure 2.9 Five-year averaged spatial map of simulated and MODIS ET. 2.3.4. Budget analysis 2.3.4.1 Linear correlation analysis of annual average results To better understand the factors influencing the hydrologic water budgets and to understand the reasons for the spatial variability in key budget components such as ET, a linear correlation analysis was performed on different water budget components against the explanatory variables including LULC and soil types. The results are summarized in Table 2.5. In the table, is the pair-wise linear correlation coefficient, and p is the probability for 30 testing the null hypothesis of no correlation against the alternative hypothesis that there is a nonzero correlation. Statistical significance is assessed at an value of 0.05. In the GR watershed, ET values are positively correlated to forested areas strongly (   0.27 ), but negatively related to impervious land use, which represents urban areas (   0.35 ). ET is also strongly correlated to wetlands (   0.26 ). In the SB watershed, ET is also strongly and positively correlated to wetlands (   0.38 ), followed by forested areas (   0.35 ). The different ET correlation coefficients in the two watersheds are due to differences in the areas of different LULC in two watersheds; for example, there is a larger fraction of forested areas in the SB watershed than in the GR watershed (in Table 3, needleleaf evergreen temperate tree is 2.44% in GR and 3.72% in SB, broadleaf deciduous temperate tree is 19.64% in GR and 25.48% in SB). We also found that ET is very sensitive to soil types. Larger percentages of sand (   0.55 in GR and   0.63 in SB) and lower percentages of clay (   0.59 in GR and   0.69 in SB) in the soil result in higher ET although some earlier work [Istanbulluoglu et al., 2012; Wang et al., 2009] found a negative dependence between mean annual ET and basin sand coverage in some of the major Nebraska Sand Hills. The correlation analysis results not only support our understanding of ET dependence on soils and land use but also explain the 5-year average ET spatial maps very well. Correlation analysis results for other hydrologic components are also given in the same table. We also note that the two major runoff generation mechanisms in the two watersheds - infiltration excess (or Hortonian runoff) and saturation excess (or Dunne runoff) are 31 strongly correlated to the urban and forested land uses respectively. The saturation excess is negatively correlated with forested land use while infiltration excess is positively correlated with urban land use. For example, in the GR watershed, the saturation excess correlation coefficient is -0.23 with forested land use. In the SB watershed, the value is -0.40. The infiltration excess correlation coefficients for urban land use are 0.39 in the GR watershed and 0.34 in the SB watershed respectively. These results partially support the view that Table 2.5 Spearman correlation coefficients Land Use / Land Cover Watershe d GR ET SB GR Inf SB GR SatE SB GR InfE SB Recharg e GR SB ρ p ρ p ρ p ρ p ρ p ρ p ρ p ρ p ρ p ρ p Forest * 0.27 0.00 0.35 0.00 -0.37 0.00 -0.37 0.00 -0.23 0.00 -0.40 0.00 -0.10 0.00 -0.26 0.00 -0.17 0.00 -0.09 0.00 Grass * -0.03 0.00 -0.11 0.00 0.03 0.00 0.13 0.00 -0.17 0.00 -0.12 0.00 0.05 0.00 0.01 0.09 0.00 0.98 0.14 0.00 Crop s -0.21 0.00 -0.28 0.00 0.56 0.00 0.52 0.00 0.23 0.00 0.40 0.00 -0.02 0.02 0.16 0.00 0.23 0.00 0.11 0.00 Urba n -0.35 0.00 -0.34 0.00 0.06 0.00 0.05 0.00 0.09 0.00 0.35 0.00 0.39 0.00 0.34 0.00 0.09 0.00 0.00 0.81 Soils Wetlan d 0.26 0.00 0.38 0.00 -0.44 0.00 -0.49 0.00 0.13 0.00 -0.13 0.00 -0.13 0.00 -0.15 0.00 -0.26 0.00 -0.29 0.00 Soils Sand 0.55 0.00 0.63 0.00 -0.04 0.00 -0.06 0.00 -0.44 0.00 -0.57 0.00 -0.57 0.00 -0.57 0.00 -0.12 0.00 -0.13 0.00 forest* : broadleaf deciduous temperate tree; grass* : c3 non-arctic grass 32 Soils Clay -0.59 0.00 -0.69 0.00 0.09 0.00 0.12 0.00 0.41 0.00 0.55 0.00 0.58 0.00 0.59 0.00 0.18 0.00 0.18 0.00 forest vegetation maintains the soil water storage and minimizes runoff. Surface runoff is also strongly correlated to soil types, negatively with soil sand percent and positively with clay percent. Infiltration is very sensitive to the land use types other than different soil type percentile coverage (small  values for soil sand and soil clay in GR and SB in Table 2.5). It is negatively correlated to wetlands and forested areas (   0.44 and 0.49 for wetlands in GR and SB respectively,   0.37 for forest in both watersheds) but strongly positively correlated to crops (   0.56 and 0.52), and no correlation was found with urban areas (   0.06 and 0.05). 2.3.4.2 Storage and partition Annual averages of different water budget components for the GR and SB watersheds are presented in Table 2.6, where SUBM is sublimation; Qoc is surface runoff or overland flow; Dperc is groundwater recharge (deep percolation) and the other terms are described earlier. The components shown in Table 2.6 represent the major hydrologic components and, as in the case of runoff, represent different processes (infiltration excess, saturation excess, melting snow and so on) contributing to the same component, therefore they are not expected to add up to a 100%. In our study, the sum of SUBM, Qoc, Qgc, Dperc and ET is about 10% to 20% larger than P, which is reasonable because there may be other inputs /outputs or storage mechanisms in the watershed. We can calculate the ratio of each hydrologic component based on the amount of precipitation; for example, the annual average ET percentages are 68.9% and 65.8% for GR and SB respectively. These ET 33 percentages are comparable to the values reported for three watersheds in Illinois by Arnold and Allen [1996]. Table 2.6 Annual budgets (unit: mm/yr) P SUBM GR 875.5 13.1 GR (% P) 100.0 1.5 SB 886.8 16.0 SB(%P) 100.0 1.8 InfE SatE Dperc Qoc Qgc Qout ET 47.3 28.8 173.6 223.4 21.3 243.2 603.1 5.4 3.3 19.8 25.5 2.4 27.8 68.9 57.7 74.5 142.4 247.9 20.6 267.8 583.2 6.5 8.4 16.1 27.9 2.3 30.2 65.8 Figure 2.10 shows the cumulative storage in the GR watershed computed using equation 1. We notice that over the period of study, storage is increasing in both watersheds, although the rate of increase appears to be higher for the GR watershed. The study of mean water balances is usually performed over long term. Over a long period, positive and negative water storage variations tend to balance and the change in storage is usually ignored in mean annual and interannual water balance calculations. However, the interannual variability of storage change can be an important component in annual water balance during dry or wet years [Wang, 2011]. Time series plots of the anomaly of GRACE equivalent water thickness product (http://geoid.colorado.edu/grace/grace.php, hereafter termed as GRACE) and the anomaly of the summation of model variables are shown in Figure 2.11 for the GR watershed. The curves represent the anomaly of water storage 34 change, S '(t )  S (t )  S , where S is calculated using equation (2) and S denotes the average over the simulation period. As expected, S accrues in cold seasons, when there is no plant growth (small ET), and declines as we proceed to spring and summer. The figure vividly illustrates the seasonal energy cycle, together with individual precipitation events. The sampling of all GRACE grids is 1 degree in both latitude and longitude (approximately 111 km at the equator). Data from two adjacent GRACE grid cells are not independent because of the smoothing applied. For the GR watershed we calculated the average of 4 GRACE grids which enclose a rectangular area around the GR watershed. The differences of amplitude variation shown in the figure are due to the coarse resolution of GRACE data, the larger rectangular area of the GRACE dataset than the area covered for the GR watershed boundary in our model and the smoothing processes used in the GRACE data. In spite of these differences both the model result and the GRACE data show a similar increasing trend for watershed storage over the five-year period considered in this work. A comparison of Figures 2.10 and 2.11 also brings out important differences in the storage computed using equations (1) and (2) highlighting the importance of the subsurface domain for this climatic region. 35 Figure 2.10 Cumulative storage changes for GR (a) and SB (b) watersheds. Figure 2.11 Comparison of water thickness anomaly between gravity-based satellite measurements (GRACE data) and PAWS model descriptions. 36 2.3.4.3 Budyko curves Figure 2.12 shows the annual ET/P versus PET/P (hereafter termed the Budyko pairs) in the Budyko framework [Wang, 2011] for 4 sub-watersheds in the GR watershed and the whole GR watershed for a 13-year simulation period (1995 to 2007) with the annual evaporation estimated by two methods: (1) the annual storage change is assumed to be negligible, and ET is computed as: ET  P  Q ; (2) the model simulated values of ET are used. PET is the annual potential ET, P is the annual precipitation and Q is the annual outflow. The potential ET is calculated as the alfalfa (mature, 40cm canopy height)-based reference ET, using the Penman-Monteith equation. In this figure the horizontal straight line indicates the arid or water-limited conditions, while the 1 to 1 line indicates the humid or energy limited limit. The x-axis, annual potential evaporation to annual precipitation ratio, can also be considered as a dryness index. The annual ET to P ratio for all the sub-watershed is also within the range of 0.6 to 0.8, which is consistent with the above discussion. Because of the effects of annual water storage carryover, the data points (black dots in Figure 2.12) tend to move downward for small values of dryness index but move upward for large values of dryness index as compared with the case when annual carryover is ignored (yellow dots in Figure 2.12). Our results confirm the previous findings by [Istanbulluoglu et al., 2012] that the Budyko curve obtained by ignoring water storage effect, i.e., method (1), can be very misleading. The yellow diamonds in the figure are calculated using method (1) and as a result, ET/P has almost no correlation with PET/P. In contrast, when we use the simulated ET, the Budyko pairs generally fall in agreement with the worldwide Budyko Curve (WBC) given in [Fu, 1981], and there is a strong positive correlation. In this climatic region, the 37 points fall below the WBC, indicating a relatively large stagger between energy input and moisture availability: although the overall PET usually falls somewhere between 1 to 1.5 times of precipitation, ET is still limited by moisture in the summer. On the other hand, although there is ample moisture in fall and winter, there is very little PET in these seasons so that most of the precipitation exits the system as discharge. Figure 2.12 Budyko curve for the GR watershed. Figure 2.12 also provides insight into the inter-annual and spatial variability of the Budyko pairs among the 4 major subbasins of the Grand River watershed, namely, Upper Grand (upperGR), lower Grand (lowerGR), Maple and Thornapple. Within the basin, annual average PET/P can vary between 1 to 1.6 and ET/P can vary between 0.6 and 0.8. [Cheng et al., 2011] found that the inter-annual variability of the Budyko pairs for each basin can 38 be described by strong linear relationships. The slopes found for each subbasin of the 2 Grand River watershed and goodness of fit (R ) values are given in the Table 2.7. The correlations found in our study are weaker than the satellite-based estimates in [Cheng et al., 2 2011], indicating more complexity in the inter-annual variability, but the R values are nevertheless above 0.72 for all subbasins. This difference reflects the difference in ET estimation approaches. For vegetated areas, the satellite-based estimates in [Cheng et al., 2011] used modified Penman-Monteith with biome-specific NDVI observed by the satellite, whereas the relevant CLM algorithms use the conductance approach which considers fractions of wet/dry canopy, atmospheric stability, stomatal conductance (coupled to CO2 assimilation and in turn to vegetation and nutrient states). The inter-annual difference in nitrogen cycling could explain the difference in inter-annual variability of the Budyko pairs between PAWS+CLM and [Cheng et al., 2011]. It is difficult to judge which method is more accurate given the systematic biases associated with the satellite observations [Kalma et al., 2008] and the uncertainties with the translation from NDVI to ET. It would be interesting to probe if the subbasins in the watershed have statistically different responses to climatic input, in terms of locations of the Budyko pairs. As it is clear from Table 7, each subbasin has different average slopes. We calculate the average Mahalanobis distance D of each subbasin using the whole Grand River basin as a reference sample. The resulting values are DupperGR = 3.08, DlowerGR=2.40, Dmaple=2.85 and Dthornapple=2.012. 39 Table 2.7 Subbasin characteristics of the Grand River watershed in relation to Budyko pairs 0.363 0.275 0.317 0.277 Whole GR 0.320 0.767 0.722 0.755 0.738 0.737 2.4 2.85 2.0 3.08 -- 33.0 22.6 48.8 1.034 23.3 11.3 76.9 0.566 23.1 20.7 62.29 1.215 9.1 20 56.2 0.784 22.3 19.6 57.9 0.904 0.6862 0.6509 0.6910 0.6772 0.6886 LowerGR Slope R2 of linear fit to inter-annual variability Mahalanobis distance from whole GR Subsurface storage capacity (m) Forest percentage Agricultural percentage Average slope (%) Annual average ET(+)/ET(-) ratio Maple Thornapple UpperGR The Mahalanobis distance [Mahalanobis, 1936] is a measure of distance in a multidimensional parameter space which is similar to the Euclidean distance but takes into account the covariance among dimensions of the reference sample. Large D indicates more dissimilarity. The upper Grand River watershed is characteristically different from the GR as a whole, and this is rather apparent from the Figure: Upper GR has consistently higher ET/P values compared to other subbasins, given similar PET/P. On the other hand Maple tends to have lowest ET. The difference in responses is explained by their differences in geologic construct, land use, and potentially other factors. We calculate the subsurface storage capacity as the subbasin-average subsurface pore space from groundwater surface to bedrock top in Table 2.7. The storage of water in the bedrock is not considered here. UpperGR has the lowest subsurface storage capacity, generally the highest ET/P and the largest D (therefore behaving most differently from other subbasins). Less storage capacity, or shallower bucket depth, exposes more moisture for ET demand and produces higher actual ET. Therefore, the basin with shallower bedrock (lower subsurface storage) 40 generates comparatively higher ET. Land use also seems to play an important role in determining the responses. Maple and Thornapple have almost identical subsurface storage capacity but Thornapple has higher forest cover, which gives it higher ET/P. More analysis will be needed to determine the importance of each factor but it is out of the scope of this paper. Figure 2.13 shows the effect of groundwater flow in shifting the Budyko pairs. For each watershed, the areas are divided into recharging (+) versus discharging (-) areas based on the direction of the long-term average infiltration. (+) areas have a long term infiltration into the ground and (-) areas have water exfiltrating from the ground to the surface. In sharp contrast, discharging areas have far higher ET than the recharging areas, due to groundwater flow that transports moisture from upland to lowland areas. For certain low precipitation years ET/P in the discharging areas can exceed 1.1. Table 2.7 shows the annual average ET(+)/ET(-) for each subbasin. The difference in this ratio is well explained by the different topography. Maple has the flattest terrain and as a result it has the lowest ET(+)/ET(-) ratio. While Thornapple has the most hilly terrain and the largest ET(+)/ET(-) ratio. This is natural because in areas with more relief, groundwater plays a more significant role of moving water from upland areas to lowland areas. This analysis indicates that the spatial variability of the Budyko pairs can be significant due to groundwater flow. Therefore the Budyko pairs obtained for smaller watersheds may not be meaningful due to uncertain groundwater divide. 41 Figure 2.13 The effect of groundwater flow in shifting the Budyko pairs. 2.4. Conclusions A process-based hydrologic model has been applied to two of the largest watersheds in the State of Michigan and the water budgets and storage were analyzed. In addition to estimating the water budgets and examining their dependence on several key factors, this paper also demonstrated the application of a process-based hydrologic model to large watersheds in the Great Lakes region. The model was able to simulate different hydrologic components and states including surface runoff, channel flow, groundwater, ET, soil moisture, soil temperature and storage. The vegetation growth dynamics were also 42 simulated by coupling PAWS with the land surface model CLM. The model results match the observation well as discussed earlier. The linear correlation analysis suggests that the variations of different hydrologic components are sensitive to LULC and soil types. The process-based model has the ability to predict the magnitude of each hydrologic component and the partitioning of precipitation. Thus the water storage of the watershed can be derived from the partitioning. Although the water balance for long-duration of large watershed tends to be zero, the inter-annual variation still shows the clear seasonal variation. The annual averages of ET range from 600mm to 700mm, approximately 60% to 70% of the annual precipitation, which is comparable to estimates from the MODIS data and other studies in the literature. GRACE anomaly equivalent water thickness and the anomaly based on model results (equation 2) are in agreement showing the seasonal water storage variations and the annual increasing trend in storage for the GR watershed. Budyko curve analysis shows that the inter-annual and spatial variability of Budyko pairs are correlated to topography, subsurface water storage and LULC. 43 Chapter 3 Modeling Solute Transport at the Watershed Scale Using a Process-Based, Distributed Hydrologic Model with Application to Bacterial Fate and Transport Modeling There is a surge of interest in the application of mechanistic or process-based integrated hydrologic models (PBHMs) in recent years. While PBHMs have been applied to address issues of water quantity and to better understand linkages between hydrology and vegetation / land surface processes, applications of these models to address water quality issues are relatively limited. PBHMs have the ability to predict fluxes of bacteria and other contaminants exported from watersheds to downstream receiving water bodies by accurately describing the fate and transport of bacteria associated with overland, channel and subsurface flows in an integrated fashion. In this chapter, we first describe a general operator splitting strategy to describe watershed-scale solute transport. We then test transport modules in different hydrologic domains using available analytical solutions. This is followed by comparison of the numerical solutions with an analytical solution for groundwater transit times with interactions between surface and subsurface flows. Finally, we demonstrate the successful application of the model to bacterial fate and transport in the Red Cedar River watershed in Michigan. The watershed bacterial fate and transport model is expected to be useful for making near real-time predictions at marine and freshwater beaches. 44 3.1. Introduction The integrated hydrology including flows in channel networks, overland flow, and subsurface flows and interactions between different domains were simulated by the distributed hydrologic model PAWS. The computational efficiency of the PAWS model allows for long-term, large-scale simulations and makes it possible to simulate transport in large watersheds. An operator-splitting strategy combined with a Lagrangian particle transport modeling approach with reactions to describe transport in different hydrologic units was applied. The algorithms were tested extensively using analytical solutions (where available), data from plot-scale experiments before applying the bacterial fate and transport model to describe monitoring data from watershed in the Great Lakes region. 3.2. Methods 3.2.1. 2D Overland Flow Transport Equation A commonly accepted approximation of the St Venant equations (Singh, 1996) for overland flow is the kinematic wave equation which can be expressed as (3.1): h (uh) (vh)   S t x y (3.1) where h is the overland flow water depth [L], and v are the x  and y  direction water velocities [L/T], S is either the source term to represent precipitation or exfiltration from the subsurface, or a sink term to represent infiltration or evaporation. 45 The solute transport equation for overland flow can be expressed in a conservative form in terms of the solute flux as (3.2) (Deng et al., 2006): (Ch) (Cuh) (Cvh)   C    C      hDx   rDm (CS  C )    hDy t x y x  x  y  y  (3.2) where C stands for the cross-sectionally averaged solute concentration or the mass of 3 solutes per unit volume of runoff with a dimension of [M/L ] ;.. represents the solute concentration in the mixing soil zone; r is a mass transfer coefficient with a dimension of [1/T]; Dm is the mixing zone thickness which is proportional to flow depth; Dx and D y are the x  and y  direction diffusion coefficients. Equation (3.2) can be derived as shown below. (2)  C Dx h h C  (uh) C  (vh) C h C  uh C  vh  t t x x y y  2C x 2  Dx h C  2C h C  Dy h  Dy  rDm (CS  C ) x x y y y 2   h  (uh)  (vh)   C C C   2C  2C  C   u v   h   h  Dx 2  D y 2    x y  x y  x y   t  t   h C h C Dx  Dy  rDm (CS  C ) x x y y (3.3) Using the kinematic wave approximation which ignores the pressure gradient terms, we have h / x  0 and h / y  0 . Substituting equation (3.1) into the above equation and dividing both sides of the above equation by the flow depth h gives: 46 CS  C C C   2C  2C rDm  u v (CS  C )   Dx 2  D y 2  h  t x y  h x y C C C  2C  2C rDm S  u v  Dx  Dy  (CS  C )  C 2 2 t x y h h x y (3.4) 3.2.2. 1D Channel Flow Equation The one-dimensional channel transport equation is given by (Gunduz, 2004):  (Cr A)  C    (vACr )   ADL r t x x  x C g  Cr * nsed Dsed wr  qL 2CL 2  0 mr  *   kCr A  qL1CL1   (3.5) * * where Cr is the solute concentration within the channel, CL1 and CL 2 are lateral seepage and overland flow associated concentrations respectively, qL1 and qL 2 are lateral seepage and overland flows per channel length (positive for inflow and negative for outflow), A is the cross-sectional area, DL is longitudinal dispersion, Dsed is vertical dispersion coefficient in the sediment layer, k is decay coefficient, nsed is porosity of the sediment layer, C g is groundwater concentration, mr is thickness of the sediment layer, wr is river wetted perimeter, CO is overland domain concentration and  Cr ,  * CL1   Cg ,  qL1  0 qL1  0 C , * CL 2   r CO , 47 qL 2  0 qL 2  0 (3.6) * * This means that the values of CL1 and CL 2 change according to the direction of the lateral flow terms such that when lateral flow is towards the channel (i.e., inflow to the channel), these values take the associated concentrations coming from the groundwater (i.e., C g ) and overland flow domains (i.e., CO ) whereas they become the concentration in the channel (i.e., Cr ) when the lateral flow is away from the channel (i.e., outflow from the channel). Initial solute concentrations need to be specified along the one-dimensional channel domain: Cr ( x,0)  Cr 0 ( x) , where Cr 0 is the initial concentration distribution along the channel network. The channel transport model can accommodate several upstream boundary conditions and a single downstream boundary condition. At any upstream boundary, a specified concentration time series can be used as the boundary condition. The concentration time series is either available from continuous measurements (i.e., specified concentration) or from simple contaminant mass loading computation (i.e., specified mass flux): Cr (0, t )  Cru (t )  M (t ) Q(t ) (3.7) where M (t ) is the mass loading rate from some upstream source, Q(t ) is the river flow at the upstream boundary, Cru is the upstream boundary concentration time series. At the downstream boundary a zero concentration gradient is generally used as the boundary condition when the boundary x  Ld is far away from the contaminant zone: 48 ADL Cr |x  Ld  0 x (3.8) which states that advection dominates at the outflow. When boundary is not far, and the outflow concentration is measured, f the total mass flux, is specified. QC  ADL C |x  Ld  f x (3.9) Two or more channels intersecting within a channel network form a junction where internal boundary conditions are specified to satisfy the contaminant mass balance. The mass balance equation at a junction can be specified as: m  Ak J k  A0 J 0  k 1 dM dt (3.10) where m is the total number of inflowing channels to the junction, J k is the total mass th flux. Ak is the area at the end of the k inflowing channel to the junction, J 0 , A0 are total mass flux and the area at the beginning of the outflowing channel from the junction respectively. The change in mass storage within a junction is assumed to be negligible compared to the change in mass within a channel: m  Ak J k  A0 J 0  0 k 1 49 (3.11) Furthermore, the continuity of concentration at junctions guarantees that all the concentrations must be equal at junctions: (Cr )k  (Cr )0 (3.12) A forward characteristic is drawn from the unknown concentration value to the know time line. The concentration value corresponding to the point where the characteristic line cuts the known time line is called the foot of the characteristics and the concentration at this point is essentially the required concentration value. In general, a linear interpolation of the two neighboring concentration values at nodes N  1 and N is considered to be sufficiently accurate. If the position of the foot of the characteristic at the known time line is represented with x foot : x foot  xN  vt (3.13) Then the concentration value at this point is calculated using a linear interpolation function: x  vt  xN 1 j j j j C foot  (Cr ) N 1  N ((Cr ) N  (Cr ) N 1) xN 1 (3.14) In the above, the superscript j denotes time. Finally the concentration at the unknown time level (j+1) is equated to the calculated concentration value at the foot of the characteristic: j j (Cr ) N1  C foot 50 (3.15) Figure 3.1 Method of characteristics at junctions for advection step For a junction shown in above figure the junction mass balance is written as: j (Cr )1,1 C3  j j 1 j j 1 (Cr ) N1 1QNC1  (Cr ) N1 2 QNC 2 ,C , ,C , j 1 j 1 QNC1  QNC 2 , , (3.16) Once this concentration is found from the mass balance, it serves as the upstream boundary condition of the outflowing channel. 51 Figure 3.2 Treatment of junctions j j (Cr ) N1 1  (Cr )1,1 ,C C3 j j (Cr ) N1 2  (Cr )1,1 ,C C3 (3.17) The treatment of junctions is slightly different in the dispersion/decay/source/sink operator. A time-weighted mass balance can be written as: n j j   M in1 i  M out1   (3.18) i 1 where the total mass flux term for an inflowing or outflowing channel is given as the summation of the dispersive mass flux, decayed mass flux and the mass flux due to source/sink terms: 52 M in,out  ADL Cg  Cr Cr x x * x * x  kCr A  qL1CL1  nsed Dsed wr  qL 2CL 2 x 2 2 z 2 2 (3.19) As an example, the discretization for the first channel is shown below: M C1  Ab ( DL )b (Cr ) N ,C1  (Cr ) N 1,C1 * (qL1 ) N ,C1 (CL1 ) N ,C1 xN 1,C1 xN 1,C1 * (qL 2 ) N ,C1 (CL 2 ) N ,C1  k (Cr ) N ,C1 AN ,C1  (nsed Dsed )i xN 1,C1 2 (C g ) N ,C1  (Cr ) N ,C1 z 2 xN 1,C1  ( wr ) N ,C1 xN 1,C1 2 2 (3.20) 3.2.3. Interactions between overland flow and channel flow 3 The exchange volume M (m ) between land and channel is computed with the procedure originally outlined in (Shen and Phanikumar, 2010) and shown in Figure 3.3. 53 Figure 3.3 The process calculating river/land exchange. Here, Z Bank is the elevation of the bank (m), h is the depth of overland flow (m),  is the average free surface elevation of the cell (m), hr is the channel flow depth and ch is the stage of the channel (m). The first condition (1) is to determine the direction of the flow. If it is from the land to the river cell (   ch ) and there is water on the land (condition 2), we first attempt to compute the contribution explicitly using the diffusive wave equation: 2 L t  M ex  c h5/3 nm x 1/2 2 L t   max(ch , zBank )  c h5/3 nm (x / 2) 1/2 in which, Lc is the length of the channel segment that overlaps with the land cell (m). 54 (3.21) However, this flux cannot exceed the amount of currently available water on the land cell, Ma  A h (3.22) Also, there is an equilibrium state, at which the river stage will be the same as the land free surface elevation. Denoting this stage as  * , the mass transfer equation is written as: (* ch ) Ab  (  * ) A (3.23) where A is again the area of the land cell and Ab is the area of the river segment that spans this land cell. Then we can find the mass transfer that leads to  * :   A  ch Ab  ( ch ) Ab M E  ( *  ch ) Ab    ch  Ab  1  Ab / A  A  Ab  (3.24) Thus the exchange mass will be the minimum of M ex , M E , M a . On the other hand, if the river stage rises higher than the land free surface elevation (and also the bank elevation), flooding would occur, which is solved using a backward Euler implicit approach to enhance stability. The exchange mass solved by the implicit method is denoted as M im . Again using the diffusive wave formulation, the mass exchange function can be written in the form of two ordinary differential equations: A M 2L   d   im   c ch (ch  z Bank )5/3 dt t nm x / 2 d d Ab ch   A dt dt 55 (3.25) The exchange mass (also the exchange volume) M ( m3 ) or M im ( m3 ) between land and channel is computed from the above steps. Thus the exchange of contaminant can be simply calculated as: * qL 2CL 2  M CO , from overland to channel Lc * qL 2CL 2  M im Cr , from channel to land Lc 3.2.4. Transport in the Vadose zone Solute transport within the soil profile below the mixing layer is controlled by both infiltration and diffusion and can be described by the following advection-diffusion equation (Dong and Wang, 2013):  CS   CS    DS  ICS  t z  z  (3.26) where CS is the solute concentration in the soil water below the mixing layer,   b k p   , b is the dry soil bulk density, and k p is a constant partition coefficient. For non-adsorbed chemicals,  is equal to the soil moisture  , t is time, z is the vertical dimension, I is the infiltration rate in the soil, DS is the dispersivity of the chemical in the soil and is taken as the sum of the molecular diffusivity and the mechanical diffusion coefficient (Ahuja, 1990; Bear and Bachmat, 1990; Bresler, 1973). 56  ' DS   '   DS (3.27) ' DS  D0 aeb a and b are empirical coefficients, with typical values of a  0.005 , b  10 (Olson and 2 Kemper, 1961). D0 is the dispersivity [L /T]of solutes in free water. e is the exponential constant (2.7183);  is a constant (1),  is the velocity of flow in soil pores, and   is the 2 dispersion rate[L /T]. The  values for loam and sandy loam were 0.35 and 0.38. The boundary condition at the underlying soil surface is related to two stages during rainfall-runoff. d  (0, t )C (0, t )  dt  J (0, t )  ICS (0, t ) 0  t  tP (3.28) d  (0, t )C (0, t )  dt  J (0, t )  ICS (0, t )  IC (0, t ) tP  t where J is the diffusion rate of the solute from the soil below the mixing layer, which is described by Fick's law: J   DS CS . t P is the time from when the rainfall started to z when runoff took place. Thus before t P the initial concentration of CS is related to the solute diffusion in the soil and change of solute concentration in the runoff, as well as the infiltration rate. While after t P the solute concentration in the runoff will also infiltrate to the soil. 57 The initial condition for the solute concentration in the runoff is: C ( z,0)  CS 0 (3.29) where CS 0 is the solute concentration in the soil when rainfall begins. The vadose zone transport equation can be derived as follows:  CS   CS    DS  ICS  t z  z   CS C       CS   DS S    ICS  t t z  z  z CS  2CS C  DS CS I   CS   DS  CS I S 2 t t z z z z z C  2CS  D  C   I    S     CS   S  I  S  DS t  t z   z  z z 2 (3.30) The equation is discretized as shown below: n n 1 n 1 C n 1  CSi  CSi 1 C   I  n 1  DS   Si     CSi    I  Si 1  t 2z  t z   z  DS n 1 n n 1 CSi 1  2CSi 1  CSi 1 let fa  z 2 D  D    I    , fc   S  I  / 2z , fd  S , then , fb    t  t z   z  z 2  58 (3.31) n n n n 1 n 1 n 1 n n 1  faCSi 1  faCSi  fbCSi 1  fcCSi 1  fcCSi 1  fdCSi 1  2 fdCSi 1  fdCSi 1 n 1 n n 1 n    fc  fd  CSi 1    fa  fb  2 fd  CSi 1   fc  fd  CSi 1   faCSi  ( b k P   )  in 1  in    t t t t     '   ' DS   '   DS   z z z   (3.32) (3.33)   D ' S z    ( D0aeb ) (eb )  '   '  D0 a z z z z     (b )   '  D0aeb  '  D0aeb b z z z z     '  D0abeb z z   ' in 1  in 1 1 1 2z  n 1  n 1   n 1 i 1 i 1   D0 ab(eb ) i 2z (3.34) 3.2.5. Two-dimensional groundwater transport equation A general partial differential equation that governs two-dimensional solute transport in the fully saturated groundwater domain, considering advection, dispersion, fluid sinks/sources is as follows (Zheng and Bennett, 2002): 59 C   C    C    C    C     Dxx     Dyy     Dyx     Dxy  t x  x  x  y  y  y  y  x     q x C   q y C  qs C s x y R   (3.35) Where R  retardation factor (dimensionless) -3 . C  .dissolved concentration (ML ) -1 qx , q y  Darcy velocity (LT ) 2 -1 Dxx , Dxy , Dyy , Dyx  dispersion coefficient tensor (L T ) -1 qs  flow rate of a fluid source or sink per unit aquifer volume (T ) -3 Cs  concentration of the fluid source or sink flux (ML )   porosity (dimensionless) 2 v2 vx y Dxx   L  T v v Dxy  Dyx   L  T  60 vx v y v (3.36) (3.37) Dyy   L v2 y v2  T x v v (3.38)  L is the longitudinal dispersivity, a property of the porous medium describing dispersive transport in the direction of flow; T is the transverse dispersivity, a similar property of the porous medium describing dispersive transport normal to the direction of flow; v is the qy q 2 magnitude of the seepage velocity, v  vx  v 2 ; vx  x , v y  . The cross terms Dxy y   and D yx in the above equation cannot in general be eliminated by a particular alignment of the coordinate axes, as is normally done with the cross terms of the hydraulic conductivity tensor in flow simulation. Unless the flow is unidirectional, all components of the dispersion coefficient tensor must be retained. The source and sink term qsCs represents the solute exchange between soil vadoze zone and the first layer of the groundwater, and the solute exchange with groundwater at the bed of river channels, as well as the solute exchange between confined and unconfined aquifer. 3.2.6. Particle tracking scheme and mass balance The operator splitting method is applied to solve all the advection-dispersion-reaction transport equations in different domain. The Lagrangian particle tracking scheme is first used to solve the advection term while the rest dispersion and reaction terms are solved by finite difference method. Most numerical methods for solving the advection-dispersion equation can be classified as Eulerian, Lagrangian, or mixed Eulerian-Lagrangian (Baptista, 61 1987; Neuman, 1984). Eulerian methods offer the advantage and convenience of a fixed grid, are generally mass conservative, and handle dispersion-dominated problems both accurately and efficiently. However, for advection-dominated problems that exist in many field situations, an Eulerian method is susceptible to excessive numerical dispersion or artificial oscillation (Anderson and Cherry, 1979; Pinder and Gray, 1977). These types of error may be limited by using a sufficiently fine spatial grid and small time steps; however, the computational effort required for a field-scale problem, especially for large watershed simulations may become prohibitive. Lagrangian methods, on the other hand, provide an accurate and efficient solution to advection-dominated problems by essentially eliminating numerical dispersion (Kinzelbach, 1986; Prickett et al., 1981; Tompson and Gelhar, 1990). However, the velocity interpolation needed in particle tracking can also result in local mass-balance errors (LaBolle et al., 1996). When Eulerian methods are used, the time-step restriction associated with the advection step tends to be fairly restrictive. This limitation can be overcome using a Lagrangian method; therefore we use a Lagrangian method for the advection step in this work since our interest is in solute transport in large watersheds. If fluid density is constant, the pathlines of solutes under advection alone coincide with the pathlines of the flow, and are governed by the following equation: dp  v  p, t  dt 62 (3.39) where the unit vector p  xi  yj  zk is the position vector, and v  vxi  v y j  vz k is the velocity vector, a function of particle position and time. The solution of this first-order ordinary differential equation for particle location at any time t can be expressed as: t p  t   p  t0    v  p, t  dt t0 (3.40) where p  t0  is the position at initial time t0 . The solution of the integration equation requires the evaluation of the velocity field ( v ) at an arbitrary point (x, y, z) and any given time t. At any point inside cell (i, j), the x component of the velocity can be estimated by linear interpolation between the components at the two cell interfaces normal to the x direction: vx  1   x  vx1   xvx 2 where  x  (3.41) x  x1 , vx1 and vx 2 are the x direction velocity on the cell interfaces, x , x1 and x x2 are the x position of the particle and the two cell interfaces respectively, x is the cell x direction size. The same calculation applies to y and z direction velocity. The Euler method to numerically solve the integration equation is: pn 1  pn  v  pn  t 63 (3.42) where p n1 is the particle location at the new time level, n  1 ; p n is the particle location at old time level, n and t is the particle tracking step size. Euler method is only accurate to the first order thus the step size t must be sufficiently small for the method to be accurate. A fourth order Runge-Kutta method was applied in the present work: pn 1  pn   vp1  2vp2  2vp3  vp4  t 6 (3.43) where vp1 , vp2 , vp3 and vp4 are the velocity components at points p1 , p 2 , p3 and p 4 , respectively and t is the time step. The velocity components at these points are interpolated from the velocities at the cell interfaces once the positions of the four points are defined. At the beginning of a particle tracking step (time = tn ), the particle is located at p n . The velocity components v n , or vp1 are first computed using the velocity interpolation scheme. The coordinates for the first trial midpoint, p 2 , are computed using the velocities at p1 : p 2  p1  vp1 t 2 (3.44) where p1  pn . Again, the velocity components at p 2 are obtained through interpolation and used to calculate the coordinates for the second trial midpoint, p3 : p3  p1  vp2 64 t 2 (3.45) And finally the coordinates for the trial endpoint, p 4 , are computed using the velocity components obtained by interpolation for p3 : p4  p1  vp3 t (3.46) The velocity components at p 4 are again obtained through interpolation. The velocity components at p1 through p 4 are substituted into the fourth Runge-Kutta scheme equation to obtain the coordinates for the particle at the beginning of the next particle tracking step ( tn 1  tn  t ). Assuming that the number of particles in cell m at time t is NPm , the concentration for the cell at time t can be estimated as Cm  t   NPm  p 1 fp (3.47) where f p is the mass fraction represented by the pth particle. Mixed Eulerian-Lagrangian methods attempt to combine the advantages of the Eulerian and Lagrangian methods by solving the advection term with a Lagrangian approach (particle tracking) and the dispersion and other terms with an Eulerian approach. Some commonly used mixed Eulerian-Lagrangian procedures, such as the method of characteristics, do not guarantee mass conservation. Depending on the use of different Lagrangian techniques to approximate the advection term, various mixed Eulerian-Lagrangian methods may be 65 loosely grouped as the forward-tracking method of characteristics (MOC), the backwardtracking modified method of characteristics (MMOC) and a combination of these two methods, which tracks mass associated with fluid volumes to conserve mass locally and globally (Zheng and Bennett, 2002). There are two approaches for placement of initial particles in the particle tracking procedure: uniform and dynamic. The uniform approach places a fixed number of particles per model cell uniformly over the entire grid while the dynamic approach places a variable number of particles per model cell according to a given criterion. Whereas the uniform approach inserts and maintains an approximately uniform particle distribution throughout the model domain, the dynamic approach adjusts the distribution of moving particles repeatedly during the simulation, adapting to the changing nature of the concentration field. In many practical modeling applications, the contaminant plumes may occupy only a small fraction of the computational domain, and strong concentration gradients exist only at sharp fronts. In these cases, the total number of particles used under dynamic particle allocation is much smaller than that required in the uniform particle allocation approach, thereby dramatically increasing the efficiency of the MOC with little loss in accuracy. The MOC has the advantage of being virtually free of numerical dispersion and at the same time leads to a concentration solution that is less discrete than that obtained from the random-walk method. The major drawback of the MOC is that it is not efficient computationally when it is necessary to track a large number of moving particles, and like some other mixed Eulerian-Lagrangian solution techniques it can also lead to large mass-balance discrepancies under certain situations. 66 The MMOC is similar to the MOC except in the treatment of the advection term. Unlike the MOC which tracks a large number of moving particles forward in time and keeps track of the concentration and position of each particle, the MMOC places one particle at the nodal point of the fixed grid at each new time level. The particle is tracked backward to find its position at the old time level. The concentration associated with that position is used to approximate the concentration at intermediate time level. Therefore the MMOC is normally more computationally efficient than the MOC, and much less computer memory. The MOC and the MMOC can be combined such that when sharp concentration fronts are present the advection term is solved by an MOC technique through the use of moving particles dynamically distributed around each front, while away from such fronts the advection term is solved by an MMOC technique. By selecting an appropriate criterion for controlling the switch between the MOC and MMOC schemes, the adaptive procedure can provide accurate solutions with virtually no numerical dispersion. In the present work a combination of MOC and MMOC methods (called the Hybrid MOC or HMOC, Zheng, C., Bennett, G.D., 2002) are used. 3.2.7. Microorganisms release from manure and loss rate One of the major sources of bacteria is animal manure. In watersheds where CAFO (Concentrated Animal Feeding Operation) centers are located, it is important to quantify the production rate of bacteria associated with animal manure. We estimated the manure production using probability distributions based on the animal types and numbers based on literature values. The probability distribution representing animal shedding intensity in log 67 number of bacterium per gram of fresh weight manure was estimated following the gamma distribution (Medema et al., 2001): P( zi )  (i , i ) (3.48) where  i and i are parameters for the gamma distribution that model fecal shedding intensity, estimated from published studies for a given census category.  i and i were estimated using the method of moments (Dougherty, 1990) as: n 1 2 S X n i  ,   n 1 2 i X S n 2 (3.49) where S 2  sample variance, X  sample mean, and n  number of sample points. The parameters of the gamma distribution for different animal types were summarized in (Dorner, 2004). Manure production rates were estimated “as excreted”. Values have been modified from rates on a volumetric basis to a mass basis based on an assumed fresh manure density of approximately 1 kg/L (Dorner, 2004). A Monte Carlo method is used to estimate the probability distribution of the total daily bacterium production. Each iteration of the simulation provided a single estimate for the total daily bacterium production in the probability distribution. The number of bacteria shed by a single animal was taken from a random sampling of a gamma distribution. Then given the number of animals in the region, using daily manure production rates, the daily pathogen production was calculated by summing the number of bacteria shed by all animals. 68 The total daily bacteria production rate was obtained by summing the daily bacterium production from all individual census categories. The solution and the number of iterations necessary to integrate the probability distributions is based on the stabilization of the mean and standard deviation of the log-transformed data. Several models have been proposed and tested to simulate the release of manure-borne microorganisms. The assumption of the exponential release is frequently used in watershed scale models such as SWAT (Soil and Water Assessment Tool) (Sadeghi and Arnold, 2002): M R  M S k1Q (3.50) and HSPF (Hydrological Simulation Program – FORTRAN) (Bicknell et al., 1997): M R  M S 1  exp(k2Q) (3.51) Here M R is the count of microorganisms released during a runoff event; M S is the count of microorganisms in the manure storage layer of soil in the beginning of the runoff event; Q is the runoff yield during the runoff event, cm; k1 and k2 are the release rate parameters, cm-1. The research at finer scales reveals more complex dependencies of the released microorganisms on time or on runoff yield than the previous two models can predict. These dependencies were successfully simulated with two-parametric models; one such model was derived from the work of (Vadas et al., 2004) as: 69  q M r  M m a( w t )b Md (3.52) and another model was modified after (Bradford and Schijven, 2002):  11/    M r  M m (1   t )      (3.53) M r is the count of microorganisms released during the runoff event; M m is the count of microorganisms in the surface-applied manure before the runoff event; M d is the mass of -2 - the applied manure per unit area, g(dry mass)cm ;  w is the water density equal to 1 g cm 3 -1 -1 ; q is the runoff rate, cm h ;  , h , and a , b ,  are dimensionless fitting parameters. The exponential release model was applied and k2 was found to be 0.69315 cm -1 assuming that a 1 cm runoff removed 50% bacteria after a 1h rainfall event. The total loss rate for bacteria can be represented as (Chapra, 1997): ' kb  kb1  kbi  kbs (3.54) i where kb  total loss rate ( d 1 ), kb1  base mortality rate ( d 1 ), kbi  loss rate due to solar radiation ( d 1 ), kbs  settling loss rate ( d 1 ). The following equation can be used to calculate a base mortality rate for total coliforms (Mancini, 1978; Thomann and Mueller, 1997): 70 kb1  (0.8  0.006Ps )1.07T  20 (3.55) where Ps  percent sea water. Thus this formulation assumes a freshwater loss rate of 0.8 d-1. The total loss is then modified to account for temperature effects by including correction factor for temperature, T .The bacteria loss due to the effects of sunlight can be represented as (Thomann and Mueller, 1997): kbi   I (3.56) where   a proportionality constant, I  average light energy (ly hr-1). Using data from (Gameson and Gould, 1974; Thomann and Mueller, 1997) concluded that  is approximately unity. The exponential model suggested by (Chick, 1908): N  N0exp(t ) (3.57) is the predominant model used to describe the survival of manure-borne microorganisms in stored manure, soil, land applied manure, streams, and groundwater. Here N is the number of indicator bacteria at the time t , N 0 is the original number of indicator bacteria, and  is ' the die-off rate constant, while we use the total loss rate, kb , above instead. The fate and transport of bacteria at the watershed-scale are modeled using the above formulations for bacterial die-off combined with mechanistic routing methods in different domains including overland flow, channel flow and the vadose zone. In what follows, we present a systematic 71 evaluation of the transport models I different domains before applying the bacterial fate and transport model to a real watershed (the Red Cedar River watershed in Michigan). 3.3. Results and Discussion 3.3.1. One-dimensional channel transport The one-dimensional channel contaminant transport model is tested against the available analytical solutions in a single channel framework. The one-dimensional advection, dispersion, reaction (ADR) equation can be written as: Cr C  2Cr  v r  DL  kCr  0 t x x 2 (3.58) which defines not only the advective-dispersive transport but also the first order decay of contaminants in a simple channel. The analytical solution with constant coefficients can be expressed as (O'Loughlin and Bowmer, 1975): vx  vx   x  vt   2 D 1   x  vt   C0  2 DL 1  Cr  x, t   e erfc  erfc  e L  2 D t   2 D t  2  L  L       where   1  4 ,   kDL v2 , with the initial and boundary condition Cr ( x,0)  0 , Cr (0, t )  C0 . 72 (3.59) To test the model with analytical solution given above, a hypothetical rectangular channel domain is assumed so that in the channel there is steady uniform flow for all times. A constant velocity of 1m/s is prescribed through a 10000 m long rectangular channel. Initially, it is assumed the contaminant concentration is zero in the channel. A constant specified unit concentration is implemented at the upstream boundary of the channel. Decay coefficient k is set to zero for simplicity. Three different tests are performed with different dispersion coefficients representing (i) an essentially pure advection flow with 2 2 very low dispersion ( DL  1.0E  8 m /s); (ii) a medium dispersion ( DL  30 m /s) flow; 2 and, (iii) a high dispersion ( DL  100 m /s) flow. The numerical simulation and analytical solution are compared in Figure 3.4 (a, b, c). It can be seen that for the pure advection case the pattern and timing of the sharp front are properly captured by the numerical method at three different time steps. For the second test, a moderate amount of dispersion, the simulated results are either identical or very close to the analytical solution. In the third test case, the simulation is performed with a significantly high dispersion value. The dispersion coefficient is so high that the plume shows a 3000m deviation from its concentration center. Similar to the previous two test cases, an excellent agreement is obtained between the numerical and analytical solutions. 73 Figure 3.4 1D channel transport simulation with constant velocity but three different values of dispersion. 74 3.3.2. Two-dimensional overland flow transport For this test case, we consider two-dimensional transport from a “patch” source in a unidirectional flow field, assuming that the x axis is aligned with the unidirectional velocity vector ( v ) and no reaction or other source/ sink terms exist. The governing transport equation can be written as: C  2C  2C C  Dx  Dy v t x x 2 y 2 (3.60) with the initial condition C ( x, y,0)  0 and the following boundary conditions Left boundary: C (t ) C (0, y, t )   0  0  y0  y  y0 otherwise Right boundary: C (, y, t )  0 Upper boundary: C ( x, , t )  0 Lower boundary: C ( x, , t )  0 The analytical solution derived for this problem is as follows (Zheng and Bennett, 2002): 75   x  v 2   exp    4 Dx    3/2 t   x C ( x, y , t )  0   4  Dx    erfc  y  y0   erfc  y  y0   d   2 D y   2 D y         C0 (t   ) 1 (3.61) where C0 (t ) is the inflow concentration at the left boundary ( x  0 ) within the patch extending from  y0 to y0 . 76 Figure 3.5 Contour plots of 2D overland transport showing analytical and numerical solution at the same time step. 2 2 The coefficients are specified as: v  1 [L/T], Dx  1 [L /T], Dy  0.1 [L /T] and C0  1 3 [M/L ]. A comparison between the numerical and analytical solutions is show in Figure 3.5a and b in the form of 2D contour plots. The 2D contour plots comparison look almost identical in terms of the shape and the color distribution of the patches of the concentration gradient. The concentration distribution is much larger on the x axis than that on the y axis because there is no velocity along y axis and the x direction dispersion is much larger 77 than that on the y direction. To better compare the two solutions, we plot the concentration profiles at two different ways: First, in Figure 3.6a, we show the breakthrough curves at different downstream locations ( d  5,10, 20 [L]) from the center of the patch source ( y  0 [L]); Second, in Figure 3.6b, the longitudinal concentration distribution at different times ( t  10, 20,30 [T]) are shown from the center of the patch source ( y  0 [L]). At the downstream location d  5 m, the concentration breakthrough curve reaches a plateau after a certain time since we assume a constant concentration C0 at the source. The same trends are noted at different downstream locations ( d  10 m and d  20 m), but the further the downstream location away from the patch source, the smaller the maximum concentration and the later the equilibrium state reached due to the dispersion and travel time of transport. In the longitudinal concentration distribution plot when t  10 [T] the concentration front should reach d  10 [L] downstream, but due to dispersion the front actually extends to d  20 [L]. Good agreement is noted between the numerical and analytical solutions. 78 Figure 3.6 2D overland transport results showing the breakthrough curves at different downstream locations (a) and the longitudinal concentration distribution at different times (b) from the center of the patch source. 79 3.3.3. One-dimensional solute transport in a soil column For non-adsorbed solutes, transport in the vadose zone follows a simpler form if the soil moisture  and infiltration rate I are constant. If, in addition, the dispersivity DS is nonzero and a constant, then the advection dispersion equation can be applied: C  2C C D V t z z 2 (3.62) With D  DS 1 and V  I 1 . For the boundary and initial conditions C ( z,0)  Ci z0 C (0, t )  C0 t 0 The solution follows (van Genuchten and Alves, 1982; Warrick, 2003): C  Ci   C0  Ci  F  z, t    z  Vt   z  Vt   Vz    exp   erfc   F  z, t   0.5 erfc      2  Dt 0.5   2  Dt 0.5   D       80 (3.63) For the case of leaching of a solute that is already present in a profile, the invading concentration is less than the initial concentration and the previous results are essentially reversed. The solution is C  C0  Ci 1  F  z, t    (3.64) The solution for a finite “pulse” can also be found using the above equations. The resulting concentration profile would be a square wave increasing from zero to C0 , extending for the length of the pulse and dropping promptly to zero. In this case, the solution is by superposition C  Ci   C0  Ci  F  z, t  0  t  t0 C  Ci   C0  Ci  F  z, t   C0 F  z , t  t0  t  t0 (3.65) where C ( z,0)  Ci C C  0, t    0  Ci z0 0  t  t0 t  t0 The calculation of breakthrough curves based on equations 3.62 to 3.65 is demonstrated in Figure 3.7. The results are presented for a steady input, a removal (leaching), and a finite -1 2 -1 pulse. Parameter values are z  L  30 cm, V  0.79 cm min and D  0.5 cm min . The 81 first result (Figure 3.7a) is a plot of C / C0 as a function of t from 0 to 100 min. The results are plotted first with time and then with pore volumes as the abscissa in Figure 3.7 a and b. These first two plots contain the same information, and there is a one-to one correspondence using the pore volume to time relationship: pvol  Vt . Note that the value z of C / C0  0.5 at one pore volume. Figure 3.7c also shows a leaching example, that is, the results for an initial concentration C0 and an incoming solution concentration of zero (C). The result is reversed from before, at small values of pvol (or t ) the outlet concentration is C / C0  1 and then the value decreases to C / C0  0 . Also shown (panel d) is the solution for a pulse of finite duration (for a time from 0 to 0.5 pore volumes). The resulting concentration begins as for the continuous input, but it increases to a finite maximum of about 0.8 and then decreases back to zero. As shown in the calculations, both the leaching and the pulse solutions are simple extensions using the same F ( z, t ) as for the original cases. The numerical results match the analytical solutions very well in each case. 82 Figure 3.7 1D soil transport showing relative concentrations for step increase as a function of time (a), step increase as a function of pore volumes (b), set decrease as a function of pore volumes (c), and pulse as a function of pore volumes (d). We also examined the effects of D , V , and L on breakthrough curves. The effect of D is 2 -1 examined first by calculating results with D  0.05, 0.5, and 5 cm min . The velocity V is 0.79 cm min -1 and the length z  L is 30 cm, as in the previous calculations. The results are presented in Figure 3.8 a. The large D results in an earlier breakthrough and the 83 breakthrough is spread out over several pore volumes. The smaller D gives a much sharper breakthrough than either the large or the intermediate values; in fact, the smaller D gives results close to that of pure advection breakthrough. -1 Results of varying V (0.079, 0.79, and 7.9 cm min ) are given in Figure 3.8 b. In this case, the breakthrough for the slow velocity comes through first and is spread out more than for -1 the faster velocities. The highest velocity (7.9 cm min ) results in a very sharp breakthrough. These curves would appear quite differently if time were plotted on the abscissa; in that case, the higher V would result in the earliest breakthrough, which would appear far ahead of those for the slower velocities. Finally, the shortest length ( z  3 cm) gives a wide breakthrough, when plotted as a function of pore volumes as in Figure 3.8 c. Conversely, the longest length (300 cm) gives a much sharper pattern. Note that the curves in Figure 3.8a correspond exactly to the three curves in Figure 3.8b. This may be shown algebraically by rewriting F ( z, t ) of equation 3.63 in terms of the pore volumes (pvol) and the Peclet number Pe , defined by Pe  VL D (3.66) For z  L we find that F ( z, t ) can be written alternatively in the form   erfc 0.5P0.5 pvol 0.5  pvol 0.5   e      F  L, t   0.5  exp  Pe  erfc 0.5Pe0.5 pvol 0.5  pvol 0.5     84        (3.67) Thus, F ( L, t ) can be plotted as a function of Pe . A tenfold increase in V has the same effect as a tenfold increase in L or an order of magnitude decrease in D . This is consistent with the results in Figure 3.7. If we used t rather than pore volumes as the abscissa, then the three sets of curves (Figure 3.8 a, b and c) would not be the same. The numerical simulations are able to reproduce the results which match the analytical solutions very well in all cases. 85 Figure 3.8 1D soil transport results showing effect of apparent diffusion, velocity, and displacement length on relative concentration plotted as a function of pore volumes. 86 3.3.4. Groundwater transit time simulation We consider a vertical section through a saturated flow system within an unconfined aquifer underlain by an impervious boundary (Chesnaux et al., 2005) showing in Figure 3.9. -1 A uniform annual average rate of infiltration, W [LT ], is applied over the entire upper boundary of the aquifer. We consider only the case where W is positive, representing net recharge. The left-hand boundary is the imaginary impermeable ground water divide and flow discharges through the right-hand fixed-head boundary. Flow is also assumed to be at steady state, which implies that the numerical model converges to the state that the downstream discharge equals the ground water recharge and that the water table position does not change with time. The hydraulic head solution follows the Dupuit-Forchheimer ellipse: h  x    W 2 2 2 L  x  hL K (3.68) The analytical solution for the traveling time from the arbitrary departure point xi at the water table to the right exit boundary follows:       1 1 1 1  t  xi   ne  L 2    xi 2    ln  KW xi  L      2 where ne is porosity,   L2  KhL / W . 87  xi  L    xi2  xi2   1     1    (3.69) Figure 3.9 Conceptual model sketch for flow in a uniformly recharged, unconfined, horizontal aquifer used for developing the general analytical solution. The left-hand boundary is impermeable, and the flow discharges through the right-hand fixed-head boundary. The travel time is calculated between the two arbitrary points: xi (the departure point at the water table) and x. For this simplified case, the 2D groundwater equation is:   h    h   K xx    K yy   W  0 x  x  y  y  (3.70) where K xx and K yy are the horizontal and vertical components of the hydraulic conductivity. For determining travel times, we will also make use of the stream function solution, which is governed by the following equation: 88   1     1      0 x  K yy x  y  K xx y    (3.71) Lines of constant stream function represent streamlines that are everywhere tangent to the velocity field. A stream tube is a region bounded by two streamlines. At the left and bottom impermeable boundaries the first-type condition is  0  0 is used, while across the upper water table, the stream function increases down gradient according to  i   0   Wdx . S Downward across the right face, the Darcy flux q (calculated previously from the head solution) is used to define decreasing stream functions using  i   0   qdy . Following s the stream function solution, the transit time t for a fluid particle to travel a distance s 2 -1 within a stream tube defined by stream function interval  [L T ] can then be calculated using: t ne p  s  ds  s (3.72) where p is the width of the stream tube [L]. The integral term in the above equation is the area of the chosen stream tube along s. This area was computed by interpolating the nodal stream function values along the finite difference mesh. For this test case, the area integrations were computed based on a stream function interval of 8 108 m2/s (equivalent to 100 stream tubes). 89 The travel time can also be computed by using particle tracking scheme in our model. We released the particles at different positions right on the top of water table. The particles will finally flow to the right-hand exit boundary. Simply by summing up the number of time steps and multiplying the time step size gives the total travel time. We compare the analytical solution (Eq. 3.69), the numerical solution using stream function (Eq. 3.72), the numerical solution using particle tracking and the analytical solution developed by (Vogel, 1968) and (Gelhar and Wilson, 1974), which can be expressed as: L h n t  xi   L e ln   W  xi  90 (3.73） Figure 3.10 Groundwater travel time results ( K  104 m/s) showing model comparison of elapsed transit time (a), the corresponding water table profile and streamlines (b), and the velocity vector field (c) as function of distance. 91 Figure 3.11 Groundwater travel time results ( K  103 m/s) showing model comparison of elapsed transit time (a), and the corresponding water table profile and streamlines as function of distance. 92 Figure 3.12 Groundwater travel time results ( K  105 m/s) showing model comparison of elapsed transit time (a), and the corresponding water table profile and streamlines as function of distance. We compared the four methods with three different conductivity K values showing in Figures 3.10 through 3.12. The analytical solution matches the two numerical solutions very well; however, the error due to the assumption in the Vogel / Gelhar and Wilson model that the water table height is everywhere assumed to remain close to the down gradient boundary head can therefore be considerable, and especially with small K value 93 the water table raises nearly 20 meters. This test case shows that the model (particle tracking scheme) can accurately estimate the groundwater transit time and at some cases it is even more accurate than the numerical method derived from the commonly used streamfunction method. It approves that the particle tracking scheme is suitable for the groundwater transport problem especially with interaction between the soil and groundwater. 3.3.5. Plot scale simulations of manure-borne fecal coliforms Overland flow experiments with the manure-borne fecal coliforms (FC) were performed in October 2003 on a two-sided lysimeter of 21.34 m long and 13.2 m wide located at the Patuxent Wildlife Research Refuge (Beltsville, MD, Guber et al., 2009). The slope of the lysimeter surface was 20%. The average simulated rainfall rate measured was 5.8  1.9 cm -1 h applied to the area. The manure slurry was uniformly applied on the top of the plots in a -2 30-cm wide strip at the rate of 11.7 L m . Irrigation was started immediately after manure application. Rainfall was simulated for 1 h after the initiation of runoff. Twelve plots of 2 m width and 6 m length were established on both sides of the lysimeter, in the bare and the vegetated areas, resulting in total of four treatments: vegetated clay loam, bare clay loam, vegetated sandy loam, and bare sandy loam. Runoff volume and fecal coliform concentrations were measured for each runoff sample from each plot. More details about the experiment can be referred to (Guber et al., 2009). The overland flow and bacteria transport model was calibrated for the bare clay loam case by minimizing the root mean square error of observation and simulated result using an optimization algorithm. The 94 overland flow component was calibrated on data of the runoff rate time series. The boundary condition of h  0 at x  0 was used. The Manning roughness coefficient n was set to 0.035. The best values of the saturated hydraulic conductivity ( K ), and other parameters related to the soil properties were found by fitting the model equation to the runoff experimental data. The transport model was calibrated using FC breakthrough curve measure at the runoff plot outlet. The hydraulic component of the model after the calibration performed satisfactorily is shown in Figure 3.13a. The performance of the calibrated transport model is shown in Figure 3.13b. Both hydraulic component and transport results are satisfactory, especially the simulations of concentrations in the first portions of runoff where measured values of these concentrations are much larger than the modeled ones in the literature result (Figure 2 in Guber et al., 2009), however, our model captured this large concentration front quite well. 95 Figure 3.13 Plot scale simulation results showing the accumulative runoff (a) and the bacteria concentration (b) as function of time. 96 3.3.6. Modeling Escherichia coli fate and transport in the Red Cedar River watershed The model was applied to simulate the fate and transport of Escherichia coli (E. coli) in the Red Cedar River watershed (RCR). The total area of the RCR watershed is 1,169km2 (Figure 3.14). The watershed has a relatively low relief with the maximum elevation recorded as 324 m and a minimum of 249 m. 30 m resolution National Elevation dataset (NED) from the USGS was used for slope and surface runoff calculations. Land use and land cover (LULC) data was based on the 30 m resolution IFMAP (Integrated Forest Monitoring Assessment and Prescription) datasets from the Michigan Department of Natural Resources (MDNR, 2010). Depending on the grid size, a cell may possess a mixture of land use / land cover types in the PAWS model. LULC data were re-classified into model classes, which are represented by several generic plant types (called the representative or functional plant types). Meteorological data including precipitation, air temperature, relative humidity and wind speed were obtained from several National Climatic Data center (NCDC, 2010) and MAWN (Michigan Automated Weather network) (Enviro-weather, 2010) stations within the watershed. The Soil Survey Geographic (SSURGO) (Soil Survey Staff) database from the U.S. Department of Agriculture, Natural Resources Conservation Service was used for soil classification and for calculating the parameters associated with soil properties. The National Hydrography dataset (NHD) from USGS was used for channel network and for routing and flow computations. The stream flow and groundwater head data for comparison were obtained from a number of USGS gaging stations within the watershed. Hydraulic conductivity data were obtained by 97 processing the information in the MDEQ Wellogic database for the groundwater model to estimate local hydraulic conductivities, groundwater heads and thicknesses of the glacial drift layer. Spatial fields were obtained by kriging after removing noise in the data (Simard, 2007). Figure 3.14 Map of RCR watershed. Elevation is shown as color gradient. NHD streams, USGS gauges, WWTP, CAFO and NCDC weather stations are shown. There are three concentrated animal feeding operation (CAFO) centers and two waste water treatment plants (WWTP) (shown in Figure 3.14) located in the watershed which are 98 considered to be the most significant sources of E. coli in the watershed. The reported animal manure in the MSU CAFO is 9,145 tons per year, while the estimated average and maximum amounts are 9,385 and 10,727 tons per year based on the animal types and amounts. The estimated and reported values are 22,958 and 32,020 tons per year for the Kubiak CAFO. The numbers are 13,560 and 5,500 tons per year for the Ma Jo Lo CAFO. The mismatches between estimated and reported values of animal manure can either due to incorrect count for manure collection, un-counted manure loss during the collection process, or inaccurate animal amount reported. The daily discharge and fecal coliform bacteria concentration time series data for East Lansing and Williamston WWTP were also provided as model inputs. For the WWTP only fecal coliform data were available (obtained from MDEQ), as E. coli were not measured. A ratio of 0.6 was used to calculate the E. coli concentrations from fecal coliform measurements based on literature values (Patrick et al. 2003). 99 Figure 3.15 Comparison of observed and simulated E. coli concentration at three monitoring sites on Red Cedar river. 100 In Figure 3.15 the red lines are the simulation results at three different observation sites (S11, S-KZ and S_WEB_A) using both WWTP and CAFO inputs as sources which accounted for a majority of the bacterial pollution in the watershed. Comparing to the E. coli background concentration (reported values which are less than a hundred CFU/100ml for rural areas, Nouho et al. 2012) the bacteria concentration released from manure and in the discharge from WWTP are dominant. The black lines are the simulation results without soil sorption. We can see that with or without soil sorption, the simulation results are not changed significantly which means the soil sorption mechanism is not important at the watershed scale. This conclusion is also supported by the conclusion from Guber et al. (2009) at plot scale experiment. Since we are unable to clarify the sources from upstream and the upstream effects were not captured by the background values for different land use and land cover; therefore, the observed time series data (at site S_WEB_A) were used as upstream boundary condition. 101 Chapter 4 Modeling Nutrient Transport in Regional Great Lakes Watersheds 4.1. Introduction Surface water eutrophication resulting from nonpoint source phosphorus and nitrogen inputs is the most common water quality problem in the United States (USEPA, 1996). Point sources of pollution, including phosphorus from industrial discharges and municipal waste water treatment plants, are easily identified and quantified. Point source phosphorus loading to surface waters has been greatly reduced since the implementation of the 1972 Federal Clean Water Act. However, control of diffuse nonpoint source pollution has been less successful and is considered the main cause of eutrophication in lakes, streams and coastal areas in the United States. In the Great Lakes region, the amount of nutrients and organic material entering the lakes has increased with intensified urbanization and agriculture. Nutrient loading increased with the advent of phosphate detergents and inorganic fertilizers. Increased nutrients in the lakes stimulate the growth of green plants, including algae. In recent years, harmful algal blooms and drinking water related issues continue to be causes for concern highlighting the need for accurate transport models. It has also been found (Shen et al., 2013) that nitrogen significantly controls transpiration, through which it influences other hydrologic fluxes. In this chapter two models (Christiane et al., 2012; P. D’Odorico et al., 2003) that describe carbon-nitrogen cycles and carbon-phosphorus cycles in the soil were adapted and coupled to the PAWS+CLM, to describe the nutrients concentrations in the soil and use the 102 transport routing scheme to simulate the nitrate and dissolved phosphorus concentration observed in the streams in the Kalamazoo watershed in Michigan. The soil moisture controls several key fluxes that regulate phosphorus availability and the mineralization and immobilization fluxes of nitrogen, including the decomposition rate of litter and humus, the mortality rate of microbes, the plant uptake, and leaching losses. The influence of soil moisture dynamics on soil carbon, nitrogen and phosphorus cycles is modeled by coupling our PAWS+CLM model to a system of nonlinear differential equations that describe the temporal evolution of phosphorus and nitrogen in different pools over a wide range of time scales (i.e., from 10 minute time step as used in the overland flow simulation in PAWS to years of the whole simulation period). The soil carbon and nitrogen dynamics of the litter pool and humus pool were simulated in CLM and the soil moisture was calculated in PAWS, while the soil phosphorus and nitrogen models include inorganic and organic pools that are coupled with the soil carbon pool. The interaction between nutrient-cycles and hydrology has been recognized as important for a long time. However, the linkages are far from being completely understood. The nitrogen cycle needs to be coupled to the water balance that depends on the type of vegetation in a region, which, in turn, determines litter composition. The soil nitrogen cycle is closely related to the soil carbon budget and to the soil water content, which in turn control microbial activity, mineralization, leaching, and plant uptake. Phosphorus availability in soils can be one of the major factors limiting plant growth. Vegetation has adapted to the natural phosphorus limitation by efficiently cycling phosphorus stored in organic forms between the soil and vegetation. The soil microbial biomass has been found 103 to be important in determining the extent of release of phosphorus and in turn, the availability of this nutrient for vegetation uptake. Moreover, phosphorus and nitrogen are co-limited in several ecosystems worldwide, however, phosphorus is vulnerable to external actions by leaching and occlusion while nitrogen is more prone to renovation because of the anthropogenic emissions. 4.2. Carbon model The model for the carbon cycle is based on CLM (Thornton et al., 2007) and Porporato et al. (2003). The carbon cycle is modeled using different carbon pools as shown in Figure 4.1. The pools represent the main components of the system, the vertically averaged values of carbon concentrations over the active soil depth. All the state variables of the soil carbon, in terms of mass per unit volume of soil (e.g. grams of carbon per m3 of soil) are needed to characterize the system: CL Carbon concentration in the litter pool; CH Carbon concentration in the humus pool; CM Carbon concentration in the biomass pool. 104 The models of carbon in litter and humus pools in Porporato 2003 are replaced by CLM variables, the total litter carbon and total soil organic matter carbon, which is dominated by humus. The biomass pool is tracked using a variable as described in a later section. 4.2.1. Carbon in litter pool The dynamics of carbon in the litter is expressed as dCL  LF  MD  DECL dt (4.1) where LF is the rate at which litter is recycled back to the forest floor, MD is the rate at which carbon returns to the litter pool due to the death of microbes and DECL is the decomposition rate of leaf litter. DECL   f d  s  kd CM  CL   (4.2) where kd is the decomposition rate of the litter, the coefficient  is a nondimensional factor which specifies a possible reduction of the decomposition rate when the litter is lack of nutrient and the immobilization is deficient needed by the bacetria to integrate the nitrogen, f d ( s) is a nondimensional factor that depicts how soil moisture affect the decoposition and is expressed as: 105  s s  fc fd  s     s fc  s  s  s fc (4.3) s  s fc where s fc is the soil moisture at field capacity. 4.2.2. Carbon in humus pool Carbon enters the humus pool when microbes decompose litter ( DECL ), whereas carbon is lost from the humus pool when microbes decompose humic substances ( DECH ) and this carbon is incorporated into the microbial biomass pool. Therefore the equilibrium state of carbon in the humus pool is expressed using dCH  rh DECL  DECH dt (4.4) where DECH is the decomposition rate of the humus pool and is modeled similarly to the decomposition of leaf litter, DECH   f d  s  khCM  CH   (4.5) where kh being the decomposition rate of humus and  being the same nondimensional reduction factor. The coefficient rh is in the range 0.15 to 0.35, depending on the composition of the plant residues. 106 4.2.3. Carbon in microbial biomass pool The inputs to the microbial pool are due to the decomposition of litter ( DECL ) and the decomposition of humus ( DECH ). Because DECL is reduced by rh , (1  rh ) is the fraction of decomposing litter that goes into the microbial pool. Both a fraction of DECH and DECL are lost by means of respiration ( rr ). The fraction, rr ( 0  rr  1  rh ) defines the fraction of decomposed organic carbon that goes into respiration ( CO2 production), which is usually estimated in the 0.6 to 0.8 interval [Brady, 1996]. The change of carbon in the microbial biomass pool is expressed as: dCM  1  rh  rr  DECL  1  rr  DECH  MD dt (4.6) where MD is the mortality rate of microbial biomass: MD  Max[kmax fm (s), kopt ]CM (4.7) where kopt is the microbial death rate that happens under optimal condition, kmax is the maximum microbial death rate and f m ( s) is the control that soil moisture, s , has on the death rate, which is expressed as a linearly decreasing function up to the soil moisture at field capacity, s fc , and a linearly increasing function from s fc to saturation (i.e., s  1 ), 107  s fc  s   s fc fm  s     s  s fc 1  s fc  s  s fc (4.8) s  s fc Figure 4.1 Schematic representation of the soil carbon budget. Carbon pools are represented by boxes with bold text while fluxes between pools are represented by arrows with italicized text. Carbon stored in living biomass (as well as the flux due to FIRE) is not considered explicitly in this model. 4.3. Nitrogen model The soil nitrogen cycles represented in Figure 4.2 including nitrogen input from atmospheric nitrogen fixation and fertilizer, mineralization, immobilization process as 108 nitrification, ammonia volatilization, ammonium fixation and leaching. The model for the nitrogen cycle is based on CLM and Porporato et al. 2003. The nitrogen cycle is highly intertwined with the carbon cycles, showing in Figure 4.3, which represents the main components of the system, the vertically averaged values of nitrogen concentrations over the active soil depth. All the state variables of the soil nitrogen, in terms of mass per unit volume of soil (e.g. grams of nitrogen per m3 of soil) are needed to characterize the system: N L Organic nitrogen concentration in the litter pool; N H Organic nitrogen concentration in the humus pool; N M Organic nitrogen concentration in the biomass pool; N  Ammonium concentration in the soil; N  Nitrate concentration in the soil. The models of nitrogen in litter and humus pools in Porporato 2003 are also replaced by CLM variables, the total litter nitrogen and total soil organic matter nitrogen which is dominated by humus. 109 Figure 4.2 Systematic sketch of nitrogen cycles in the soil 110 Figure 4.3 Schematic representation of the components of the nitrogen cycles. 4.3.1. Nitrogen in litter pool The nitrogen balance in the litter pool is similar to the litter carbon equation, with each term divided by the C/N ratio of its respective pool, i.e. dN L DECL LF MD    dt (C / N ) LF (C / N ) M (C / N ) L (4.9) where (C / N ) L and (C / N ) M are the carbon nitrogen ratio in litter and microbial biomass pools respectively, (C / N ) LF is the C/N ratio of the added plant residues or fresh litter, which has been found to vary widely, ranging from 10 in legumes and young green leaves to more than 200 in sawdust [Brady, 1996]. This large variability may produce pronounced 111 changes in the C/N ratio of the litter pool, which has a very important role in regulating decomposition, immobilization, and mineralization. 4.3.2. Nitrogen in humus pool The nitrogen balance equation may be simply obtained by dividing Eq. (4.4) by the carbon nitrogen ratio in humus pool, (C / N ) H , i.e. dN H DECL DECH  rh  dt (C / N ) H (C / N ) H (4.10) This implies the assumption that the products of the humification process from litter have the same characteristics, and thus also the same C/N ratio, as the soil humus. As a consequence, the value of (C / N ) H remains constant in time, making the above equation redundant. 4.3.3. Nitrogen in microbial biomass pool Consequently the balance of nitrogen component in the biomass may be expressed as dN M  (C / N ) L  1  rh dt (C / N ) H   DECL DECH MD      (C / N ) L (C / N ) H (C / N ) M (4.11) The first two terms on the right hand side of equation above represent the incoming nitrogen from decomposition and do not contain rr because the respiration process only involves the carbon component. The third term is the output of nitrogen due to microbial 112 death, while the fourth term,  , takes into account the contribution due to either the net mineralization or to the immobilization. We will discuss this term in detail in the next section. 4.3.4. Nitrogen mineralization and immobilization rates The term  in Eq. (4.11) attains positive or negative values in relation to the difference  between the rate of gross mineralization and the total rate of immobilization of NH 4 ,    IMM N , and NO3 , IMM N , i.e.   MIN N  IMM N (4.12)   where IMM N  IMM N  IMM N . We only modeled the net mineralization and immobilization amounts as if they were mutually exclusive processes. Thus, when   0 , we assume that  MIN N  ,   IMM N  0 (4.13)  MIN N  0   IMM N   (4.14) while, when   0 , The switch between the two states is determined by the condition that (C / N ) M is constant in time, which is implemented analytically as 113 d (C / N ) M dCM dN M   (C / N ) M  0 dt dt dt (4.15) which, using Eqs. (4.6) and (4.11), yields  1  rr 1   DECH    (C / N ) H (C / N ) M     f d ( s)CM khCH     rh 1  rh  rr  1      DECL    (C / N ) L (C / N ) H (C N ) M   1  rr 1    (C / N ) H (C / N ) M   rh 1  rh  rr   1      kl CL     (C / N ) L (C / N ) H (C / N ) M    (4.16) Note that the above equation, and not Eq. (4.11), is the real dynamical equation to be associated to Eq. (4.6) for the biomass evolution. The previous condition makes Eq. (4.11) redundant, as the biomass C/N ratio is set constant. When the term in curly brackets of Eq. (4.16) is positive, net mineralization takes place, while no net immobilization occurs, as indicated by Eq. (4.13). In such conditions, humus and litter decomposition proceeds unrestricted and the parameter  is equal to 1. In the opposite case, when the term in curly brackets of Eq. (4.16) is negative, net mineralization is is halted and immobilization sets in Eq. (4.14). The latter is partitioned proportionally between ammonium and nitrate on the basis of their concentrations and according to two suitable coefficients, ki and ki , i.e. 114  ki N   IMM N ,  IMM N    ki N  ki N    ki N    IMM N  IMM N .  ki N   ki N   (4.17) The rate of immobilization may be limited by environmental factors, biomass concentration, and especially by insufficient mineral nitrogen. For this reason, we assume the existence of an upper bound for the rate of mineralization, i.e. IMM N  IMM max which after Eq. (4.16) becomes   1  rr  1   kh C H     (C / N ) H (C / N ) M   f d ( s )CM  rh 1  rh  rr 1 k C   l L  (C / N )  (C / N )  (C / N )  L H M       IMM max    (4.18) The maximum immobilization rate is assumed to depend on the concentration of the microbial biomass and soil moisture in a similar way to the decomposition process, IMM max  (ki N   ki N  ) f d (s)CM (4.19) If the immobilization rate was lower than IMM max , this means that the bacteria can meet their nitrogen requirement and decompose the organic matter at a potential rate. At this case the coefficient  is equal to 1. Otherwise the model reduces  to a value lower than 1, so that the immobilization rate is equal to IMM max . By imposing the equality in Eq. (4.18), 115 i.e. IMM N  IMM max , the value of  by which the decomposition rates must be reduced is obtained,   1  rr  1   kh C H    (C / N ) H (C / N ) M      ki N   ki N  /  rh 1  rh  rr 1  kl CL       (C / N ) L (C / N ) H (C / N ) M           (4.20) 4.3.5. Ammonium and Nitrate The balance of ammonium and nitrate in the soil can be modeled respectively as dN    MIN N  IMM N  NIT  LE   UP  dt (4.21) dN    NIT  IMM N  LE   UP  dt (4.22) The nitrification rate can be modeled as first order kinetics, i.e. NIT  f n (s)knCM N  (4.23) where the constant kn defines the rate of nitrification, and f n ( s) accounts for the soil moisture effects on nitrification. 116  s  s  fc f n ( s)    1 s 1  s fc  s  s fc (4.24) s  s fc For both ammonia and nitrate, leaching occurs when the nitrogen in the soil solution percolates below the root zone. It is thus simply proportional to the leakage in the vadose zone modeled in PAWS, LE   a  L( s )  N snZ r (4.25) The nondimensional coefficients a  , 0  a   1 , are the fractions of the dissolved ammonium and nitrate, respectively, and are related to the corresponding solubility coefficients. Since nitrate is a mobile ion, a  can be taken to be equal to one, while a  is much lower, as a large fraction of ammonium is absorbed by the soil matrix. The last terms of Eqs. (4.21) and (4.22), UP  and UP  , are the ammonium and nitrate plant uptake rates. For both ammonium and nitrate, passive and active uptake can be regarded as additive processes, i.e.   UP  UPP  UPa (4.26) The passive uptake can simply be modeled as proportional to the transpiration rate, T ( s) , which is modeled in PAWS, and to the nitrogen concentration in the soil solution, i.e. 117  UPp  a  T (s)  N snZ r (4.27) Three possible cases thus occur in the representation of the active component,  0    UPa   ku N      DEM  UPp  if  DEM   UPp  0, if  ku N   DEM   UPp , if  ku N   DEM   UPp  0 (4.28) where DEM  is a given plant demand, than which only if the passive uptake is lower we assume that the plant tries to compensate for the deficit with the active mechanism of uptake. The parameter ku ( d 1 ) describes the dependence of the diffusion process on the soil moisture level and can be modeled as ku  a Fs d snZr (4.29) where the term F is a rescaled diffusion coefficient, and d expresses the nonlinear dependence of the diffusion process on soil moisture. 4.4. Phosphorus model Phosphorus availability in soil plays an import role in agricultural system as well as in the natural ecosystem. Almost one third of the global arable land was evaluated that is 118 occupied by phosphorus limited soil. The rock phosphate can create phosphorus fertilizer to stimulate crop yields, but only available for short term phosphorus limitations. Vegetation could expeditiously promote phosphorus circulation in organic forms to reduce the limitation because of the biological and physical effect of soil moisture on phosphorus availability. In biological processes, soil moisture influences phosphorus availability by controlling the oxygen supply and affecting microbial activity and death via soil water. In physical processes, soil moisture leads to phosphorus losses because the soil erosion and runoff preferentially remove the phosphorus rich fractions and the plant uptake influenced by the level of water potential at the rooting zone and foliage section. Soil moisture is the hydrological variable that synthesizes the effects of the hydro-climatic conditions, soil properties and vegetation type on the phosphorus cycle. Any change of these factors could affect soil moisture dynamics, and in turn control phosphorus availability. We linked the models of soil moisture dynamics, soil carbon and soil phosphorus cycles over a wide range of time scales (i.e., from daily to seasonal to annual). The model for the phosphorus cycle is based on Christiane et al. (2012). The phosphorus cycle is also highly intertwined with the carbon cycle, as shown in Figure 4.4, which represents the main components of the system, the vertically averaged values of phosphorus concentrations over the active soil depth. All the state variables of the soil phosphorus, in terms of mass per unit volume of soil (e.g. grams of phosphorus per m3 of soil) are needed to characterize the system: P Organic phosphorus concentration in the plant; V 119 PL Phosphorus concentration in the litter pool; PH Phosphorus concentration in the humus pool; PM Phosphorus concentration in the biomass pool; PI Inorganic phosphorus concentration in the soil; POCC Occluded phosphorus concentration in the soil. Phosphorus is found in two major pools, organic and inorganic. The total inorganic phosphorus pool, which includes the available fraction and less available, occluded fraction can account for 30-70% of the total phosphorus found in the soil. External phosphorus inputs to the system (Figure 4.4) are due to atmospheric deposition (ATM); whereas, external phosphorus outputs are due to leaching (LE), erosion and runoff (ER), and fireinduced losses as smoke, air-borne ash particulate, and the removal of ashes and other plant residues by overland flow (FIRE). 120 Figure 4.4 Schematic representation of the soil and biomass phosphorus budget. Organic and Inorganic phosphorus pools are separated in two groups by boxes with bold text while fluxes between pools are represented by arrows with italicized text. 4.4.1. Organic phosphorus 4.4.1.1 Organic phosphorus stored in the living biomass pool ( P ) V Phosphorus stored in the living biomass pool includes live aboveground vegetation, belowground biomass such as roots, and mycorrhizae. We account for the changes in the amount of phosphorus that is stored in the plant as, 121 dP V  UP  LF (t ) dt (C / P) LF (4.30) where (C / P) LF is the average C/P ratio of fresh litter, which has been found to vary widely, ranging from 100 to 4100 [McGroddy, 2004; Stark 1972]. The resorption of nutrients prior to leaf senescence is one of the main factors that affect the C/P ratio of fresh litter. Another factor that affects the C/P ratio of fresh litter is the falling leaves. When the plant phosphorus demand can be met by passive uptake ( UPp ), this term is expressed as, PI   k pT ( S ) nZ s r  UPp   * k T ( S ) P  p nZ r s  PI  P* (4.31) otherwise where k p is a parameter accounting for the role of mycorrhizal hyphae in the movement of phosphate ions to the plant and PI is the concentration of PI in soil water for a soil with nZ r s a root depth, Z r , porosity, n , and relative soil moisture content, s , P* is the amount of PI required before it becomes limiting. UP not only consists of passive uptake, but also active uptake due to the role that enzymes play in increasing phosphorus availability. 122  DEM * ,  UP   UPp  UPa ,  where DEM *  k pT ( s) PI  P* (4.32) PI  P* P* . nZ r s There are two main strategies by which soil phosphorus in more recalcitrant forms is made available for active uptake through enzymatic release ( UPa ), (1) the dissolution of phosphorus containing minerals through a combination of soil acidification and the release of metal complexing agents (predominantly organic acid anions) ( ENZ I ) and (2) the enzymatic breakdown of organic P by releasing extracellular enzymes ( ENZO ).  0, DEM *  UPp  0   UPa   DEM *  UPp , 0  DEM *  UPp  ENZO,max  ENZ I ,max  ENZO,max  ENZ I ,max  DEM *  UPp  ENZO  ENZ I ,  (i.e., PI P* ) (4.33) where  0, if   ENZO   DEM *  UPp , if   ENZO,max , if  DEM *  UPp  0 0  DEM *  UPp  ENZO,max ENZO,max  DEM *  UPp and 123 (4.34)  0, if    ENZ I   DEM *  UPP  ENZO,max , if   ENZ I ,max , if   DEM *  UPP  ENZO,max ENZO,max  DEM *  UPP  ENZO,max  ENZ I ,max (4.35) ENZO,max  ENZ I ,max  DEM *  UPP The maximum rates, ENZO, max and ENZ I , max , are modeled as a first order kinetics as ENZ I , max  kenzi P OCC (4.36) ENZO, max  kenzo PH (4.37) with kenzi and kenzo dependent on microbial activity. 4.4.1.2 Phosphorus stored in plant litter ( PL ) The amount of phosphorus that is stored in leaf litter is expressed as dPL DECL LF MD    dt (C / P) LF  C / P M (C / P) L (4.38) where (C / P) LF is taken as a constant value, (C / P) M is the carbon to phosphorus ratio of the microbial biomass and (C / P) L is the carbon to phosphorus ratio of the decomposing leaf litter. 124 4.4.1.3 Phosphorus in microbial biomass pool ( PM ) In general, phosphorus is immobilized by microbes from the inorganic phosphorus pool when decomposing litter is nutrient poor, whereas phosphorus is mineralized by microbes when the decomposing litter is nutrient rich. The changing of phosphorus in microbial biomass pool is expressed by dividing the respective flux by its C/P ratio and taking into account net immobilization from the decomposing litter and humus pools, rh dPM  1  rr    C / P  dt C / P  HC L    DECH MD   IMM P (4.39)  DECL  1  rr    C / P  H  C / P M  where (C / P) H is the carbon phosphorus ratio in humus pool, (C / P) HC is a relatively constant C/P ratio for the recalcitrant organic matter which microbes are unable to decompose (i.e., the fraction going into the humus pool). In the model immobilization of phosphorus has two components: a term that directly depends on soil organic matter decomposition and another that depends on (and is proportional to) soil available inorganic phosphorus.   DEC  L  DECH IMM P  1     rr    kimm PI CM   C / P  C / P H   L        (4.40) where kimm is the rate at which microbes immobilize phosphorus from the PI pool and  is a term that accounts for the relative degree to which microbes are phosphorus limited. Similarly, the mineralization rate is modeled as, 125  DEC  L  DECH MIN   rr    C / P L C / P H    (4.41) In order to control the range of the C/P ratio for microbes, we use (C / P)THR  as the lower bound for (C / P) M to define the point at which the microbial biomass becomes limited by carbon and (C / P)THR  as the upper bound for (C / P) M , which defines the point where the microbial biomass becomes nutrient limited. Selected values for these ratios should be site specific. For microbial C/P ratio that are less than (C / P)THR  , all of the phosphorus associated with respired C will be mineralized. In contrast, for microbial C/P ratios that are greater than (C / P)THR  , yet less than (C / P)THR  , the residual phosphorus associated with respired carbon that is not required to satisfy microbial demand (1   ) will be mineralized (  ). Finally, for microbial C/P ratios that are greater than or equal to (C / P)THR  , all of the phosphorus associated with respired carbon will be immobilized. 1, if     (C / P)THR   (C / P) M      if   (C / P)THR   (C / P)THR    0, if  (C / P) M  (C / P)THR  (C / P)THR   (C / P) M  (C / P)THR  (C / P) M  (C / P)THR  (4.42) 4.4.1.4 Phosphorus in humus pool ( PH ) The change of phosphorus in humus pool can be expressed as 126 dPH DECL DECH  rh   ENZO dt (C / P) HC (C / P) H (4.43) where (C / P) HC is a constant value reflecting the inability of microbes to decompose more recalcitrant compounds. The value of  C / P  HC must also meet the following condition: (1  rr ) (C / P) HC  rh to satisfy the conservation of mass. Because the C/P of leaf litter (C / P) LF (C / P) LF is a constant value, the condition is based on the most conservative C/P ratio of litter that could occur. This condition was obtained because the amount of phosphorus going into the microbial pool during decomposition has to be greater than or equal to zero. This condition is expressed mathematically as (1  rr ) DECL DECL  rh  0. (C / P) LF (C / P) HC 4.4.2. Inorganic phosphorus Available soil phosphorus ( PI ) is typically found in low quantities because, (i) it is rapidly taken up and stored in vegetation (UP), (ii) can be immobilized by microbes (IMMP), (iii) can be fixed (FIX) to more occluded forms, which slowly becomes available over longer time periods, and (iv) is slowly weathered from mineral forms (WE). dPI  MIN  ATM  WE  kimm PI CM  FIX  UP  ERI  LE dt (4.44) where ERI accounts for losses of PI due to erosion, and soil particles that are transported during runoff. ERI is expressed as proportional to the amount of PI (i.e., ke PI ). ATM is 127 the flux of PI deposited atmospherically and WE is the input due to weathering from mineral forms, which is expressed proportionally to POCC (i.e., WE  kw P OCC ). Factors controlling FIX include soil PH and the presence of soil minerals such as Fe and Al that bind phosphorus compounds. Ignore the release of secondary minerals during soil formation and express FIX as proportional to PI : FIX  k f PI with k f dependent on the average soil PH. LE is expressed as a function of deep leaching losses as in Porporato et al., LE  L( S )kl PI nZ r s (4.45) where kl describes the solubility of PI in the rooting zone soil moisture, L( S ) is the rate of leakage of deep infiltration. Changes in the occluded phosphorus pool ( POCC ) are expressed as, dP OCC  FIX  WE  ENZ  ER I O dt (4.46) where ERO is proportional to POCC while using the same time constant as ke P OCC . 4.5. Model results To test the nutrient models we first applied them to the test cases described in the reference [P. D’Odorico 2003; Christiane 2012] using the same parameter sets to test the carbon, nitrogen and phosphorus evolution respectively. We then coupled the CLM-carbon and the 128 nitrogen and phosphorus models described above, and use the coupled model to simulate the observed nitrate and dissolved phosphorus concentrations in the stream in Kalamazoo watershed. 4.5.1. Nitrogen test case We tested the carbon-nitrogen model using the same stochastic approach to simulate nitrogen and carbon cycles in the broad-leafed savanna at Nylsvley (South Africa) as described in [P. D’Odorico 2003]. The infiltration was modeled as a Poisson process, and the soil moisture, transpiration and vertical percolation are modeled in the same way following this reference. The parameters representative of the climate, soil, vegetation and nutrient pools can be referred to [P. D’Odorico 2003] and are listed in table 4.1 also. Table 4.1 Parameters representative of the climate, soil, vegetation and nitrogen pool. Soil parameters Porosity n 0.4 Field capacity s fc 0.3 Hygroscopic point sh 0.02 -1 Saturated hydraulic conductivity K s (m d ) 1.1 Soil depth Z r (m) 0.80 Average storm frequency  (d-1) 0.23 Average storm depth  (mm) 11 Rainfall parameters Vegetation parameters Max. evapotranspiration Emax (mm d-1) 4.5 Litter fall (average rate) LF (gC m-2 d-1) 1.5 129 Table 4.1 cont’d Soil-vegetation parameters Point of incipient stress s* 0.17 Permanent wilting point sw 0.065 Litter fall (C / N ) L F 58 Microbial biomass (C / N ) M 11.5 Humus (C / N ) H 22 C/N ratios Average reported C, N amount Humus pool CH (gC m-3) 8500 Biomass pool CM (gC m-3) 12.5 – 125 Litter pool CL (gC m-3) 960 – 1400 Ammonium pool N  (gN m-3) 0 Nitrate pool N  (gC m-3)  1.25 Average soil moisture s 0.11 Rate of uptake (gN m-3 d-1) 0.02 – 0.03 Rate of mineralization (gN m-3 d-1) 0.021 Rate of litter decomposition (gN m-3 d-1) 1.2 Fraction of dissolved N+ a 0.05 Fraction of dissolved N- a 1 Plant demand of N+ DEM  0.2 Plant demand of N- 0.5 Nonlinear diffusion coefficient DEM  d Rescaled diffusion coefficient F 0.1 Litter decomposition rate kl 6.5 105 130 3 Table 4.1 cont’d N+ partition coefficient ki 1 N- partition coefficient ki 1 Humus decomposition rate kh 2.5 106 Nitrification rate coefficient kn 0.6 Isohumic coefficient rh  (C / N ) H  min  0.25,  (C / N ) L   Respiration fraction coefficient rr 0.6 Figure 4.5 shows examples of model generated time series for some of the state variables relevant to the carbon and nitrogen cycles, and Figure 4.6 shows the corresponding decomposition, net mineralization, uptake, and leaching rates. In both cases, the average values observed at the site (broken lines in both figures) are in good agreement with the simulations. The nitrogen and carbon dynamics are driven by fluctuations of soil water content, which in turn depend on the intermittent and unpredictable input of precipitation. As shown in Figure 4.5, depending on the inertia of the various pools and on the degree of dependence on soil moisture, the random fluctuations imposed by precipitation are filtered by the temporal dynamics of the state variables in a very different manner. Some variables  (e.g. NO3 ) preserve much of the high frequency variability imposed by the random forcing of precipitation, while some others ( CH , CL , and CM ) show smoother fluctuations. The  high frequency component of NO3 fluctuations (period of days to weeks) can be linked to the direct dependence of mineralization and nitrification on soil moisture, which transfers 131  the random fluctuations of the rainfall forcing to the budget of NO3 . On the other hand, the low frequency variability (period of seasons to years) resembles that of organic matter (in particular of the microbial biomass) and depends on the inertia imposed on the dynamics by the dimension of the soil carbon and nitrogen pools, which is very large compared to the fluxes. 132 Figure 4.5 Temporal dynamics of soil moisture (a), carbon litter (b), carbon humus (c), carbon biomass (d), ammonium (e), and nitrate (f) simulated for the case of the broadleafed savanna at Nylsvley. The broken lines represent the average values or the range of values observed at Nylsvley. The average values are reported whenever observations are available. 133 Figure 4.6 Simulated rates of litter decomposition (a), net mineralization (b), nitrate uptake (c) and nitrate leaching (d) in the broad-leafed savanna at Nylsvley. The broken lines represent the average values observed at Nylsvley. The nitrate leaching (d) is reported to be negligible at Nylsvley. 4.5.2. Phosphorus test case The modeling framework integrates several key processes to examine the controls on phosphorus availability over a wide range of time scales, which makes it applicable in any phosphorus-limited ecosystem. The carbon-phosphorus model was applied to a well 134 documented Cerrado ecosystem using data from a Long Term ecological Research site located in central Brazil as described in [Christane 2012]. The model could correctly represent all the key driving factors that impact the long and short term phosphorus, soil moisture and carbon dynamics. Again the precipitation was modeled as a Poisson process and the soil moisture, transpiration and vertical percolation are modeled in the same way following the reference. The litter fall data used a prescribed litter fall rate from Fig. 4 on [Christane 2012]. The parameters related to soil, climate, vegetation, carbon and phosphorus cycles and the source of these values that are representative of the simulation area can be referred to [Christane 2012] and also list in table 4.2. 135 Table 4.2 Parameters values for soil, climate, vegetation, carbon and phosphorus cycles that are representative of a Brazilian Cerrado ecosystem. Porosity n 0.625 Soil depth Z r (cm) 100 Saturated hydraulic conductivity K s (cm day-1) 452 Soil moisture when plants close stomata 0.27 Soil moisture at field capacity s* s fc Soil moisture at wilting point sW 0.17 Evaporation rate at wilting point EW (cm day-1) 0.01 Hygroscopic soil moisture sh 0.14 Average wet season storm depth  w (cm) 1.26 Average dry season storm depth  d (cm) 0.4 Average wet season storm frequency w (day-1) 0.46 Average dry season storm frequency d (day-1) 0.05 Role of AM in movement of phosphorus ions kp 10 Maximum transpiration rate Emax (cm day-1) 0.15 Maximum evaporation rate Tmax (cm day-1) 0.13 Litter decomposition rate kd (ha Mg C-1 day-1) 4 104 Leaching loss constant kl (day-1) 0.035 Maximum enzymatic dissolution rate from PH kenzo (day-1) 0.04 Fraction of litter that undergoes humification rh 0.18 Fraction of carbon lost to respiration rr 0.45 0.6 - Microbial C/P ratio where carbon is limting (C / P)THR  (kg C kg P 1 ) (C / P)THR  (kg C kg P-1) C/P ratio of leaf litter (C / P) LF (kg C kg P-1) 1300 Humus decomposition rate kh (ha Mg C-1 day-1) 3 106 Microbial C/P ratio where nutrients are limiting 136 65 56 Table 4.2 cont’d Microbial death rate km (day-1) 0.001 Optimal microbial death rate kopt (day-1) 8 106 Maximum enzymatic dissolve rate from POCC kenzi (day-1) 0.014 Rate at which microbes immobilize PI kimm (ha Mg C-1 day-1) 0.015 -1 Rate at which PI is fixed to POCC k f (day ) 0.0005 Amount of PI required before limiting P* (kg PI ha-1) ATM (kg PI ha-1 day-1) 0.054 PI deposited atmospherically 1.6 105 Results shown in Figures 4.7 and 4.8 indicate that the model captures well the mean behavior of the key state variables and fluxes. The solutions are close to those presented in the original paper (Christiane et al., 2012). The observed average values are shown in red broken lines in the figures. For instance, Resende et al. [2010] found that CM in the top -1 100cm is 2.6 Mg C ha , CL is 6.7 Mg C ha -1 -1 and PI is 1.2 kg P ha . The simulation results are fluctuated around the observed values and the mean values are comparable. 137 Figure 4.7 Simulated long term dynamics of net mineralization (a), phosphorus uptake (b), litter decomposition (c), and leaching of PI beneath the rooting zone (d). The dashed lines represent the mean observed fluxes for this system. 138 Figure 4.8 Long term temporal dynamics for simulated state variables in a Brazilian Cerrado ecosystem: (a) soil moisture; (b) carbon in the litter pool; (c) carbon in the microbial biomass; and, (d) available inorganic phosphorus in soil. 139 4.5.3. Nitrogen and Phosphorus model application to Kalamazoo watershed 4.5.3.1 Study area The Kalamazoo River drains into Lake Michigan and the watershed is one of the important watersheds in Michigan known for its natural resources. The Kalamazoo River watershed is located in the southwest portion of Michigan Lower peninsula. It drains approximately 2,020 square miles in ten counties. Geographically (Figure 4.9), the watershed is about 162 miles long and varies in width from 11 to 29 miles. Some of the larger tributaries of the Kalamazoo are Rice Creek, Wilder Creek, Wabascon Creek, Battle Creek River, Augusta Creek, Portage Creek, Gun River, Swan Creek, and Rabbit River. There are approximately 400,000 residents. The land use and land cover map of the watershed is shown in Figure 4.10. The hydrograph comparison at USGS site 04108660 located close to the outlet of the watershed is shown in Figure 4.11. The river and its watershed are notable for historical industrial contamination and the more recent oil spill in the river. Point sources of pollution from sewage and industrial activity are treated and their discharges are regulated. Although point sources of pollution are under control, excessive nutrients and other contaminants in the watershed are almost entirely due to non-point sources. Efforts are currently underway to develop nutrient (nitrogen and phosphorus) and contaminant transport models for the watershed. 140 Figure 4.9 Map of Kalamazoo watershed showing the location of USGS gauging, NCDC stations and nitrogen sampling sites et. al. 141 Figure 4.10 Land use and land cover map for Kalamazoo watershed. 142 Figure 4.11 Streamflow comparison for USGS 04108660 in Kalamazoo watershed 4.5.3.2 Nitrogen results The model of carbon-nitrogen cycles described in sections 4.2 and 4.3 was applied to the Kalamazoo watershed for simulating the soil nitrogen concentration in different pools. The time evolution of carbon and nitrogen variables were simulated in each computational grid. The litter carbon and nitrogen, and soil organic matter carbon and nitrogen concentrations from CLM were used and combined with carbon-nitrogen ratios in different pools to calculate the initial values for all other carbon and nitrogen state variables in the model. The nitrate concentration in soil simulated by the model in each time step was used as a source term which was washed off by surface runoff calculated in PAWS model and routed in the transport model. The nitrogen fertilizer application was simulated by using records 143 from United States Department of Agriculture (USDA) website. The fertilizer application timing should be pre-determined of the year for certain crops and there is no such records available. Thus we can only simply estimate the timing based on the precipitation records (assuming the fertilizer is only applied during the period without significant rainfall). The total amount of nitrogen fertilizer used is averaged over the computation grids based on the different type of crop land. After the plant removal, which is reported as percentage of fertilizer usage on USDA website, the residue was calculated as nitrogen input for the agricultural land. The nitrate concentration in the groundwater in the area is considered to be a significant source, but no detailed time series data were recorded. We use a value of xxx averaged from oberservation at different location and time. The nitrate observations in several streams in the Kalamazoo watershed were reported in [Baas, 2009]. Figure 4.12 shows the plot of simulated results comparing with observed data from early April to end September in 2006 at several sites shown in Figure 4.9. The grey bands in the each subplots of Figure 4.12 show the 50% of nitrogen fertilizer application. The model generally describes the trend of the change of nitrate concentration in the river as well as the peak concentrations. Note that since these sites are along the same river the change of simulated concentration profiles have the same pattern but different magnitudes. It has been found that groundwater is a major source contributing to the nitrate observation in the streams. Knowing the exact timing for fertilizer application and with the help of parameters estimation will further improve the model performance. 144 Figure 4.12 Nitrate concentration at four sites in Kalamazoo River. 145 4.5.3.3 Dissolved phosphorus results The model of carbon-phosphorus cycles described in section 4.2 and 4.4 was applied to Kalamazoo watershed for simulating the soil phosphorus concentration in different pools. The litter and soil organic matter carbon concentrations from CLM were used and combined with carbon-phosphorus ratios in different pools to calculate the initial values for all other carbon and phosphorus state variables in the model. The inorganic phosphorus concentration in soil simulated by the model in each time step was used as a source term which was washed off by surface runoff calculated in PAWS model and routed in the transport model. The phosphorus fertilizer application was simulated by using records from USDA website. The fertilizer application timing was simply estimated as described earlier and the total amount of phosphorus fertilizer used is averaged over the computation grids based on the different type of crop land. Similarly, the fertilizer residue and groundwater are also calculated as sources as described earlier but with different values. The dissolved phosphorus observations in several streams in Kalamazoo watershed were provided in [Baas, 2009]. Figure 4.13 shows the plot of simulated results comparing with observed data from early April to end September in 2006 at several sites shown in Figure 4.9. The bands in each subplots of Figure 4.13 are also due to the 50% fertilizer amount applied. Note that the much wider bands in Figure 4.13 than those in Figure 4.12 due to the different calculation of fertilizer residue of nitrogen and phosphorus fertilizer. The model generally describes the trend of the change of dissolved phosphorus concentration in the river as well as the peak concentrations. Almost all observation data can be covered by the grey bands in the Figure 4.13, which implies that the fertilizer application is the most significant source 146 for dissolved phosphorus concentraton in the streams. Note that the change of simulated concentration profiles have the same pattern but in different magnitude along the same river like Kalamazoo River. It would be noted that the dissolved phosphorus concentration shown in Figure 4.13 is a little portion in the total phosphorus observation. During resuspension process the occluded phosphorus or particulate phosphorus which occupy a large amount of phosphorus, the re-suspended phosphorus can be dissolved and change the phosphorus concentration. Also, the exact fertilizer application timing and parameters estimation will further improve the comparisons. 147 Figure 4.13 Dissolved phosphorus concentration at four sites in Kalamazoo River. 148 4.6. Conclusions The models presented in this chapter address the role of hydrological processes on the soil carbon and nitrogen, phosphorus cycles from the mathematical modeling point of view. The resulting schemes are systems of nonlinear ordinary differential equations for the carbon and nitrogen and phosphorus dynamics in the soil, driven by the water balance simulated by PAWS. Carbon and nitrogen in litter and soil organic matter simulated in CLM were adopted to the models' framework. The influence of rainfall, vegetation and soil characteristics is directly incorporated, through soil moisture dynamics, into the process of mineralization and immobilization, plant uptake and leaching. Two models were first tested against the observation in a broad-leafed savanna at Nylsvley (South Africa) and a well-documented Cerrado ecosystem located in central Brazil. The models were able to capture the mean behavior of key fluxes and state variables. Then the models were applied to Kalamazoo watershed and compared with the observed nitrate and dissolved phosphorus concentration in streams. The models generally described the trend of the change of the concentration profile and captured the peak concentrations at different observation sites. Although CLM carbon-nitrogen cycles were considered accurate and the carbon and nitrogen in litter and soil organic matter pools from CLM were used to estimate the initial condition for other state variables, the models are sensitive to initial conditions and many parameters used in the models are site-specific. Thus identifying the values of those sitespecific parameters as well as the initial conditions for different state variables for different 149 land use and land cover and soils types is still needed to improve the models performance. The wash off of fertilizer application on crop land can cause huge nitrogen and phosphorus concentration peak in the stream during storm event. Although the model also took the fertilizer application into consideration, the accurate fertilizer amount applied and the timing are vague in the current model estimation. Moreover, inorganic phosphorus only occupies a little portion in the total soil phosphorus amount, so does the dissolved phosphorus in total phosphorus observation. A large amount of phosphorus is in the occluded phosphorus form or particulate phosphorus, which is strongly dependent on soil particles and sediment fate. Thus a sediment transport module is needed to fulfill the phosphorus transport modeling framework, which is considered as one of the most important future work. 150 Chapter 5 Conclusion A process-based hydrologic model has been applied to four different watersheds in the State of Michigan. The water budget and storage were analyzed for two of the largest watersheds, Grand River watershed and Saginaw Bay watershed. The model was able to simulate different hydrologic components and states and the model results match the observation well. Based on the linear correlation analysis it has been found that the variations of different hydrologic components are sensitive to LULC and soil types. The process-based hydrologic model addresses the complexity of the flow in both subsurface and surface domains and explicitly connects the different hydrologic components. Thus it is suitable for large scale watershed transport simulation. A transport model was developed using an operator-splitting strategy combined with a Lagrangian particle transport modeling approach with reactions. The transport model was tested extensively using analytical solutions and data from plot-scale experiments. The hydrologic model together with the transport model were applied to Red Cedar River watershed to simulate the bacterial fate and transport. The statistic methods were used to estimate the release of the manure borne bacteria and fecal coliform discharged from WWTP. The nutrients transport model addresses the role of hydrological processes on the soil carbon and nitrogen, phosphorus cycles. Two models simulating the carbon-nitrogen and carbon-phosphorus cycles were tested against the observation using parameters values from literature. The hydrologic model combined with nutrients transport model were applied to Kalamazoo watershed and compared with the observed nitrate and dissolved phosphorus concentration in streams. 151 The carbon-nitrogen and carbon-phosphorus models are sensitive to initial conditions and many parameters used in the models are site-specific. In order to improve the model performance for the nutrients transport, further work is needed to identify the values of sitespecific parameters as well as the initial conditions of different state variables on different land use and land cover and soils types. To identify more unknown bacterial contamination sources accurately is also important for the bacteria transport model. The accurate timing and amount of fertilizer application can greatly improve the model results, however, these records are not accessible in this thesis research. Moreover, a large amount of phosphorus is in the occluded phosphorus form or particulate phosphorus, which is strongly dependent on soil particles and sediment fate and can be described by a sediment transport module. 152 REFERENCES 153 REFERENCES Adam, W., et al. (1995), The Great Lakes: an environmental atlas and resource book, 3rd ed, edited by K. Fuller, H. shear and J. Wittig, Government of Canada, Toronto, Ontario, Canada, United States Environmental Protection Agency, Great Lakes National Program Office, Chicago, Illinois, U.S.A., ISBNs: 0662234413, 9780662234418. Ahuja, L.R., 1990. Modeling soluble chemical-transfer to runoff with rainfall impact as a diffusion process. Soil Science Society of America Journal, 54(2): 312-321. Anderson, M.P., Cherry, J.A., 1979. Using models to simulate the movement of contaminants through groundwater flow systems. C R C Critical Reviews in Environmental Control, 9(2): 97-156, doi: 10.1080/10643387909381669. Arnold, J. G., and P. M. Allen (1996), Estimating hydrologic budgets for three Illinois watersheds, Journal of Hydrology, 176(1-4), 57-77, doi: 10.1016/0022-1694(95)02782-3. Berbery, E. H., Y. Luo, K. E. Mitchell, and A. K. Betts (2003), Eta model estimated land surface processes and the hydrologic cycle of the Mississippi basin, Journal of Geophysical Research: Atmospheres, 108(D22), 8852, doi: 10.1029/2002JD003192. Baptista, A.E.d.M., 1987. Solution of advection-dominated transport by EulerianLagrangian methods using the backwards method of characteristics, Massachusetts Institute of Technology, Massachusetts Institute of Technology. Bear, J., Bachmat, Y., 1990. Introduction to Modeling of Transport Phenomena in Porous Media. Theory and Applications of Transport in Porous Media (Book 4). Springer, 580 pp. Beven, K. (2006), Searching for the Holy Grail of scientific hydrology: Qt   S , R, t  A as closure, Hydrol. Earth Syst. Sci., 10(5), 609-618, wos: 000241763700001. Bicknell, B.R., Imhoff, J.C., Kittle, J.L.J., Donigian, A.S., Jr., Johanson, R.C., 1997. Hydrological simulation program-FORTRAN: User's manual for verison 11. U.S. Environmental Protection Agency National Exposure Research Laboratory, Athens, Ga. Black, P. E. (1997), Watershed functions, Journal of the American Water Resources Association, 33(1), 1-11, doi: 10.1111/j.1752-1688.1997.tb04077.x. Bradford, S.A., Schijven, J., 2002. Release of Cryptosporidium and Giardia from dairy calf manure: Impact of solution salinity. Environmental Science & Technology, 36(18): 39163923, doi: 10.1021/es025573l. 154 Bresler, E., 1973. SIMULTANEOUS TRANSPORT OF SOLUTES AND WATER UNDER TRANSIENT UNSATURATED FLOW CONDITIONS. Water Resources Research, 9(4): 975-986, doi: 10.1029/WR009i004p00975. Canale, R.P., Auer, M.T., Owens, E.M., Heidtke, T.M., Effler, S.W., 1993. Modeling fecal coliform bacteria—II. Model development and application. Water Research, 27(4): 703-714, doi: http://dx.doi.org/10.1016/0043-1354(93)90180-P. Chapra, S.C., 1997. Surface water-quality modeling. McGraw-Hill Companies, 844 pp. Cheng, L., Z. Xu, D. Wang, and X. Cai (2011), Assessing interannual variability of evapotranspiration at the catchment scale using satellite-based evapotranspiration data sets, Water Resources Research, 47(9), W09509, doi: 10.1029/2011WR010636. Cherkauer, K. A., and T. Sinha (2010), Hydrologic impacts of projected future climate change in the Lake Michigan region, J. Gt. Lakes Res., 36, 33-50, doi: 10.1016/j.jglr.2009.11.012. Chesnaux, R., Molson, J.W., Chapuis, R.P., 2005. An Analytical Solution for Ground Water Transit Time through Unconfined Aquifers. Ground Water, 43(4): 511-517. Chick, H., 1908. An investigation of the laws of disinfection. Journal of Hygiene, 8(1): 92158. Cho, K.H., Pachepsky, Y.A., Kim, J.H., Kim, J.W., Park, M.H., 2012. The modified SWAT model for predicting fecal coliforms in the Wachusett Reservoir Watershed, USA. Water Research, 46(15): 4750-4760, doi: 10.1016/j.watres.2012.05.057. Christensen, J. H., et al. (2007), Regional climate projections, in Climate Change 2007: The physical science basis. Contribution of working group I to the fourth assessment report of the intergovernmental panel on climate change, edited by S. Solomon, D. Qin, M. Manning, Z. Chen, M. Marquis, K. B. Averyt, M. Tignor and H. L. Miller, Cambridge University Press, Cambridge, United Kingdom, New York, NY, USA, ISBNs: 978-0-521-88009-1, 978-0-521-70596-7. Collatz, G. J., J. T. Ball, C. Grivet, and J. A. Berry (1991), Physiological and environmental regulation of stomatal conductance, photosynthesis and transpiration: a model that includes a laminar boundary layer, Agricultural and Forest Meteorology, 54(2–4), 107-136, doi: 10.1016/0168-1923(91)90002-8. Dai, Y. J., R. E. Dickinson, and Y. P. Wang (2004), A two-big-leaf model for canopy temperature, photosynthesis, and stomatal conductance, Journal of Climate, 17(12), 22812299, doi: 10.1175/1520-0442(2004)017<2281:ATMFCT>2.0.CO;2. 155 Deng, Z.Q., de Lima, J., de Lima, M.I.P., Singh, V.P., 2006. A fractional dispersion model for overland solute transport. Water Resources Research, 42(3), doi: 10.1029/2005wr004146. Dickinson, R. E. (1983), Land surface processes and climate surface albedos and energybalance, Adv Geophys, 25, 305-353, wos: A1983QS81300007. Dickinson, R. E., et al. (2002), Nitrogen controls on climate model evapotranspiration, Journal of Climate, 15(3), 278-295, doi: 10.1175/15200442(2002)015<0278:NCOCME>2.0.CO;2. Dong, W.C., Wang, Q.J., 2013. Modeling Soil Solute Release into Runoff and Transport with Runoff on a Loess Slope. Journal of Hydrologic Engineering, 18(5): 527-535, doi: 10.1061/(asce)he.1943-5584.0000622. Dorfman, M., Rosselot, K.S., 2011. Testing the Waters: A Guide to Water Quality at Vacation Beaches, Twenty-first Annual Report, Natural Resources Defense Council, Washington D.C. Dorner, S.M., 2004. Waterborne pathogens: sources, fate, and transport in a watershed used for drinking water supply. Dissertation Thesis, University of Waterloo, 395 Wellington Street, Ottawa ON, K1A 0N4, Canada, 366 pp. Dougherty, E.R., 1990. Probability and Statistics for the Engineering, Computing and Physical Sciences. Prentice Hall, New Jersey, 784 pp. Enviro-weather, 2010. Enviro-weather Automated Weather Station Network, (formerly known as MAWN) East Lansing, Michigan, 48823, Available at: http://www.agweather.geo.msu.edu/mawn/ Farquhar, G. D., S. V. Caemmerer, and J. A. Berry (1980), A biochemical-model of photosynthetic CO2 assimilation in leaves of C-3 species, Planta, 149(1), 78-90, doi: 10.1007/BF00386231. Fraser, R.H., Barten, P.K., Pinney, D.A.K., 1998. Predicting stream pathogen loading from livestock using a geographical information system-based delivery model. Journal of Environmental Quality, 27(4): 935-945. Fu, B. P. (1981), On the calculation of the evaporation from land surface, Chinese Journal of Atmospheric Sciences (in Chinese), 5(1), 23-31, doi:10.3878/j.issn.10069895.1981.01.03. Gao, H., Q. Tang, C. R. Ferguson, E. F. Wood, and D. P. Lettenmaier (2010), Estimating the water budget of major US river basins via remote sensing, Int J Remote Sens, 31(14), 3955-3978, doi: 10.1080/01431161.2010.483488. 156 Gaumnitz, L., T. Asplund, and M. R. Matthews (2004), A growing thirst for groundwater, Wisconsin Natural Resources magazine. Gameson, A.L.H., Gould, D.J., 1974. Effects of Solar Radiation on the Mortality of Some Terrestrial Bacteria in Sea Water, Proc. Int. Symp. on Discharge of Sewage from Sea Outfalls. Pergamon Press, London, Great Britain. Gelhar, L.W., Wilson, J.L., 1974. - Ground-Water Quality Modelinga. - 12(- 6): - 408 Gleeson, T., Y. Wada, M. F. P. Bierkens, and L. P. H. van Beek (2012), Water balance of global aquifers revealed by groundwater footprint, Nature, 488(7410), 197-200, doi: 10.1038/nature11295. Guber, A.K., Yakirevich, A.M., Sadeghi, A.M., Pachepsky, Y.A., Shelton, D.R., 2009. Uncertainty Evaluation of Coliform Bacteria Removal from Vegetated Filter Strip under Overland Flow Condition. Journal of Environmental Quality, 38(4): 1636-1644, doi: 10.2134/jeq2008.0328. Gunduz, O., 2004. Coupled flow and contaminant transport modeling in large watersheds. Dissertation Thesis, Georgia Institute of Technology, Georgia Institute of Technology. Gupta, R. S. (2008), Hydrology and Hydraulic Systems, 3rd ed., 896 pp., Waveland Press, Incorporated, Long Grove, IL., ISBN: 1577664558. Istanbulluoglu, E., T. Wang, O. M. Wright, and J. D. Lenters (2012), Interpretation of hydrologic trends from a water balance perspective: The role of groundwater storage in the Budyko hypothesis, Water Resources Research, 48(3), doi: 10.1029/2010WR010100. Jamieson, R., Gordon, R., Joy, D., Lee, H., 2004. Assessing microbial pollution of rural surface waters: A review of current watershed scale modeling approaches. Agricultural Water Management, 70(1): 1-17, doi: http://dx.doi.org/10.1016/j.agwat.2004.05.006. Jia, Y. W., G. H. Ni, Y. Kawahara, and T. Suetsugi (2001), Development of WEP model and its application to an urban watershed, Hydrological Processes, 15(11), 2175-2194, doi: 10.1002/hyp.275. Kalma, J., T. McVicar, and M. McCabe (2008), Estimating land surface Evaporation: A review of methods using remotely sensed surface temperature data, Surv Geophys, 29(4-5), 421-469, doi: 10.1007/s10712-008-9037-z. Kinzelbach, W., 1986. Groundwater Modeling, An Introduction with Sample Programs in BASIC. Elsevier, Amsterdam, Netherlands, 333 pp. 157 Kling, G. W., et al. (2003), Confronting climate change in the Great Lakes region: Impacts on our communities and ecosystems, 104 pp, Union of Concerned Scientists, Cambridge, Massachusetts, Ecological Society of America, Washington, D.C. Kondo, J., and S. Ishida (1997), Sensible heat flux from the earth’s surface under natural convective conditions, Journal of the Atmospheric Sciences, 54(4), 498-509, doi: 10.1175/1520-0469(1997)054<0498:SHFFTE>2.0.CO;2. LaBolle, E.M., Fogg, G.E., Tompson, A.F.B., 1996. Random-walk simulation of transport in heterogeneous porous media: Local mass-conservation problem and implementation methods. Water Resources Research, 32(3): 583-593, doi: 10.1029/95wr03528. Li, K. Y., J. B. Boisvert, and R. D. Jong (1999), An exponential root-water-uptake model, Canadian Journal of Soil Science, 79(2), 333-343, WOS:000081215100012. Lu, D., and Q. Weng (2006), Use of impervious surface in urban land-use classification, Remote Sensing of Environment, 102(1–2), 146-160, doi:10.1016/j.rse.2006.02.010. Luo, Y., E. H. Berbery, and K. E. Mitchell (2005), The operational eta model precipitation and surface hydrologic cycle of the Columbia and Colorado basins, Journal of Hydrometeorology, 6(4), 341-370, doi: 10.1175/JHM435.1. Mahalanobis, P. C. (1936), On the generalised distance in statistics, paper presented at Proceedings National Institute of Science, India, 2(1), 49-55. Mancini, J.L., 1978. NUMERICAL ESTIMATES OF COLIFORM MORTALITY-RATES UNDER VARIOUS CONDITIONS. Journal Water Pollution Control Federation, 50(11): 2477-2484. MDNR, 2010. 2001 IFMAP/GAP Lower Peninsula Land Cover, available at http://www.mcgi.state.mi.us/mgdl/?rel=thext&action=thmname&cid=5&cat=Land+Cover+ 2001, retrieved 2009-Nov-28. Medema, G.J., Ketelaars, H.A.M., Hoogenboezem, W., 2001. Cryptosporidium and giardia : occurrence in sewage, manure, and surface water / editors, G.J. Medema, H.A.M. Ketelaars, W. Hoogenboezem ; authors: W. Hoogenboezem ... [et al.]. RIWA, Association of River Waterworks. Mitchell, K. E., et al. (2004), The multi-institution North American Land Data Assimilation System (NLDAS): Utilizing multiple GCIP products and partners in a continental distributed hydrological modeling system, Journal of Geophysical Research: Atmospheres, 109(D7), D07S90, doi: 10.1029/2003JD003823. 158 Moore, J.A., Smyth, J.D., Baker, E.S., Miner, J.R., Moffitt, D.C., 1989. MODELING BACTERIA MOVEMENT IN LIVESTOCK MANURE SYSTEMS. Transactions of the Asae, 32(3): 1049-1053. Mu, Q., M. Zhao, and S. W. Running (2011), Improvements to a MODIS global terrestrial evapotranspiration algorithm, Remote Sensing of Environment, 115(8), 1781-1800, doi: 10.1016/j.rse.2011.02.019. Mu, Q., F. A. Heinsch, M. Zhao, and S. W. Running (2007), Development of a global evapotranspiration algorithm based on MODIS and global meteorology data, Remote Sensing of Environment, 111, 519-536, doi: 10.1016/j.rse.2007.04.015. NCDC, 2010. National Climatic Data Center, http://www.ncdc.noaa.gov/oa/climate/climatedata.html#daily available at Neuman, S.P., 1984. ADAPTIVE EULERIAN LAGRANGIAN FINITE-ELEMENT METHOD FOR ADVECTION DISPERSION. International Journal for Numerical Methods in Engineering, 20(2): 321-337, doi: 10.1002/nme.1620200211. Niu, G. Y., and Z. L. Yang (2006), Effects of frozen soil on snowmelt runoff and soil water storage at a continental scale, Journal of Hydrometeorology, 7(5), 937-952, doi: 10.1175/JHM538.1. Nouho Koffi Quattara, Anouk de Brauwere, Gilles Billen, Pierre Servais, 2012. Modelling faecal contamination in the Scheldt drainage network. Journal of Marine Systems, doi:10.1016/j.jmarsys.2012.05.004. O'Loughlin, E.M., Bowmer, K.H., 1975. Dilution and decay of aquatic herbicides in flowing channels. Journal of Hydrology, 26(3–4): 217-235, doi: http://dx.doi.org/10.1016/0022-1694(75)90004-9. Oleson, K. W., et al. (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, CO., ISSNs: 2153-2397, 2153-2400. Olson, S.R., Kemper, W.D., 1961. Movement of nutrients to plant roots. Adv.Agron., 20: 91-151. Pan, M., and E. F. Wood (2006), Data assimilation for estimating the terrestrial water budget using a constrained ensemble Kalman filter, Journal of Hydrometeorology, 7(3), 534-547, doi: 10.1175/JHM495.1. Patrick P. Pasmussen and Andrew C. Ziegler, 2003. Comparions and Continuous Estimates of Fecal Coliform and Escherichia Coli Bacteria in Selected Kansas Streams, May 1999 Through April 2002, USGS, Water-Resources Investigations Report, 03-4056. 159 Pinder, G.F., Gray, W.G., 1977. finite elements simulaion in Surface and subsurface hydrology. Academic Press, San Diego, 295 pp. Prickett, T.A., Naymik, T.G., Lonnquist, C.G., 1981. A "Random-Walk" Solute Transport Model for Selected Groundwater Quality Evaluations. In: Resources, I.D.o.E.a.N. (Ed.), State of Illinois. Reed M. Maxwell, et al. (2013), Surface-subsurface model intercomparison: A first set of benchmark results to diagnose integrated hydrology and feedbacks, Water Resources Research. Roads, J., et al. (2003), GCIP water and energy budget synthesis (WEBS), Journal of Geophysical Research: Atmospheres, 108(D16), 8609, doi: 10.1029/2002JD002583. Roads, J. O., S. C. Chen, M. Kanamitsu, and H. Juang (1999), Surface water characteristics in NCEP global spectral model and reanalysis, Journal of Geophysical Research: Atmospheres, 104(D16), 19307-19327, doi: 10.1029/98JD01166. Rodell, M., J. Chen, H. Kato, J. Famiglietti, J. Nigro, and C. Wilson (2007), Estimating groundwater storage changes in the Mississippi River basin (USA) using GRACE, Hydrogeol J, 15(1), 159-166, doi: 10.1007/s10040-006-0103-7. Rodriguez-Iturbe, I., and J. M. Mejí (1974), On the transformation of point rainfall to areal a rainfall, Water Resources Research, 10(4), 729-735, doi: 10.1029/WR010i004p00729. Rolle, K., Gitau, M.W., Chen, G., Chauhan, A., 2012. Assessing fecal coliform fate and transport in a coastal watershed using HSPF. Water Science and Technology, 66(5): 10961102, doi: 10.2166/wst.2012.282. Sadeghi, A.M., Arnold, J.G., 2002. A SWAT/Microbial Sub-Model for Predicting Pathogen Loadings in Surface and Groundwater at Watershed and Basin Scales, Total Maximum Daily Load (TMDL) Environmental Regulatoins: Proceedings of the March 11-13, 2002 Conference, Fort Worth, Texas, USA, pp. 56-63. Sawicz, K., T. Wagener, M. Sivapalan, P. A. Troch, and G. Carrillo (2011), Catchment classification: empirical analysis of hydrologic similarity based on catchment function in the eastern USA, Hydrology and Earth System Sciences, 15(9), 2895-2911, doi: 10.5194/hess-15-2895-2011. Sayama, T., J. J. McDonnell, A. Dhakal, and K. Sullivan (2011), How much water can a watershed store?, Hydrological Processes, 25(25), 3899-3908, doi: 10.1002/hyp.8288. Seilheimer, T. S., P. L. Zimmerman, K. M. Stueve, and C. H. Perry (2013), Landscapescale modeling of water quality in Lake Superior and Lake Michigan watersheds: How 160 useful are forest-based indicators?, J. Gt. Lakes Res., 39(2), 211-223, ISSN: 0380-1330, 10.1016/j.jglr.2013.03.012. Sellers, P. J. (1985), Canopy reflectance, photosynthesis and transpiration, Int J Remote Sens, 6(8), 1335-1372, WOS:A1985AQY5100003. Shen, C., J. Niu, and M. S. Phanikumar (2013), Evaluating controls on coupled hydrologic and vegetation dynamics in a humid continental climate watershed using a subsurface-land surface processes model, Water Resources Research, 49, 1-21, doi:10.1002/wrcr.20189. Shen, C. P., and M. S. Phanikumar (2010), A process-based, distributed hydrologic model based on a large-scale method for surface-subsurface coupling, Advances in Water Resources, 33(12), 1524-1541, doi: 10.1016/j.advwatres.2010.09.002. Simard, A., 2007. Predicting groundwater flow and transport using Michigan's statewide wellogic database, Michigan State University, East Lansing, 109 pp. Singh, V.P., 1996. Kinematic Wave Modeling in Water Resources, Surface-Water Hydrology. Wiley-Interscience, 1424 pp. Sivapalan, M. (2006), Pattern, process and function: Elements of a unified theory of hydrology at the catchment scale, in Encyclopedia of Hydrological Sciences, edited, John Wiley & Sons, Ltd., doi: 10.1002/0470848944.hsa012. Soil Survey Staff, Survey Geographic (SSURGO) Database for Michigan. Natural Resources Conservation Service, United States Department of Agriculture. Available online at http://soildatamart.nrcs.usda.gov Accessed 06/01/2010. Su, F., and D. P. Lettenmaier (2009), Estimation of the surface water budget of the La Plata Basin, Journal of Hydrometeorology, 10(4), 981-998, doi: 10.1175/2009JHM1100.1. Tague, C., G. Grant, M. Farrell, J. Choate, and A. Jefferson (2008), Deep groundwater mediates streamflow response to climate warming in the Oregon Cascades, Climatic Change, 86(1-2), 189-210, doi: 10.1007/s10584-007-9294-8. Tavener, B., and M. Iqbal (2003), The development of a hydrologic budget to determine the nitrogen and phosphorus loads of the Cedar River watershed in Iowa, Environ Geol, 43(4), 400-407, doi: 10.1007/s00254-002-0657-1. Thomann, R.V., Mueller, J.A., 1997. Principles of Surface Water Quality Modeling and Control. Prentice Hall, New York, 656 pp. Thornton, P. E., and N. A. Rosenbloom (2005), Ecosystem model spin-up: Estimating steady state conditions in a coupled terrestrial carbon and nitrogen cycle model, Ecological Modelling, 189(1-2), 25-48, doi: 10.1016/j.ecolmodel.2005.04.008. 161 Thornton, P. E., and N. E. Zimmermann (2007), An improved canopy integration scheme for a land surface model with prognostic canopy structure, Journal of Climate, 20(15), 3902-3923, doi: 10.1175/JCLI4222.1. Tian, Y.Q., Gong, P., Radke, J.D., Scarborough, J., 2002. Spatial and temporal modeling of microbial contaminants on grazing farmlands. Journal of Environmental Quality, 31(3): 860-869. Tompson, A.F.B., Gelhar, L.W., 1990. Numerical simulation of solute transport in threedimensional, randomly heterogeneous porous media. Water Resources Research, 26(10): 2541-2562, doi: 10.1029/WR026i010p02541. Toon, O. B., C. P. Mckay, T. P. Ackerman, and K. Santhanam (1989), Rapid Calculation of Radiative Heating Rates and Photodissociation Rates in Inhomogeneous MultipleScattering Atmospheres, J Geophys Res-Atmos, 94(D13), 16287-16301, doi: 10.1029/JD094iD13p16287. Trapp, R. J., N. S. Diffenbaugh, H. E. Brooks, M. E. Baldwin, E. D. Robinson, and J. S. Pal (2007), Changes in severe thunderstorm environment frequency during the 21st century caused by anthropogenically enhanced global radiative forcing, Proceedings of the National Academy of Sciences, 104(50), 19719-19723, doi: 10.1073/pnas.0705494104. USEPA, 2001. Protocol for Developing Pathogen TMDLs. EPA 841R00002. In: Office of Water (4503F), U. (Ed.), Washington D.C., pp. 132. Vadas, P.A., Kleinman, P.J.A., Sharpley, A.N., 2004. A simple method to predict dissolved phosphorus in runoff from surface-applied manures. Journal of Environmental Quality, 33(2): 749-756. van Genuchten, M.T., Alves, W.J., 1982. Analytical solutions of the one-dimensional convective-dispersion solute transport equation. In: Agriculture, U.S.D.o. (Ed.), pp. 151. Vogel, J.C., 1968. INVESTIGATION OF GROUNDWATER FLOW WITH RADIOCARBON. Journal Name: pp 355-69 of Isotopes in Hydrology. Vienna, International Atomic Energy Agency, 1967.; Other Information: From Conference on Isotopes in Hydrology, Vienna. See STI/ PUB--141; CONF- 661133. Orig. Receipt Date: 31-DEC-68: Medium: X. Wagener, T., M. Sivapalan, P. Troch, and R. Woods (2007), Catchment Classification and Hydrologic Similarity, Geography Compass, 1(4), 901-931, doi: 10.5194/hess-15-28952011. 162 Walker, S.E., Mostaghimi, S., Dillaha, T.A., Woeste, F.E., 1990. MODELING ANIMAL WASTE MANAGEMENT-PRACTICES - IMPACTS ON BACTERIA LEVELS IN RUNOFF FROM AGRICULTURAL LANDS. Transactions of the Asae, 33(3): 807-817. Wang, D. (2011), Evaluating interannual water storage changes at watersheds in Illinois based on long-term soil moisture and groundwater level data, Water Resources Research, doi:10.1029/2011WR010759. Wang, T., E. Istanbulluoglu, J. Lenters, and D. Scott (2009), On the role of groundwater and soil texture in the regional water balance: An investigation of the Nebraska Sand Hills, USA, Water Resources Research, 45(10), W10413, doi: 10.1029/2009WR007733. Warrick, A.W., 2003. Soil Water Dynamics. Oxford University Press, USA, 198 Madison Avenue, New York, New York, 10016. Wilkinson, J., Jenkins, A., Wyer, M., Kay, D., 1995. MODELING FECAL-COLIFORM DYNAMICS IN STREAMS AND RIVERS. Water Research, 29(3): 847-855, doi: 10.1016/0043-1354(94)00211-o. Wuebbles, D., and K. Hayhoe (2004), Climate change projections for the United States Midwest, Mitigation and Adaptation Strategies for Global Change, 9(4), 335-363, doi: 10.1023/B:MITI.0000038843.73424.de. Zheng, C., Bennett, G.D., 2002. Applied Contaminant Transport Modeling. Wiley Interscience, 621 pp. 163