MODELING THE MOVEMENT OF WATER，BACTERIA AND NUTRIENTS
ACROSS HETEROGENEOUS LANDSCAPES IN THE GREAT LAKES REGION
USING A PROCESSBASED 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 PROCESSBASED HYDROLOGIC MODEL
By
Jie Niu
The development and application of processbased 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 largescale hydrology remains a
topic that is relatively unexplored. Understanding controls on largescale 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
remotelysensed datasets for evapotranspiration (ET) and land water thickness equivalent
(GRACE). We then report a budget analysis of major hydrologic fluxes and compute
annualaverage 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 interannual 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 operatorsplitting 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 watershedscale 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 processbased 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. ShuGuang 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. ShuGuang
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 ProcessBased,
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. Twodimensional 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. Onedimensional channel transport ····································72
···································
3.3.2. Twodimensional overland flow transport······························ 75
······························
3.3.3. Onedimensional solute transport in a soil column ·······················80
······················
vi
3.3.4. Groundwater transit time simulation ···································87
··································
3.3.5. Plot scale simulations of manureborne 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 Fiveyear 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 gravitybased 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 lefthand
boundary is impermeable, and the flow discharges through the righthand fixedhead
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 104 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 103 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 105 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 broadleafed 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 processbased 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 intermodel variability in water
budget estimates to projected future changes due to climate change.
The
hydrologic
water
budget
represents
the
balance
between
precipitation,
evapotranspiration, surface runoff, 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 realtime 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 satellitebased 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]. Spatiallydistributed hydrologic fluxes from model reanalysis have become more popular (e.g. NLDAS [Mitchell et al., 2004]). Welltested,
processbased 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, processbased 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
storagedischarge relations due to complex subsurface heterogeneity [Beven, 2006; Sayama
et al., 2011]. The application of physicsbased watershed models to large river basins is a
topic of increasing interest. Combining gravitybased satellite measurements [Rodell et al.,
2007], MODIS products for evapotranspiration and leaf area index and groundbased
3
measurements for other key hydrologic variables (e.g., soil moisture, streamflows,
groundwater heads, temperature) opens up the possibility to test physicsbased 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 physicsbased watershed model and independent
(satellite and groundbased) 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 welltested
processbased models? What is the longterm relation between different components, as
expressed, for example, by the Budyko curve? How do storage and groundwater flow
influence the longterm 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 nonpoint 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 swimmingassociated 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 manureborne
bacteria releasing process and the transport of bacteria through runoff, however, the dieoff
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 longterm
observations of flow and bacterial concentrations. Models that incorporate both landscape
and instream 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 dieoff, 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 instream 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
longterm, daily simulation of nutrient load in large catchments [Arnold and Fohrer, 2005;
Viney et al., 2000], whereas others are eventbased and are clearly unsuitable for longterm
continual predictions. Although the longterm 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 processbased, watershedscale 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, gridbased models is expected to facilitate seamless
integration of watershed models with similar (i.e., gridbased) lake and ocean models to
make realtime (or near realtime) 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 carbonnitrogen cycles and carbonphosphorus 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 remotelysensed data and a processbased hydrologic model that includes detailed representations of subsurface and land
surface processes. Model performance is evaluated using groundbased observations
(streamflows, groundwater heads, soil moisture, soil temperature) as well as satellitebased
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 annualaverage 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 physicsbased watershed
models that explicitly describe surface (e.g., St. Venant equations for channel flow) and
subsurface (e.g., Richards and/or Darcy equationbased) 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
processbased models to large watersheds. The model has been tested extensively using
available analytical solutions, laboratory data, streamflow, sitemeasured 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 (20032007). 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 modelbased
ET estimates and comparisons are shown as both timeseries of spatial averages and 2D
spatial fields. PAWS model descriptions of storage changes are then compared with
gravitybased 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 longterm
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 realworld, 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 subdomains  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 subgrid heterogeneity
within each land cell. Several generic plant functional types (PFT) [Thornton and
10
Zimmermann, 2007] represent model classes and are used to reclassify 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
Quasi3D
Recharge / Discharge (Vadose zone /
Groundwater interaction)
Noniterative 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, processbased 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 freezethaw 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 5year 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 overland 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 processbased 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 EnviroWeather
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 nonarctic 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 subwatershed 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 massbalance 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 5year 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 5year 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 underestimation 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 EnviroWeather stations
within the two watersheds. We compared soil moisture and soil temperature at 10 cm depth
for EnviroWeather 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 twoyear
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 twobigleaf model
[Dai et al., 2004], while the MODIS product is based on the PenmanMonteith 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 fiveyearaverage 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 Fiveyear 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 pairwise 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 5year 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 nonarctic 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 fiveyear 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 gravitybased 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 subwatersheds in the GR watershed and the whole
GR watershed for a 13year 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
PenmanMonteith equation. In this figure the horizontal straight line indicates the arid or
waterlimited conditions, while the 1 to 1 line indicates the humid or energy limited limit.
The xaxis, annual potential evaporation to annual precipitation ratio, can also be
considered as a dryness index. The annual ET to P ratio for all the subwatershed 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 interannual 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 interannual 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 satellitebased estimates in [Cheng et al.,
2
2011], indicating more complexity in the interannual 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 satellitebased estimates in [Cheng et al.,
2011] used modified PenmanMonteith with biomespecific 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 interannual difference in
nitrogen cycling could explain the difference in interannual 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 interannual
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 subbasinaverage 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 longterm 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 processbased 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 processbased 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
processbased 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 longduration of large watershed
tends to be zero, the interannual 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 interannual 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
ProcessBased, Distributed Hydrologic Model with Application to
Bacterial Fate and Transport Modeling
There is a surge of interest in the application of mechanistic or processbased 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 watershedscale 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 realtime 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 longterm, largescale simulations and makes it possible to simulate transport in
large watersheds. An operatorsplitting 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 plotscale 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 crosssectionally 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 onedimensional 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 crosssectional 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 onedimensional 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 vt
(3.13)
Then the concentration value at this point is calculated using a linear interpolation function:
x vt 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 ) N1 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 ) N1 1QNC1 (Cr ) N1 2 QNC 2
,C
,
,C
,
j 1
j 1
QNC1 QNC 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 ) N1 1 (Cr )1,1
,C
C3
j
j
(Cr ) N1 2 (Cr )1,1
,C
C3
(3.17)
The treatment of junctions is slightly different in the dispersion/decay/source/sink operator.
A timeweighted mass balance can be written as:
n
j
j
M in1 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 advectiondiffusion
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 nonadsorbed 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
rainfallrunoff.
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
2z
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 / 2z , 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
2z
n 1 n 1 n 1
i 1
i 1
D0 ab(eb ) i
2z
(3.34)
3.2.5. Twodimensional groundwater transport equation
A general partial differential equation that governs twodimensional 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 advectiondispersionreaction
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 advectiondispersion
equation can be classified as Eulerian, Lagrangian, or mixed EulerianLagrangian (Baptista,
61
1987; Neuman, 1984). Eulerian methods offer the advantage and convenience of a fixed
grid, are generally mass conservative, and handle dispersiondominated problems both
accurately and efficiently. However, for advectiondominated 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 fieldscale problem, especially for large watershed
simulations may become prohibitive.
Lagrangian methods, on the other hand, provide an accurate and efficient solution to
advectiondominated 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 massbalance errors
(LaBolle et al., 1996). When Eulerian methods are used, the timestep 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 firstorder
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 n1 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 RungeKutta 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 RungeKutta 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 EulerianLagrangian 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 EulerianLagrangian 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 EulerianLagrangian methods may be
65
loosely grouped as the forwardtracking 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 randomwalk
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
EulerianLagrangian solution techniques it can also lead to large massbalance
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 logtransformed data.
Several models have been proposed and tested to simulate the release of manureborne
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 k1Q
(3.50)
and HSPF (Hydrological Simulation Program – FORTRAN) (Bicknell et al., 1997):
M R M S 1 exp(k2Q)
(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, cm1.
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 twoparametric 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):
11/
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 surfaceapplied 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
d1. 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 hr1). 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 manureborne 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 dieoff rate constant, while we use the total loss rate, kb , above instead. The fate and
transport of bacteria at the watershedscale are modeled using the above formulations for
bacterial dieoff 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. Onedimensional channel transport
The onedimensional channel contaminant transport model is tested against the available
analytical solutions in a single channel framework. The onedimensional 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 advectivedispersive 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. Twodimensional overland flow transport
For this test case, we consider twodimensional 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. Onedimensional solute transport in a soil column
For nonadsorbed 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
z0
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
z0
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 oneto 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 lefthand boundary is the imaginary impermeable ground water divide and
flow discharges through the righthand fixedhead 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 DupuitForchheimer
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 lefthand
boundary is impermeable, and the flow discharges through the righthand fixedhead
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 firsttype 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 108 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 righthand 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 104 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 103 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 105 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 manureborne fecal coliforms
Overland flow experiments with the manureborne fecal coliforms (FC) were performed in
October 2003 on a twosided 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
30cm 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 reclassified
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)
(Enviroweather, 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, uncounted 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, SKZ 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
carbonnitrogen cycles and carbonphosphorus 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 nutrientcycles 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
colimited 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 hydroclimatic 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 3070% 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, airborne 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 CLMcarbon 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 carbonnitrogen model using the same stochastic approach to simulate
nitrogen and carbon cycles in the broadleafed 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
(d1)
0.23
Average storm depth
(mm)
11
Rainfall parameters
Vegetation parameters
Max. evapotranspiration
Emax (mm d1) 4.5
Litter fall (average rate)
LF (gC m2 d1) 1.5
129
Table 4.1 cont’d
Soilvegetation 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 m3) 8500
Biomass pool
CM (gC m3) 12.5 – 125
Litter pool
CL (gC m3)
960 – 1400
Ammonium pool
N (gN m3)
0
Nitrate pool
N (gC m3)
1.25
Average soil moisture
s
0.11
Rate of uptake
(gN m3 d1)
0.02 – 0.03
Rate of mineralization
(gN m3 d1)
0.021
Rate of litter decomposition
(gN m3 d1)
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 105
130
3
Table 4.1 cont’d
N+ partition coefficient
ki
1
N partition coefficient
ki
1
Humus decomposition rate
kh
2.5 106
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 broadleafed 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
phosphoruslimited ecosystem. The carbonphosphorus 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 day1)
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 day1)
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 (day1)
0.46
Average dry season storm frequency
d (day1)
0.05
Role of AM in movement of phosphorus ions
kp
10
Maximum transpiration rate
Emax (cm day1)
0.15
Maximum evaporation rate
Tmax (cm day1)
0.13
Litter decomposition rate
kd (ha Mg C1 day1)
4 104
Leaching loss constant
kl (day1)
0.035
Maximum enzymatic dissolution rate from PH
kenzo (day1)
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 P1)
C/P ratio of leaf litter
(C / P) LF (kg C kg P1)
1300
Humus decomposition rate
kh (ha Mg C1 day1)
3 106
Microbial C/P ratio where nutrients are limiting
136
65
56
Table 4.2 cont’d
Microbial death rate
km (day1)
0.001
Optimal microbial death rate
kopt (day1)
8 106
Maximum enzymatic dissolve rate from POCC
kenzi (day1)
0.014
Rate at which microbes immobilize PI
kimm (ha Mg C1 day1)
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 ha1)
ATM (kg PI ha1 day1)
0.054
PI deposited atmospherically
1.6 105
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 nonpoint 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 carbonnitrogen 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 carbonnitrogen 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 predetermined 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 carbonphosphorus 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 carbonphosphorus 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 resuspended 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 broadleafed savanna at Nylsvley
(South Africa) and a welldocumented 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 carbonnitrogen 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 sitespecific. 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 processbased 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
processbased 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 operatorsplitting strategy combined with a Lagrangian particle transport modeling
approach with reactions. The transport model was tested extensively using analytical
solutions and data from plotscale 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 carbonnitrogen and carbonphosphorus 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 carbonnitrogen and carbonphosphorus models are sensitive to initial conditions and
many parameters used in the models are sitespecific. 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 chemicaltransfer to runoff with rainfall impact as a
diffusion process. Soil Science Society of America Journal, 54(2): 312321.
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): 97156, doi: 10.1080/10643387909381669.
Arnold, J. G., and P. M. Allen (1996), Estimating hydrologic budgets for three Illinois
watersheds, Journal of Hydrology, 176(14), 5777, doi: 10.1016/00221694(95)027823.
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 advectiondominated 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), 609618, wos: 000241763700001.
Bicknell, B.R., Imhoff, J.C., Kittle, J.L.J., Donigian, A.S., Jr., Johanson, R.C., 1997.
Hydrological simulation programFORTRAN: 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), 111, doi: 10.1111/j.17521688.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): 975986, 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): 703714,
doi: http://dx.doi.org/10.1016/00431354(93)90180P.
Chapra, S.C., 1997. Surface waterquality modeling. McGrawHill Companies, 844 pp.
Cheng, L., Z. Xu, D. Wang, and X. Cai (2011), Assessing interannual variability of
evapotranspiration at the catchment scale using satellitebased 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, 3350, 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): 511517.
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): 47504760, 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: 9780521880091,
9780521705967.
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), 107136, doi:
10.1016/01681923(91)900028.
Dai, Y. J., R. E. Dickinson, and Y. P. Wang (2004), A twobigleaf model for canopy
temperature, photosynthesis, and stomatal conductance, Journal of Climate, 17(12), 22812299, doi: 10.1175/15200442(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, 305353, wos: A1983QS81300007.
Dickinson, R. E., et al. (2002), Nitrogen controls on climate model evapotranspiration,
Journal
of
Climate,
15(3),
278295,
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): 527535, doi:
10.1061/(asce)he.19435584.0000622.
Dorfman, M., Rosselot, K.S., 2011. Testing the Waters: A Guide to Water Quality at
Vacation Beaches, Twentyfirst 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.
Enviroweather, 2010. Enviroweather 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 biochemicalmodel of
photosynthetic CO2 assimilation in leaves of C3 species, Planta, 149(1), 7890, 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 systembased delivery model. Journal of
Environmental Quality, 27(4): 935945.
Fu, B. P. (1981), On the calculation of the evaporation from land surface, Chinese Journal
of Atmospheric Sciences (in Chinese), 5(1), 2331, 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),
39553978, 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.  GroundWater 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), 197200, 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): 16361644, 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): 117, 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), 21752194, 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(45),
421469, doi: 10.1007/s107120089037z.
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), 498509, doi:
10.1175/15200469(1997)054<0498:SHFFTE>2.0.CO;2.
LaBolle, E.M., Fogg, G.E., Tompson, A.F.B., 1996. Randomwalk simulation of transport
in heterogeneous porous media: Local massconservation problem and implementation
methods. Water Resources Research, 32(3): 583593, doi: 10.1029/95wr03528.
Li, K. Y., J. B. Boisvert, and R. D. Jong (1999), An exponential rootwateruptake model,
Canadian Journal of Soil Science, 79(2), 333343, WOS:000081215100012.
Lu, D., and Q. Weng (2006), Use of impervious surface in urban landuse classification,
Remote Sensing of Environment, 102(1–2), 146160, 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), 341370, 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), 4955.
Mancini, J.L., 1978. NUMERICAL ESTIMATES OF COLIFORM MORTALITYRATES
UNDER VARIOUS CONDITIONS. Journal Water Pollution Control Federation, 50(11):
24772484.
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 2009Nov28.
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 multiinstitution 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): 10491053.
Mu, Q., M. Zhao, and S. W. Running (2011), Improvements to a MODIS global terrestrial
evapotranspiration algorithm, Remote Sensing of Environment, 115(8), 17811800, 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, 519536, 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 FINITEELEMENT
METHOD FOR ADVECTION DISPERSION. International Journal for Numerical
Methods in Engineering, 20(2): 321337, 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), 937952,
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):
217235,
doi:
http://dx.doi.org/10.1016/00221694(75)900049.
Oleson, K. W., et al. (2010), Technical Description of version 4.0 of the Community Land
Model (CLM). NCAR Technical Note Rep. NCAR/TN478+STR, National Center for
Atmospheric Research, Boulder, CO., ISSNs: 21532397, 21532400.
Olson, S.R., Kemper, W.D., 1961. Movement of nutrients to plant roots. Adv.Agron., 20:
91151.
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),
534547, 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, WaterResources Investigations Report, 034056.
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 "RandomWalk" 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), Surfacesubsurface 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), 1930719327, 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), 159166, doi: 10.1007/s1004000601037.
RodriguezIturbe, I., and J. M. Mejí (1974), On the transformation of point rainfall to areal
a
rainfall, Water Resources Research, 10(4), 729735, 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 SubModel for Predicting Pathogen
Loadings in Surface and Groundwater at Watershed and Basin Scales, Total Maximum
Daily Load (TMDL) Environmental Regulatoins: Proceedings of the March 1113, 2002
Conference, Fort Worth, Texas, USA, pp. 5663.
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), 28952911, doi:
10.5194/hess1528952011.
Sayama, T., J. J. McDonnell, A. Dhakal, and K. Sullivan (2011), How much water can a
watershed store?, Hydrological Processes, 25(25), 38993908, 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 forestbased indicators?, J. Gt. Lakes Res., 39(2), 211223, ISSN: 03801330,
10.1016/j.jglr.2013.03.012.
Sellers, P. J. (1985), Canopy reflectance, photosynthesis and transpiration, Int J Remote
Sens, 6(8), 13351372, 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 subsurfaceland
surface processes model, Water Resources Research, 49, 121, doi:10.1002/wrcr.20189.
Shen, C. P., and M. S. Phanikumar (2010), A processbased, distributed hydrologic model
based on a largescale method for surfacesubsurface coupling, Advances in Water
Resources, 33(12), 15241541, 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, SurfaceWater
Hydrology. WileyInterscience, 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), 981998, 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(12), 189210, doi: 10.1007/s1058400792948.
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),
400407, doi: 10.1007/s0025400206571.
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 spinup: Estimating
steady state conditions in a coupled terrestrial carbon and nitrogen cycle model, Ecological
Modelling, 189(12), 2548, 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),
39023923, 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):
860869.
Tompson, A.F.B., Gelhar, L.W., 1990. Numerical simulation of solute transport in threedimensional, randomly heterogeneous porous media. Water Resources Research, 26(10):
25412562, 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 ResAtmos, 94(D13), 1628716301, 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), 1971919723, 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 surfaceapplied manures. Journal of Environmental Quality,
33(2): 749756.
van Genuchten, M.T., Alves, W.J., 1982. Analytical solutions of the onedimensional
convectivedispersion 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 35569 of Isotopes in Hydrology. Vienna,
International Atomic Energy Agency, 1967.; Other Information: From Conference on
Isotopes in Hydrology, Vienna. See STI/ PUB141; CONF 661133. Orig. Receipt Date:
31DEC68: Medium: X.
Wagener, T., M. Sivapalan, P. Troch, and R. Woods (2007), Catchment Classification and
Hydrologic Similarity, Geography Compass, 1(4), 901931, doi: 10.5194/hess1528952011.
162
Walker, S.E., Mostaghimi, S., Dillaha, T.A., Woeste, F.E., 1990. MODELING ANIMAL
WASTE MANAGEMENTPRACTICES  IMPACTS ON BACTERIA LEVELS IN
RUNOFF FROM AGRICULTURAL LANDS. Transactions of the Asae, 33(3): 807817.
Wang, D. (2011), Evaluating interannual water storage changes at watersheds in Illinois
based on longterm 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 FECALCOLIFORM
DYNAMICS IN STREAMS AND RIVERS. Water Research, 29(3): 847855, doi:
10.1016/00431354(94)00211o.
Wuebbles, D., and K. Hayhoe (2004), Climate change projections for the United States
Midwest, Mitigation and Adaptation Strategies for Global Change, 9(4), 335363, doi:
10.1023/B:MITI.0000038843.73424.de.
Zheng, C., Bennett, G.D., 2002. Applied Contaminant Transport Modeling. Wiley
Interscience, 621 pp.
163