A .. .. \fi fir, . . . .mewwwhmrx u. 3. .(WK. . unri— n‘uil l: I .139? : .Wm. I... 3.. .. . . up. afimunfiwf. J .flhfltalazyr} . at a» 1kg. .634, $509». . K9...l§ 24v .‘ 1 :19. 3:: i. 2. Jaguar. .- ‘ H {E w; ‘ Hm ,uirvbwtbré. i. I TH‘Z-‘STK- . ‘2 200% This is to certify that the thesis entitled EXPLORING RELATIONSHIPS BETWEEN LAND USE AND ECOHYDROLOGY USING MULTIVARIATE STATISTICS AND PROCESS-BASED MODELS presented by o 23 >4 :5 .3? g 5 § Nicklaus R. Welty a it; 4 .2 o 2 has been accepted towards fulfillment of the requirements for the Master of degree in Geological Sciences Science '._ ,. H. i /; wl’mwxv Major Péfessor’s Signature iZ/é 1/2005, Date MSU is an Affirmative Action/Equal Opportunity Institution PLACE IN RETURN Box to remove this checkout from your record. To AVOID FINES return on or before date due. MAY BE RECALLED with earlier due date if requested. DATE DUE DATE DUE DATE DUE 2/05 a/CIRC/DateDueindd-p. 1 5 EXPLORING RELATIONSHIPS BETWEEN LAND USE AND ECOHYDROLOGY USING MULTIVARIATE STATISTICS AND PROCESS-BASED MODELS By Nicklaus R. Welty A THESIS Submitted to Michigan State University in partial fulfillment of the requirements for the degree of MASTER OF SCIENCE Department of Geological Sciences 2005 ABSTRACT EXPLORING RELATIONSHIPS BETWEEN LAND USE AND ECOHYDROLOGY USING MULTIVARIATE STATISTICS AND PROCESS-BASED MODELS By Nicklaus R. Welty Although land use is known to affect stream ecosystem health, it can be difficult to directly relate stream health indicators, like dissolved oxygen, to the landscape. In this study we examined land use and stream ecosystem health indicators using statistical relationships and process-based models. Statistics were performed on synoptically collected baseflow water samples from 125 stream sites within four Michigan watersheds. We measured parameters like nitrate that are suspected to have land use influences. Principal components and correlation analyses both indicate strong associations between agricultural areas and increased stressor levels, while forested regions generally have lower stressor levels. Urban areas correlate with increased stream temperature and chloride concentrations. The land use-stream ecosystem relationship was explored with a coupled groundwater-surface water flow and nitrate transport and fate model of Cedar Creek, in southwestern Michigan. The flow model used process- based ecohydrologic principles by using NEXRAD precipitation and PET to calculate recharge. Modeling results indicate stream water nitrate has urban and agricultural sources. The majority of dissolved oxygen levels in Cedar Creek can be explained through this simulation of nitrate transport, highlighting the importance of understanding the sources and pathways of nitrate. ACKNOWLEDGMENTS I thank my advisor, David Hyndman, for his guidance, and my two other committee members, R. Jan Stevenson and Phanikumar Mantha for their contributions and suggestions. Anthony Kendall, Dushmantha Jayawickreme, and Amy Lansdale provided helpful insight. This work would not have occurred without the generous financial support of the Environmental Protection Agency. Supplemental funding for field sampling was provided by the Great Lakes Fisheries Trust and the National Science Foundation. Water level data were provided by Mike Wiley’s group at University of Michigan, and the most of the nutrient analyses were conducted by Lara Panayotoff. We also greatly appreciate the contributions of Bryan Pij anowski, Mike Wiley, and Lara Panayotoff. Finally, I thank my partner Katie for her patience and support the past two years. iii TABLE OF CONTENTS LIST OF TABLES ............................................................................................................. vi LIST OF FIGURES .......................................................................................................... vii CHAPTER 1: MULTIVARIATE STATISTICAL RELATIONSHIPS BETWEEN LAND USE AND STREAM ECOHYDROLOGY ........................................................................ 1 INTRODUCTION .............................................................................................................. 2 MATERIALS AND METHODS......................., ................................................................ 4 Study Sites .................................................................................................................. 4 Muskegon River Watershed .................................................................................... 6 Sycamore Creek ...................................................................................................... 8 Black Creek ............................................................................................................. 8 Grand Traverse Bay region .................................................................................... 9 Synoptic Sampling and Laboratory Analysis ........................................................... 10 GIS Analysis ............................................................................................................. 14 Land use ................................................................................................................ 14 Sourceshed generation .......................................................................................... 15 Riparian buffer generation ................................................................................... 16 Statistical Procedures ................................................................................................ 16 Order of Analysis ...................................................................................................... 19 RESULTS AND DISCUSSION ....................................................................................... 19 Total Sourceshed Analysis ........................................................................................ 19 Chemistry and total sourceshed land use ............................................................. 19 Dissolved oxygen, temperature, chemistry, and total sourceshed land use .......... 25 Additional variables .............................................................................................. 31 Riparian Buffer Analysis .......................................................................................... 33 Exploratory Analysis of Mass Fluxes Through Differential Sourcesheds ................ 38 CONCLUSIONS ............................................................................................................... 38 CHAPTER 2: SIMULATIN G THE ECOHYDROLOGY OF A SMALL MICHIGAN WATERSHED WITH MIXED LAND USE .................................................................... 40 INTRODUCTION ............................................................................................................ 41 METHODS ....................................................................................................................... 43 Study Region: Cedar Creek Watershed .................................................................... 43 Climatological, GIS, and Ecohydrological Data ....................................................... 45 Climate data .......................................................................................................... 47 Discharge data ...................................................................................................... 48 Water chemistry .................................................................................................... 49 Groundwater and Surface Water Models .................................................................. 51 Modeling ................................................................................................................... 52 General modeling approach ................................................................................. 52 Groundwater flow model and parameter optimization ......................................... 53 Transport model .................................................................................................... 61 iv Stream ecohydrology model .................................................................................. 64 RESULTS AND DISCUSSION ....................................................................................... 65 Groundwater Flow Optimization .............................................................................. 65 Nitrate Transport Model ........................................................................................... 72 Stream Ecohydrology Model .................................................................................... 77 Relative influence of land use practices ................................................................ 88 CONCLUSIONS ............................................................................................................... 90 REFERENCES CITED ..................................................................................................... 92 LIST OF TABLES Table 1: Land use proportions of the study regions ............................................................ 7 Table 2: Range and mean of variables examined. ............................................................ 21 Table 3: Correlation coefficients of land use and water chemistry. Correlations with a p- test value >99% in bold. The p-test value is the probability that the correlation occurred randomly if the true correlation is zero. ............................................................................ 22 Table 4: First three principal components of land use and chemistry. For clarity, loadings <-0.3 and >03 bolded. ...................................................................................................... 23 Table 5: First three principal components for DO, temperature, and total sourceshed land use. For clarity, loadings <—O.3 and >03 bolded. ............................................................ 27 Table 6: First three principal components for DO, temperature, and riparian buffer land use. For clarity, loadings <-O.3 and >03 bolded. ............................................................ 37 Table 7: Parameters for calculation of PET. ..................................................................... 57 Table 8: Parameter values for recharge model. ................................................................ 60 Table 9: Optimized recharge and conductivity parameters. Till conductivity was tied to outwash conductivity with a scaling factor of 0.6. ........................................................... 66 Table 10: Values for nitrate recharge concentrations used in MT3DMS. ........................ 74 vi LIST OF FIGURES Figure 1: Regions of study in Michigan. A: Grand Traverse watershed, B: Muskegon River watershed, C: Black Creek watershed, D: Sycamore Creek watershed. ................... 5 Figure 2: DO concentration patterns in the Sycamore Creek watershed are consistent for four different dates. ........................................................................................................... 1 1 Figure 3: Dissolved oxygen correlates positively with openland and forest and negatively with agriculture. ................................................................................................................ 28 Figure 4: Dissolved oxygen tends to decrease with increasing agricultural land use. Points with stars represent a 2.5 km reach of the Sycamore Creek system with suspected localized impacts. Point with circle has very low flow (0.016 fl3/s (0.00045 m3/s)). The solid line is a linear regression including the aberrant points and the dashed line excludes the points. .......................................................................................................................... 29 Figure 5: Stream temperature generally increases with urban land use. Crosses represent all sites, while open squares are sites with urban land use >10%. Solid line is a linear regression calculated from all sites, dashed line is a linear regression calculated from sites with urban land use >10% ................................................................................................. 30 Figure 6: A linear decreasing trend exists between total sourceshed forest land use and conductivity ....................................................................................................................... 32 Figure 7: Buffer agricultural land use correlates more strongly with SRP and TP than total sourceshed agricultural land use. .............................................................................. 34 Figure 8: Map of land use [MDNR, 2001] and hydrography with hillshaded topography within the Cedar Creek groundwatershed. The portion of Cedar Creek displayed in bold is the tributary simulated in QUALZK. The inset map shows the watershed location as a small black region in the southern Muskegon River watershed (gray outline) relative to the state of Michigan. The expanded surface watershed of the regional Muskegon River watershed is displayed in bold on the state map. .............................................................. 44 Figure 9: Conceptual model of recharge model. NEX: original NEXRAD precipitation, IC: infiltration capacity, PRE: precipitation that becomes groundwater recharge, PT: potential transpiration, PE: potential evaporation, PM: Penman-Monteith Equation, SWE: snow water equivalent, 20: previous day snow depth, Z: current day snow depth, p: current day snow density, p0: previous day day snow density, pf: fresh snow density, Q0: previous day stream discharge, Q: current day stream discharge, SUB: sublimation, SMO: previous day soil moisture, SM: current day soil moisture, FC: field capacity, RCH: current day recharge, RCHO: previous day recharge ......................................................... 54 Figure 10: Simulated vs. observed heads for the flow optimization. Open circles are head observations from the period of the model simulation (1988-2004), closed circles from prior to the start of the model (1969-1987) ....................................................................... 69 vii Figure 11: Simulated baseflow (dotted line), monthly extracted baseflow (points) and actual streamflows [Wiley and Richards, unpublished data] (line) for the upstream transducer on Cedar Creek. The optimization routine used the monthly extracted baseflow and the simulated baseflow for calculation of the objective function. Locations of this transducer is shown in Figure 1. ............................................................................ 70 Figure 12: Simulated baseflow (dotted line), monthly extracted baseflow (points) and actual streamflows [Wiley and Richards, unpublished data] (line) for the downstream transducer on Cedar Creek. Locations of this transducer is shown in Figure 1 ............... 71 Figure 13: Simulated vs. kriged observed groundwater nitrate concentrations (r2 = 0.22) in residential wells. Dashed envelope is i 2 mg L". Observed concentrations are from MDCH/MDEQ database. .................................................................................................. 75 Figure 14: Simulated groundwater nitrate concentrations by MT3DMS in the Cedar Creek watershed after 42 years of land use practices. ...................................................... 76 Figure 15: Measured DO concentrations in Cedar Creek on two separate sampling trips. QUAL2Kw simulated the conditions of 8/16/2005, as it was the trip water chemistry was collected. ........................................................................................................................... 78 Figure 16: QUALZKw simulated and observed streamflow for the 8/16/2004 synoptic sampling. ........................................................................................................................... 79 Figure 17: Simulated and observed water temperatures for the 8/16/2004 synoptic sampling. Error bars represent the accuracy of the YSI temperature measurement. An incoming groundwater temperature of 10°C ..................................................................... 80 Figure 18: Simulated and observed nitrate concentrations for the 8/16/2004 synoptic sampling. ........................................................................................................................... 82 Figure 19: Simulated and observed dissolved oxygen concentrations, and observed chlorophyll a concentrations for the 8/16/2004 synoptic sampling. Simulated dissolved oxygen saturation is dashed, and increases in the upstream portion of Cedar Creek where the temperature is dropping, as colder water can hold more dissolved oxygen. DO error bars represent the accuracy of the DO probe. ................................................................... 84 Figure 20: DO levels are much closer to observed values if reaeration is included in QUALZKw ........................................................................................................................ 85 Figure 21: Simulated algae areal concentrations and observed chlorophyll a concentration for the 8/16/2004 synoptic sampling ................................................................................. 87 Figure 22: The influence of different land use practices can be examined by simulating different land use scenarios in MT3D and running QUALZKw with the results. ............ 89 viii CHAPTER 1: MULTIVARIATE STATISTICAL RELATIONSHIPS BETWEEN LAND USE AND STREAM ECOHYDROLOGY INTRODUCTION Land use threatens the integrity of surface waters [Malmqvist and Rundle, 2000], reduces biodiversity [Sala et al., 2000], and disrupts key hydrologic processes [DeFries and Eshleman, 2004]. The transport of nutrients from agricultural and urban land uses is one way land use can degrade stream ecosystems. Hydrologic investigations provide critical information to help link regional land uses to stream integrity. A study by the USEPA [2000] determined that almost 40% of 5.79 million kilometers of US streams and rivers were too polluted to support fishing and swimming, mainly due to agricultural inputs. The Michigan Department of Natural Resources considers “Absence of Land Use Planning” a severe environmental risk [Rustem et al., 1992]. Understanding the causes of water quality impairment allows informed land management decisions that will preserve the aesthetics and biotic integrity of our hydrologic resources for future generations. Frequently, these decisions deal with the mitigation of common stream stressors, including low DO levels, increased temperature, and elevated nutrient concentrations, all of which can negatively impact stream biota. In this paper, we refer to the stressors measured, which include DO, temperature, nutrients, and chlorophyll-a, as ecohydrological variables. Ecohydrology (or hydroecology) is the interdisciplinary study of the interaction between biota and hydrology at multiple spatial and temporal scales [see Hannah et al., 2004; Porporato and Rodriguez-[turbe, 2002]. Controlling these ecohydrological variables is vital, as they harm stream organisms (e. g., low oxygen levels killing fish) as well as humans (e.g., nitrate-induced methemoglobinemia). Land use proximal and distal to a stream can also damage physical habitats by altering runoff and sediment yields [Jacobson et al., 2001]. One approach for linking land use and ecohydrology is synoptic baseflow sampling. Synoptic studies, where samples are collected from many sites in a short period of time, allow researchers to examine the spatial patterns of surface water chemistry [Grayson et al., 1997]. Baseflow conditions are when groundwater is the dominant source of stream water, and are often used in harmony with synoptic sampling. During baseflow the stream chemistry is largely a function of the incoming groundwater chemistry, and the effects of runoff are negligible, so synoptic sampling reveals land use- water quality patterns. Land use proximal and distal to a stream have unique impacts on stream ecohydrology. A common approach in land use-water quality studies is to compare the water quality at a point to the land use of the watershed. We used a more refined technique of generating a watershed for each sampling point, termed sourcesheds [Wayland et al., 2003]. Land use within the riparian buffer surrounding the stream was also analyzed in order to determine the influence of spatial distribution of land uses. The overall goal of this study is to link ecohydrological variables to land uses, using different land use area descriptions and ecohydrological formulations. First, we compare sourceshed land use and ecohydrological concentrations. Second, we analyze riparian buffer land use and ecohydrological concentrations. Finally, we examine sourcesheds and ecohydrological mass fluxes, rather than concentrations. MATERIALS AND METHODS Study Sites For this investigation, sites were selected across Michigan with a range of land uses that are thought to impact streams, via overland flow and shallow groundwater pathways (Figure 1). We evaluate relationships between land use and stream ecosystem stressors for these sites using an approach based on synoptic water quality data, high- resolution 1997-2000 land use data, and multivariate statistics. While most studies have analyzed land uses and major ion chemistry, we consider the influence of land uses on ecohydrological stressors. Diurnally fluctuating variables were monitored pre-dawn and the latest land use data were used. Land use proportions are evaluated along with water quality variables using multivariate statistics. Soils in these regions are well-drained, so tile drains are not used. o 20 40 so N '— Kilometers 0 Sampling s'tes -Srearns - Urban/Roads Agri cultu re .. E 3i Fuest/grasslmd/Wetlaids mwi River watershed, C: Black Creek watershed, D: Sycamore Creek watershed. Muskegon River Watershed The Muskegon River watershed (Figure 1) is a critical fishery for the Great Lakes, and has already lost two native species, the Arctic grayling (Thymallus arcticus) and the Great Lakes muskellunge (Esox masquinongy) [0'Neal, 1997]. Land use impacts like logging and dam construction are cited as a cause for this disappearance. Other fish may be in danger as well: in 1990, DO levels fell below the Michigan Water Quality Standards level of 7 mg/L 50% of the time in summer downstream of Croton Dam. Furthermore, the watershed also has a large number of state-designated trout streams, and two of the state’s 52 Blue Ribbon Trout Streams [0'Neal, 1997]. A variety of glacial sediments drape the 7057 km2 watershed, with most areas covered by glacial tills, glacial outwash, and lacustrine sands and gravels [Farrand and Bell, 1982]. The primary bedrock underlying the glacial sediments of the watershed are Mississippian and Pennsylvanian sandstone and shale [MDEQ, 1987]. Land use across the watershed is dominated by agricultural fields, forests, and a few significant urban areas with populations between 10,000 and 40,000: Muskegon, Cadillac, and Big Rapids (Figure 1 and Table 1). Several researchers have studied the result of land use practices in the watershed. Tang et a1. [2005] predicted that this watershed will see increases in runoff, nutrients, oil and grease, and heavy metal concentrations due to projected urban expansion of runoff and non-point source pollution, based on a land use change model and an empirical environmental impact model. Watershed % Agriculture % Forest % Urban % Open/ % Wetlands % Open bare water Muskegon River 18.3 40.6 3.6 16.0 17.4 4.1 Grand Traverse 18.0 41.5 4.0 19.9 9.7 6.9 Sycamore Creek 51.3 12.5 16.7 8.5 10.6 0.4 Black Creek 70.7 9.6 3.6 5.4 9.0 1.7 Table 1: Land use proportions of the study regions. Sycamore Creek Sycamore Creek is well suited for our analysis as it flows through a mélange of agricultural, forested, and urban land (Figure 1 and Table 1). Corn and soybeans are the dominant crops in the region [Lombardo et al., 2001]. Previous studies of Sycamore Creek concluded that sediment, phosphorus, nitrogen, and agricultural pesticides are the major pollutants [Lombardo et al., 2001; Suppnick, 1996]. These studies also suggest that Sycamore Creek may serve as a proxy for many southern Michigan streams. The upstream portions of this 275 km2 watershed are dominated by agricultural land which usually contains small forested buffers around streams, while the downstream portion flows through the suburban outskirts of southern Lansing before discharging into the Red Cedar River (Figure 1). The dominant surficial geology} surrounding Sycamore Creek is medium-textured glacial till [Farrand and Bell, 1982], while bedrock mainly consists of Pennsylvanian interbedded sandstone and shales of the Saginaw Formation [MDEQ, 1987]. Black Creek The Black Creek watershed (Figure 1) was chosen as a study basin due to the expected stressed ecosystem conditions related to a confined animal feeding operation (CAFO). Stream habitats in the Black Creek catchment are impaired relative to other streams in the region [Roth et al., 1996]. Healthy streams in southeastern Michigan may have the highest biological diversity in the state [Allan et al., 1997], so it is important to link impaired water bodies with practices causing impairment. The surface sediments of this 102.3 km2 region studied are dominated by fine textured glacial till [Farrand and Bell, 1982], which overlie the Mississippian Coldwater Shale [MDEQ, 1987]. Approximately 70% of the land is agricultural, with row and forage crops (Table 1). Grand Traverse Bay region The Grand Traverse Bay watershed, which covers 3275 km2 of the northwestern Southern Peninsula of Michigan (Figure 1), was chosen because of its rapid urbanization. Over 75,000 people live in Grand Traverse County, centered around Traverse City, the largest city in the watershed. This county has one of the fastest grth rates in the state, with population projected to increase by over 60% by 2025 [Grand Traverse County, 2002]. The land use changes associated with this influx should be examined in the context of water quality and ecosystem integrity to prevent potential detrimental effects. For example, endangered species like the Bald Eagle (Haliaeetus leucocephalus) reside in the area, but land use change is threatening their further survival [Grand Traverse County, 2002]. Protecting the Boardman River, the largest river in the watershed, is particularly important as it is designated to be managed as a “Natural River” by the Michigan Department of Natural Resources and is one of Michigan’s top ten trout streams [MDNR, 2002]. Surficial glacial sediments in the watershed include till, dune sand, glacial outwash, and lacustrine sand and gravel [Farrand and Bell, 1982]. These overlie bedrock that is mainly Devonian and Mississippian shale [MDEQ, 1987]. Forests are the dominant land use followed by agriculture (Table l). Synoptic Sampling and Laboratory Analysis A strategy of sampling streams at multiple points over a short time interval with relatively constant hydrologic conditions was adopted during the summer of 2004. This synoptic sampling approach is advantageous because it provides a ‘snapshot’ of water quality at a particular instant in time [Grayson et al., 1997]. Baseflow conditions, where groundwater inputs dominate streamflow, exist during the late summer in most Michigan streams. These conditions were used for synoptic sampling because the water quality is indicative of regional land use impacts [Wayland et al., 2003], as the solutes are carried by groundwater to the streams. Regional hydrographs were used to ensure baseflow conditions were present. Although the data analysis presented here represents a single summer synoptic sample set, the relationships established here are not limited by the temporal range of the dataset. Ecohydrological variables were measured for more than one synoptic event on four streams, and only minor differences were identified from one sampling event to the next (Figure 2). 10 10 9 u G"“. N. ._ Q ‘~ , O ‘ T _ _ 7‘ T‘ 8’ .:T\TP.HD7TT;\ ‘ . 0“. I g -- s 7 -- , ‘ *~ .. - " E. 6 __ Mt .~« .A, L +-7/16/05 ; o - 2;. i + 8/5/05 - a 5 l (i. » + 8/10/05 1 4 «g l __-+ 3/13/05 ] 3 , . . 0 1o 20 30 4° 50 Distance downstream, km Figure 2: DO concentration patterns in the Sycamore Creek watershed are consistent for four different dates. 11 13‘" Ecohydrologic data for multivariate analysis must be collected with attention to the temporal variability in these variables. For example, DO varies diumally with maxima late in the afiemoon and minima early in the morning. Early morning DO provides a temporally-refined indicator of DO stress. If DO measurements from morning to afiemoon are combined in analyses, correlations between land use and DO stress will be weakened. Correlations between land use and low D0 are vital to establish land use- stressor relationships, as oxygen levels are critical to the overall ecosystem health of streams. Likewise, water temperature was monitored near dawn to capture the stream temperature independent of variation in factors such as local shading and temporal fluctuations in solar radiation. Only early morning DO and water temperature measurements were used in this study. The ecohydrological stressors we measured were selected because they provide information linking land use practices to ecosystem health. Water chemistry was tested for chemicals that can provide land use and ecosystem stress information. Nitrogen and phosphorus, and their associated species, are the major nutrients responsible for eutrophication of surface waters, and have land use sources [Carpenter et al., 1998; Jenkinson, 2001]. Silica was measured as to provide information on diatom activity. Chloride was also measured, as it is primarily derived from deicing agents [Amrhein et al., 1992, NRC, 1991], and has been correlated to urban land use in Michigan [Boutt et al., 2001, Wayland et al., 2003]. DO levels were monitored, since they determine the overall health of an lotic ecosystem, as it is necessary for all invertebrate and fish life [Diaz, 2001]. For example, fish grow less [Bejda et al., 2001], consume less food [Bejda et al., 2001; Chabot and Dutil, 1999], and have increased mortality [Plante et al., 1998] 12 at low DO levels. Water temperature was also monitored, as it has dramatic effects on stream biota like fish, whose body temperature is largely a fimction of the surrounding water temperature [Coutant, 197 6]. Temperature can be both directly and indirectly lethal to fish [Coutant, 1976], as well as disrupt key processes like migration, spawning, and resistance to disease and pollutants [Armour, 1991]. Ecosystem stressors were evaluated at the upstream side of bridges (when possible) using water quality probes and discrete samples collected and analyzed in the laboratory. Dissolved oxygen, water temperature, pH, and conductivity were measured using Yellow Springs Instruments (Yellow Springs, OH) multi—parameter water quality probes. Each probe was calibrated for DC with a water-saturated air technique before sampling. Dissolved oxygen sampling occurred within the period either two hours before dawn or two hours after dawn to identify locations with oxygen stress. Water samples were collected in acid-washed 75 mL polyethylene bottles and flash-frozen on dry ice in the field to prevent chemical transformations and biological degradation. Sample bottles were rinsed three times with stream water before collecting a sample. These samples were analyzed to determine the concentrations of total nitrogen, NO3-N02 (hereafter referred to as NOx), NH3, total phosphorus (TP), soluble reactive phosphorus (SRP), SiOz, and Cl’. Water samples were analyzed according to standard methods (APHA, 1998) to determine the concentrations of total nitrogen (TN), NO3-N02 (N Ox), NH3, total phosphorus (TP), soluble reactive phosphorus (SRP), SiOz, and Cl' [Panayotoffi unpublished data]. Water was filtered prior to analyses for SRP, NH3, and Cl'. SRP was assayed with the ascorbic acid method. TP was assayed as SRP following persulfate digestion. TN was assayed as NOx by second derivative UV spectroscopy 13 following persulfate digestion. Four assays were run on a Scalar Sans-Serif segmented- flow autoanalyzer: NOx by automated cadmium reduction, NH3 by the automated phenate method, SiOz by the molybdosilicate method, and C1' by the automated ferricyanide method. 100 mL of stream water was filtered for chlorophyll assay through 0.7 micron Whatman GF/F filters in the field. The filter was packed in aluminum foil and frozen, also in the field. Chlorophyll-a concentrations were determined by extracting pigments for 24 hours in 90% ethanol, sonicating the sample for 15 minutes during the first hour of extraction, fluorometric assay with a Turner Designs TD700, and correcting for phaeophytin with acidification. A Marsh-McBirney Flo-Mate 2000 (Frederick, MD) was used in tandem with a top-setting wading rod to measure stream velocities and calculate total discharge. Logistics prevented measurement of all variables at all 125 sites. Water chemistry was analyzed for all 125 stream sites, while early morning DO and water temperature were measured at a subset of 74 sites. Suspended chlorophyll was collected at 41 sites, while conductivity and pH were collected at 47 sites. Stream discharge was collected at 93 of the 125 original sites, 28 of which also had early morning DO samples. GIS Analysis Land use Recent technological advances in remote sensing and geographic information systems (GIS) have facilitated the study of the land use and hydrology of large watersheds. Researchers in Michigan frequently use the Michigan Resource Information System land use dataset [e. g. Johnson et al., 1997; Wayland et al., 2003], which is based l4 on aerial photography from the late 1970’s. However, this study uses a new high- resolution land use dataset from the Michigan Department of Natural Resources, which is based on 1997-2000 Landsat data [MDNR, 2001]. The Integrated Forest Monitoring Assessment and Prescription (IFMAP) is a detailed digital land use map of Michigan that was developed from a composite of Landsat imagery from spring, summer, and fall of 1997, 1999, and 2000 [MDNR, 2001]. Land uses are broken down into 30 m cells categorized using a hierarchical classification scheme. Although IFMAP defines 33 distinct land use categories, these were aggregated into six categories for analysis: (1) agriculture, (2) forest, (3) urban, (4) openland and barren, (5) wetlands, and (6) open water. Sourceshed generation To evaluate land use proportions, sourcesheds were generated using ArcInfo [ESRL 2003]. A sourceshed is the land area that drains to a particular point in a stream. Sourcesheds are important to consider in ecohydrological studies, because the water quality of a sampling point is a function of the upstream land uses. Rather than use pre- existing watershed boundaries, a surface water sourceshed was generated for each of the 125 sampling points visited for this study using the United States Geological Survey 30 m National Elevation Dataset (NED) Digital Elevation Model (DEM). Two types of sourcesheds were used: total and differential. A total sourceshed represents all the land draining to a point, while a differential sourceshed represents the land draining to a point since the last sampling point. These sourcesheds were intersected with IFMAP land use data and the land use percentages within each sourceshed were calculated for use in statistical analyses. 15 Riparian bufler generation Riparian buffers, the region of vegetation surrounding a stream, can reduce the flux of nutrients like nitrate [Hill, 1996], as well as moderate stream temperature [Osborne and Kovacic, 1993]. Burkat et al. [2004] suggest riparian buffers are especially important around small streams, where they have a greater potential for runoff and groundwater interception. Therefore, we examined the land uses proximal to streams to explore these localized land use impacts on water quality. The MDEQ recommends a riparian buffer of vegetated land of at least 100 fi on each side of a stream [MDEQ, 1997]. As per the recommendation, the land use proportions within a 30 m buffer on each side of the stream were calculated within each differential sourceshed. Statistical Procedures Multivariate statistics are commonly used to examine the factors that influence water quality [Brown, 1998]. Specifically, Johnson and Gage [1997] recognized that the combination of GIS and multivariate statistics yields insight into complex relationships between land use and stream ecology. Principal components analysis (PCA), one multivariate approach, has been widely used in water quality studies. Mendiguchia et al. [2004] characterized water quality in the Guadalquivir River (Spain) through PCA and suggested land use practices explained the variation in water quality reflected in the principal components. However, these authors did not explicitly include land use proportions in the PCA data matrix. Other authors have documented success when including land use as a variable in PCA. In Michigan, Wayland et al. [2003] conducted a statistical study of water quality and land use in the Grand Traverse Bay region through 16 factor analysis, a technique similar to PCA. Key correlations between land use and water quality were identified, including increased agriculture with elevated Ca+2, Mg”, and alkalinity, and higher Na+, KL, and Cl' with urban areas. The agricultural association was likely due to dissolution of soil minerals during soil cultivation and the urban signature was likely due to road salt application. Wayland et al. [2003] suggested stronger correlations would have been determined if only early moming DO measurements were included; high-resolution land use data may also have improved correlations. Wayland et al. [2003] concentrated on characterizing the major cations and anions associated with land uses, and examined relatively few nutrients. In contrast, we focus on the nutrient chemistry and other stressors such as DO and temperature, as they have more direct effects on stream ecosystems. Johnson et al. [1997] also examined land use-water quality relationships in Michigan using multivariate methods. They found that agricultural land proportion within the watershed strongly correlated with nitrate. Some variables, like total phosphorus correlated best with land use within a 100 m stream buffer, while total nitrogen and nitrate were associated with total watershed land use. Land use-water quality relationships in the Raisin River watershed, in southeastern Michigan, were the target of a study by Allan et al. [1997]. Subcatchment land uses, particularly the proportion of agricultural land, were strongly associated with stream quality variables and buffered land uses were weakly correlated in their study. Ecohydrological variables such as DO, stream temperature, and suspended chlorophyll were not considered by Johnson et al. [1997] or Allan et al. [1997]. 17 In this study, relationships between land use and ecohydrologic variables were evaluated using both Pearson Product Moment Correlation (hereafter referred to as correlation) and principal components analysis (PCA). The former technique measures the strength of the linear individual relationships between water quality variables and land uses, while the latter identifies the major, orthogonal variable groupings in large datasets with multiple underlying forces [Brown, 1998]. Values of the correlation coefficient near 5:1 indicate strong linear relationships, while values near 0 indicate no linear relationship. The goal of PCA is to reduce the dimensionality of a dataset and extract a minimum number of factors that explain the maximum amount of the data variance [Brown, 1998]. Here, PCA was used to evaluate the influence of land use on stream ecosystem stressors. To perform PCA, land use percentages and water quality observations (including water temperature, suspended chlorophyll, DO, and nutrient concentrations) were placed in a data matrix. The matrix was standardized by subtracting the mean of each dataset from the observation and dividing by the standard deviation of each observation set (i.e., the mean water temperature was subtracted from each individual water temperature measurement, and then each individual water temperature measurement was divided by the standard deviation of all the water temperatures). The principal components of this standardized matrix were then calculated in MATLAB [The Mathworks, 2002] using the statistics toolbox. Traditionally, stream chemistry is analyzed in a concentration form. Here, this convention was used, but we also examined solute mass fluxes, calculated as the product of stream discharge and solute concentration. In some cases, it should be possible to narrow the range of processes that affect solute concentrations by using available mass 18 flux data. For example, concentrations frequently decrease downstream, but it is difficult to determine whether dilution via tributaries/ groundwater or a bio geochemical process is responsible. If the mass flux decreases, some biogeochemical processes are responsible. For a conservative constituent like Cl‘, the total mass flux should not decrease in gaining stream environments. Order of Analysis First, stream chemistry and total sourceshed land use were examined for the entire 125 site dataset using correlation analysis and PCA. Second, these analyses were conducted on the data from 74 sites where DO and temperature could be added to the chemistry and land use data. Third, suspended chlorophyll and conductivity were evaluated relative to total sourceshed land use for the 41 available sites. Fourth, the influence of 30 m riparian buffered land uses on stream ecohydrology were evaluated. Finally, differential sourceshed land use and ecohydrologic mass fluxes were investigated. RESULTS AND DISCUSSION Total Sourceshed Analysis Chemistry and total sourceshed land use Correlation coefficients and PCA loadings were calculated using data from all 125 sample sites within the four watersheds (Tables 2, 3 and 4). Urban, agricultural, and forested land use proportions all were distinctly correlated with various ecohydrological variables. Cl' concentrations strongly correlated with urban land use areas, likely due to 19 the deicing agent NaCl and septic tanks in suburban regions. Agriculture and forested land uses presented inverse signatures: agriculture is associated with higher P and N, while forests correlate oppositely. The main source of both N and P is likely fertilizer from agricultural regions, and additional sources of N include fixation by legumes and mineralization of soil organic N [Jenkinson, 2001]. Petery'ohn and Corell [1984] suggested the considerable nutrient reduction capability of forests, and similar processes may be influencing such correlations, including uptake and denitrification. However, forested regions may simply have lower nutrient source concentrations relative to the anomalously high concentrations of agricultural areas. Silica concentrations are negatively correlated with open water, perhaps due to diatom utilization of Si. Regions with open water generally correspond to lakes or broad reaches in the river network, enabling diatoms to flourish and deplete Si concentrations. While the mechanics of Si depletion in lakes are well known [e.g., Schelske and Stoermer, 1971], there has been far less study of this phenomenon in rivers and wetlands. Recently, Stevenson et al. [In press] correlated Si depletion in rivers in Kentucky with increased benthic chlorophyll a resulting from N and P enrichment. 20 Variable Range Mean DO (% sat.) 2.00 - 118 69.8 Water temperature (°C) 10.7 - 20.0 14.9 Conductivity (p8) 218 - 938 604 pH (-log [H+]) 6.95 — 9.07 7.96 Total nitrogen (pg L") 57.0 - 5390 1220 NO3-N02 (pg L“) 11.0 - 5670 752 NH; (pg L") 10.0 - 469 34.1 Total phosphorus (pg L']) 2.62 - 181 33.0 Soluble reactive phosphorus (11g L'I) 1.00 - 84.6 8.07 Cl' (mg L“) 1.50 - 82.4 20.9 SiOz (mg L") 2.17 - 16.5 9.13 Suspended chlorophyll (pg If) 2.59 - 1280 38.8 Table 2: Range and mean of variables examined. 21 .83 m. scum—oboe on... on. b 382...... 3.5.8.. cote—8.50 on. .3. bzfimnoa on. mm 2:? $8-.. on? .20.. E .xéoA 2:...» .8... m 5.3 82.28.80 .bmmfioco .033 new 8: new. we 35.0508 22.2880 ”m o3...- ...... o...- 2...- mo..- :...- 3...- 2...- .5...- mm...- 2... cm... 3.... 8...- .m o... :.... .2... 3.... «m... :... 2...- 8... ms..- mu...- me... 8... .6 8.. an... E... on... we... 8...- 2.... S...- 3.... :.... om... 7F 8.. .2... :.... an... .8... mm... 8.... 8...- 2... :... £2 8.. 2... :.... .5... 2.... R...- wo...- mm... E... .52 8.. E... 8...- 3... 3....- ...o. R... o..... n:- ...... . 2.- 2... 3....- _ _...- :.... 8...- $5 8.. 8... :.... N... no..- 8...- 85:63 2... 8...- S...- N.....- .3... .283 8.. a... 2....- mm...- .85. 8.. S... z... 22.5.5 2... 2... 23.3.3. 8.. :35 2.3 .m -6 29 £2 .52 A: ham 85.83 .283 .85... 22.0 6.3.3:»... :85 22 Component Variable 1 2 3 Urban -0.20 0.68 -0.09 Agriculture -0.44 -0.29 0.21 Open/bare 0.22 0.52 -0.06 Forest 0.46 0.03 -0.09 Water 0.46 0.03 -0.09 Wetlands -0.07 -0.07 -0.74 TP -0.25 0.01 -0.10 TN -0.30 0.14 0.07 Cl' -0.33 0.36 0.01 SiOz 0.20 0.16 0.60 % Explained 40. 9 54. 6 67. 4 Table 4: First three principal components of land use and chemistry. For clarity, loadings <-0.3 and >03 bolded. 23 PCA identified complex relationships between land use proportions and ecohydrologic variables. Over 40% of the dataset variance is explained by TN, Cl", and agriculture loading opposite forest and water (Table 4). This component can be interpreted again as fertilizer inputs, fixation, or mineralization. The second component exhibits a strong urban, openland/bare and Cl“ association, inferred again as either road salt or septic signature. An additional 13% of the variance is explained by an inverse relationship between wetlands and Si, which is similar to the open water-Si relationship. The PCA results displayed some of the same relationships present in the simple correlation matrix, but illuminates situations where a dominant process appears to affect multiple correlations. The first component is interpreted as the signature of agriculture, which causes correlations among nutrients and land uses, while the second component is representative of urban land uses, which causes correlation with Cl'. Boutt et al. [2001] and Wayland et al. [2003] attributed surface water Cl' concentrations in the Grand Traverse Bay region mainly to road salt, with septic systems and oil field brines as secondary sources. Here, Cl' concentrations correlated strongly with urban regions and moderately with agricultural regions. The source of urban Cl' can be further explored by evaluating roads separately (urban-roads) from other high to low intensity urban areas (urban-excluding roads), a valuable feature of the IFMAP land use database. A correlation with only urban-roads would likely indicate road salt the main Cl' source; however we found moderate correlations between Cl' and both types of urban land, with R=0.6 for urban-roads and R=0.5 for urban-excluding roads. This suggests both road salt and septic systems provide significant Cl' inputs to the study watersheds. While strong correlation exists between urban-roads and Cl’ concentrations, several sites 24 with low urban proportions had high Cl' concentrations. Some of these sites occur in regions with oil and gas fields. These oil and gas fields commonly occur in agricultural areas, which could cause the agriculture-CT correlation. Although seemingly innocuous, Cl' can be important to control in some regions as it can increase the mobilization of heavy metals [Amrhein et al., 1992], harm roadside vegetation, and be linked in high concentration to incidence of hypertension in humans [NRC, 1991]. Dissolved oxygen, temperature, chemistry, and total sourceshed land use Prior to analysis with land use, correlations among ecohydrological variables were examined. DO was found to have a weak negative correlation with SRP and NH3, which may be a manifestation of SRP stimulation of autotrophic and heterotrophic microbial activity and nitrification. Water temperature was not significantly correlated with other water chemistry variables. Dissolved oxygen was better related to land use than nutrient concentrations in our sample sourcesheds. Correlation and PCA indicate land use linkages between DO and temperature (Figures 3 and 5, Table 5). Dissolved oxygen, openland/bare, and forest proportions loaded opposite agriculture in the first PCA component, which accounted for nearly 39% of the variance (Table 5). Dissolved oxygen correlated with agriculture and combined openland/bare and forest (r = -0.51 and 0.53, respectively) (Figures 3 and 4). DO levels near saturation only occur at sites with very low proportions of agriculture (<5%) while sites with intensive agriculture have DO saturations in the range of 20—80% (Figure 4). Only two land use categories display a significant influence on water temperature: urban (Figure 5) and wetlands/water. The correlation between urban land 25 use proportions and stream temperature is much stronger when the evaluation is made with sites that have greater than 10% urban land use proportions (Figure 5). Baseflow reduction, increased reception of shortwave solar radiation due to removal of vegetation, and channel widening have all been suggested as causes of increased water temperature in urban areas [LeBlanc et al.,1997]. Increased residence times of stream water in wetlands provides greater thermal inertia, a greater surface area for radiation absorption, and a proportionately smaller input of colder groundwater, which warms water naturally in these habitats. 26 Component Variable 1 2 3 DO -0.44 0.1 1 -0.25 Temperature 0.08 0.63 0. 18 Urban -0. 12 0.57 -0.40 Agriculture 0.57 —0.1 1 -0. 13 Open/bare -0.49 -0.01 -0. 14 Forest -0.45 -0. 13 0.35 Water 0.03 0.34 0.73 Wetlands -0. 15 -0.34 0.21 % Explained 38. 9 61. 2 76. 8 Table 5: First three principal components for DO, temperature, and total sourceshed land use. For clarity, loadings <-0.3 and >03 bolded. 27 0.8~***** 777777 7~~7 77777 ., Open/bare Forest Urban - , , Temperature Wetlands Water 1 l l l C Agriculture Figure 3: Dissolved oxygen correlates positively with openland and forest and negatively with agriculture. 28 raw O 20 4O 60 80 100 o . 't r _ — — _ - y = -0.384x + 81.752 I'M'w" " ° y = -O.3792x + 84.885 R2 = 0.2565 R2 = 0.3884 Figure 4: Dissolved oxygen tends to decrease with increasing agricultural land use. Points with circles represent a 2.5 km reach of the Sycamore Creek system with suspected localized impacts. Point with square has very low flow (0.016 ft3/s (0.00045 m3/s)). The solid line is a linear regression including the aberrant points and the dashed line excludes the points. 29 9 1 é i 3 . 2 ‘ o i g + .2 3:9 + + All Urban 0 Urban >10% 4; y = 0.0252x + 14.4 y = o_059x + 14.0 i + = 0.153 2 = l 1 0 r r 0.595 ] 0 20 40 60 80 “/6 Urban Figure 5: Stream temperature generally increases with urban land use. Crosses represent all sites, while open squares are sites with urban land use >10%. Solid line is a linear regression calculated from all sites, dashed line is a linear regression calculated from sites with urban land use >10%. 30 Additional variables Chlorophyll was negatively correlated with D0, which is likely a result of high decomposition rates of accumulating algae and DO depletion by bacterial respiration. A negative correlation between DO and chlorophyll was also recognized by Stevenson and White (1995). In PCA, chlorophyll loaded positively with temperature and water/wetlands. Sabater et al. [2000] examined chlorophyll and water quality variables and identified a similar positive linear relationship between chlorophyll and temperature. Conductivity was strongly correlated with agriculture (R = 0.69) and forest (R = - 0.78) (Figure 6). Conductivity also correlated positively with Cl', a relationship in Michigan streams identified previously by Stevenson et al. [In press], and had a negative correlation with DO. The positive correlation between conductivity, a proxy for TDS, and agriculture may be due to soil disturbance and dissolution of newly exposed minerals in agricultural fields. 31 1000 w -—— —~ _ ”I- ,I- in»; J :——-—--—— -_—-_____-_.-__ ___- l i : y=-11.821x+ 889.63 1 1 ° ‘9’; R2 = 0.6141 1 m 7504| ° 1 :1 l 9 l a 1 l §5004fi l :1 l l 2 = l O l i 0 250 ~: ] F l I l 0 444-4474444444444—444—4 0 25 50 75 "/6 Forest Figure 6: A linear decreasing trend exists between total sourceshed forest land use and conductivity. 32 Riparian Buffer Analysis Agricultural land use in the riparian buffers was more strongly correlated with SRP and TP than agricultural land use proportions over the total sourcesheds (Figure 7). This result is consistent with the conclusion Johnson et al. [1997] reached when examining a 100 m buffer and phosphorus levels in Michigan. Both studies demonstrate the effectiveness of riparian buffers on water quality. A sourceshed may have high agricultural land use, but if a forested buffer is present, some of the P may be intercepted. 33 0.50 lTotal Sourceshed l ”flgRiparian Buffer _.._._—_..._____ Ma____._ _.._.__. _l n: 0.25 - 0.00 - SRP Figure 7: Buffer agricultural land use correlates more strongly with SRP and TP than total sourceshed agricultural land use. 34 Only a weak negative correlation was observed between forest proportions in stream buffers and nutrient concentrations. This is likely due to the dominance of groundwater in Michigan streams during baseflow conditions. Riparian buffers, which tend to be forested, can reduce nitrate in shallow groundwater up to 90% [Osborne and Kovacic, 1993]. Runoff can carry high levels of nutrients that are intercepted by buffers. During baseflow, however, there are minimal runoff contributions to streams. It is likely that a strong negative correlation will exist between forest buffers and N concentrations during wet periods, because overland flow is more likely to transport nutrients to a stream than groundwater. The weak negative correlation is likely due to the reduction of nutrients by buffers in shallow groundwater like Osborne and Kovacic [1993] reported. Such a negative correlation may not exist for P because it is commonly believed to be introduced with sediments and then recycled, a conclusion also reached by Johnson et al. [1997] Analysis of buffered land use proportions also demonstrates the importance of proximal land uses to stream temperature. Earlier PCA of DO, temperature, and sourceshed land uses identified urban areas and water/wetlands as the only significant controls on stream temperature. However, PCA of DO, temperature, and riparian buffer land use suggests forested buffers affect temperature more than urban or water/wetland land uses. PCA indicates that an inverse relationship between forested land use buffer proportion and temperature account for nearly 20% of the variance (Table 6). Water/wetlands and urban proportions still exhibit moderate negative loadings with temperature (Table 6). Thus, the spatial distribution of forests along buffers can be more important than the ratio of forest to non-forest land use. These results suggest that stream 35 temperature will remain at lower levels provided a dense canopy surrounds the stream, results similar to conclusions reached by Osborne and Kovacz'c [1993]. Thus, researchers interested in stream temperature and land use should focus on the land use immediately surrounding the stream as well as examining land uses across the watershed. 36 Component Variable 1 2 3 DO -0.57 -0. 12 0.18 Temperature 0.06 -0.64 -0.24 Urban -0.21 -0.45 -0. 10 Agriculture 0.60 -0.05 0.36 Open/bare -0.34 -0.24 0.26 Forest -0.37 0.30 0.32 Water 0.1 1 -0.41 0.02 Wetlands -0.1 1 0.24 -0.78 % Explained 25.0 44.9 62. 6 Table 6: First three principal components for DO, temperature, and riparian buffer land use. For clarity, loadings <-0.3 and >03 bolded. 37 Exploratory Analysis of Mass Fluxes Through Differential Sourcesheds Mass fluxes of SRP, TP and C1' are all moderately correlated with urban land uses (r = 0.53, 0.65, 0.64, respectively). SRP and TP mass fluxes correlate with urban areas likely due to both higher runoff and phosphorus inputs from suburban septic systems [Carpenter et al., 1998]. We also examined the mass fluxes of D0 to evaluate stressor-land use relationships. The DO differential mass flux is the product of the stream discharge and the DO concentration, instead of the DO percent saturation. A correlation coefficient of 0.73 was calculated between agricultural land pr0portion and DO mass flux, which is stronger than the relationship between DO percent saturation and agriculture for the same dataset (R = 0.47). Thus, some land use impacts may be better related to DO using a mass flux analysis. CONCLUSIONS A variety of relationships between land uses and stream ecohydrology were evaluated using multivariate analysis. Our approach uses Landsat-based land use proportions, early morning DO measurements, and baseflow chemistry conditions as a variables in PCA to investigate relationships between land uses and ecosystem stressors. Multiple spatial scales (sourcesheds and riparian buffers), multiple representations of stream chemistry (concentrations and mass fluxes), and detailed land use data (e. g., urban divisible into road/non-road) were combined to examine correlations between land use 38 and stream ecosystem health indicators. This approach yielded correlations and principal components that possess distinct land use signatures. Across our Michigan study sites, streams in agricultural regions were generally associated with higher TN and conductivity, and lower DO. Although TP correlates with agricultural land proportions, the land use within the riparian buffer appears to more strongly influence the P levels in these streams. Agricultural regions also generally had the lowest DO mass fluxes. The strongest correlations in this study were between Cl' and urban areas. Our analysis indicates road salt and septic systems are the dominant sources of Cl', and isolated samples appear to have been influenced by oil/gas field-derived brines. Mass fluxes of TP and SRP are elevated in urban areas, possibly due to increased runoff. Stream temperature is frequently higher in urban areas, but the presence of a healthy vegetated buffer can alleviate the situation. In fact, the presence of forested land immediately surrounding a stream appears to moderate the water temperature more than any other variable examined in this study. Forests serve another beneficial role as well: in our four watersheds, TN levels are rarely elevated where forests are abundant. Management decisions ofien concentrate on direct runoff to streams as the controlling factor in stream water quality. In Michigan and other humid regions, this neglects the fact that groundwater provides the overwhelming majority of streamflow during the summer. As a result, groundwater plays a critical role in determining levels of stream ecosystem stressors during the baseflow season. 39 CHAPTER 2: SIMULATING THE ECOHYDROLOGY OF A SMALL MICHIGAN WATERSHED WITH MIXED LAND USE 40 INTRODUCTION Land uses can have significant influences on the health of stream ecosystems. Nearly 40% of rivers in the United States are impaired, mainly due to agricultural practices [USEPA, 2000]. In Michigan, elevated nutrient levels have been correlated with the proportion of agricultural and urban land within the source areas of samples [Johnson et al., 1997; Wayland etal., 2003; Welty et al., submitted]. Nitrate is one of the most common groundwater contaminants, detected in over 70% of 2,000 wells with varying land use in the US [Nolan and Stoner, 2000]. The highest median concentrations were in shallow groundwater beneath agricultural areas. However in some areas, urban land uses rival agricultural land uses as a source of nitrate to groundwater due to the density of possible sources in a small region [Wakida and Lerner, 2005]. In Michigan, understanding the transport of nitrate from different land uses requires sound knowledge of groundwater flow patterns, because of the dominance of groundwater in Michigan stream flow. For example, groundwater contributes over 90% of the annual discharge of the Manistee River in the northwestern portion of Michigan’s Southern Peninsula [Hendrickson and Doonan, 1972]. The Muskegon River watershed, where our study site is located, also has groundwater-dominated streamflow [0 ’Neal, 1997]. The detrimental effects of high nitrate concentrations can be partly mitigated by managing a watershed using knowledge gained from simulations of flow and solute transport. Groundwater models can simulate subsurface transport of solutes from various land use sources, providing insight to the pathways and fate of harmful substances. For example, Molenat and Gascuel-Odoux [2002] modeled groundwater flow and nitrate transport in France, demonstrating the impact of different land use strategies on the spatial distribution of groundwater nitrate concentrations. Puckett and Cowdeljy [2002] 41 combined groundwater modeling and age dating to link land use and groundwater nitrate in Minnesota. The influence of fertilizer application rates on groundwater nitrate concentrations in France was demonstrated by Conan et al. [2003]. In Michigan, Boutt et al. [2001] and Wayland et al. [2002] documented a link between land use and water quality using groundwater flow and solute transport models. Such models have even more utility when they examine the effect of different management practices [Almasri and Kaluarachchi, 2005]. Surface water quality models have also been used to examine the complex ecohydrological relationships between land uses and nutrients, oxygen, and temperature in streams. However, stream water quality models are generally used to examine the impacts of point sources rather than non-point sources, which are best examined with groundwater models. For example, Chaudhury et al. [1998] demonstrated the potential for low dissolved oxygen beneath a wastewater treatment plant due to nitrification using the QUALZE. Such surface water quality models have the means to simulate non-point sources, but this ability is rarely used, as it requires the existence of a calibrated groundwater model for the region contributing water to the stream. Despite both the evolution of hydrologic models and the recognition of the significance of non-point nitrate sources to streams, researchers rarely simulate flow and reactive transport in coupled surface and ground water systems. van Lanen and Dijksma [1999], modeled two-dimensional groundwater transport of nitrate to a groundwater-fed stream, but did not simulate in-stream reactions. A robust simulation of groundwater and surface water flow, transport, and reactions must quantify spatially variable groundwater and nutrient fluxes to streams. These non-point fluxes are difficult to estimate without a 42 regional groundwater model, which rarely exists for a region. Here, we describe an ecohydrological model developed by coupling widely used codes for groundwater flow, solute transport, and stream water quality. This coupled model was used to explore process-based linkages between land uses and stream ecosystem health. METHODS Study Region: Cedar Creek Watershed Cedar Creek, in southwestern Michigan (Figure 8), is well-suited to exploring relationships between land use and water quality. It flows through the lower half of the Muskegon River watershed (7,052 kmz), which has reduced biodiversity due to land use [0 ’Neal, 1997]. Land use changes in this region are projected to further increase runoff volumes, nutrient loads, and heavy metal concentrations over the next 35 years [Tang et al., 2005]. Small watersheds, such as Cedar Creek (124 kmz), can exercise considerable influence over nutrient exports to large streams and lakes [Peterson et al., 2001]. The spatial distribution of land uses within the Cedar Creek watershed facilitates simulation of land use-derived solutes such as nitrate since the upstream portion of this watershed is dominated by agriculture, including five concentrated animal feeding operations (CAFOs), while the downstream portions are predominantly forested. Overall, 60% of the watershed is forested/openlands/wetlands, while 36% is agricultural and the remaining 4% is urban. Medium and coarse-textured glacial tills drape the northern watershed and transition to glacial outwash and lacustrine sand and gravel in the central and southern watershed. 43 ® Transducer . _. J [:] CAFO ,4," “ g .3“. = = Streams and lakes {if} I. L14 ' . Residential wells 11/ I. - Urban. roads - Agriculture . Forest. openland. wetlands ‘ 0 1.25 2.5 S .. _ 3’ Muskegon river surface watershed Muskegon River expanded model Figure 8: Map of land use [MDNR, 2001] and hydrography with hillshaded topography within the Cedar Creek groundwatershed. The portion of Cedar Creek displayed in bold is the tributary simulated in QUAL2K. The inset map shows the watershed location as a small black region in the southern Muskegon River watershed (gray outline) relative to the state of Michigan. The expanded surface watershed of the regional Muskegon River watershed is displayed in bold on the state map. 44 The groundwater source area of Cedar Creek was delineated using a multi-layer model of the region encompassing the Muskegon River watershed [Kendall et al., in preparation] (Figure 8). The groundwatershed was used for this study rather than the surface watershed because regional modeling of the Grand Traverse Bay watershed in Michigan by Boutt et al. [2001] indicated that surface watersheds can have significant differences from groundwater source areas. The regional Muskegon River groundwater model was developed by expanding boundaries to significant hydrologic features (i.e., the next large stream beyond the surface watershed to avoid this issue at regional scales) (Figure 8) [Kendall et al. , in preparation]. Automated parameter estimation routines were applied to this 18,800 km2 model to calculate optimal hydraulic conductivity values for aquifer sediments in geologic zones parameterized using a digital map created from Farrand and Bell [1982]. Climatological, GIS, and Ecohydrological Data Before constructing the groundwater flow and transport models and the stream ecohydrology model, we collected and analyzed the necessary GIS, climate, discharge, and water chemistry data. We established a GIS database for the Cedar Creek region with hydrography, geology, topography, and land use characteristics. These datasets were compiled from the Michigan Geographic Data Library, established by the Michigan Department of Environmental Quality. Land surface elevations, which were used for the groundwater model surface boundary and gradient calculations in the stream ecohydrology model, were defined based on the National Elevation Dataset 26 m digital elevation model (DEM). The DEM was also used to calculate a drainage network in the 45 watershed, which we used for delineating the location of Cedar Creek. The Michigan Department of Natural Resources Michigan Resource Information System (MiRIS) provided locations of lakes for groundwater modeling and confined animal feeding operations (CAF Os) for transport modeling. Hydrogeologic zones were parameterized according to a surficial geology coverage derived from Farrand and Bell [1982]. The geometry of the aquifer base was interpolated between measured bedrock elevations by kriging the elevations of this contact from water well logs that were deep enough to intersect this boundary. The distribution of forested and agricultural lands were determined from the Integrated Forest Monitoring Assessment and Prescription (IFMAP), which is a statewide digital land use map with 30 m resolution derived from 1997-2000 Landsat data [MDNR, 2001]. We used the IFMAP to determine the location and distribution of agricultural lands for transport modeling, as well as for estimating riparian shade in the stream ecohydrology model. We obtained nitrate concentrations (discussed in Water Chemistry section) from the Michigan Department of Community Health (MDCH) and Michigan Department of Environmental Quality (MDEQ) residential well databases. Groundwater nitrate concentrations and addresses of wells in the two counties of the watershed were geocoded with a 60% spelling match, which resulted in 161 observations spanning 1980- 2003. Observations were from 124 wells, and the average was used for sites with more than one observation. Wells with nitrate concentrations below the detection limit (either 0.5 or 0.2 mg/L) were assigned a concentration of 50% of the detection limit. Since we are simulating regional behavior with average inputs, we compared our simulated concentrations to a kriged field of measured nitrate concentrations. 46 GIS grids of leaf area index (LAI), the ratio of one-sided green leaf area to ground area (Myneni et al., 2002), were used in calculations of potential evapotranspiration (PET) (discussed in Modeling section). For this study, we chose to use the remotely- sensed eight-day LAI averages from NASA’s Moderate Resolution Imaging Spectroradiometer (MODIS) instrument. Climate data The primary precipitation data used in this study were derived from the Next- Generation Radar (NEXRAD) network. However, 6% of the hourly NEXRAD data were missing, and if the gap was larger than one day, these periods used point observations from the closest National Oceanic and Atmospheric Administration (N OAA) precipitation gage, in Muskegon, Michigan (30 kilometers from the watershed). NEXRAD gaps less than one day long were filled by linear interpolation between the gaps. The use of NEXRAD data minimizes the error associated with using precipitation data that was recorded at a point some distance from the study area. Also, NEXRAD data integrates a precipitation signal across the region of interest, which has advantages over point precipitation which could have significant bias in convective storms. We resampled 4 km resolution NEXRAD data to 50 m, and calculated the average precipitation in the grid cells covering the watershed. Although NEXRAD and gage precipitation both have errors that alter the precipitation recorded, monthly observed NEXRAD data closely correlate to point precipitation data for this portion of Michigan, with r'2 values of 0.75 and 0.92 for 2003 and 2004, respectively [Jayawickreme and Hyndman, submitted, 2005]. Hourly NEXRAD data were summed into daily totals to 47 match the daily time steps used in our MODF LOW model. We used a variety of additional NOAA data for our linked model, including snow depths, air temperatures, dew point temperatures, wind speeds, and percentages of cloud cover. For PET calculations (described below in the Modeling section), we used solar radiation, relative humidity, and air temperature data from the Ludington, Michigan station of the Michigan Automated Weather Network (MAWN), operated by the Michigan State University Extension, the Michigan Agricultural Experiment Station and the Michigan Department of Agriculture. Discharge data Two pressure transducers recorded stream stage at hourly to sub-hourly intervals [Wiley and Richards, unpublished data]. One transducer was installed in the northern, agricultural-dominated portion of the watershed, while the other was installed in forested land near the watershed outlet (See Figure 8). Stream discharge measurements [Wiley and Richards, unpublished data] were used to construct rating curves between stage and discharge. For the upstream transducer 20 measured streamflows were plotted against transducer-recorded stage to create a logarithmic stage-discharge relationship. For the downstream transducer, only five measured streamflows were available to create a rating curve. A flood event at this site caused a shift in base level, requiring us to use a different rating curve before and afier the flood. Four measured flows collected before the flood were used to create a logarithmic stage-discharge curve, which was shifted after the flood to correspond to the remaining measured streamflow. 48 Baseflow was estimated from the discharge records using the USGS standard sliding interval method of Pettyjohn and Henning [1979]. This technique is founded on the relationship between drainage area and surface runoff. The area draining to each pressure transducer (148 km2 and 31 km2 for the downstream and upstream site, respectively) was calculated using ArcGIS. Next, this drainage area is used to calculate the length of time for surface runoff to cease following precipitation. This characteristic time is halved and reduced by one to produce a baseflow window. Then, for each day, the technique looks to past and firture days by the number of days defined by the baseflow window, and the lowest discharge in that period is considered baseflow. Water chemistry Water chemistry was sampled using a synoptic strategy, where samples were collected at multiple locations over a one to three day period while streams were at baseflow conditions. Synoptic baseflow chemistry data is ideal for this study, as it provides a ‘snapshot’ of water quality at essentially one point in time [Grayson et al., 1997]. Although QUAL2Kw was only used to simulate the concentrations of a single synoptic event, we measured DO concentrations that closely matched the synoptic event on an additional sampling trip. The ecohydrological variables that were sampled included DO, water temperature, conductivity, total nitrogen (TN), NO3-N02 (N Ox), NH4, total phosphorus (TP), soluble reactive phosphorus (SRP), SiOz, Cl', and suspended chlorophyll-a. Our models used inputs of nitrogen, phosphorus, temperature, and DO; SiOz and Cl' were measured for use in a different study. 49 Dissolved oxygen, temperature, and conductivity were measured pre-dawn using Yellow Springs Instruments multi-parameter water quality probes that were calibrated for DO with a water-saturated air technique prior to sampling. Water chemistry samples were collected in acid-washed 75 mL polyethylene bottles and flash-frozen on dry ice in the field to prevent chemical transformations and biological degradation. Water samples were analyzed according to standard methods (APHA, 1998) to determine the concentrations of total nitrogen (TN), NO3-NOz (N Ox), NH3, total phosphorus (TP), soluble reactive phosphorus (SRP), SiOz, and Cl' [Panayotoffi unpublished]. Water was filtered prior to analyses for SRP, NH3, and Cl'. SRP was assayed with the ascorbic acid method. TP was assayed as SRP following persulfate digestion. TN was assayed as NOx by second derivative UV spectroscopy following persulfate digestion. Four assays were run on a Scalar Sans-Serif segrnented-flow autoanalyzer: NOx by automated cadmium reduction, NH3 by the automated phenate method, SiOz by the molybdosilicate method, and C1' by the automated ferricyanide method. 100 mL of stream water was filtered for chlorophyll assay through 0.7 micron Whatman GF/F filters in the field. The filter was packed in aluminum foil and frozen, also in the field. Chlorophyll-a was determined by extracting pigments for 24 hours in 90% ethanol, by sonicating the sample for 15 minutes during the first hour of extraction, by fluorometric assay with a Turner Designs TD700, and by correcting for phaeophytin with acidification. 50 Groundwater and Surface Water Models Our ecohydrological model was developed by linking commonly used codes for groundwater flow, groundwater transport, and stream water quality. We use this linked set of publicly available codes to examine the influence of land use on the ecohydrology of Cedar Creek. We simulated the three-dimensional transient groundwater flow in this region using MODFLOW-2000 [Harbaugh et al., 2000]. The flow field calculated by the groundwater flow model provided the primary input for our simulation of nitrate transport through groundwater to streams using the MT3DMS code [Zheng and Wang, 1999]. Finally, we simulated stream water quality using QUAL2Kw [Pelletier and Chapra, 2005] based on the simulated nutrient inputs from the groundwater transport code, assuming steady state stream flow that is well mixed laterally and vertically. Stream velocities were calculated using the Manning Equation. A comprehensive heat balance accounts for loadings from sources including the atmosphere, tributary inflow, and groundwater inflow. QUAL2Kw simulated the conditions of the synoptic sampling event that were generated by the transient groundwater flow and solute transport models. It is reasonable to use a steady-state stream model coupled to transient groundwater flow and transport models, because the transient models are simulating the conditions leading up to the synoptic sampling event. Then, we simulated the conditions of the day of the synoptic sampling, which a steady-state model is suitable for. 51 Modeling General modeling approach We used three modeling codes in this study: MODFLOW, MT3DMS, and QUALZKW. MODF LOW was used to simulate groundwater fluxes fi'om the water table to surface water bodies, and to estimate hydraulic conductivity and unsaturated zone delay (discussed below) over the period of 1988-2004. Precipitation for 2003-2004 were from NEXRAD, while years before this were assigned the precipitation from 2004, which appeared to be a representative average year. MATLAB scripts were used to create the MODFLOW input files and to optimize model parameters [Kendall et al, unpublished]. Model Optimization was based on observations from 2003 and 2004. Baseflow-extracted discharges from the stream pressure transducers were used as our flow observations, which were compared to the output from the Stream Gage Package. In the optimization routine, we weighted flows 100 times higher than heads to account for the lower relative confidence in head measurements. We used the optimization toolbox of MATLAB to calculate the best values of unsaturated zone delay and hydraulic conductivity. To ensure computational efficiency, the MATLAB function frnincon used a linear search gradient to minimize the objective function, which was the sum of head and flow absolute residuals, which were normalized by the number of observations of each type. Cell-by—cell fluxes generated by MODFLOW were used by MT3DMS to simulate transport of nitrate transport through groundwater to streams from 1962 to 2004. The fluxes of nitrate to Cedar Creek on August 16, 2004 (the date of our synoptic sampling) were exported from MT3D to QUAL2Kw, where in-stream reactions were simulated. 52 Groundwater flow model and parameter optimization We developed a novel approach to evaluate recharge rates based on a simple soil water balance model. Unfortunately, MODFLOW does not account for unsaturated zone processes unless the new VSF package [Ihoms et al., in review] is implemented, which is rather complex. When precipitation falls on the surface evapotranspiration reduces the water that infiltrates and recharges groundwater after some delay period. MODFLOW instantly applies groundwater recharge directly to the water table, which means that precipitation data cannot be directly applied as recharge, since the quantity would be too large (due to ET) and recharge would be applied too quickly (due to unsaturated zone flow). Our recharge model has four components: calculation of precipitation, PET, snow water equivalent (SWE), and recharge using a soil moisture mass balance (Figure 9). 53 Precipitation Yes( \No l Rompute PT, PE from PM? i. Snow model Yes 1/ \, No Yes ( \NO [WWW] Yes I \No . Yes .{ \No [SWE=po(Z-Zol(SUB) I 1 ISMO+PRE+SWE>PT?j Yes .{ \ No [SM = 5M0 + pas + SWE - PT - RCHOJ SM =0 l l Yes ( \NO Yes .{ \. No ECH=SM-FC-PE HRCH=SM~FC| Soil moisture model Figure 9: Conceptual model of recharge model. NEX: original NEXRAD precipitation, IC: infiltration capacity, PRE: precipitation that becomes groundwater recharge, PT: potential transpiration, PE: potential evaporation, PM: Pemnan-Monteith Equation, SWE: snow water equivalent, 20: previous day snow depth, Z: current day snow depth, p: current day snow density, p0: previous day day snow density, pf: fresh snow density, Q0: previous day stream discharge, Q: current day stream discharge, SUB: sublimation, SMO: previous day soil moisture, SM: current day soil moisture, FC: field capacity, RCH: current day recharge, RCHO: previous day recharge. 54 A simple infiltration capacity filter was applied to daily NEXRAD precipitation to separate precipitation that would become runoff rather than groundwater recharge. A daily precipitation total above 30 mm was considered to exceed the infiltration capacity of the soil in the watershed. When precipitation was above 30 mm, the first 30 mm of precipitation was passed to the soil moisture mass balance, but any excess precipitation was considered to become runoff. We used MATLAB to separately calculate daily potential transpiration (PT) and daily potential evaporation (PE), which combined represent potential evapotranspiration (PET). This separation is possible because the equation for PT and PE is the same except for a canopy term. We calculated PET using a version of the Penman-Monteith equation (Monteith, 1965) adapted by Chen et al. (2005) for dependence on leaf area index and stomatal conductance (Table 7). Potential transpiration was calculated by . eS—e ARni+pcp E. : rai 2, [3+7 1+5?— 1" ai (1) E-, is the water transpired by vegetation Rn, is the net radiation of the vegetation (MAWN) p is the density of moist air (Calculated via ideal gas law) cp is the specific heat of air at a constant pressure (1005 J/kg K) es is the saturated water vapor pressure [T etens, 193 0] e is the actual water vapor pressure (product of e, and hourly humidity/ 100) r,“ is the aerodynamic resistance to vapor transport, set to 15, based on Chen et al. [2005] 55 2., is the latent heat of vaporization of water [Harrison, 1963] A is the slope the saturated vapor pressure-temperature curve [T etens, 1930] y is the psychrometric constant [ant, 1952] re, is the canopy resistance to vapor transport. Radiation fluxes were provided by the MAWN station in Ludington, Michigan, which is the closest available radiation flux tower. We assumed that the meteorological conditions across the watershed were similar to those at the PET station, allowing us to use the PET estimate for our watershed. Canopy resistance to vapor transport was calculated by Ci (2) g si L i gs, is stomatal conductance, set to 0.0122 [Schulze et al., 1994] L, is the leaf area index, varied by season according to the input grids of the watershed Potential evaporation is calculated the same as PT, except rci is set to one: eS—e Aan-‘FPCPT ai E . W' 2, (A + r) ‘3’ 56 Parameter Value/range Solar radiation 0-2950 ly Air temperature -23-15 °C Humidity 26.9-101.3 % Specific heat of air at a constant pressure 1005 j/kg K Atmospheric pressure 100 kPa Aerodynamic resistance to vapor 15 s/m transport Density of moist air 1.18 kg/m3 LAI 1.5 - 6 Stomatal conductance 0.0122 m/s Table 7: Parameters for calculation of PET. 57 Since NEXRAD data is limited to liquid precipitation, a simple snowmelt model was used to calculate snow water equivalent (SWE) from NOAA snow depth data in Muskegon (Figure 9). Fresh snow was assumed to have a density of 0.1 g/mL, based on the commonly assumed 10:1 ratio of snow thickness to liquid water equivalent. This simplifying assumption [Roebber et al., 2003; Judson and Doesken, 2000] was necessary because no snow density data exists for the study area. The temperature and snow depth record were compared to the stream discharge at the two transducers. If snow depth decreased but discharge did not increase, compaction of the snowpack was assumed and calculated by the product of the initial snow density and the relative change in snow depth. This provided a new snow density for the day. If snow depth increased, the new snow density is now a depth weighted average of fresh snow and the older snow densities. When snow depth decreased and discharge increased, snow water equivalent was calculated by the product of the current snow density and decrease in snow depth. A simple soil mass balance model was developed to compute recharge from precipitation, PET, and SWE (Figure 9, Table 8). We chose to develop the most parsimonious model that reasonably describes unsaturated zone flow for the Cedar Creek system. More parameters could be easily added for more complex systems. In this case, initial soil moisture was set to zero, and incoming precipitation, SWE, and PT were used to estimate the next day’s soil moisture. For each day, PT was subtracted from the sum of soil moisture, SWE, and precipitation, and if PT exceeded the sum of soil moisture and precipitation, soil moisture was set to zero. Any recharge from the previous day was removed from this quantity to maintain a proper mass balance. When soil moisture exceeded the mean field capacity of a medium-coarse soil, 0.75 cm/m [USDA, 1998], 58 potential recharge was calculated. Finally, potential evaporation was subtracted from this potential recharge if precipitation occurred on the day of calculation. There is a significant delay between precipitation and groundwater recharge, which depends on the water table depth Therefore, we optimized a delay parameter which represents unsaturated flow. The unsaturated delay is the time for a pulse of water to move one meter vertically through the unsaturated zone. A map of depth to water was used in the optimization routine for calculation of the delay. Mathematically, the delay is the slope of the line between the water table depth and infiltration time, assuming the delay is a linear function of the unsaturated zone thickness. A term could be added to this function to account for the hydraulic properties. With incoming recharge calculated, delay and hydraulic] conductivity were optimized to observations of head and discharge with monthly stress periods and daily timesteps. We used head observations that were recorded upon installation of 37 residential water wells in the region. However, only 11 of these wells were installed during the period of the MODF LOW optimization, 1988-2004. So, in order to have an adequate number of head observations, all head observations were used in the years 2003 and 2004 on the day the observation was recorded to account for seasonal variability. This technique assumes that long-term changes in the groundwater levels in the watershed are small. So, a historical head level was repeated in 2003 and 2004 on same calendar day it first occurred. 59 Parameter Value Initial snow density 0.1 g/mL Sublimation of snow during melt 5% Field capacity 0.75 cm/m Table 8: Parameter values for recharge model. 60 Regional groundwater flow patterns were simulated in MODF LOW, using a one layer model with 12,371 active cells, each 100 m by 100 m. The lower model boundary was defined using the kriged bedrock surface from residential water wells, while the top of the model domain was developed by interpolating the surface elevations from the DEM to the model cells. The edges of the groundwater source area, determined from regional modeling as discussed earlier, were considered to be no-flow boundaries. The exchange of water between the stream network and the aquifer was simulated using the Streamflow-Routing Package (SFR) of MODFLOW [Prudic et al, 2004]. The SF R used stream stages interpolated between stream crossing contours from a digital USGS map, and Manning’s equation was used to calculate stream depth for each cell’s flow, which was the summation of upstream groundwater inputs. A sub-routine of the SFR package, the Stream Gaging Station package, was used to calculate stream discharge at discrete locations. Transport model MT3DMS (version 5) [Zheng, 2005] was used to simulate nitrate transport from different sources to Cedar Creek based on fluxes calculated in the MODF LOW model. We simulated the effect of 42 years of constant nitrate input, based on land use data from 1998 (MiRIS) and 1997-2000 (IF MAP) [MDNR, 2001], on the regional groundwater nitrate concentrations using monthly stress periods to incorporate varying groundwater recharge rates and daily transport timesteps. Four sources of nitrate were included as recharge concentrations: CAF Os, agricultural inputs, atmospheric deposition, and septic systems. CAFO and agricultural nitrate sources were manually calibrated to a field of 61 kriged observed values, which represent the average spatial trend in nitrate concentrations. Urban and atmospheric nitrate were not allowed to vary in the calibration routine, because these sources are well-defined. The goal of the transport model is to generate fluxes of nitrate for simulation in QUAL2Kw. Thus, it is not critical to precisely match the concentration observed in the wells, provided the stream model simulates the observed levels well. CAF Os produce very large quantities of animal waste, which is then generally stored in lagoons before application as fertilizer to the surrounding land [Wilde et al., 2000]. Although the waste contains little nitrate, it contains high amounts of ammonia, which bacteria can oxidize to nitrite and nitrate [Chapra, 1997]. Leaking lagoons and infiltration from land application are thus potential sources of nitrate in the Cedar Creek watershed. CAFO nitrate recharge concentrations can vary greatly. For instance, Krapac et al. [2002] observed groundwater nitrate concentrations from <1 to >20 mg/L near Illinois swine facilities, while Harter et al. [2002] recorded a mean concentration of 64 mg/L in shallow groundwater below a California dairy and Fraters et al. [1997] noted mean shallow groundwater concentrations around intense Netherlands dairy farms from 100 to >200 mg/L. Nitrate recharge concentrations from CAFOs were calibrated based on simulated/observed groundwater nitrate concentrations from the MDCI-I/MDEQ database. Recharge concentrations were assigned to locations of CAF Os according to their location in the MiRIS land use map. Agricultural nitrate, representing nitrogen-based fertilizer and manure, was allowed to vary from 25-75% of the calibrated CAFO nitrate concentration. The spatial distribution of agricultural lands was defined from the IFMAP land use map [MDNR, 62 2001]. Recharge concentrations fi'om agriculture were also calibrated based on simulated/observed groundwater nitrate concentrations from the MDCH/MDEQ database. Although atmospheric deposition of nitrogen is usually small, it can significantly influence surface water bodies, as demonstrated by Jaworski et al. [1997]. As a result, we included atmospheric nitrate from the National Atmospheric Deposition Program station in Wellston, Michigan (approximately 100 km from the study watershed) in our transport model. The 20 year average atmospheric deposition value of 0. 1 8 mg/L was added to the applied recharge nitrate concentrations in all land use types. Urban-derived nitrate can be more significant than agricultural nitrate due to the focus of possible sources in a small region [Wakida and Lerner, 2005]. In this rural watershed septic systems provide the most likely source of urban nitrate inputs. Therefore, we considered the influence of septic nitrate on the watershed using the locations of residential wells in the MDEQ database. We used a representative mass flux to simulate the influence of urban nitrate. The mass flux of nitrate from a septic system is the product of the flow and nitrate concentration of the system. To approximate the flow from a septic system in the area, we used the estimate of 100 gallons of water per day used by the average American [USEPA, 2004]. The census-designated communities of Holton and Twin Lakes in the watershed have an average household size of 2.78 people [USCB, 2002], so on average septic systems contribute 278 gallons/day of water. We used a value of 1 mg/L of nitrate for the nitrate concentration exiting a septic system. Literature values for groundwater nitrate concentrations due to septic systems range from 10 mg/L and above [Wilhelm et al., 1994; Gold et al., 1990] to < 1 mg/L [Arnade, 1999; 63 Whelan and Barrow, 1984]. We decided 1 mg/L is a realistic value, as nitrate levels exiting a septic system are quickly reduced by bacterial denitrification to near this level. So, the mass flux of nitrate from a septic system in the watershed was simulated as the product of 278 gallons/day and 1 mg/L. Septic system nitrate fluxes were added to mass fluxes from agriculture, CAFOs, and atmospheric deposition, which were calculated by dividing the assigned concentration by the recharge. Then, these mass fluxes were divided by the recharge to yield nitrate concentrations. Stream ecohydrology model Transient groundwater flow and solute transport were simulated, however the QUAL2Kw code can only simulate steady state stream transport and reactions. Therefore, QUAL2Kw was set up to simulate the conditions of August 16, 2004, when the nutrient data were collected during a part of our regional synoptic sampling. QUAL2KW required input of hydraulic, climatic, and water chemistry data. For the hydraulic representation of stream flow, Cedar Creek was split into 6 reaches, each with unique hydraulic characteristics. Manning’s Equation was used to calculate stream discharge in the reaches, using channel cross sections measured in the field, and Manning coefficients of 0.03. The model uses inputs for a range of regional climatic data, including air temperature, dew point temperature, wind speed, and cloud cover. The fraction of blocked solar radiation was estimated for each stream reach using data from field visits, aerial photographs, and from the IFMAP land use map. These three sources allowed us to characterize the riparian area surrounding the stream as either forested or agricultural. Riparian forested land was estimated to block 75% of incoming solar 64 radiation while riparian agricultural land was estimated to block 50%. For water chemistry and reactions, non-point fluxes from MT3DMS and one point source was simulated. Groundwater discharging to Cedar Creek was assigned a temperature of 10°C, based on our measurements of shallow groundwater temperatures in the nearby Grand Traverse Bay Watershed, Michigan. Our only available groundwater temperature data in the watershed was too shallow, and could not be used for this model as it showed clear indications of stream temperature fluctuations. One tributary, which enters Cedar Creek in the upstream portion of the watershed, approximately 2.2 km from the head, was simulated as a point source in QUAL2Kw based on streamflows calculated by the SFR package (Figure 8). Nutrients, DO, and temperature were all measured for this tributary and the values were assigned in QUAL2K. The recommended values for rate constants were used. RESULTS AND DISCUSSION Groundwater Flow Optimization Hydraulic conductivity and unsaturated zone delay were optimized to water levels in wells and streamflow observations in a 16 year MODFLOW simulation (1988-2004) (Table 9). A paucity of wells in till regions required tying till conductivity to outwash conductivity with a scaling factor of 60%, resulting in values of 5.01 and 8.35 m/d, respectively. This scaling factor was based on estimates by other researchers [Spansky, unpublished data, 2004]. 65 Parameter Optimal Value Delay 2.47 d/m Outwash hydraulic conductivity 8.35 m/d Till hydraulic conductivity 5.01 m/d Table 9: Optimized recharge and conductivity parameters. Till conductivity was tied to outwash conductivity with a scaling factor of 0.6. 66 The delay parameter for the Cedar Creek watershed was estimated to be 2.47 d/m. Since no transient groundwater level data were available for the Cedar Creek watershed, we evaluated this delay parameter using transient head data from 7 of our groundwater transducers in the Grand Traverse Bay Watershed (GTBW). The date of the midpoint between the initial water-table rise and the peak of water table elevation due to snowmelt in 2004 (the center of mass of the water table rise) was identified in the seven GTBW wells. The center of mass for each well was plotted on a graph of center of mass arrival date vs. depth to water at center of mass arrival. Assuming horizontal flow around the well is negligible, a linear regression of these points yields an estimate of the unsaturated zone delay. The calculated delay for the GTBW 3.93 d/m, which is in a reasonable range of the 2.47 d/m optimized for the Cedar Creek Watershed given the differences in watershed characteristics. The wells in the GTBW that were used for the delay calculation were dominated by glacial till, while our watershed is dominated by higher- conductivity outwash. If enough data were available, delay parameters could be estimated for each geologic unit. Our model provides reasonable estimates of heads and transient streamflow given its simplicity and ease of implementation. Static head observations compared well to simulated levels, considering the head observations were recorded when wells were installed from 1969-2000 while the model simulated a period from 1988-2004 (Figure 10). As expected, the model provides a better match to the head levels that were recorded during the model simulation period (1988-2004) than to those recorded prior to the start of the model (1969-1987). This is partly due to longer term trends in water levels as well as the assumption that the seasonal variations during years prior to our transient data are 67 the same as those recorded in 2004. The simulation provided a reasonable match to the baseflow estimated from the transducers at the two sites on Cedar Creek (Figures 11, 12). Measured and simulated baseflows agreed slightly better at the upstream transducer. The downstream simulated baseflow has several peaks that are absent in the observed baseflow, which is likely due to a few periods of time when the precipitation exceeded the capability of the infiltration capacity filter to separate runoff from infiltration. A portion of the mismatch between simulated and observed baseflows is also likely due to our use of solar radiation and humidity data from a station outside our watershed. 68 230 -- R2: 0 All values: 0.79 g . 22° * Pre-1988: 0.77 'l a .. S 1988-2004: 0.87 0’ co '3' o g 210 . O . a E {6. g, 200 » (I) .D o . _- _ ._ ___ _-_ O . . O Pre-1988 Observations 190 4 l l ‘0 : 0 1988-2004 Observations 180 C _ __ _ _ 180 190 200 210 220 230 Simulated Head. m Figure 10: Simulated vs. observed heads for the flow optimization. Open circles are head observations from the period of the model simulation (1988—2004), closed circles from prior to the start of the model (1969-1987). 69 80000 4— Actual Streamflow l l l 0 Extracted Baseflow .- 1 60000 L-ijgmulated Baseflow ] 2 M E 40000 1 20000 - i i [A l ’I it A l, ‘ -l - ‘- W 11 \ “YR” 0 - 1/1/2003 1/1/2004 Figure 11: Simulated baseflow (dotted line), monthly extracted baseflow (points) and actual streamflows [Wiley and Richards, unpublished data] (line) for the upstream transducer on Cedar Creek. The optimization routine used the monthly extracted baseflow and the simulated baseflow for calculation of the objective function. Locations of this transducer is shown in Figure 1. 70 Q, m3ld 200000 i ii ii— Actual Streamflow 150000 g 1 il l . EXtraCted BaseflOW l. d 1::- Simulated Baseflow 100000~«/l”91 T1’r .11 l I :1, . N & .IHJ . {1’ ‘1 1‘. irifil {‘1‘ ti rilii 1 i ii i inky.“ f‘ (SN. /\ l- .44....1 “11.3.7.3“; ‘ 45.4 11. . 50000 0 - 1/1/2003 1/1/2004 Figure 12: Simulated baseflow (dotted line), monthly extracted baseflow (points) and actual streamflows [Wiley and Richards, unpublished data] (line) for the downstream transducer on Cedar Creek. Locations of this transducer is shown in Figure 1. 71 Nitrate Transport Model MT3DMS was designed to generate nitrate fluxes to Cedar Creek for use in QUAL2Kw. The optimized recharge and conductivity values were used to generate a flux field in MODFLOW that was used by MT3DMS for a 42 year simulation of nitrate transport. Recharge concentrations of nitrate from agriculture, CAFOs, and septic systems were calibrated to best match the observed nitrate levels in wells across the region. The calibrated values of nitrate source concentrations indicate CAFOs have the highest input concentrations (Table 10). The majority of simulated nitrate concentrations are within 2 mg/L from the kriged nitrate observations in the residential wells (Figure 13). We are comparing simulated nitrate for 2004 to a kriged field of nitrate observations from 1980-2000, so a low r2 value (0.22) is expected. Furthermore, the purpose of the transport model is to provide nitrate fluxes to the stream model, so simulated levels in the stream model are more critical. As expected, the agricultural regions have elevated nitrate levels relative to forested regions, with moderate levels in urban areas (Figure 14). For simplicity, the nitrate sources were considered to be active for the entire simulation. In reality, there are temporal variations in the flux of nitrate from agricultural lands, CAFOs, and septic systems. In the past, the regional population was lower, and fewer septic systems would have been in use. Thus, our septic nitrate concentration is likely slightly high for the early period of the model. Furthermore, researchers have observed significant temporal variations in septic nitrate inputs [Arnade, 1999], which we did not consider due to lack of data. These variations in septic nitrate do not exert much influence over the spatial distribution of groundwater nitrate, though, as they are contributing a relatively small mass of nitrate to the system. CAFOs, however, have the 72 highest nitrate source concentrations in the model, and thus nitrate fluxes are more important to simulate during the correct period. We based our CAF O locations on the 1998 MiRIS land use database, and inspection of aerial photos from 1979 also reveals the presence of the CAFOs used in the model. The presence of CAFOs in 1979 supports our use of them at the simulation start, in 1962. However, we have not accounted for any new CAF Os that may be operating in the watershed since 1998. 73 Nitrate Source Recharge Concentration, mg/L Atmospheric 0. 1 8 Agriculture 20 CAFO 40 - Septic 1.0 Table 10: Values for nitrate recharge concentrations used in MT3DMS. 74 Observed N0 , mg L'1 12 10 L- 1. ” l l 6 Simulated N03, mg L'1 10 Figure 13: Simulated vs. kriged observed groundwater nitrate concentrations (r2 = 0.22) in residential wells. Dashed envelope is 3: 2 mg L”. Observed concentrations are from MDCH/MDEQ database. 75 >10 mglL 0 mg/L Figure 14: Simulated groundwater nitrate concentrations by MT3DMS in the Cedar Creek watershed afier 42 years of land use practices. 76 Stream Ecohydrology Model Based on inputs from our groundwater flow and transport modeling, QUAL2Kw simulated streamflow, stream water temperature, DO, and nitrate concentrations on August 16, 2004. An additional sampling survey indicates the conditions of August 16 were representative of baseflow (Figure 15). While QUAL2Kw has the capability to model many other parameters, we only had data to support a robust simulation of streamflow, stream water temperature, DO, and nitrate. The QUAL2Kw simulation was based on flow estimates from groundwater to the Cedar Creek stream system from the MODF LOW simulation, which overestimates stream flow in the lower reaches of Cedar Creek (Figures 12 and 16). Four measured streamflows compare well to simulated streamflow. Two of these streamflows were used to calculate the objective function in the groundwater flow model optimization. A better match between these simulated and observed flows on this day could have be obtained if more than two of these measurements had been used in the optimization and if the importance of these measurements were weighted more heavily than those on other days, however this did not seem justified. QUAL2Kw simulated stream temperature to within 1 C° of observed levels (Figure 17). Water temperatures are highest in the agricultural headwaters of Cedar Creek, and decrease over 2 C° by the furthest downstream measurement point. We used the model to evaluate the influence of shading and incoming groundwater temperature on the longitudinal temperatures in Cedar Creek. This sensitivity analysis indicates the temperature decrease is caused by a shift from little to no riparian shade in the headwaters to high shade levels downstream. The gradient of the stream temperature decrease is too sharp to be caused by a change in groundwater temperature. 77 12 10 - X :1 8 X i x X x U) E, 6,- M“... g 4 1 +8/16/2005 1 2 x8/23/2005 ‘ 0 ~ . — VP CW __;_--.. O 5 10 15 20 25 30 Distance Downstream, km Figure 15: Measured DO concentrations in Cedar Creek on two separate sampling trips. QUAL2Kw simulated the conditions of 8/16/2005, as it was the trip water chemistry was collected. 78 1.25 1 laminar-.2555 * 7 “ _, — I] 1 00 “i I Observed streamflow (not used in objective function) i ,5 “A Obseryejlflstreamflow (usedjn objective function) - _ l «‘5 E. 0.75 4- E r 5 0.50 -- 9 a 0.25 -- 0 5 10 15 20 Distance Downstream (km) Figure 16: QUAL2Kw simulated and observed streamflow for the 8/16/2004 synoptic sampling. 79 14 .__._.---. --__-..-..,-_. _L OD .3 N 1 Water Temperature (°C) a : 9 -_ 3 —Simulated Water Temperature 1 3 I Observed Water Temperature__] 8 l 1 F _,___.____ l 0 5 10 15 20 Distance Downstream (km) Figure 17: Simulated and observed water temperatures for the 8/16/2004 synoptic sampling. Error bars represent the accuracy of the YSI temperature measurement. An incoming groundwater temperature of 10°C. 80 The simulated stream nitrate concentrations are similar to those we measured in mid-August 2004 (Figure 18). Non-point source nitrate concentrations were input to each of the six reaches based on the MT3D concentrations and the MODFLOW fluxes, and point-source nitrate was simulated from the one upstream tributary (Figure 8). Both observed and simulated nitrate reach their maximum approximately 2.5 km from the beginning of Cedar Creek (kilometer O). The source of this nitrate at this point is a mixture of agricultural and CAFO sources along with septic nitrate from the nearby community of Holton. Nitrate concentrations gradually decrease downstream, and are at approximately half of the peak value by 10 km downstream. The processes responsible for the reduction in nitrate could be dilution via incoming groundwater or a lower mass flux of nitrate entering the stream. Our analysis indicates that a combination of the two processes is the cause. The mass flux of nitrate to the stream is highest from 4.5 to 7 km, and is over 50% lower from kilometer 7 to the furthest downstream point. Dilution operates on top of this process, as stream discharge increases, serving to further reduce nitrate concentrations. 81 I l l l l l l 3 I Observed Nitrate | ' l l 1 3 . Q 1500 __ —Simulated Nitrate a ____ 3 .93 £5 1000 - z + £3 E a: 500 1* I I 0 l l i l 0 5 10 15 20 Distance Downstream (km) Figure 18: Simulated and observed nitrate concentrations for the 8/16/2004 synoptic sampling. 82 Nitrate is relatively simple to model in a stream since the longitudinal concentration is largely a firnction of the incoming groundwater and tributary nitrate sources, but D0 is more difficult to simulate as it is affected by in-stream reaeration, photosynthesis and groundwater concentrations. QUAL2Kw generally underpredicted DO concentrations, but simulates the overall DO concentration pattern (Figure 19). Concentrations of D0 are very low in the headwaters but increase to near saturation by 2.5 km downstream. The DO concentration decreases by approximately 1 mg/L over the next 7 km. After this point DO levels increase, likely due to an increase in stream slope, which produces higher velocities, encouraging reaeration along the reach. The influence of reaeration was considered by prescribing an extremely low reaeration rate, 0.001 d'l. As Figure 20 demonstrates, if atmospheric reaeration of oxygen is not accounted for in a stream ecohydrology model the simulated oxygen levels will be too low. 83 12 ._.___-__--___- 160 1 Simulated Dissolved Oxygen : I Observed Dissolved Oxygen 3 - - - Simulated Dissolved Oxygen Saturation l on O Chlorophyll (ug/L) l.__ _. __._..___. Dissolved Oxygen (mg/L) oo 6 _ - 2 Observed Chlorophyll ________ 40 4 ' l ' 1' . l l T 0 0 5 10 15 20 Distance Downstream (km) Figure 19: Simulated and observed dissolved oxygen concentrations, and observed chlor0phyll a concentrations for the 8/16/2004 synoptic sampling. Simulated dissolved oxygen saturation is dashed, and increases in the upstream portion of Cedar Creek where the temperature is dropping, as colder water can hold more dissolved oxygen. DO error bars represent the accuracy of the DO probe. 84 12 DO, mglL O) V \- 3 Reaeration 3— — No Reaeration ‘ ml Observed DO ———_ _——-_——_———_ 1O 15 20 Distance Downstream, km Figure 20: DO levels are much closer to observed values if reaeration is included in QUAL2Kw. 85 Although hydrologic studies commonly ignore biotic interactions with the hydrosphere we simulated algal concentrations in QUAL2Kw (Figure 21). Areal algae concentrations are highest in the upstream agricultural headwaters and decrease rapidly once nitrate inputs are reduced. The algae concentration pattern looks similar to another proxy of stream biota, suspended chlorophyll-a (Figure 21). We do not have algae data for our watershed, however suspended chlorophyll-a is a closely related quantity. Again, concentrations of chlorophyll-a are highest in the headwaters and decrease rapidly. The spike in chlorophyll-a in the headwaters of Cedar Creek is likely caused by the increased nitrate concentrations there, a condition also reported by Stevenson et al. [In press]. Additionally, suspended chlorophyll-a displays a concentration pattern opposite that of dissolved oxygen, with high concentrations upstream and decreasing downstream (Figure 21). In multivariate analysis of land use and ecohydrological variables, Welty et al. [Submitted] and Stevenson and White [1995] identified a correlation between increased chlorophyll-a and low DO. Our model is simulating the process responsible for the correlation noted by the authors. 86 l O.) C a r Chlorophyll, ug/L Bottom Algae mg/m2 N o _L o 1 1 0 '1 ' 1‘ .l 1 To 0 5 1O 15 20 Distance Downstream (km) Figure 21: Simulated algae areal concentrations and observed chlorophyll a concentration for the 8/16/2004 synoptic sampling. 87 Relative influence of land use practices MT3DMS and QUAL2Kw can be used to examine the influence of agricultural, urban, and CAFO nitrate sources. A MT3DMS simulation was run for four different land use scenarios: (1) CAFO, agriculture, urban, and atmospheric set to the optimal values; (2) CAFO nitrate was set to zero; (3) agricultural nitrate set to zero; and (4) urban nitrate set to zero. For the latter three cases the active nitrate sources were set to the optimal values. The simulated nitrate concentration in Cedar Creek corresponds best to the observed nitrate concentrations for Case 1 (Figure 22). When CAFO nitrate is turned off (Case 2), the peak nitrate concentration in the upstream portion of Cedar Creek is not reached, but the nitrate pattern downstream is similar to the observed values (Figure 22). Altemately, when agricultural nitrate is turned off, the simulated nitrate is too low for the entire stream (Figure 22). Finally, urban nitrate is not simulated, which causes the nitrate concentrations in the stream to be slightly lower than the observed values throughout the stream (Figure 22). This comparison has several implications. First, while CAFO nitrate is likely responsible for the peak nitrate concentration in Cedar Creek, the overall concentration pattern is controlled by agricultural nitrate inputs. Urban nitrate sources are present throughout the watershed, and thus adjusting the urban nitrate concentration can vertically shifi the simulated nitrate curve. Due to the low magnitude of the source concentration (0.18 mg/L), the relative influence of atmospheric nitrate was not simulated. 88 1 800 3 i 662E36— 3 1500 I I 3 —Optimized 3 412004 3 ""NOCAFO 3 T33: 1 —— No Agriculture 3 «5‘ 9001 3 --~-——NoUrban 3 o ,' _ ,-” .___.___.. - . .. z 600 1 x \ I ‘44- emu...“ 300 ~ ‘ ~ ~~~~~~~~ I 0 . 4 0 5 10 15 20 Distance Downstream, km Figure 22: The influence of different land use practices can be examined by simulating different land use scenarios in MT3D and running QUAL2Kw with the results. 89 CONCLUSIONS The relationship between land use and stream ecohydrology was simulated using process-based groundwater flow, nitrate transport, and stream ecohydrology models of Cedar Creek. This technique yielded several significant conclusions. First, the model results demonstrate the importance of groundwater flow to low order streams in Michigan. Although runoff during wet periods is responsible for some nutrient transport, groundwater nutrient fluxes explain the main spatial concentration patterns during baseflow in this watershed. Our approach has generated valuable nitrate flux estimates which were not previously available for CAFOs, agriculture, or septic systems in our watershed. CAFOs in the watershed contribute an average recharge concentration of 40 mg/L of nitrate, while agricultural lands appear to contribute roughly half of this level, and 1 mg/L is added per septic system. This provides input for future nutrient transport models in the region, which could eventually be used for management of land use practices. There are a number of assumptions and simplifications in our approach that are important to consider for application to other regions. Our model simulated the conditions of a single synoptic event, so additional simulations would strengthen the results. However, our synoptic event appears representative of this baseflow period with a survey collected roughly one week later. The use of PET data calculated from a region outside the watershed adds uncertainty to the results. In general, our approach provides realistic estimates of recharge, however there are periods of larger discrepancy between simulated and observed flows partly due to the difference in meteorological conditions between the MAWN site and Cedar Creek. Despite these potential drawbacks, our 90 process-driven approach is easily transferable to other watersheds. Climate, land use, elevation, hydrogeologic and geochemical data are widely available for most regions. Our model simulates ecohydrological processes including PET, saturated zone flow, nitrogen cycle and stream biota reactions, rather than relying on empirical relationships that may be unique to a particular watershed. This approach extends the traditional use of hydrologic models to represent ecosystem processes that exert considerable influence over variables such as DO. It also allows for the evaluation of the relative contribution of CAP Os, agriculture, and septic systems to the nitrate mass balance of the watershed. Through comprehensive ecohydrological watershed models we can examine the influence of land use practices on the health of aquatic ecosystem. 91 REFERENCES CITED Allan, J. D., D. L. Erickson, and J. Fay. 1997. The influence of catchment land use on stream integrity across multiple spatial scales. Freshwater Biology 37: 149-161. Almasri, M. N., and J. J. Kaluarachchi. 2005. Multi-criteria decision analysis for the optimal management of nitrate contamination of aquifers. J oumal of Environmental Management 74:365-381. Amrhein, C., J. E. Strong, and P. A. Mosher. 1992. Effect of Deicing Salts on Metal and Organic-Matter Mobilization in Roadside Soils. Environmental Science & Technology 26:703-709. APHA. 1998. Standard methods for the evaluation of water and wastewater, 20th edition. American Public Health Association, Washington, DC. Armour, C. L. 1991. Guidance for evaluating and recommending temperature regimes to protect fish. Biological Report 90(22), US. Fish and Wildlife Service, Fort Collins. Amade, L. J. 1999. Seasonal correlation of well contamination and septic tank distance. Ground Water 37:920-923. Bejda, A. J ., B. A. Phelan, and A. L. Studholme. 1992. The Effect of Dissolved-Oxygen on the Growth of Young-of-the-Year Winter Flounder, Pseudopleuronectes- Americanus. Environmental Biology of Fishes 34:321-327. Boutt, D. F ., D. W. Hyndman, B. C. Pijanowski, and D. T. Long. 2001. Identifying potential land use-derived solute sources to stream baseflow using ground water models and GIS. Ground Water 39:24-34. Brown, C. E. 1998. Applied multivariate statistics in geohydrology and related sciences. Springer, Berlin. Brunt, D. 1952. Physical and dynamical meteorology, 2nd edition. Univ. Press, Cambridge. Burkart, M. R., D. E. James, and M. D. Tomer. 2004. Hydrologic and terrain variables to aid strategic location of riparian buffers. Journal of Soil and Water Conservation 59:21 6-223. Carpenter, S. R., N. F. Caraco, D. L. Correll, R. W. Howarth, A. N. Sharpley, and V. H. Smith. 1998. Nonpoint pollution of surface waters with phosphorus and nitrogen. Ecological Applications 8:559-568. 92 Chabot, D., and J. D. Dutil. 1999. Reduced grth of Atlantic cod in non-lethal hypoxic conditions. Journal of Fish Biology 55:472-491. Chapra, S. C. 1997. Surface Water-Quality Modeling. WCB/McGraw-Hill, Boston. Chaudhury, R. R., J. A. H. Sobrinho, R. M. Wright, and M. Sreenivas. 1998. Dissolved oxygen modeling of the Blackstone River (Northeastern United States). Water Research 32:2400-2412. Chen, J. M., X. Y. Chen, W. M. Ju, and X. Y. Geng. 2005. Distributed hydrological model for mapping evapotranspiration using remote sensing inputs. Journal of Hydrology 305: 15-39. Conan, C., F. Bouraoui, N. Turpin, G. de Marsily, and G. Bidoglio. 2003. Modeling flow and nitrate fate at catchment scale in Brittany (France). Journal of Environmental Quality 32:2026-2032. Coutant, C. 1976. Thermal effects on Fish Ecology. Pages 891-896 in J. R. Pffafflin and E. N. Ziegler, editors. Encyclopedia of Environmental Science and Engineering. Gordon and Breach, New York. de Azevedo, L. G. T., T. K. Gates, D. G. F ontane, J. W. Labadie, and R. L. Porto. 2000. Integration of water quantity and quality in strategic river basin planning. Journal of Water Resources Planning and Management-Asce 126285-97. DeFries, R., and N. K. Eshleman. 2004. Land-use change and hydrologic processes: a major focus for the future. Hydrological Processes 18:2183-2186. Diaz, R. J. 2001. Overview of hypoxia around the world. Journal of Environmental Quality 30:275-281. ESRI. 2003. ArcInfo Workstation. Environmental Research Systems Institute, Redlands, CA. Farrand, W. R., and D. L. Bell. 1982. Quaternary Geology of Southern Michigan. The University of Michigan, Ann Arbor, MI. Fraters, B., H. A. Vissenberg, L. J. M. Boumans, T. de Haan, and D. W. de Hoop. 1997. Results of the monitoring programme: quality of upper groundwater on farms in the sandy areas of the Netherlands 1992-1995. RIVM-rapport 714801014, RIVM, Bilthoven, the Netherlands. Gold, A. J ., W. R. Deragon, W. M. Sullivan, and J. L. Lemunyon. 1990. Nitrate-Nitrogen Losses to Groundwater from Rural and Suburban Land Uses. Journal of Soil and Water Conservation 45:305-310. 93 Grand Traverse County. 2002. Comprehensive Plan...A Vision for the Future. Williams and Works, Traverse City, MI. Grayson, R. B., C. J. Gippel, B. L. F inlayson, and B. T. Hart. 1997. Catchment-wide impacts on water quality: the use of 'snapshot' sampling during stable flow. Journal of Hydrology 199: 121-134. Hannah, D. M., P. J. Wood, and J. P. Sadler. 2004. Ecohydrology and hydroecology: A 'new paradigm'? Hydrological Processes 18:3439-3445. Harbaugh, A. W., E. R. Banta, M. C. Hill, and M. G. McDonald. 2000. MODFLOW- 2000, The US. Geological Survey modular ground-water model- user guide to modularization concepts and the ground water flow processes. Open-File Report 00-92, USGS, Reston, VA. Harrison, L. P. 1963. Fundamentals concepts and definitions relating to humidity. in A. Wexler, editor. Humidity and moisture. Reinhold Publishing Co., NY. Harter, T., H. Davis, M. C. Mathews, and R. D. Meyer. 2002. Shallow groundwater quality on dairy farms with irrigated forage crops. Journal of Contaminant Hydrology 55:287-315. Hendrickson, G. E., and C. J. Doonan. 1972. Hydrology and recreation on the cold-water resources of Michigan's Southern Peninsula. Water Information Series Report 3, US. Geological Survey/Michigan Geological Survey, Lansing, MI. Hill, A. R. 1996. Nitrate removal in stream riparian zones. Journal of Environmental Quality 25:743-755. Jacobson, R. B., S. R. Femmer, and R. A. McKenney. 2001. Land-Use Changes and the Physical Habitat of Streams: A Review with Emphasis on Studies within the US. Geological Survey Federal-State Cooperative Program. Circular 1175, USGS, Reston, VA. J aworski, N. A., R. W. Howarth, and L. I. Hetling. 1997. Atmospheric deposition of nitrogen oxides onto the landscape contributes to coastal eutrophication in the northeast United States. Environmental Science & Technology 31 : 1995-2004. J ayawickreme, D. H., and D. W. Hyndman. Submitted. Evaluating the influence of land cover on seasonal water budgets. Water Resources Research. J enkinson, D. S. 2001. The impact of humans on the nitrogen cycle, with focus on temperate arable agriculture. Plant and Soil 228:3-15. 94 Johnson, L. B., C. Richards, G. E. Host, and J. W. Arthur. 1997. Landscape influences on water chemistry in Midwestern stream ecosystems. Freshwater Biology 37 :193- 208. Johnson, L. B., and S. H. Gage. 1997. Landscape approaches to the analysis of aquatic ecosystems. Freshwater Biology 37:1 13-132. Judson, A., and N. Doesken. 2000. Density of freshly fallen snow in the Central Rocky Mountains. Bulletin of the American Meteorological Society 81 : 1 577-15 87. Krapac, I. G., W. S. Dey, W. R. Roy, C. A. Smyth, E. Storment, S. L. Sargent, and J. D. Steele. 2002. Impacts of swine manure pits on groundwater quality. Environmental Pollution 120:475—492. LeBlanc, R. T., R. D. Brown, and J. E. FitzGibbon. 1997. Modeling the effects of land use change on the water temperature in unregulated urban streams. Journal of Environmental Management 49:445-469. Lombardo, L. A., G. L. Grabow, D. E. Line, J. Spooner, and D. L. Osmond. 2001. 2001 Summary Report: Section 319 National Monitoring Program Projects. National Nonpoint Source Watershed Project Studies, NCSU Water Quality Group, Biological and Agricultural Engineering Department, North Carolina State University, Raleigh, NC. Malmqvist, B., and S. Rundle. 2002. Threats to the running water ecosystems of the world. Environmental Conservation 29: 1 34-1 53. Mathworks, T. 2002. MATLAB 6.5. Natick, MA. MDEQ. 1997. Guidebood of Best Management Practices for Michigan Watersheds. Michigan Department of Environmental Quality, Surface Water Quality Division, Lansing, MI. MDEQ. 1987. Bedrock Geology of Northern Michigan. Lansing, MI. MDNR. 2001. Integrated Forest Monitoring Assessment and Prescription. Michigan Department of Natural Resources, Lansing, MI. MDNR. 2002. Boardman River Natural River Plan. Michigan Department of Natural Resources, Lansing, MI. Mendiguchia, C., C. Moreno, M. D. Galindo-Riano, and M. Garcia-Vargas. 2004. Using chemometric tools to assess anthropogenic effects in river water - A case study: Guadalquivir River (Spain). Analytica Chimica Acta 515: 143-149. 95 Molenat, J ., and C. Gascuel-Odoux. 2002. Modelling flow and nitrate transport in groundwater for the prediction of water travel times and of consequences of land use evolution on water quality. Hydrological Processes 16:479-492. Monteith, J. L. 1965. Evaporation and the environment. Pages 205-234 in Proceedings of the Symposium on Experimental Biology 19. Myneni, R. B., S. Hoffman, Y. Knyazikhin, J. L. Privette, J. Glassy, Y. Tian, Y. Wang, X. Song, Y. Zhang, G. R. Smith, A. Lotsch, M. Friedl, J. T. Morisette, P. Votava, R. R. Nemani, and S. W. Running. 2002. Global products of vegetation leaf area and fraction absorbed PAR from year one of MODIS data. Remote Sensing of Environment 83:214-231. National Research Council. 1991. Highway Deicing: Comparing Salt and Calcium Magnesium Acetate. Special Report 235, Transportation Research Board, Washington DC. Nolan, B. T., and J. D. Stoner. 2000. Nutrients in groundwaters of the conterrninous United States 1992-1995. Environmental Science & Technology 34:1156-1165. O'Neal, R. P. 1997. Muskegon River Watershed Assessment. Fisheries Special Report 19, Michigan Department of Natural Resources Fisheries Division, Lansing, MI. Osborne, L. L., and D. A. Kovacic. 1993. Riparian Vegetated Buffer Strips in Water- Quality Restoration and Stream Management. Freshwater Biology 29:243-258. Pan, Y. D., A. Herlihy, P. Kaufrnann, J. Wigington, J. van Sickle, and T. Moser. 2004. Linkages among land-use, water quality, physical habitat conditions and lotic diatom assemblages: A multi-spatial scale assessment. Hydrobiologia 515:59-73. Pelletier, G., and S. C. Chapra. 2005. QUAL2Kw Theory and Documentation (version 5.1). Washington State Department of Ecology, Olympia, WA. Peterjohn, W. T., and D. L. Correll. 1984. Nutrient Dynamics in an Agricultural Watershed - Observations on the Role of a Riparian Forest. Ecology 6521466- 1475. Peterson, B. J ., W. M. Wollheim, P. J. Mulholland, J. R. Webster, J. L. Meyer, J. L. Tank, E. Marti, W. B. Bowden, H. M. Valett, A. E. Hershey, W. H. McDowell, W. K. Dodds, S. K. Hamilton, S. Gregory, and D. D. Morrall. 2001. Control of nitrogen export from watersheds by headwater streams. Science 292:86-90. Pettyjohn, W. A., and R. Henning. 1979. Preliminary estimate of ground-water recharge rates, related streamflow and water quality in Ohio. Project Completion Report Number 552, Ohio State University Water Resources Center. 96 Plante, S., D. Chabot, and J. D. Dutil. 1998. Hypoxia tolerance in Atlantic cod. Journal of Fish Biology 53:1342-1356. Porporato, A., and I. Rodriguez-Iturbe. 2002. Ecohydrology - a challenging multidisciplinary research perspective. Hydrological Sciences J oumal-J oumal Des Sciences Hydrologiques 47:811-821. Prudic, D. E., L. F. Konikow, and E. R. Banta. 2004. A new stream-flow routing (SFRl) package to simulate stream-aquifer interaction with MODFLOW-2000. Open-File Report 2004-1042, USGS, Reston, VA. Puckett, L. J ., and T. K. Cowdery. 2002. Transport and fate of nitrate in a glacial outwash aquifer in relation to ground water age, land use practices, and redox processes. Journal of Environmental Quality 31:782-796. Roebber, P. J ., S. L. Bruening, D. M. Schultz, and J. V. Cortinas. 2003. Improving snowfall forecasting by diagnosing snow density. Weather and Forecasting 18:264-287. Roth, N. E., J. D. Allan, and D. L. Erickson. 1996. Landscape influences on stream biotic integrity assessed at multiple spatial scales. Landscape Ecology 11:141-156. Rustem, W. R., W. E. Cooper, S. Harrington, and A. S. Armoudlian. 1992. Michigan's Environment and Relative Risk. Michigan Department of Natural Resources, Lansing, MI. Sabater, S., J. Armengol, E. Comas, F. Sabater, I. Urrizalqui, and I. Urrutia. 2000. Algal biomass in a disturbed Atlantic river: water quality relationships and environmental implications. Science of the Total Environment 263: 185—195. Sala, O. E., F. S. Chapin, J. J. Armesto, E. Berlow, J. Bloomfield, R. Dirzo, E. Huber- Sanwald, L. F. Huenneke, R. B. Jackson, A. Kinzig, R. Leemans, D. M. Lodge, H. A. Mooney, M. Oesterheld, N. L. Poff, M. T. Sykes, B. H. Walker, M. Walker, and D. H. Wall. 2000. Biodiversity - Global biodiversity scenarios for the year 2100. Science 287:1770-1774. Schelske, C. L., and E. F. Stoermer. 1971. Eutrophication, Silica Depletion, and Predicted Changes in Algal Quality in Lake Michigan. Science 173:423-424. Schulze, E. D., F. M. Kelliher, C. Komer, J. Lloyd, and R. Leuning. 1994. Relationships among Maximum Stomatal Conductance, Ecosystem Surface Conductance, Carbon Assimilation Rate, and Plant Nitrogen Nutrition - a Global Ecology Scaling Exercise. Annual Review of Ecology and Systematics 25:629-660. 97 Stevenson, R. J ., and K. D. White. 1995. A Comparison of Natural and Human Determinants of Phytoplankton Communities in the Kentucky River Basin, USA. Hydrobiologia 297 :201-216. Stevenson, R. J ., S. T. Rier, C. M. Riseng, R. E. Schultz, and M. J. Wiley. In press. Comparing effects of nutrients on algal biomass in streams in 2 regions with different disturbance regimes and with applications for developing nutrient criteria. Hydrobiologia. Suppnick, J. 1996. Informational total maximum daily load for dissolved oxygen in Sycamore Creek. MDEQ Surface Water Quality Division, Lansing. Tang, 2., B. A. Engel, B. C. Pijanowski, and K. J. Lim. 2005. Forecasting land use change and its environmental impact at a watershed scale. Journal of Environmental Management 76:35-45. Tetens, O. 1930. Uber einige meteorologische Begriffe. z. Geophys. 6:297-309. Thoms, R. B., R. L. Johnson, and R. W. Healy. In review. The US. Geological Survey Modular ground-water model--user guide to the variably saturated flow capability (VSF). Reston, VA. USCB. 2002. 2000 Census of population and housing, summary population and housing characteristics. PHC-1-24, US. Census Bureau, Washington, DC. USDA. 1998. Estimating soil moisture by feel and appearance. Program Aid Number 1619, United States Department of Agriculture Natural Resources Conservation Service, Washington DC. USEPA. 2000. National Water Quality Inventory 2000 Report. Report EPA-841-F-02- 003, United States Environmental Protection Agency Office of Water, Washington, DC. USEPA. 2004. Water Facts. Report 816-F-04-036, United States Environmental Protection Agency Office of Water, Washington DC van Lanen, H. A. J ., and R. Dijksma. 1999. Water flow and nitrate transport to a groundwater-fed stream in the Bel gian-Dutch chalk region. Hydrological Processes 13:295-307. Wakida, F. T., and D. N. Lerner. 2005. Non-agricultural sources of groundwater nitrate: a review and case study. Water Research 39:3-16. 98 Wayland, K. G., D. W. Hyndman, D. Boutt, B. C. Pijanowski, and D. T. Long. 2002. Modelling the impact of historical land uses on surface-water quality using groundwater flow and solute-transport models. Lakes Reserv Res Manage 7: 1 89- 199. Wayland, K. G., D. T. Long, D. W. Hyndman, B. C. Pijanowski, S. M. Woodhams, and S. K. Haack. 2003. Identifying relationships between baseflow geochemistry and land use with synoptic sampling and R-mode factor analysis. Journal of Environmental Quality 32: 1 80-190. Welty, N. R., D. W. Hyndman, and R. J. Stevenson. Submitted. Investigating linkages between land use and ecosystem stressors. Journal of Environmental Quality. Whelan, B. R., and N. J. Barrow. 1984. The Movement of Septic-Tank Effluent through Sandy Soils near Perth .1. Movement of Nitrogen. Australian Journal of Soil Research 22:283-292. Wilde, F. D., L. J. Britton, C. V. Miller, and D. W. Kolpin. 2000. Effects of animal feeding operations on water resources and the environment. Open-File Report 00- 204, USGS, Reston, VA. Wilhelm, S. R., S. L. Schiff, and W. D. Robertson. 1994. Chemical Fate and Transport in a Domestic Septic System - Unsaturated and Saturated Zone Geochemistry. Environmental Toxicology and Chemistry 13:193-203. Zheng, C., and P. P. Wang. 1999. MT3DMS, A modular three-dimensional multi-species transport model for simulation of advection, dispersion and chemical reactions of contaminants in groundwater systems; documentation and user's guide. Contract Report SERDP-99-1, US. Army Engineer Research and Development Center, Vicksburg, MS. Zheng, C. 2005. MT3DMS v5 Supplemental User's Guide. Technical Report to the US. Army Engineer Research and Development Center, Department of Geological Sciences, University of Alabama. 99 l33331333333133313333333331331