m x n m m 3 Eva. L. .2 1 r .1? 71mm ”Human. “fl!” 3mm -G’¢s;4w§z T! mun nun A‘IW‘MWWN'NF-HHW (J2? . LIBRARY ‘ Michigan State University This is to certify that the dissertation entitled Predicting the Impacts of Land Use and Climate on Regional- . Scale Hydrologic Fluxes presented by Anthony D Kendall has been accepted towards fulfillment of the requirements for the Ph.D. degree in Environmental Geosciences 0am) {/27 2%“ 4 Major Professor’s Signature 51/13/200 {1 Date MSU is an Affirmative Action/Equal Opportunity Employer 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 MAY 2 8 20115 ;3.T- " 1‘3 5/08 K‘lProi/AchiPres/CIRC/DateDue indd PREDICTING THE IMPACTS OF LAND USE AND CLIMATE ON REGIONAL-SCALE HYDROLOGIC FLUXES By Anthony D Kendall A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of DOCTOR OF PHILOSOPHY Environmental Geosciences 2009 ABSTRACT PREDICTING THE IMPACTS OF LAND USE AND CLIMATE ON REGIONAL-SCALE HYDROLOGIC FLUXES By Anthony D Kendall Climate and land use change are altering the terrestrial hydrologic cycle locally, regionally, and globally. Stream discharge and runoff from local to continental scales are increasing in much of the world, even as peak runoff timing shifts because of earlier snow melt. Regional precipitation patterns are shifting, changing the quantity and distribution of moisture. Land use and land cover changes due to both human and natural activities further alter hydrologic fluxes, changing the balance between runoff and evapotranspiration. Managing water resources in the face of these changes requires a predictive under- standing of their impacts at all scales. Most decisions about land use occur at local scales, smaller than individual watersheds. On the other hand, water withdrawal regulations are typically implemented at the statewide scale. These decisions then impact regional water budget and fluxes in large basins, often influencing the flows entering other states or nations. Little is known about the quantitative influence of land use and cover types on terrestrial hydrologic fluxes, including groundwater recharge and evapotranspiration. At the plot scale, the processes governing biosphere-hydrosphere interactions have been intensively studied. But the tools necessary to examine fluxes at larger scales are not yet fully developed. This dissertation presents the Integrated Landscape Hydrology Model (ILHM) de- SiE’Sned to simulate fine-resolution terrestrial hydrologic fluxes over large domains. It is a process-based, fully-distributed model that is pararneterized solely using physically- measurable parameters. ILHM can be used with little calibration to determine the water budgets of regional watersheds. An application to the Muskegon River Water- shed (MRW) in lower Michigan showed uncertainties in water budget estimates on the order of 10%. ILHM also provides hourly values of a host of hydrologic fluxes and storages, and here is used to understand the relatively influences of regional climate variability, land use, and soil texture on groundwater recharge and evapotranspiration (ET) in the MRW. Off-growing season lake-effect precipitation patterns strongly influence groundwater recharge, but not ET, resulting in nearly 50% more recharge near the lake shore than further inland. Soil textures play a very important role in governing both recharge and ET, with fine textures reducing recharge relative to sandy soils by nearly 40% and increasing upland ET by 30%. Variability between land use classes in sandy soils is almost 20% of upland ET and 45% of recharge. This work is dedicated to my wife Cheryl, whose tireless support included countless hours of assistance in the field and expertise in the lab, and to my daughter Lydia, whose life began after Chapter 3 was written and has since helped her Daddy make it all the way through the end of Chapter 5. iv ACKNOWLEDGMENTS I would first like to thank my advisor, David Hyndrnan. His unwavering support during the development of ILHM helped me slog through some very difficult problems. His understanding of the physical processes that drive terrestrial hydrologic fluxes proved far more valuable than any textbook. Most importantly, however, has been his constant encouragement of my education and career. Had it not been for his suggestion now seven years ago that I pursue a graduate degree in the Geosciences, this Dissertation would be very different indeed. Though I have already dedicated this work to my wife Cheryl, she deserves a place here as well. Her skills both in the field and the lab allowed me to accomplish far more in the early part of my degree than I would have otherwise. There are a number of graduate students, and some undergraduates, who have lended their talents as well. In particular, Dush Jayawickreme and Nick Welty pro- vided both field assistance and helped to develop early portions of ILHM code. Oth- ers, including Chris May, Matt Spansky, Jeff Eggleston, and Jason Bernstein lent numerous hours to collect data in the field. Finally, my funding has come from a number of sources, including grants from the National Science Foundation, Great Lakes Fisheries Trust, and MSU’S Center for Water Sciences. Additionally, MSU and the Department of Geological Sciences provided support through teaching assistantships, fellowships, and scholarships. TABLE OF CONTENTS LIST OF TABLES viii LIST OF FIGURES xi 1 Introduction 1 1.1 Integrated Hydrologic Modeling ..................... 2 1.2 Description of Chapters ......................... 5 2 Examining Watershed Processes Using Spectral Analysis Methods Including the Scaled-Windowed Fourier Transform 7 2.1 Abstract .................................. 7 2.2 Introduction ................................ 8 2.3 Datasets and Study Area ......................... 10 2.4 Methods .................................. 12 2.5 Results and Discussion .......................... 31 2.6 Conclusions ................................ 46 Evaluating Temporal and Spatial Variations in Recharge and Streamflow Using the Integrated Landscape Hydrology Model (ILHM) 49 3.1 Abstract .................................. 49 3.2 Introduction ................................ 50 3.3 Methods .................................. 53 3.4 Results ................................... 70 3.5 Discussion and Conclusions ....................... 79 ILHM Development, Validation, and Simulation of Regional-Scale Stream Discharge 84 4.1 Introduction ................................ 84 4.2 The Integrated Landscape Hydrology Model .............. 88 4.3 Model Domain and Input Data Preparation .............. 103 4.4 Uncalibrated Model Results and Discussion ............... 122 4.5 Conclusions ................................ 131 Influence of Land Use, Climate, and Soils on Recharge and Evapo- transpiration Across a Regional Watershed 133 5.1 Introduction ................................ 133 5.2 Methods and Study Site ......................... 136 vi 5.3 Model Results and Discussion ...................... 5.4 Conclusions . . APPENDICES oooooooooooooooooooooooooooooo A ILHM Model Development A.1 Evaporation and Transpiration ...................... A.2 Infiltration, Percolation, Throughflow, and Exfiltration ........ A.3 Runoff Routing .............................. B ILHM Regional Upscaling B.1 Surface Domain .............................. B.2 Deep Unsaturated Zone ......................... 83 Saturated Zone BIBLIOGRAPHY vii 170 170 170 177 180 182 182 195 196 197 2.1 2.2 3.1 3.2 4.1 4.2 4.3 4.4 4.5 5.1 LIST OF TABLES List of the nominal time periods used in each seasonal analysis, and of the years omitted from the average unit responses. Data were from 1990-2001 .................................. Slope break locations and 6 values for Figure 2.4. ........... List of adjustable parameters in the UEB Snowmelt Model and their assumed values .............................. List of calibrated parameters for the unsaturated and saturated groundwater models. ........................... Input climate and LAI data sources. Start and end dates are general- ized. Locations of gage data are shown in Figure 4.4. ......... Table of observations and validation data. Location abbreviations “MR” refers to “Muskegon River”. The “Area” column refers to the surface watershed area ........................... Wetland depths assigned based on NWI classifications. The depths are based on a sinusoidal water table fluctuation of ~2.5 m ......... Table of optimized constants for the NDVI ——> LAI mapping function 4.12, and coefficient of determination R2. Open water and barren are assumed to have LAI = 0 in ILHM .................... Total Nash-Sutcliffe efficiency and average annual flux error values for each of the 6 USGS gages. Efficiency values are calculated using log- transformed simulated and gaged values. Average annual flux error is Sim—obs divided by the catchment area of each gage. The final column is average annual error normalized by average annual precipitation for that catchment ............................... Input climate and LAI data sources with the start and end dates of each source. Locations of gage data are shown in Figure 5.1 ...... viii 19 34 64 71 105 106 109 117 130 140 5.2 5.3 5.4 5.5 5.6 5.7 5.8 5.9 Tables of multiple regression models for annual watershed-average A) Upland ET, and B) Upland DP. Columns include the root-mean square residual between the multiple regression model and ILHM, along with the correlation coefficient R2, and the model parameters. Note, the final five years of the simulation were excluded from this analysis be- cause DP seems to behave exceptionally during this period compared the rest of the simulation. ........................ Cross-correlation values between monthly climate indices ........ VVatershed-average annual mean, standard deviation, and coefficient of variation, cu, of hydrologic fluxes ..................... Watershed-averaged mean, standard deviation, and coefficient of vari- ation, 0v, of upland water budget components .............. Average annual upland ET and DP by soil texture in agricultural ar- eas between 75-150 km from the lake shore along the prevailing wind direction. ................................. Average annual upland ET and DP by land use in sandy soils 75-150 km from the lake shore along the prevailing wind direction ....... Nested linear multiple-regression models for upland ET. A) displays the model intercepts, coefficients, correlation coefficients, and root- mean square residuals (RMSR) between the ILHM results and the linear models, listed in B) are the ranges of values for each linear model parameter, C) shows the influence of each parameter relative to the mean field (the intercept value, abbreviated Intcpt.), calculated as the product of the range of each parameter and the model coefficients from part A. PCS: PNGS: TG S and TNGS are growing season and non- growing season precipitation and temperature, respectively. F C is soil field capacity, and all land use classes are abbreviated. ........ Nested linear multiple—regression models for upland DP. A) displays the model intercepts, coefficients, correlation coefficients, and root- mean square residuals between the ILHM results and the linear models, listed in B) are the ranges of values for each linear model parameter, C) shows the absolute value of the relative influence of each parameter, calculated as the product of the range of each parameter and the model coefficients from part A. Column headers as in Table 5.8 ........ ix 150 151 154 154 159 162 165 166 Al A2 A3 31 B2 B3 B4 B5 B6 B7 B8 B9 Parameters used in calculating PET ................... 172 Land-Use parameters ........................... 174 Soil properties .................... i ........... 178 Land use dependent parameters in ILHM ................ 186 Depression capacity, SDEpmax [m x 103]. Values are from Manfreda et al. (2005), assuming no soil-texture dependence, and taking values for clay soil textures. ........................... 187 Soil-texture dependent parameters. Texture name abbreviates S, L, Si, C, and M are sand, loam, silt, clay, and muck, respectively. If multiple sources are listed for a particular parameter, values obtained from multiple sources were averaged ................... 188 Table listing physical constants ..................... 189 Table listing albedo and emissivity parameters, and sources ...... 189 Table of domain-indepedent physical parameters and sources ..... 190 Table of domain-dependent physical parameters, values are explained in Chapter 4 ................................ 190 Lookup table of unsaturated zone wetting front velocities. ...... 195 Table of preliminarily-calibrated saturated hydraulic conductivities for deep sediments .............................. 196 2.1 2.2 2.3 2.4 LIST OF FIGURES Images in this dissertation are presented in color. Map of the Muskegon River and Grand vaerse Bay watersheds, with an inset map of Michigan showing their locations. Groundwater wells with transducers are marked with triangles, and the Evart gauge sub- basin is shaded in grey ........................... SWFT results for a test function, f (:r) a) Curves representing the successive superposition of the three sinusoids with different frequen- cies, and b) the SWFT scalogram. The shaded contours represent amplitudes between 0.5 and 1.0 as shown in the colorbar. The dashed line represents the boundary beyond which the SWFT scalogram is undefined (where half the window width is greater than the available number of points within the dataset) ................... a) Comparison of the heuristic snowmelt model to the UEB model, b) the difference between the two (heuristic-UEB) and the cumulative difference. ................................. Composite plot of stream discharge, well-head relaxation, and pre- cipitation spectra. Grey lines plot the raw FT spectra, solid lines are the binned-mean amplitudes, and dashed lines give the upper 95% confidence intervals. a) An overlay of the spectra in parts c—e and their linear fits, b) NEXRAD precipitation spectrum for 4/ 1 / 2004-11 / 30 / 2004, not plotted on part a, c) stream discharge spec- tra for 10/1/1997—9/30/2004, d) well-head relaxation for well B13 from 7/ 1/ 2003-7/ 1 / 2005, and e) watershed-available precipitation from 1 / 1 / 1990-1 / 1 / 2001. Note that, in part c the limited resolution makes it appear as if much of the stream discharge spectrum is above the confidence limit, however fewer than 5% of points exceed the limit in a spectrum with over 130,000 points. .................. xi 11 26 32 2.5 2.6 2.7 2.8 2.9 2.10 SWFT of stream discharge at the USGS gage in Evart, MI. a) Stream discharge time-series divided by the area of the watershed above this gage, b) Scale-averaged amplitudes, Ap, c) Filled contour scalogram of the amplitude (A) vs. both period (P) and time (T), along with the location of the amplitude-weighted mean period, Pt (white line), d) Time—averaged amplitude spectrum, At. ................ SWFT of watershed-available precipitation at Big Rapids, MI. a) watershed-available precipitation time-series, b) Scale-averaged ampli- tudes, Ap, c) Filled contour scalogram of the amplitudes (A) vs. both period (P) and time (T), along with the location of the amplitude- weighted mean period, Pt (white line), Note, white areas have ampli- tudes < log(amplitude):-3.5, d) T ime-averaged amplitude spectrum, At. ..................................... SW FT of stream discharge at Evart from 9/1/1991 to 8/ 31 / 2004. a) Stream discharge time-series, b) Filled contour scalogram of the am- plitude (A) vs. both period (P) and time (T). The year labels are centered at January 1. .......................... SWFT of groundwater-head fluctuations in Well BIO in the Grand Traverse Bay Watershed. a) Well-head fluctuation time-series, b) Filled contour scalogram of the amplitudes (A) vs. both period (P) and time (T) ...................................... Stream unit response hydrographs. a) Unit response functions obtained using data from 9/ 1 / 1999 to 8/31/2000, and b) Seasonal unit response functions using selected years (see Table 2.1) from 1990-1999 ...... Plots of convolved stream discharge from 09/01/1999 to 8 / 31 / 2000 us- ing the annually- and seasonally-derived unit hydrographs (initial mod- eled discharged matched to stream discharge) on top of the measured stream discharge (a), and on top of the measured stream discharge. b) Plot of the percentage difference between each model and the measured discharge (modeled-measured) / (measured) *100. ............ xii 38 40 41 43 44 45 3.1 3.2 3.3 3.4 3.5 3.6 3.7 On the left is a map for the Cedar Creek watershed (shaded) along with the groundwater contributing area to Cedar Creek (dashed outline). Also displayed on this map are locations of two stream gages, nine discharge measurement cross sections, and residential drinking water wells located within the watershed. On the right, a map of the lower peninsula of the state of Michigan shows the Cedar Creek watershed within the greater Muskegon River watershed. The boundary of a regional groundwater model of the Muskegon River is shown in bold on the state map. The three precipitation and climate gage locations are Hesperia (H), and Fremont (F). ................... Static GIS datasets that were integrated into the ILHM framework. From left to right, the datasets are land use/ cover from IFMAP, qua- ternary geology from (Farrand and Bell, 1982), SSURGO soil textures, and land surface elevation from the NED 26.5 m DEM. ....... Time series data inputs to ILHM for 2003 and 2004: A) LAI of forest and agricultural land covers; B) monthly rain and snow along with daily average temperature, a horizontal line at 0°C is included as a visual aid; C) weekly averaged windspeed and solar flux; D) inset of hourly windspeed and solar flux data from 6/21/04 through 6/ 24/04; E) weekly averaged relative humidity (RH) and temperature; F) inset of hourly RH. and temperature from 6/ 21 / 04 through 6/ 24/04. Conceptual model and simplified box diagram of the ILHM. ..... Plot of measured discharge versus hydraulic radius in streams in the greater Muskegon River Watershed and adjacent watersheds. The hy- draulic radius calculation assumes rectangular cross-section geometry. The power-law fit to the data produced a correlation coefficient, R2, of 0.77 .................................... Uncalibrated simulated vs. observed groundwater heads plotted on top of a 1:1 line. Observations are from 96 residential wells installed in the Cedar Creek Watershed between 1990 and 2004 (see locations in Figure 3.1) ................................ A & C) Upper and B & D) Lower Cedar Creek stream discharges with simulated values shown in white and gage values in black calculated based on stage discharge relationships. Manually measured discharge values are shown as black circles. .................... xiii 54 55 60 62 70 72 3.8 3.9 3.10 3.11 3.12 4.1 4.2 4.3 4.4 4.5 Uncalibrated simulated vs. observed flows for the Upper and Lower Cedar Creek gages plotted on top of a 1:1 line from approximately 20 manual discharge measurements taken at each of the upper and lower Cedar Creek gage locations. ....................... ILHM-simulated average annual recharge for the Cedar Creek water- shed. Annual total precipitation averaged approximately 83 cm during the period of simulation. ......................... ILHM-simulated average annual precipitation excess runoff for the Cedar Creek watershed. Annual total precipitation averaged approxi- mately 83 cm during the period of simulation. ............. A) Monthly average deep percolation (bars) in the Cedar Creek wa- tershed plotted with the ratio of percolation to total precipitation for each month (stars). B) Monthly difference in deep percolation between forested and agricultural land covers. .................. Average monthly evaporation (E), transpiration (T), and evapotran- spiration (E+T) plotted for both forested and agricultural landscape as stacked bars along with the difference (Forest—Agriculture) between the two land-use types plotted as circles. ................ Graphic depiction of the spatial discretization scheme used in ILHM. ILHM generalized process flow diagram, focusing on surface processes. Dashed outlines indicate optional items. ................ Detail of hydrologic fluxes calculated by ILHM for upland (left) and wetland (right) portions of cells ...................... Map of the Muskegon River Watershed and surrounding model bound- ary. In part A, major municipalities along the drainage network are highlighted. Part B shows the locations of input gages and observation data locations ................................ Static GIS inputs for ILHM expanded model boundary ........ xiv 75 76 77 78 80 90 94 104 108 4.6 4.7 4.8 4.9 4.10 4.11 4.12 Plot of annual model-averaged solar radiation, from 1980 through 2007. There are two primary periods, 1) only cloudiness observations are available for shortwave radiation modeling, 2) cloudiness observations along with limited solar radiation data are available ........... Plot of monthly model-averaged precipitation from daily interpolated NCDC precipitation gages, and NEXRAD radar from 1996 - 2008. The quality of the N EXRAD data can be differentiated into three periods, 1) all monthly observations are a poor fit to gage data, 2) NEXRAD provides a good match to months with liquid precipitation only, and 3) NEXRAD data and gage data match closely in almost all months. Plot of monthly model-averaged LAI from three sources: 1) 8 km GIMMS AVHRR data, converted from NDVI to LAI, 2) 1 km AVHRR data, uncorrected, converted from NDVI to LAI, and 3) MODIS LAI data ..................................... Simulated and observed (gage) hydrographs for three USGS stream gages, A) Muskegon River at Evart, B) Clam River at Vogel Center, and C) Bear Creek near Muskegon. ................... Plots of simulated and observed streamflows and groundwater head, observation data are described in Table 4.2. A) Observed daily dis- charge at the Evart, MI gage vs. unbiased residual (sim — obs) /obs, B) observed discharge from 105 points across the MRW during low-flow conditions vs. simulated baseflow , C) observed groundwater head vs. simulated head, and D) annual median and +/- 1 standard deviation of the head residuals sim — obs. All horizontal and diagonal lines are reference lines, not linear fits. ...................... Plot of water-year (Oct - Sep) Nash-Sutcliffe efficiency of the log- transformed simulated and gage streamflow values for all six USGS gages. “MR” refers to “Muskegon River”. ................ Plots of stream discharge error Sim—obs for all six USGS gages normal- ized by the area of each gage catchment. Plotted in A) are water-year error values, and B) are average monthly error values across the 27 water years in the simulation. ...................... XV 113 115 120 124 126 128 129 5.1 5.2 5.3 5.4 5.5 5.6 5.7 Map of the Muskegon River Watershed and surrounding model domain boundary, with an inset showing the location of the model domain within lower Michigan, USA. Filled symbols on the lower map are climate gage locations for data listed in Table 5.1 ............ Maps of select ILHM input datasets. Sub-figures A and B show distri- butions of major and minor land use classes, C illustrates the dominant soil textural classes, and D plots the average annual precipitation across the model domain. ............................ Monthly values of select climate and leaf area index (LAI) inputs. The heavy solid line is the median monthly value during the simulation pe- riod, which is surrounding by a shaded range indicating the 5% and 95% values. Plotted in A) monthly temperatures, B) monthly precip- itation, C) monthly average solar radiation, and D) monthly average LAI for four land use classes ........................ Maps of total cell mean annual modeled A)evapotranspiration (ET), B) mean annual modeled deep percolation (DP), C) coefficient of variation (ratio of standard deviation to the mean) of annual ET, D) coefficient of variation, 0., of annual DP. Note the color scales differ for each sub- figure. ................................... Plots of annual watershed-averaged fluxes for the MRW. On the left, watershed-average ET, storage, and runoff for the MRW are plotted, while on the right upland ET, overland flow, root-zone storage, and DP within upland landunits are shown. Parts A 81, B plot annual fluxes in cm. Those quantities, divided by annual precipitation, are plotted in parts C & D. Note, for storage, negative values indicate release of water from storages in the model. .................... Average monthly fluxes of A) watershed ET, runoff (overland + ground- water discharge), and storage, and B) upland ET, overland flow, root- zone storage, and DP within upland landunits. ............ Filled area plot of watershed-average monthly evaporation and tran- spiration components for the MRW. Shown are upland transpiration, wetland transpiration, soil evaporation, canopy evaporation, snow and ice sublimation, and open water evaporation ............... xvi 137 138 139 146 148 152 1 5 3 5.8 5.9 5.10 5.11 Filled value plots of average monthly A) precipitation, B) upland evap- otranspiration (ET), and C) upland deep percolation (DP) vs. distance from Lake Michigan along the prevailing wind direction. All values are averages of only cells within the model domain containing greater than 75% deciduous forests, with sand as the dominant texture. Values are averaged over the entire model period. ................. Plots of monthly average upland ET and DP across generalized soil texture classes as a percent of water-year (Oct - Sep) precipitation. Values are from cells within the model domain containing greater than 75% agricultural land use between 75-150 km from the lake shore along the prevailing wind direction. ...................... Plots of monthly average upland ET and DP across land use classes as a percent of water-year precipitation. Values are from cells within the model domain located between 75—150 km from the lake shore along the prevailing wind direction with either sands (A & B) or loams (C & D) as the dominant soil texture class both vertically and horizontally. Maps of linear regression modeled annual average A) Upland ET and B) Upland deep percolation (DP) for 1980 - 2007. Parts C) and D) Show mean root mean square residuals (RMSR) between the ILHM annual fluxes and the linear regression modeled values of ET and DP respectively. ................................ xvii 161 167 CHAPTER 1 Introduction Climate and land use change are altering the terrestrial water cycle worldwide (Labat et al., 2004). Effects of these changes vary regionally in response to heterogeneous forc- ings, landscape properties, and anthropogenic influences. Impacts of these changes are felt across a variety of scales, from hectares to continents. These impacts include changes to soil moisture, groundwater recharge, evaporation and transpiration, and stream discharge. Changes to the soil moisture regime may be the single most economically-disruptive impact of climate warming. Increased evapotranspiratve demand coupled with a longer growing season and shifting seasonality and quantity of precipitation will alter the productivity of non-irrigated agricultural. In many regions, including the mid- western US, growing season precipitation is projected to decline, even while annual totals increase (Christensen et al., 2007). Increased demand for irrigation will stress highly-allocated water resources. In much of the western US, net water demand greatly exceeds annual supplies ( Vo'ros- marty et al., 2000). Elsewhere, much of what is not allocated to human use provides stream discharge and surface water levels to ecosystems adapted to the status quo. Irrigation alters the water balance of a watershed by increasing evapotranspiration, and thereby reducing groundwater levels and stream discharge. Agricultural practices are not the only land use changes that disrupt terrestrial wa— ter cycle. Reforestation and deforestation, urbanization, draining and channelization of wetlands, and impoundment of streams all affect the balance between ET, runoff, and storage. Many of these land use changes occur simultaneously, particular within larger basins, though the character of the aggregate varies regionally. Our ability to quantitatively forecast the impacts of climate and land use change is limited due primarily to the complexity of the interactions among landscape pro- cesses, and their characteristic spatial and temporal scales. Precipitation, tempera- ture, solar radiation, and other climatic factors in the hydrologic cycle vary globally, regionally, and more locally due to orographic and lake effects. Discharge within a large river system is the aggregate of the response of each individual square meter within that basin, whose infiltration capacities, land cover, and relief vary consid- erably over hectare scales. Beneath the surface, groundwater systems of varying complexity and heterogeneity transmit fluxes to surface water bodies, with varying degrees of influence depending on the hydrogeologic setting. 1.1 Integrated Hydrologic Modeling To develop predictive capabilities, we must first understand how climate and land use impact terrestrial hydrology. One of the primary tools to understand environmen- tal processes is long-term monitoring. Unfortunately, except for stream discharge, very few hydrologic time series are available for more than a few decades. Critical information, such as plot-scale soil moisture, groundwater levels, wetland inundation extents, shortwave solar radiation, plant phenological development, soil temperature, evapotranspiration, and many others are collected in very few locations and gener- ally over short periods only. As a result, direct water-balance calculations generally require estimation of one or several terms of the water budget. The result is that critical landscape and climate input heterogeneties are not resolved. Process-based models can be used both to understand the behavior of complex systems and as a predictive tool to forecast impacts of changes. Recent advances in processing capabilities, and the general avaiability of GIS data, makes the task of simulating the terrestrial water cycle at the relevant spatial and temporal scales feasible. There are dozens of widely-available hydrologic models, but many do not directly simulate physical proceses such as evaporation, transpiration, snow accu- mulation and melt, soil moisture vertical redistribution, groundwater recharge, and streamflow. Many do, however, simulate physical processes, though few meet the following set of requirements: 1. The code must be open-source, so that process descriptions in the model can be modified and examined 2. The model must simulate key hydrologic fluxes, including overland flow runoff and routing, groundwater recharge, surface and groundwater evapotranspira- tion, snow processes, and soil moisture redistribution. 3. Surface, near-surface, and ground water processes must all be explicitly de- scribed, not simply by linear reservoir routing. 4. Fine-resolution, large-domain modeling must be possible and feasible, both computationally and practically. This includes being able to run on modest equipment (desktop machines or small clusters), and incorporate GIS inputs to facilitate rapid development. Historically, hydrologic models have been developed for application at a narrow range of scales. Land surface models beneath global and regional climate models run at hourly timesteps over hundreds of square kilometer cells, but typically do not include comprehensive groundwater modules. Watershed models operate on catchments at a variety of scales while lumping the processes within those catchments into a collection of calibratable parameters. Plot-scale and one-dimensional root zone models simulate a host of physical processes with vertical layers as small as a few centimeters. Fully- coupled hydrologic models describe small catchments with dozens of vertical layers and cell or mesh sizes as small as a few square meters. None of these classes of models are well suited for simultaneously describing the locale, regional, and continental impacts of land use and climate change. More re- cently, a new class of integrated hydrologic model has arisen, coupling one-dimensional surface process and three dimensional saturated zone models. These have mostly adapated existing watershed models, some fully- or semi-distributed, others lumped parameter, and coupled them with existing groundwater modeling codes. The her- itage of the separate models varies widely, but all are an attempt to address the issue of fine-resolution, large-scale terrestrial hydrologic modeling. When my research began in 2004, few of these models were available. Those that were had been created prior to the general adoption of GIS technologies, and thus ill-suited for fine-resolution, large-domain modeling. In order to satisfy the four cri- teria above, a combination of existing and newly-development models were brought together into an integrated suite called the Integrated Landscape Hydrology Model (ILHM). ILHM is similar to this new class of integrated hydrologic models in that it ex- plicitly simulates the processes of the surface routing, canopy and root-zone, deep unsaturated-zone, and saturated zone domains. It was written in a modern, high- level programming language, MATLAB, that facilitates rapid development and test- ing. The code is modular, and developing and implementing new capabilities requires days, rather than weeks or months common to older codes written in compiled lan- guages. Eventually, the core code itself may be re-written in a compiled language to enhance performance and interoperability with other models. The development effort was focused along three pathways: 1) writing the core code and process capabilities, 2) writing an interface to the code that was fast, powerful, flexible, and scalable to eventually enable much larger simulations, and 3) developing a series of methods and practices for building large-domain, fine-resolution hourly climate inputs over a three decade period. In total, the code is now 70,000 lines, not including pre-compiled libraries or external interface tools. From the beginning, my intention has been that the code may eventually be adopted by others, thus significant effort was invested in extensive comments within the code, documenting the modules themselves, and creating a support infrastructure to handle multiple contributors and developers. ILHM is intended to be both a fully-explicit water— and energy-balance code, with run-time coupling between all domain modules (surface routing, canopy and root zone, unsaturated zone, and saturated groundwater). Currently, however, coupling is one way though full-coupling is under active development. Additionally, while the entire model preserves water-balance, only portions such as the snOWpack module and lake temperature module, preserve energy-balance. A new set of physical process modules have been developed that will complete the energy-balance capabilities of ILHM, but are not yet fully complete and tested. 1.2 Description of Chapters All of the chapters in this Thesis address the same common theme: how do climate and land use impact terrestrial hydrology? In the next Chapter, very simple mod- els are used in conjunction with time- and frequency-domain analytical techniques. The following two Chapters describe the development of ILHM, and its application first to a small catchment in central Lower Michigan, and then finally to a regional watershed and surrounding domain nearly of 19,000 km2. The Appendices present detailed equations and tables describing the structure of and processes simulated by ILHM. The final Chapter focuses on interpreting the results from a 28-year ILHM simulation to understand how local and regional climate, soil texture, and land use impact hydrologic fluxes. Chapters 2 and 3 were published in a slightly different form in 2007 (Kendall and Hyndman, 2007; Hyndman et al., 2007). They have been modified to fit the formatting of this document. CHAPTER 2 Examining Watershed Processes Using Spectral Analysis Methods Including the Scaled-Windowed Fourier Transform This chapter was published in 2007 in a slightly modified form in an AGU Geophysi- cal monograph, Subsurface Hydrology: Data Integration for Properties and Processes (Kendall and Hyndman, 2007). 2.1 Abstract Important characteristics of watershed processes can be extracted from hydrologic data using spectral methods. We extract quantitative information from precipita- tion, stream discharge, and groundwater head data from watersheds in northern- lower Michigan using Fourier Transform (FT) methods. By comparing the spectra of these data using similar units, we graphically illustrate the hydrologic processes that link precipitation to stream discharge and groundwater levels including evapo- transpiration. We also demonstrate how unit hydrographs can be efficiently and non- parametrically derived using the FT in a manner that allows for a quantitative sea- sonal comparison of precipitation and the resulting stream discharge response. This analysis clearly illustrates the reduction in summer discharge levels due to canopy interception and evapotranspiration. We also develop a systematic application of the FT we call the Scaled-Windowed Fourier Transform (SWFT), which extracts time— varying spectral content using a similar approach to the wavelet transform. While computationally less efficient than the wavelet transform, the SWFT allows for em- bedded detrending and tapering. Application of this method clearly illustrates the non-stationarities of spectral content within the three chosen data types, leading to a greater understanding of discharge-generating processes. 2.2 Introduction Spectral analysis (SA) provides a powerful means of extracting information from hy- drologic data. This type of analysis can reveal processes that may be obscured in direct time—series analysis by providing data not just on temporal fluctuations, but also on the spectrum of frequencies within those fluctuations. Furthermore, integrat- ing multiple datatypes in a quantitative and meaningful way is relatively simple using spectral analysis. Here we apply both existing and novel SA techniques to hydrologic data from two watersheds in northern lower-Michigan. We use these methods to il- luminate the processes that link precipitation to stream discharge and groundwater levels. Spectral analysis has been applied within the hydrologic sciences for decades (e.g., Bras and Rodriguez-[tame (1985); Hameed (1984)), most commonly as a means of estimating coefficients for linear autoregressive or stochastic models (e.g., Naff and Gatjahr (1983); Jukic and Denic-Jukic (2004); Zhang and Schilling (2004)). During the last decade, SA has been applied to examine fractal and multi-fractal behav- ior of systems (e.g. Tessier et al. (1996); Pelletier and Turcotte (1997); Kirchner et al. (2000)), and to examine linkages between hydrologic and climatic processes (eg. Tessier et al. (1996); Pellctier and Turcotte (1997); Coulibaly and Burn (2004)). Here we apply SA techniques to explore linkages between hydrologic processes and to provide a deeper understanding of those processes. Previous SA process comparison studies have generally not assured the similarity of measurement units, nor have they (with the exception of more recent wavelet studies including, e.g. Gaucherel (2002)and Coalibaly and Burn (2004)) considered the time—variant nature of process spectra. We demonstrate the quantitative utility of comparing data and spectra with similar physical units for inference of process relationships. Also, using an application of the common Fourier transform we call the Scaled-Windowed Fourier Transform (SWFT), we illustrate the non-stationary behavior of precipitation, stream discharge, and groundwater head spectra. Comparing Fourier spectra to the SWF T output, we illustrate how ignoring such non-stationarity can result in misinterpretations of spectra and thus of inferred process details. Finally, using a spectral derivation of seasonal unit hydrographs, we demonstrate how simple stream discharge models can be improved by considering the seasonality of watershed processes. Recognizing that spectral analysis is not a standard tool in the hydrologic science toolbox, Fleming et al. (2002)published a practical introduction to SA that focused on applications of the Fourier Transform (FT) including direct frequency domain inves- tigation, spectral filtering, and spectral simulation model validation. Here we apply SA in a similar manner, explaining our assumptions and the common complications, and demonstrating how these methods can be used across the hydrological sciences. Additionally, parts of the methods section are intended to contribute to a set of "best practices" of SA in the hydrologic sciences. 2.3 Datasets and Study Area Six data types from two northern lower-Michigan watersheds were used in this study: daily temperature, precipitation, and snowfall, both 15-minute and hourly stream discharge, and bi—hourly groundwater heads. Figure 2.1 shows the locations of both Michigan watersheds, as well as the locations of our groundwater transducers. The 3711 km2 Evart sub-basin of the Muskegon River Watershed (MRW) was selected for this study because there are no actively controlled flow structures on the Muskegon River or its tributaries upstream of this station. Because the Evart sub-basin lacks a network of monitored groundwater wells, we used water level data from our pressure transducers in the nearly adjacent Grand Traverse Bay Watershed (GTBW), which has a similar hydrogeological setting. Un- consolidated sediments within the Evart sub-basin and the GTBW were deposited by the same set of glacial episodes. These sediments are characterized primarily by coarse to fine sands and gravels with a small percentage of fine-grained material (Far- rand and Bell, 1982). The similarity in depositional history and topographic variation between these two basins suggests that GTBW groundwater head fluctuations can be used as a proxy for the nearby upper MRW. Preliminary 15- and 60-minute discharge data from the United States Geological Survey (USGS) stream gauge on the Muskegon River at Evart (Station ID: 4121500) were used in this study (USGS, 2005). Hourly data for this gauge were used from October 1989 until 15-minute discharge data became available in 1997. Combining the hourly data with resampled 15-minute data, the total data record extends from 10/1/1989 through 9/ 30/ 2004. For finer temporal-scale resolution, we processed hourly 4-km N EXRAD data (An- dresen, 2004) for 2004 and extracted the mean rainfall over the Evart sub-basin at each hourly interval. NEXRAD data Show a high degree of correlation to the cor- responding gauge values in this region (Jayawickreme and Hyndman, 2007). These 10 Grand Traverse Bay Wat rshed Muskegon River Watershed Watersheds 7 I: Evartsub-basln 7x72 Lake Michigan 0 10 20 40km Figure 2.1: Map of the Muskegon River and Grand Traverse Bay watersheds, with an inset map of Michigan showing their locations. Groundwater wells with transducers are marked with triangles, and the Evart gauge sub-basin is shaded in grey. 11 radar-based precipitation data are only available for approximately 9 months out of the year in this region (from 4/1/2004 to 11 / 30/ 2004) due to errors in snowfall estimates, so they are only suitable for evaluating shorter-period system behavior. We installed a network of 17 pressure transducers in USGS groundwater wells across the GTBW. Water table elevations were recorded every two hours beginning either in June 2003 (9 transducers) or June 2004 (8 transducers). Data through 7/ 1 / 2005 were used in this analysis. 2.4 Methods 2.4.1 Fourier Transform Fourier’s Theorem states that any complex periodic function can be decomposed into a set of periodic basis functions of varying amplitude, period, and phase shift. The most common such decomposition technique is the Fourier transform (FT), which uses sinusoidal basis functions. Fleming et al. (2002) present the basic theory, and a full mathematical treatment of this technique can be found in textbooks (e.g. Bras and Rodriguez-Iturbe (1985); Percival and Walden (1993)). A common implementation of the discrete FT is the discrete Fast Fourier Transform (FFT) popularized by Cooley and Tukey (1965). In this study, we use the FFTW libraries that are integrated into the MATLAB computing environment (Frigo and Johnson, 1998). This particular set of general-radix algorithms does not constrain the user to the requirement in Cooley and Tukey (1965) that the number of samples (N) be a power of 2. The output of the FFT algorithm is a complex array representing the magnitude and phase shifts of the Fourier coefficients. For instance, the F FT of a sinusoid of unit period and amplitude is an array with a single non-zero value corresponding to a period of 1. The FFT of a complex summation of sinusoids would produce non-zero 12 peaks in the spectrum. If instead the time-series input to the F F T were a broad, single-peaked curve such as a gaussian, exponential, or gamma function, spectral amplitudes would be non-zero across a broad range of periods. Fourier spectra are generally plotted as power vs. frequency, however we find that . plots of amplitude vs. period provide a more natural basis for viewing spectra of interrelated physical phenomena. The spectral amplitude of a time-domain input corresponds directly to hydrologic flux quantities such as precipitation and stream discharge. Note that the log-log slope of the amplitude spectrum is l; 6, where (3 is the slope of spectral power in log-log coordinates. A system is typically considered fractal (or multi-fractal) if its power spectrum roughly follows behavior (Avnir et al., 1998) If one assures the similarity of the units of each time-series dataset, the ampli- tudes of multiple spectra can be compared in a physically meaningful manner. To accomplish this, a suitable unit for comparison must first be chosen. For this study, we chose length/ time (L/ T) units because precipitation is the most common hydro- logic forcing mechanism. In order to change the volumetric discharge units (L3 / T) of stream discharge to L/ T units, the discharge was divided by the drainage area of the watershed upstream of the gauge. Groundwater head data, measured in L units can be differentiated to yield L / T units. Using the differentiated data, an approximation of the rate of head relaxation at the well location can be obtained by selecting only the negative values. Note, this technique assumes constant lateral inflow. This was then multiplied by an estimate of drainable porosity for the sediments (0.2), to correct for water level differences in porous media vs. open water. 2.4.2 Assuring Periodicity The FFT algorithm assumes that the data are periodic, namely the starting and ending points of the dataset are identical.Violating this condition results in spurious 13 features in the output spectrum (Bach and Meigen, 1999). However, time-series data from environmental systems rarely satisfy this criterion. There are four primary means of assuring periodicity: 1) filtering, not considered here, refer to Fleming et al. (2002) for a discussion, 2) taper function multiplication (tapering), 3) data subset selection (discussed in "Reducing Aliasing and Leakage" below), and 4) trend subtraction. Tapering can be a valid means of forcing periodicity if one considers how it affects the resultant spectra. Tapering refers to the multiplication of a mean-removed sig- nal by a “tapering function” (sometimes referred to as a windowing function) that smoothly tapers from a peak of 1 at its center to 0 at the edges. There are a variety of pre—defined tapering functions (Blackman and Tukey, 1958; Harris, 1978), and each affects the spectra differently. Tapering has the side effect of reducing the amount of information in the signal, thus limiting its applicability for short data records (Flem- ing et al., 2002). The shape of the tapering function, and how gradually it tapers near the edges, controls how much information is lost in the process. There is a tradeoff since the tapering functions that reduce information loss have more severe spectral “side lobes”, which distort the transformed spectra by shifting power from primary frequencies to harmonics (integer multiples or factors) of those frequencies. To recover the amplitude of an isolated peak within a spectrum (but not the entire spectrum itself), spectra from each tapering function can be corrected to account for the effect of side lobes. The multiplicative correction factor for each tapering function is given by the ratio of unaltered- to tapered-amplitude (Harris, 1978). In this study, tapering is primarily applied within the SWFT method where distortion is minimal since only a single amplitude is selected from each FT. Trend removal is often necessary for environmental data series to assure periodicity while minimizing loss or side-lobe distortion from tapering. The simplest method is to linearly remove the trend from the data, which can be effective if the non- 14 periodicity of the data is associated with a nearly linear trend. If there is some roughly sinusoidal long period fluctuation, a more valid means of trend removal may be to subtract a half-period sinusoid from the signal. In this case, the peak and trough of the subtracted function occur at each end of the original signal. A variety of additional trend-removal techniques have also been developed (Mann, 2004). We used the sinusoidal trend removal method in this study due to its simplicity and physical basis. 2.4.3 Reducing Aliasing and Leakage If the FT is blindly applied to data without thought to dominant system processes, sampling rates, or sampling interval, aliasing and leakage may occur. Both aliasing and leakage act to shift spectral power (or amplitude) from “true” frequencies to harmonics of those frequencies, although each acts differently. Aliasing results from under sampling high-frequency fluctuations (Bras and Rodriguez-Iturbe, 1985), while leakage or overspill is caused by both the non-periodicity of the system as well as non-periodicity of the processes within that system (Bach and Meigen, 1999). It is important to note that leakage will occur if the endpoints do not match but the inverse is not always true. Even if the time-series endpoints match, leakage may occur due to variability of the processes that contribute to the sampled time-series. In theory, aliasing can be avoided by merely increasing the sampling rate until the time series is fully resolved. Specifically, one cannot resolve spectral peaks with frequencies greater than half the sampling rate (the Nyquist frequency) (Bras and Rodriguez-Itarbe, 1985). If a time-series has significant power in frequencies above half the sampling rate, aliased power will be present in the empirical spectrum. In the case of many environmental datasets, aliased peaks may not be significant because as is shown below, these datasets are typically strongly damped at high frequencies. Our analysis indicates that some data do require sampling frequencies on the order 15 of once per hour, but most of those discussed here are sufficiently sampled with daily sampling rates. Spectral leakage can be more persistent and troubling than aliasing because peri- odic processes within a system are rarely sampled over an integer number of cycles. Leakage is commonly reduced by applying a tapering function to the time-series, de- trending the time-series, or both (Fleming et al., 2002). However, when the entire FT spectrum is of interest rather than specific spectral peaks, a more effective means of decreasing leakage in environmental datasets may be to carefully select subsets of the data that correspond to natural breaks in processes, thus assuring near-integer sampling without distortion. Data subset selection is also important because most environmental processes are non-stationary, thus each occurrence of a process may vary in both period and am- plitude. If one includes multiple cycles of a periodic process in the FT, the true peak location, shape, and amplitude can be obscured by diflerences in system states across cycles. Additionally, selecting a subset of data in which some system processes are inactive can also greatly simplify spectral analyses and reveal the spectra of weaker processes in portions of datasets that would otherwise be obscured by dominant, but intermittent, processes. For example, the diurnal fluctuation of stream discharge is often clearly visible during baseflow conditions but can be obscured by runoff and near-stream groundwater discharge during late spring and early summer. Alternately, if one were only interested in the spectral behavior of runoff, for instance, then select- ing data surrounding an isolated moderate-precipitation event during a dry season yields a stream discharge response primarily to direct precipitation and runoff. 2.4.4 Unit Response Functions While there are a variety of techniques to derive unit hydrographs from discharge and precipitation data, most are either ill-posed or require assumptions about system 16 behavior (Yang and Han, 2006). But direct FT deconvolution can produce unit hydrographs quickly and deterministically. The total time-series response of a linear system to a forcing input can derived by the convolution of the system unit response and the input time-series as follows (Smith, 1997): 00 q(t) = [h(T)p(t — T)dT (2.1) 0 where q is total time-series response (i.e. stream discharge), h is the unit response (for the case of stream discharge, this unit response has the special name of “unit hy— drograph”), and p is the input precipitation time series. According to the convolution theorem, Q(k) = H(k)P(/€) (2-2) where Q(k), H (k), and P(k) are the FTS of q, h, and p, respectively. The unit hydrogeraph time-series is then h(t) = F-1 [25%] (t) (2.3) where F “1[...](t)denotes the inverse Fourier transform. Thus the unit-response hy- drograph is given by the inverse FT of the ratio of the discharge and precipitation spectra. Though not strictly necessary, the analysis is simplified if the sampling frequency and units of the precipitation and discharge time-series are identical. The data should be resampled so that the number of samples, n, and the sampling frequency, f, are equivalent. The unit-response function, h(t), has length 11, however, in this case the unit response function is only valid up to the point where it becomes negative, since precipitation can not directly produce a decrease in stream discharge (Yang and Han, 2006). The negative response is therefore the signature of some other watershed process. If baseflow separation is used, this issue can be avoided. Though some small uncertainties will remain due to the data themselves, perhaps producing 17 rri negative responses. However, we chose not to use baseflow separation as it produces a synthetic dataset not directly tied to watershed process, thus conflicting with the investigative nature of these spectral analysis techniques. The resultant unit-response hydrograph is not an invariant property of the water- shed, as it is sensitive to variations in runoff-generating processes. These processes can be studied by directly comparing different unit-response hydrographs. Differences in response curve timing, peak, and shape can all be used to infer the activity and relative influence of various watershed processes. Applying data subset selection with these process differences in mind can allow for a quantitative sensitivity analysis of system sub-processes. In this study, we compare seasonally derived unit-response hydrographs for the Evart sub—basin averaged over 10 consecutive years. The derived unit hydrograph will be incorrect if stream discharge is still responding to precipitation inputs that occurred prior to the start of the data period (Smith, 1997). However, applied over entire seasons this error, as well as any error resulting from noisy data, is greatly reduced. Nevertheless, the derived unit-response function remains highly sensitive to edge conditions of the time-series inputs, thus the nominal time period (given by Table 2.1) of each season was adjusted to remove precipitation events or sudden increases in discharge near the edges of the time-series. The nominal time period for each season does not correspond to starting and ending dates of each season because they were chosen to assure similarity of hydrologic response within a season based on assumptions of process activity, also listed in Table 2.1. Also note that the fall and winter seasons overlap because the minimum length of the time-series subset must be longer than the watershed response time, which in this case was on the order of 60 days. To assure periodicity for the FT of each season’s data, a Tukey tapering function (Blackman and Tukey, 1958)was applied with relatively steep taper (coefficient of 0.1). 18 Table 2.1: List of the nominal time periods used in each seasonal analysis, and of the years omitted from the average unit responses. Data were from 1990-2001. Season Nominal Period Justification Years Omitted Winter 12/1-3/ 15 Mainly frozen precipitation 1991, 1992 Spring 3/15—5/ 15 Little evapotranspiration 1990, 1993-991 Summer 6/15-10/1 Maximum evapotranspiration 1990-91, 1993, 1995-96 Fall 10/15—12/25 Mostly liquid precipitation 1991, 1994-95 Note: I Only four years were included in the spring average discharge response, 1991, 19.92, 2000, 2001. Tapering was chosen over trend-removal because the magnitude of the discharge re- sponse to precipitation was of primary importance. After making these adjustments, some seasons continued to produce non-physical results (characterized primarily by sinusoidal unit-response behavior) and were thus omitted. These omissions are jus- tified on the grounds that the non-physical results reveal only the sensitivities and limitations of the method, and nothing about watershed process. 2.4.5 Scaled-Windowed Fourier Transform A key difficulty in applying the FT to environmental datasets is that non-stationarities in the data introduce artifacts in the Fourier spectrum. Here there are two types of non-stationarity to consider. The first is non-stationarity of process period, where subsequent cycles of a process within the system have slightly different spectra due to changes in system properties. The second is best described as intermittency, where a process may be active during only a portion of the sampled data. The first results in spreading of spectral power about a central period (if the process is sinusoidal), while the second results in spurious spectral power at harmonics of the primary period. To avoid these effects and more clearly illuminate important changes in the spec- tra, we developed a method that we call the Scaled-Windowed Fourier Transform (SWFT). The SWFT is very similar to a sinusoidal wavelet transform (WT), but it differs in a number of respects. The mathematical development of the SWFT pre- 19 sented here is fundamentally different from that commonly presented for the WT, in a way we feel increases the utility of this method of hydrologic scientists. In particular we feel that there is great value in demonstrating that the SWF T produces Fourier coefficients localized in both the time- and frequency-domains, as opposed to the sim- ilarly localized WT coefficients that are wavelet-dependent. Additionally, the SWF T is capable of embedded detrending rather than relying on tapering alone to reduce leakage, potentially producing better spectral estimates. Finally as developed, the periods and times queried by the SWF T are more flexible than typical WT schemes, with the tradeoff of decreased computational efficiency.The SWFT also differs from the traditional windowed FT that only transforms data within a specific subset of the overall time series called a data window (here the window is different from a tapering function). This data window is then slid along the time series to produce a map of spectral power varying in both period and time. Unfortunately, the windowed FT forces a choice between severe aliasing of low-frequency components of the signal or poor resolution of high-frequency non-stationary processes ( Torrence and Compo, 1997). By contrast to the windowed FT, the SWF T method scales the width of the window over successive passes along the time-series. At each window position, the data are detrended, multiplied by a tapering function, and Fourier transformed. A single amplitude corresponding to the Fourier coefficient of a single frequency is selected from the complete FT at each window position. The window is then slid along the data, producing a time-varying series of amplitudes for that frequency. When the end of the dataset is reached the window width is rescaled and the process is repeated for the next frequency. Note, hereafter we use the word "period" solely to indicate the spectral period, or the inverse of frequency. We feel that the period of spectral content is generally more applicable to the hydrologic sciences than the frequency. The SWFT produces the same type of scalogram as would be generated with the 20 WT. These scalograrns are well defined mathematically and physically (if proper units are used), and can be examined using standard statistical techniques. Aside from simple examination of the scalogram, comparisons of related scalograrns such as the cross-scalogram and the coherence phase map (see Terrence and Compo (1997)) can be calculated as well. Conceptually, the SW FT scalogram is produced via the following: Fm 2 0' FUJI (jstart ijend) — 0] 'le (2.4) where p and q are indices defined mathematically below, which correspond to the periods and times at which the SWFT transform is applied. F ( . . )1: indicates a single value from the discrete Fourier transform spectrum, a: is the time-series dataset, .lstart and jend are the beginning and ending indices of the current data window. T is a tapering function (optional), and D is a detrending function (optional), each defined over the current data window. C is a multiplicative correction factor unique to each type of tapering function that is computed as C = IF (sin (—-oo : +oo))k| / |F(sin(—oo : +00) 'lel where k is the index that corresponds to the period 27r. Note that a reasonable approximation of C can be obtained using just three or four cycles of the sine function. If a tapering function is not used, C = 1. For our analysis, we chose a Tukey window with a gradual taper (Tukey window coefficient of 0.5) as a tradeoff between frequency resolution and side lobe distortion, resulting in C m 1.33. Thkey coefficients closer to 1 produce greater side-lobe distortion, while those nearer to 0 reduce side-lobe distortion but increase leakage. The complex definition of the discrete Fourier transform (DFT) modified from Press et al. (1992) is Tl . 27a . -_ Fk .__ :Ije-TM-UU 1) (2.5) i=1 where :1: is the time-series dataset, and k and j are indices running from k = [1, . . . ,n] and j = [1, . . . , n] and n is the total number of data points to be transformed. The 21 DF T spectrum F has n points of which the first is the "gain" term, and the next n/2 points correspond to the periods n/f - [1,1/2,1/3,...,2/n]. To compute the SWFT, we would like to extract a single Fourier coefficient cor- responding to a certain period, Pp, from the entire DFT spectrum. Also, we need the coefficient not for the entire dataset x, but for a windowed subset xj with inde- ces j = [nstart, . . . ,nend] (which relate to times t = j / f) where the total number of points in the window is Np. To minimize leakage, the window width Np should be chosen such that an integer number of sinusoidal cycles fit within it, which can expressed as N, = P, - f - w. (2.6) Here w is an integer multiple with possible values [1,2, . . . ,floor (n/l’p - f )]where "floor" is a function that rounds the argument toward zero. In general, higher values of a: increase the frequency resolution of the SWFT while decreasing the time resolution. Note that the entire dataset :r is used, and the standard DFT Fourier coefficient is produced when w is set equal to 72/ Pp - f. When Equation 2.5 is applied to the data window Np, we note that the indices k in the DFT output then correspond to periods Np/f- [1, 1/2, 1/3, . . . ,2/Np], and the first point in the DF T output is the DC gain term. The SWF T requires only the value 1",, corresponding to Pp. So P = Np/ f - 1/w, and therefore (replacing k with p) we get a=a=an an Substituting Equations 2.6 and 2.7 into 2.5, and recalling the definition of N , we get "end —2'l7f'lll(j_1) Pp: Z :L‘je 1if) . j=nstart The window of width Np is slid along the dataset producing an array PM where q (2.8) corresponds to the center-window time qu with indices p and q = [1,2, . . . ,qma.:rp] 22 . The maximum resolution of P is such that the product Pmaa: - f = [2,3, . . . ,n], since the shortest period allowed is given by the Nyquist critical frequency (2/ f )_1. However, for computational efficiency when Pmax spans several orders of magnitude, the user can specify any subset P of Pmar. For instance, one could choose P - f = [2,5,12,14,...,n], or any other arbitrary subset of Pmaz. The index p = [1. . . . ,length(P)] then refers to the periods at which the SWF T will be calculated, where “length” is a function that calculates the size of the 1St dimension of the array P. Substituting Equation 2.8 into Equation 2.4 yields the complete definition of the SWFT: "end —2i7rw (1n— 1) j=nstart where m is an index [1,2, . . . , Pp - f3 - w], or m = (j + 1) — n-start . ”start and new are functions of p and q: "start = (q - 1) - Pp - f - w +1 = (q — 1) - Np +1, and (2.10) ”end = "start + Np - 1- (2.11) Because the Fourier coefficients can only be time-localized to a window of width Np, we include the capability for each data window to overlap the previous one in order to increase the temporal resolution of the transform. As an example, if Np = 10 and no overlap is allowed, the minimum temporal resolution at this value of Pp would be 10. Instead overlap is allowed such that the minimum resolution can be as low as 1. This could either be done with a fixed overlap (i.e., each window overlaps half of the other across all periods), or the overlap can scaled in any other fashion, such as linearly with P. For this study, we define the overlap value, op to be given by a simple linear scaling as p 0P : floor 0min + “fix-E755 (Omar — 0min) (2'12) 23 where 0mm,- and 0min are specified by the user, and max(P) is the maximum value within the array P. The minimum value of 0mm = 1, and the maximum value of 0mm = Np. With this modification, Equation 2.10 becomes ”start = ((1 — 1) ' Np/Op + 1 (2'13) and q(max)p is given by q(max)p = floor (n — Np) - [of +1 . (2.14) ’p The amplitude or power spectra can be extracted from the full SWFT spectrum in the same way it would be for the standard FT. Here we use the amplitude spectrum that is calculated via qu = 2 - [qul /N,, (2.15) where the vertical bars indicate the magnitude of the complex value qu and the factor of 2 arises because the DFT spectrum is symmetric about the vertical axis and thus distributes half of the spectral power at a period to the each of the positive and negative instances of that period. The center-window times T associated with the arrays F' and A are given by qu=f'[(q—1)1:—:+%[=f'[Np(q—;)—l+%)[- (2-16) Note that F, A, and T are, in general non-rectangular arrays (the exception is when w = Np). MATLAB’S cell array capability was used to store these values. For convenience of both visualization and further processing, such as contouring or com- puting cross-scalograms and phase-coherence maps, these arrays can be interpolated to rectangular grids. The SWFT array F can either be computed directly via Equation 2.9, or it can be computed via Equation 2.4 using the values of jstart and jend from Equations 2~13 and 2.11. In either case, the tapering and detrending functions T and D are calculated using the data a:(jstm~t I km)- 24 The SWFT was developed primarily for flexibly visualizing the non-stationary spec- tral content of a time-domain signal. Including the ability for 0,, to scale with pe- riod greatly reduces the computational demand by calculating qu less frequently for longer periods. Allowing w to vary enables flexibility between time— and frequency- localization, as the needs of the user demand. Though not used here, w could also be a function of p allowing the frequency resolution to also scale with the period. Including explicit tapering and detrending further improves the ability of the SWFT to represent dynamic spectral content when P spans several orders of magnitude. In order to illustrate the interpretation of the SWFT scalogram, we examine a simple test case using a summation of three separate sinusoids given by: f (:11) = f1(:1:) + f2(a:) + f3(:z:), where f1(:r) = sin(4a:) over 0 S :1: S 127r, and f2(:1:) = sin(10:r) over 0 S :r g 127r, and 0 0S$£4n,87r_<_:r£127r f3(l‘)= . sin(:1:) 4 < (I: < 87r Figure 2.2a illustrates the successive superposition of the three sinusoids, the darkest curve plots f3(.r), the mid-tone curve plots f3(.’17)+ f2(a:), and the lightest curve plots f (:17)- The SWF T scalogram (Figure 2.2b) of the function f (3:) reveals the periods, ampli- tudes, and ranges of activity of each of the simple sinusoids. The method extracts the periods of the three sinusoids and reconstructs the amplitudes accurately, except for the f3(:r). This is simply because only two cycles of the sinusoid were used in f3(:v), and the window width was chosen as twice the period, thus only the very center point should reach an amplitude of 1. The shorter-period sinusoids both have peak ampli- tudes very near 1, although the middle has amplitudes >1 in some locations because 0f the interaction between the two longer-period sinusoids. Importantly the longest Period sinusoid, which is defined only for 47r < :r < 87r, only has large amplitudes in this range, thus revealing the time-varying spectral behavior of the input signal. 25 o is 1‘0 1'5 2'0. 25 x-Locatlon 0.5 0.6 0.7_ 943.. 0.9 1.0 FourieriAmplitude Figure 2.2: SWFT results for a test function, f(x). a) Curves representing the successive superposition of the three sinusoids with different frequencies, and b) the SWFT scalogram. The shaded contours represent amplitudes between 0.5 and 1.0 as shown in the colorbar. The dashed line represents the boundary beyond which the SWFT scalogram is undefined (where half the window width is greater than the available number of points within the dataset). 26 In cases where the input is something other than a summation of sinusoids, the scalogram output will exhibit large amplitudes across a range of periods, as expected from Fourier theory. If the broad-spectrum behavior of the input data is relatively time-localized, such as the runoff response to a brief storm event, the scalogram will exhibit lineations corresponding to that event. Thus, temporally continuous but spectrally limited high amplitude regions indicate periodic processes within the input data while temporally limited but spectrally broad lineations indicate broad-spectrum processes. 2.4.6 Scalogram Averaging Averaging the scalogram amplitudes across periods or time gives the period-averaged or time-averaged amplitude spectra, respectively, of the SWFT scalogram ( Torrence and Compo, 1997). All of these three averages are computed using a rectangularly interpolated grid, Asr, at user-selected periods 3 and times r, from qu. Since the period spacing in the SWFT scalogram is not necessarily uniform, we use the period- weighted mean amplitude Ap, calculated using PA Apr: :5 =1 3 ————i‘", (2.17) 23:13? where L is the dimension of the rectangular grid in the period direction. This provides aggregate information on the temporal variation of spectral amplitude. The time-averaged amplitude At, given by 1 M Where M is the dimension of A in the time direction, displays average amplitudes across time at a single period. Referred to as the global SWFT spectrum, this time- aVeraged spectrum is qualitatively similar to the FT spectrum, but differs in physical meaning. 27 ‘l t“ Bulk changes in the relative influence of short- vs. long-period fluctuations across time can be visualized using the amplitude-weighted average period Pt, = :3le PSAST' 25:11:13? All three of these scalogram averages are demonstrated below. Pi. (2.19) 2.4.7 Watershed Available Precipitation: Snowmelt Model- ing Although this study focuses on revealing and exploring watershed process using a purely data-driven SA, one model is required for the analyses. Because precipitation falls as snow during most of the winter months in northern lower-Michigan, a snow storage-and-release model is needed. The data required for a full energy / water bal- ance model were either unavailable or unreliable for the entire data period, thus we implemented a simple heuristic snowmelt model. Because data are available for both fresh snow totals and observed snow depth, the model can simply identify when a decrease in snow pack thickness corresponds to melting event or densification of the snow pack. We then compared one year of the output from the heuristic model to the full energy-balance UEB snow model (Tarboton and Luce, 1996) in order to assure that it acceptably predicts both snowmelt timing and volume. This heuristic model tracks snow water equivalent and releases snowmelt based on a three-part rule structure: 1. The density of newly-fallen snow is calculated from daily precipitation and snow fall totals. Since precipitation can be mixed frozen and liquid, a maximum new- snow density cutoff, newmax, was determined from the data. Any precipitation in excess of this cutoff is considered equivalent to rainfall and immediately released. 2. The total snowpack density is updated based on the new snow depth and the 28 accumulated water content of the pack. 3. If the average snowpack density exceeds the maximum pack density, packmax, it is assumed that some of the snow has melted. Melt is then generated equivalent to the depth of the pack multiplied by the difference between the calculated pack density and the maximum density. There are three important assumptions in this model. First, we assume that drifting does not affect the observed snow depths at the measurement location (typical of snow models). Second, we assume that the snow pack properties are uniform throughout, which is realistic in the MRW region as total pack thicknesses are typically less than a third of a meter. Finally, this model assumes that the density at which the snow pack releases water remains the same throughout the season. The maximum densities of the new snowfall and the snowpack are physical quanti- ties that are generally not equal. Maximum snowpack density, poo/imam, is determined by both its composition and water holding capacity, which are functions of the ther- mal history of the pack and new snowfall conditions. Direct comparison of snow depth and stream discharge in our study area suggests that water is released from the snowpack when the combined snow / water density reaches approximately 0.35 g/cm3, according to this model, thus this value was chosen for packmax. The maximum new snowfall, newmax density of 0.23 g/cm3 (determined as the ratio of new snow water content to new snow depth) was extracted from the data, as this was the relatively abrupt limit above which higher-density new snowfall events were obvious outliers. The model was tested relative to the UEB snow model using data from the win- ter of 1999/ 2000. The solid line in Figure 2.3a is the heuristic snowmelt model, and the dashed line is the output from the UEB snow model. Both models output the combined snowmelt and liquid precipitation, referred to hereafter as watershed- aVailable precipitation. Watershed-available precipitation is that which can be acted Upon by physical or biological processes. Note that the models provide very similar 29 a . 4 - —Heuristic Model i -Energy Balance Model; W.A. precipitation (cm) E A 3 5 E Q) g -2.5 2 9 2' 0 e if; -1 . Difference ~ -252, § -2 . . . . -5 E as Dec Jan Feb Mar Apr May a g 1999 <-|—> 2000 a Figure 2.3: a) Comparison of the heuristic snowmelt model to the UEB model, b) the difference between the two (heuristic-UEB) and the cumulative difference. results, with the UEB model predicting slightly more melt early in the season and the heuristic model predicting more late in the season. The heuristic model predicts ap- proximately 1.5 cm more snowmelt during the season than the UEB model. However, the heuristic model does not allow for sublimation, which entirely accounts for the ~1.5 cm difference at the end of the modeled period. The heuristic model was chosen for our analysis since it performs acceptably, despite minimal data requirements and a simple structure. 30 2.5 Results and Discussion 2.5.1 Examining Watershed Processes: Spectral Comparison The spectra of stream discharge, watershed-available precipitation, and relaxation of the water table elevation in well B13 can be directly compared to infer interaction timescales between the data types and to examine details of processes within each type. Well 313 (data not plotted) was chosen because the water table is deep enough (>30 meters) to preclude direct evapotranspiration effects. Four key features of the spectra will be compared to integrate the hydrologic datasets and reveal process de- tails: 1) spectral peaks, 2) log-log linear SIOpes (i.e., fractal scaling), 3) locations of slope-breaks, and 4) relative spectral amplitudes. Since watershed-available precip- itation is the primary forcing function for natural watershed processes in our study region, the spectra of well-head relaxation and stream discharge can be viewed as modifications, or fractal filters (Kirchner et al., 2000), of the watershed-available pre- cipitation spectrum. The same is true to some degree for the well-head relaxation and stream discharge spectra, as groundwater inputs to the Muskegon River account for a majority of its annual discharge (Jayawickreme and Hyndman, 2007). The spectra of these three data types are plotted in Figure 2.4 along with log-log linear fits and 95% confidence intervals. Also plotted is fixed period-width binned spectrum to aid visualization. Slopes for selected linear portions of the spectra were calculated using least-squares regression between user selected bounds. These bounds were chosen to match portions of the spectrum that exhibited a linear slope. The Slope breakpoints are then calculated at the intercepts of the separate linear fits. The 95% confidence interval (dotted line) is determined by multiplying the x-square value for a system with 1 degree of freedom by the average amplitude given by a noise (or scaling) model (Torrence and Compo, 1997). In this case, the noise model was assumed to be given by the linear fits. This enables the flexibility of applying 31 a , a E ....... ‘.‘"" E -3 .g ....... g .2- , r » .3 4 '6. """""""" '1 Precip. 3. l E ; - — Well i , 3 3 ' ; .-. Stream I 01 1 l o 1 2 ‘i s log(period) (d) yé r . . . . .l _. , '. - \J . . ................... log(amplitude) (m/d) I I I I I I : ‘l ' ; t: yéasr . log(penod) (d) Figure 2.4: Composite plot of stream discharge, well-head relaxation, and precipi- tation spectra. Grey lines plot the raw FT spectra, solid lines are the binned-mean amplitudes, and dashed lines give the upper 95% confidence intervals. a) An overlay of the spectra in parts c-e and their linear fits, b) NEXRAD precipitation spec- trum for 4/ 1 / 2004-11 / 30/ 2004, not plotted on part a, c) stream discharge spectra for 10/1/1997—9/30/2004, d) well-head relaxation for well B13 from 7/ 1 / 2003-7/ 1 / 2005, and e) watershed-available precipitation from 1/1/1990-1/1/2001. Note that, in part c the limited resolution makes it appear as if much of the stream discharge spectrum is above the confidence limit, however fewer than 5% of points exceed the limit in a spectrum with over 130,000 points. 32 statistical confidence intervals to datasets without assuming a priori a particular type of noise. This is useful when working with spectra that exhibit multi-scaling behavior ( Tessier et al., 1996; Dahlstedt and Jensen, 2004), where a single-scaling noise model, and therefore confidence test, would be inadequate. All three process spectra have annual-cycle peaks at or near 365 days, although the peak in the head relaxation data is significantly below the 95% CI boundary. This is probably due to the relatively short data record available for the head relaxation. The discharge spectrum has a series of peaks at the harmonics of the 1 day peak that are artifacts, as later discussion will demonstrate. The head relaxation spectrum also has a peak near 117 days, along with weaker peaks near 70 and 50 days (Figure 2.4d). These are near the integer factors 3, 5, and 7 of the 365 day peak, again suggesting a spectral artifact rather than processes acting at these periods. The watershed- available precipitation spectrum has two additional peaks at "174 days, and another at ~65 days. These may not indicate sinusoidal processes active at those periods, but may instead be the spectral signature of a non-sinusoidal process characterized primarily by a longer-period oscillation. In particular these two peaks are near the integer factors 2 and 6 of the primary 365 day peak. The slopes of the spectra and SIOpe breaks (Table 2.2) provide provocative evidence of linkages between watershed processes. If a quantity, such as stream discharge, is being forced directly by another, such as precipitation, then fractal scaling active in the forcing input should exhibit itself directly in the response variable (Tessier et al., 1996). However, a non-linear system response behaves like a fractal filter, modifying the input scaling behavior (Kirchner et al., 2000). There are two exam- ples of this in Figure 2.4 and Table 2.2. The first is the segment of the discharge spectrum between one and three hours with a 6 of 1.5. This is roughly similar to the slope of the NEXRAD hourly precipitation spectrum (data provided by Andresen (2004)). The Slope break locations and 6 values for Figure 2.4. 6 = 0.9 slope in the 33 Table 2.2: Slope break locations and 6 values for Figure 2.4. Discharge Watershed-Available Precip Head Relaxation Head Fluctuation Break 063 Break 1 61 Break l3 Break2 62 1.0 hour 15 0.90 0.03 3.0 hours 2.9 18.2 hours 0.30 -0.08 14.0 hours 1.9 15.4 days 1.1 9.6 days 0.01 16.5 days 1.2 19.3 days 4.4 Note: I These bold italicized entries correspond to results from the NEXRAD dataset, normal entries refer to the watershed-available precipitation from gage data. 2 These italicized entries correspond to slopes from Well BIS, whose binned-mean spec- trum is unplotted. NEXRAD spectrum is an underestimate of the true spectral slope, given that the hourly NEXRAD spectrum represents just a single year of data that undersamples precipitation, and thus suffers from high-frequency aliasing and slope-flattening as described in Kirchner (2005). If, as indicated in these empirical spectra, the “true” spectral slopes of streamflow and precipitation match in this range, we interpret the similarities of slope up to a period of three hours to indicate that direct precipita- tion is the dominant flow-generation process. Beyond three hours, processes with a different scaling relationship dominate flow. Process linkages are also apparent between long-period stream discharge and head relaxation spectra. Both spectra have a slope break at approximately 16 days and follow a ,0 = 1.1 — 1.2 scaling relationship to longer periods. This suggests that varia- tions in groundwater discharge control stream discharge variability at periods longer than approximately 16 days. The exact value of this slope break is approximate be- cause of the gradual transition between linear segments in the discharge spectrum between approximately 10 and 30 days. Zhang and Schilling (2004) observed slope breaks in Iowa streams at approximately 30 days. The similarity in slopes between head-relaxation and discharge in Figure 2.4 suggests that groundwater inputs dorn- inate streamflow in these Midwestern streams for periods longer than 10-30 days while in-stream and near-surface watershed processes appear to dominate at shorter 34 periods. The spectral slope values also provide useful information. Particularly interesting is the 6 = 2.9 slope seen at periods between about 3 hours and 16 days in the discharge spectrum. The uniformity of scaling in this portion of the spectrum means that any watershed processes active in this period range also exhibits similar scaling. Fundamentally, this is because linear processes that are additive in the time-domain also add in the spectral domain (Smith, 1997). If the scaling of any hydrologically significant watershed process differed from the others, it would preclude the observed uniform scaling. This uniformity is interesting considering the variety of watershed processes active in this period range, including bank storage and release, precipitation runoff, near-surface “interflow”, near-strearn saturated groundwater response, and evapotranspiration. Further investigation of this uniformity could be undertaken with a more concentrated set of data designed to explore these processes, or using a detailed process-based hydrological model. Analysis of the relative amplitudes of the three spectra in Figure 2.4a provides a hydrologic response time to precipitation. The amplitudes of watershed-available precipitation and discharge converge at a period of approximately 2.5 years. Thus any long-period fluctuation in precipitation will be followed by an equal magnitude fluctuation in stream discharge. Therefore this watershed and its unsaturated and saturated groundwater processes do not appear to control hydrologic fluxes with pe- riodicities greater than 2.5 years. This observation could be used in autoregressive models of basins with short discharge records but more extensive precipitation data. Another interpretation of the similarity in magnitude between the precipitation and discharge spectra is that the combined response time of both surface and groundwa- ter systems is approximately 2.5 years under current conditions, although extended drought periods could certainly affect this value. Therefore, in order to assure in- sensitivity to initial conditions, a model of this watershed must be “spun up” with 35 realistic meteorological inputs for at least 2.5 years. Beyond 10-30 days, the well-head relaxation amplitudes are approximately 2-3 times greater than those of stream discharge. This disparity is likely a combination of both physical processes as well as our assumptions. First, the well-head relaxation amplitudes would exceed those of stream discharge because evapotranspirative losses. Also, the head in an individual well may fluctuate more than the average of all wells in the watershed, particularly if the water table at that well is deep. Thus a more representative comparison would be between stream discharge and the average of spectra from wells distributed across the watershed. However, in this case many of our wells showed signs of periodic anthropogenic disturbance that violates the assumptions of our simple differencing technique. Other factors that contribute to the disparity between head relaxation and stream discharge amplitudes may include overestimation of porosity or average recharge rate because of the short data record in the well, as well as differences in annual recharge between the Evart watershed and the location of the well in the Grand Traverse Bay region. The head relaxation results presented in Table 2.2 are taken from the spectrum of smoothed-differentiated head fluctuations shown in Figure 2.4d while the head fluc- tuation data in Table 2.21 were taken from the spectrum of actual head fluctuations (not plotted). Unlike head fluctuations, head relaxation (and therefore a decrease in storage) is directly related to groundwater discharge, allowing direct comparison of their spectra. Other studies have reported head fluctuation spectra, (Zhang and Schilling, 2004; Lee and Lee, 2000; Nafir and Cutjahr, 1983), but because the basic units of head fluctuation and stream discharge differed in these studies, their reported groundwater head and baseflow scaling laws can not directly be related. 36 2.5.2 Examining Intra-Annual Spectral Variability Using the SWFT Physical interpretations of Fourier transform spectra assume that system processes are both non-intermittent and stationary. However, many watershed processes are either inactive for parts of the year (intermittent) or possess different spectral charac- teristics over successive cycles (non-stationary). Thus, interpreting spectra from many occurrences of a given process is problematic. To overcome this, we use the Scaled- Windowed Fourier Transform (SWFT), which does not require either stationarity or non-intermittency. As Figures 2.5-2.8 demonstrate, all three watershed processes examined in this study exhibit both non-stationary and intermit processes. The SWFT spectrum of stream discharge (Figure 2.5) displays its time-varying Fourier spectral content, revealing details about how discharge responds to hydro- logic inputs and watershed processes. As described previously, the vertical lineations in the SWFT scalogram are caused by the broad spectrum of time-series stream dis- charge peaks with each such lineation corresponding to a pulse increase in stream discharge. Clearly separated lineations, which primarily occur during the summer months, extend from very short periods and tend to reach maximum amplitudes at periods between 20-40 days. This range of periods corresponds to the width of the time-series discharge peak, and thus to the time-scale of surface and near-surface wa- tershed response to precipitation inputs. As will be shown in the next section, this 20-40 day time scale of surface and near-surface hydrologic response is also similar to the width of the unit hydrographs developed for this watershed. Spring months typically exhibit increased Fourier amplitudes relative to other sea- sons across all periods shorter than several hundred days. The spring of 1998 has a very different spectrum than that of 1996 or 1997. Those years had several spring discharge peaks followed by relatively high summer baseflow levels, whereas the single discharge peak of 1998 is followed by low baseflow and weakening of the 365 day am- 37 ._.../\ ~..>/\ s2 2%: ‘U .9 . 00: 0. 51,5 <..- 01x >A2. $1: 5 2g“ 8 v (0 Period (d) 10'4 10'3 Time-Averaged Amplitude (m/d) “(m/rd) ’l . Figure 2.5: SWFT of stream discharge at the USGS gage in Evart, MI. 8.) Stream discharge time-series divided by the area of the watershed above this gage, b) Scale- averaged amplitudes, Ap, c) Filled contour scalogram of the amplitude (A) vs. both period (P) and time (T), along with the location of the amplitude-weighted mean period, Pt (white line), d) Time-averaged amplitude spectrum, At. 38 plitude. This summer baseflow portion of the time-series discharge is accompanied by decreased amplitudes across nearly all periods, except near 365 days. Thus, stream discharge during this portion of the year is dominated by long-period fluctuation, most likely from groundwater inputs, as indicated by the amplitude—weighted mean period curve (Pt). Seasonal differences in spectral character between the summer and spring/ fall are not evident in the watershed-available precipitation SWFT scalogram (Figure 2.6c), but show up very prominently in discharge (Figure 2.5c). This suggests that although there is spectral power in that period range in watershed-available precipitation, sum- mer evapotranspiration and canopy interception decrease the magnitude of discharge response and thus the spectral amplitudes. Another important difference between Figures 2.5c and 2.6c is that the dominant power the in precipitation spectrum occurs at shorter periods while the reverse is true for stream discharge. This corresponds to the behavior seen in Figure 2.4c, however the time-averaged amplitude (At) spectrum of watershed-available precipitation dif- fers from the Fourier spectrum in Figure 2.4e. Such differences are artifacts produced by the violations of the assumptions of stationarity and non-intermittency inherent in the standard FT. The SWFT discharge scalogram for the water years 10/1/1991-9/30/2004 (Figure 2.7) displays the spectral content of stream discharge at the Evart gage over periods of 0.3-1000 days. Although the longer-period spectrum is not plotted, 2-3 year and 6-8 year cycles are evident in periods between 100-400 days, perhaps related to climate cycles as seen in Coulibaly and Burn (2004). For 1991-93, long-period fluctuations are most active near a period of 180 days, which then cease during 1994-95 before resuming from 1996-98 at 365 days. Again, this long-period activity switches off for two years and then resumes centered about 180 days. Another interesting set of features in this scalogram are the summer-fall periods of 1998, 2000, and 2002 in 39 3: O O1 Precip. (m/d) ---~-0.025 ' Lil 0 .W.A d 10 5’2 =2 >< 8 ‘ <12 E 5 .. ~ 4 m .- 5; 1 'O o ’I: 0 o. 101 3 5 3 2 5 Time-Averaged log(amplitude, I .. ._.,, . _,. 1., Amplltude Figure 2.6: SWFT of watershed-available precipitation at Big Rapids, MI. a) watershed-available precipitation time-series, b) Scale-averaged amplitudes, Ap, c) Filled contour scalogram of the amplitudes (A) vs. both period (P) and time (T), along with the location of the amplitude-weighted mean period, Pt (white line), Note, white-areas have amplitudes < log(amplitude)=-3.5, d) Time—averaged ampli- tude spectrum, At. 40 Discharge (m 3/s) Period (d) L, i it- ? Pi E ii F 1994 1996 1998 2000 2002 2004 -5.5 -5 . -4 -3.5 1‘: 1992 log(amplitude)<'6 (m) l I In Figure 2.7: SWFT of stream discharge at Evart from 9/ 1 / 1991 to 8/31/2004. 21) Stream discharge time-series, b) Filled contour scalogram of the amplitude (A) vs. both period (P) and time (T). The year labels are centered at January 1. which amplitudes are drastically reduced up to periods on the order of 100-200 days. These observations from Figure 2.7 enable a deeper understanding of time-averaged spectra such as the FT spectrum or the SWFT time-averaged amplitude (At). The FT spectrum of discharge (Figure 2.40) contains a weak spectral peak at 180 days along with a peak at 365 days. Prior to examining the spectrum as a scalogram, we were unable to distinguish between dominant spectral peak harmonics and peaks from independent processes at those periods. The scalogram shows that there are processes generating true peaks at both the 365-day and 180—day periods. Additionally, we can use the information from the scalogram to indicate when pro- cesses that generate particular spectral peaks are most active. A prominent example of this is the 1 day peak in the FT spectrum of discharge. Plausible interpretations of 41 this 1 day peak include diurnal fluctuation in evapotranspiration or streambed con- ductance during the summer months. However, Figure 2.7 shows that the dominant 1 day amplitudes occur in the winter and early spring months. A close examination of the time-series reveals two important details: 1) the diurnal fluctuation is strongest during periods where daily maximum temperatures are subfreezing, and 2) discharge peaks during the coldest mid-morning hours of each day. These observations suggest that the diurnal signal may be related to icing effects at the instrument. The 1-day peak seen in Figure 2.4 is a “true” spectral peak, but it does not indicate cyclic sys- tem processes so much as inaccuracies in measurements. Also, there are no significant amplitude peaks at the 0.5-day period, thus confirming that the corresponding peak in Figure 2.4c was indeed a harmonic. This information could then be used to fil- ter the discharge time-series and remove this apparently erroneous periodic discharge behavior. The SWF T of head fluctuations from Well 810 (Figure 2.8) provides another ex- ample of the importance of viewing spectra as a function of both period and time. Well B10 was chosen because it exhibited the greatest amplitude of fluctuation in the period 1-30 days of the 17 wells in our study. Figure 2.8 displays both the head fluctuation data (mean removed) as well as the SWFT scalogram for periods between 0.3 and 100 days. The largest time-series amplitude head fluctuations occur from late October to early April. The scalogram reveals that most of the time series fluc- tuation is caused by larger Fourier amplitudes of the 1-30 day periods, which are greatest during the winter. Note that the dominant period in this range is not con- stant throughout the year. During the summer, early fall, and spring, the dominant period is near the 7-10 day range, that then shifts to the shorter 2-7 day range in mid-November. The SWFT scalogram clearly reveals this information that would be difficult to directly extract from the time-series data. 42 Mean-Removed Head (m) Period (d) b .n- . mi _ “.i'iwwi ll; 1* Nov Jan Mar 2004+[—> 2005 log(amplitude) '3-‘5 4? '52-'75 '2 (m) i . Figure 2.8: SWFT of groundwater-head fluctuations in Well B10 in the Grand Tra- verse Bay Watershed. a) Well—head fluctuation time-series, b) Filled contour scalo- gram of the amplitudes (A) vs. both period (P) and time (T). 2.5.3 Spectrally-Derived Watershed Precipitation Response Functions The watershed annual unit response functions for the portion of the Muskegon River Watershed above Evart gauge were calculated using discharge and Big Rapids pre- cipitation data between 9/1/1999 and 8/ 31 / 2000 (chosen arbitrarily). Two different unit response functions are shown in Figure 2.9a, that of the discharge response to watershed-available precipitation as well as to raw precipitation. Including the snow storage-and-release model creates higher peak discharge responses with a more phys- ically realistic long tail due to groundwater discharge. The higher peak is expected because without a snow model this method treats the precipitation the same in Jan- uary as it would in July, even though the January precipitation fell as snow and was stored until later, resulting in no significant short-term discharge response. A convolution of the solid-line unit response function in 2.9a with watershed- 43 9.. ..... E ........... ....... §.----RawPrecip.. ‘ ; : — W.A. Precrp. ‘ ) 3/s/cm ........................................ (m 8 6' (I) § Stream Discharge Response to Precipitation 0 1‘0 20 30 40 50 Days Since Precipitation Event Figure 2.9: Stream unit response hydrographs. a) Unit response functions obtained using data from 9/1/1999 to 8 / 31 / 2000, and b) Seasonal unit response functions using selected years (see Table 2.1) from 1990-1999. available precipitation according to Equation 2.1, produces the dashed modeled dis- charge and residual curves in Figures 2.10a and b. Because the unit-response was truncated as described in the methods, the convolution is not a perfect reconstruc- tion. The resultant discharge is an overestimate in the summer and fall months, but an underestimate during the spring. This is to be expected as the annually-calculated response curve effectively averages the system behavior throughout the year. If the unit response is calculated seasonally rather than annually, the resultant set of unit response curves reveals the seasonal differences among runoff responses in this watershed (Figure 2%). Discharge response during the spring is much higher than either the annual curve or those of the other three seasons, perhaps due to a combination of frozen soils and higher average soil saturation prior to watershed- available precipitation events (which can be either rainfall or snowmelt). Summer discharge response, on the other hand, is highly damped due to canopy interception, lOwer average soil moisture, and evapotranspiration. The fall and winter responses 44 Z? 80 ' - -‘ Annual . m l g —- Seasonal , I' a, . .. E’ (D .C 0 .9 0 § 200; b N} o . [E 100 ‘fll\§“ . ““15, _, T) CT WW ‘._ a A J N- w M '8 _50 v v ‘Nh VI'II : J Sep Nov Jan ar May Jul Sep 1999 «F» 2000 Figure 2.10: Plots of convolved stream discharge from 09 / 01 / 1999 to 8 / 31 / 2000 using the annually— and seasonally-derived unit hydrographs (initial modeled discharged matched to stream discharge) on top of the measured stream discharge (a), and on top of the measured stream discharge. b) Plot of the percentage difference between each model and the measured discharge (modeled-measured)/(measured)*100. appear to be very similar, suggesting that there are not large differences in discharge response between these two seasons. This result was somewhat unexpected since frozen soils are generally expected during the winter months due to long periods of sub-freezing temperatures. Significant areas of frozen soils would tend to increase discharge response, as infiltration capacity is greatly reduced. The lack of a response difference suggests that the soils are not homogeneously frozen throughout the winter season, or that this freezing is not important for runoff generation in this watershed during this time period. Analysis of data from the wells in the GTBW indicates that the delays between peak spring recharge and peak saturated water table response scales approximately as 2-4 days / meter of unsaturated zone depth. Thus, for depths on the order of 30 meters, the delay between full groundwater response to a precipitation (or snowmelt) event can be as much as 120 days. As the data lengths included in the seasonal response 45 calculations are only on the order of 60-90 days the seasonal unit response curves can not properly represent the groundwater response. The tail in the annually calculated unit response curve of Figure 2.9a is likely a more reasonable representation of the average groundwater response. In addition to providing more insight into watershed processes than the annually- derived unit response curve, the seasonally-derived curves produce a much more ac- curate estimate of stream discharge when convolved with watershed-available pre- cipitation (Figure 2.10a). The residuals between the two convolutions and stream discharge (Figure 2.10b) quantitatively demonstrate the improvement in discharge estimation gained by seasonal convolution and consideration of non-stationary system behaviors. Figure 2.10a illustrates the utility of the seasonally-derived watershed unit response curves for providing a very simple means of forecasting discharge response to precipitation events. Because the seasonal curves average watershed responses across significant variability in watershed state properties (such as soil moisture), the seasonally-derived convolution underestimates the largest peaks in the discharge data by less than 50%, while predicting measured flows to an accuracy of +/- 25% in most other cases. Much of the remaining residual is because the 50-day unit response curves fail to capture much of the groundwater response to spring and fall precipitation. An analysis that accounts for the different seasonal responses between wet and dry years, and explicitly incorporates the full groundwater response, might further improve the accuracy of forecasting with this technique. 2.6 Conclusions We present the application of three spectral analysis techniques, direct spectral com- parison, the Scaled Windowed Fourier Transform (SWFT), and the derivation of the unit hydrograph via FT deconvolution. We have discussed and demonstrated how 46 each technique requires consideration of limitations and possible pitfalls in order to be applied successfully, and we elucidated a set of best practices in their application. Most importantly, we have demonstrated that these spectral analysis methods can be used to integrate hydrologic data in order to evaluate watershed processes. The spectra of related data types were directly compared to infer process linkages between hydrologic inputs, watershed processes, and stream discharge. Similarities in amplitude peaks, log-log linear fractal scaling behaviors, and breaks in scaling slopes among datasets indicate the nature of linkages. For example, fractal scaling in precipi- tation may be matched in the stream spectrum at periods shorter than approximately 3 hours for the Evart, MI sub-basin of the Muskegon River Watershed. From peri- ods of 3 hours to approximately 10-30 days, stream scaling follows a ,0 = 2.9 slope. This single scaling relationship is notable considering the variety of processes active in this period range, suggesting mathematical similarities among these processes. Be- yond 30 days, the scaling apparent in the stream spectrum appears to be controlled by groundwater inputs. But, past a period of approximately 2.5 years, fluctuations in the precipitation spectrum control the stream discharge spectrum. This overall watershed response time should be considered when developing transient predictive simulations of watershed behavior. We introduced the Scaled Windowed Fourier Transform (SWF T) technique to ex- amine the time-varying content of fundamentally non-stationary hydrologic datasets. The SWFT scalograms revealed both the non-stationarity and intermittency of stream discharge, precipitation, and groundwater head fluctuation spectra. The effect of evapotranspiration and canopy interception is evident in a comparison of the SWFT scalogram of summer precipitation events to the highly damped discharge scalogram for those seasons. Also, the 1-day peak evident in the FT spectrum of stream discharge was shown to be due largely to measurement error rather than diurnal hydrologic pro- cesses. Importantly, these are a subset of many possible observations from the rich 47 t-. “W-w set of information contained within the SWF T scalogram. Using direct FT deconvolution, spectral analysis can be also be used to estimate stream unit hydrographs. Because of the simplicity of this method, temporally- and seasonally-varying hydrographs can be quickly derived to better understand non- stationary watershed processes. For the Muskegon River above Evart, MI, the ground- water dominance of the stream discharge spectrum beyond approximately 15-20 days is confirmed by visual inspection of the main unit response peaks. These peaks Show a 15—20 day primary stream response period followed by a long-tailed groundwater response that continues out to at least 50 days. Also, the seasonally-derived unit hydrographs quantitatively reveal decreased discharge responses due to evapotranspi- ration during the summer months and augmented responses during spring snowmelt. Acknowledgements We are grateful for the financial support from the National Science Foundation grant (EAR-0233648), and secondary support from the Great Lakes Fisheries Trust. We thank Dr. Jeffrey Andresen and Dushmantha Jayawickreme for providing and pro- cessing the NEXRAD precipitation data used in this research. We also thank the hydrogeology research group at Michigan State University for their suggestions and contributions. Finally, we thank Dr. James W. Kirchner and an anonymous reviewer for their comments. 48 CHAPTER 3 Evaluating Temporal and Spatial Variations in Recharge and Streamflow Using the Integrated Landscape Hydrology Model (ILHM) This chapter was published in 2007 in a slightly modified form in an AGU Geophysi- cal monograph, Subsurface Hydrology: Data Integration for Properties and Processes {Hyndman, Kendall, and Weltt, 2007) 3.1 Abstract Projections of climate and land use changes suggest significant alterations to the hy- drology of the Upper Midwest. Forecasting those changes at regional scales requires new modeling tools that take advantage of increased computational power and the latest GIS and remote-sensing datasets. Because of the need to resolve fine-scale pro- cesses, fully-coupled numerical simulations of regional watersheds are still prohibitive. Although semi-distributed lumped-parameter models are an alternative, they are of- ten not able to accurately forecast across a broad range of hydrologic conditions such as those associated with climate and land use changes. We have developed a loosely-coupled suite of hydrologic codes called the Integrated 49 Landscape Hydrology Model (ILHM), which combines readily-available numerical and energy- and mass-balance modeling codes with novel routines. This code is used to predict hydrologic fluxes through a 130 km2 portion of the Muskegon River Water- shed in northern-lower Michigan. We combine GIS maps of the land cover, soils, and sediments with a variety of gaged and remotely-sensed data for this watershed to simu- late evapotranspiration, groundwater recharge, and stream discharge from 1990-2004. These estimates are compared to measured stream discharge data to demonstrate the capability of the ILHM to provide reasonable predictions of groundwater recharge with minimal calibration. The results begin to illustrate critical differences in hydrologic processes due to land cover and climate variability, including a demonstration that approximately 75% of precipitation becomes recharge during leaf-off periods while almost no recharge occurs during the growing season. 3.2 Introduction Land use and climate changes are expected to alter the spatial and temporal distri- bution of groundwater recharge over the next century (Bouraoui et al., 1999; IPCC, 2001). These changes could have far reaching consequences since recharge maintains groundwater supplies that are used as primary drinking water sources, and is critical to stream ecosystem health since groundwater is the main source of streamflow during dry periods. Despite the clear importance of groundwater recharge, its spatial and temporal distribution is generally poorly understood in humid regions. Many hydro- logic modeling studies ignore both spatial and temporal variations in recharge rates, either because limited measurements of critical parameters are available, or because existing modeling methods are not adequate to accurately evaluate these variations at the scales of interest. Integration of available hydrologic and landscape data can help improve estimates of historic recharge rates, and can then provide the basis for 50 evaluating the range of impacts of anthropogenic alterations of the landscape and climate on future hydrological and ecological conditions. A range of approaches have been develOped to estimate recharge rates based on relatively simple analysis of flows and levels in surface water, the unsaturated zone, and the saturated zone, as reviewed by Scanlon et al. (2002). These methods include analyzing baseflows or tracer concentrations, developing estimates based on changes in groundwater levels (reviewed by Healy and Cook (2002)), and evaluating recharge through the unsaturated zone with lysimeters or well-instrumented field sites. A variety of empirical models have also been developed to estimate recharge across a range of scales (e.g., Bogena et al. (2005)), which can provide estimates with varying degrees of reliability and spatial extent depending on the types, quality, and density of the input data (Scanlon et al., 2002). Numerical models provide a powerful framework to integrate different data types for recharge estimation. Such models can be categorized as lumped parameter models or process-based models. Lumped-parameter semi-distributed models, such as SWAT (Arnold et al., 1993) and TOPMODEL (Beuen and Kirkby, 1979) have parameters that can be adjusted to fit measurements but can not be independently measured. As a result, such models tend to have difficulty either predicting flow in a new system without independent calibration, or projecting likely changes in a currently modeled system due to changes in factors including climate and land cover. Process-based codes such as MODF LOW for groundwater flow are based on fully distributed parameters such as hydraulic conductivity, which can be independently measured based on laboratory analyses or using field evaluations such as pump or slug tests. Unfortunately most groundwater codes are not designed to estimate recharge rates because they do not incorporate important landscape and unsaturated zone processes that are critical to redistribution of precipitation from the soil surface and the vegetation canopy. 51 To address this limitation, several codes have previously been developed to link MODFLOW or other groundwater codes to landscape or watershed codes that in- corporate aspects of the hydrologic cycle beyond groundwater flow. For example, MODFLOW has been linked with SWAT (Sophocleous and Perkins, 2000) and HSPF (Said et al., 2005), which are both lumped-parameter codes. A new Variably Satu- rated Flow (VSF) package was developed as a MODF LOW module (Thorns et al., 2006) to add unsaturated zone and overland flow processes to groundwater flow simu- lations, but the data and computational requirements appear to be too great for large watershed simulations based on our analysis of this code for the watershed presented in this paper. A variety of process-based models have also been developed to simulate fully cou- pled surface water, and variably saturated subsurface flow. Such codes, which in- clude SUTRA3D (V033 and Provost, 2002), Mike-SHE (Refsgaard and Storm, 1995), WASH123D (Yeh and Huang, 2003), MODHMS (HydroGeoLogic, 2002; Panday and Huyakorn, 2004), and InHM (Vandeeraak and Loague, 2001; VanderKuiaak, 1999), provide powerful tools to examine complex interactions between flow and transport across the range of natural conditions observed in the surface and subsurface. Unfor- tunately, the data requirements and significant computational demands have generally limited the use of these codes to simulate flows through fairly small domains. In this research, we develop a new Integrated Landscape Hydrology Model (ILHM) to integrate widely-available hydrologic and landscape data in a synergistic and com- putationally efficient manner to assess temporal and spatial changes in important hydrologic processes. Since the focus of this monograph is data integration in hydrol- ogy, we begin by describing the watershed that we chose for testing and development of the code along with the available hydrologic and landscape data used in this sim- ulation. This is followed by a detailed description of the model development and results. 52 3.3 Methods 3.3.1 Study Region: Cedar Creek Watershed The Cedar Creek Watershed, in southwestern Michigan (Figure 3.1), was chosen as a site to test the ILHM since it is one of our main field sites in an ongoing eco- hydrological monitoring and modeling study. Cedar Creek flows through the lower half of the Muskegon River watershed (7,052 km2), where urbanization of previously agricultural and forested landscapes is projected to increase runoff volumes and the associated solute transport over the next 35 years based on an empirical model ( Tang et al., 2005). The spatial distribution of land uses within the Cedar Creek Watershed facilitates evaluation of differences in recharge associated with land cover types since the upstream portion of this area is dominated by agriculture while the downstream portion is predominantly forested (Figure 2a). The quaternary geology ranges from medium and coarse-textured glacial tills that drape the northern watershed, to glacial outwash and lacustrine sand and gravel in the central and southern watershed (Figure 3.2B). The groundwater source area, which we call a groundwatershed, of Cedar Creek was delineated using a two-layer groundwater model of the region encompassing the Muskegon River Watershed (Figure 3.1). The groundwatershed (~130 km2) was used in addition to the surface watershed (""100 km2) for this study because regional model- ing of the Grand Traverse Bay Watershed in Michigan by (Boutt et al., 2001) indicated that surface-and groundwatersheds can differ significantly. The regional Muskegon River groundwater model was developed by expanding the watershed boundaries to significant hydrologic features (i.e., the next large stream or lake beyond the surface watershed) to avoid this issue at regional scales (Figure 3.1). The groundwatershed boundary does fluctuate somewhat with both seasonal and long term climatic varia- tions, but for simplicity in this study we have defined the groundwatershed using the 53 A Precipitation/Climate sites 22.: Surface Watershed EZI Groundwatershed + Flow Sampling Sites a Residential Wells .-"""" <> Stream Gages — StreamsNVetlands,g" Regional model expanded boundary Cedar Creek Muske on River Watershed g watershed _ S 01 3 5km 04080 160km Figure 3.1: On the left is a map for the Cedar Creek watershed (shaded) along with the groundwater contributing area to Cedar Creek (dashed outline). Also displayed on this map are locations of two stream gages, nine discharge measurement cross sections, and residential drinking water wells located within the watershed. On the right, a map of the lower peninsula of the state of Michigan shows the Cedar Creek watershed within the greater Muskegon River watershed. The boundary of a regional groundwater model of the Muskegon River is shown in bold on the state map. The three precipitation and climate gage locations are Hesperia (H), and Fremont (F). 54 A: Land Use - Forest 'A n'cuIture -U an/Grassland Water/Wetlands B: Quaternary C: Soil Texture Geology Sand [:3 Med-textured Tlll - Loamy Sand Sand/Gravel - Silty Loam Outwash " Loam - Coarse TIII " " Muck - Coarse Moraines Sandy Loam Dune Sand D: Elevation 270 H 175 (m) Figure 3.2: Static GIS datasets that were integrated into the ILHM framework. Ftom left to right, the datasets are land use/cover from IFMAP, quaternary geology from (Farrand and Bell, 1982), SSURGO soil textures, and land surface elevation from the NED 26.5 m DEM. 55 steady-state model. 3.3.2 Data Collection and Analysis Before constructing the groundwater model for the expanded MRW region and the ILHM to define the Cedar Creek groundwatershed, we assembled the available land- scape, hydrology, and climate data for the region into a geodatabase. In many parts of the world, the types of data used for this analysis are commonly available as a free download from internet sites. However, supplementary data such as flows and water levels beyond those available from the US Geological Survey will often need to be collected for model calibration or optimization. 3.3.2.1 Hydrologic Data Two pressure transducers installed in Cedar Creek recorded stream stage at hourly to sub-hourly intervals (Wiley and Richards, 2006) from mid to late 2002 through 2004. These surface water levels provide critical information for this study. One transducer was installed in the northern, agricultural portion of the watershed, while the other was installed in forested land near the watershed outlet (See Figure 3.1). Stream discharge measurements (Wiley and Richards, 2006) were used to construct rating curves between stage and discharge. The stage discharge relationships were developed using pairs of measured streamflow paired with water levels from the transducers. For the upper watershed site 19 measurement pairs were used, while 23 measurement pairs were used for the lower watershed site. Groundwater levels for this region were collected from the Michigan Department of Environmental Quality (MDEQ) residential well database, as no monitoring well data are available in this region except at our surface water-groundwater interaction site adjacent to Cedar Creek. Unfortunately, the wells at this site are too close to the stream to provide useful information about groundwater levels for this watershed- 56 scale model testing. Observations were available from 99 wells installed across the watershed during the simulation period in the MDEQ database. For each well, one static water level measurement was taken by the well driller at the time of installation. The static water level measurements were used in a preliminary calibration of an early version of the Cedar Creek groundwater model, but there is a significant amount of error associated with these water level measurements since a variety of methods were used by well drillers to identify the location of the well, the elevation of the ground surface, and the depth to water. Residential wells were not included in the model as extraction wells since most of the extracted water is assumed to return via septic systems and the remainder is assumed to be a very small component of the water budget for this region. There are no known irrigated agricultural areas within the Cedar Creek watershed. 3.3.2.2 Geologic, Landscape, and Remote Sensing Data We established a GIS database for the Cedar Creek region with topography, land use, hydrography, and hydrogeology characteristics. These GIS datasets were compiled from the Michigan Geographic Data Library, established by the Michigan Department of Environmental Quality. These datasets are assumed to be static for the purposes of this analysis. Land surface elevations (see Figure 3.2D), which were defined based on the National Elevation Dataset 26.5 m digital elevation model (DEM), were used for a range of model inputs. Flow direction and gradients for the overland flow and near surface soil moisture redistribution modules were calculated from the full-resolution DEM and then upscaled. This 4x upscaled DEM was used to set the land surface elevation for the soil water balance model. The drainage network and lake boundaries were de- fined based on the Michigan Framework GIS dataset, with stream crossing elevations manually extracted from a digitized version of the USGS 1:24000 topography quad- 57 tangles. The watershed does not contain any large lakes that are connected to the stream system, thus lakes are not separately considered in the groundwater portion of the integrated model. The land cover distribution across the Cedar Creek model area (Figure 3.2A) is taken from the Integrated Forest Monitoring Assessment and Prescription (IF MAP) coverage, which is a statewide digital land cover map with 30 m resolution derived from 1997—2000 LANDSAT data (MDNR, 2001). This watershed is dominated by forested/openland/wetland land covers (60%), while 36% of the area is agricultural and the remaining 4% is urban. The IFMAP coverage provides inputs that are used in calculating evapotranspiration and overland flow. Land cover types are associ- ated with transpiration estimates through a variety of terms including stomatal con- ductance, canopy height, and root depth, and with evaporation through changes in canopy interception, wind speed and interception of incoming radiation. Overland flow is also associated with land cover through Manning’s roughness coefficients. The details of these connections are included in the modeling sub-sections below. Hydrogeologic zones were parameterized according to a Quaternary Geology cover- age of Farrand and Bell (1982) for the groundwater model along with the SSURGO soils database (see Figure 3.2C), which was then mapped into saturated and un- saturated zone parameters according to lookup tables based on literature values. The geometry of the aquifer base was interpolated between measured bedrock ele- vations by de—clustered and polynomial-detrended simple kriging of the elevations of the drift / bedrock contact from oil and gas wells across the entire expanded Muskegon River model domain. Initial hydraulic conductivity values were assigned to the ge- ologic zones based on an optimization of these parameters for the nearby Grand Traverse Bay watershed that has the same geologic zones (Boutt et al., 2001). GIS grids of leaf area index (LAI), the ratio of one-sided green leaf area to ground area (Myeni et al., 2002), were also used in calculations of potential evapotranspira- 58 tion (PET), canopy interception, and solar radiation interception. For this study, we used remotely-sensed LAI measurements from NASA’S Moderate Resolution Imaging Spectroradiometer (MODIS) eight-day averaged product. Spatially-averaged LAI for both forest and agricultural land-use types are plotted for 2003 and 2004 in Figure 3.3A. In cases where these data were not available (i.e., prior to 2000), we average all LAI grids for each Julian day and apply these multi-year averages to the earlier periods. This will have little effect on our results because we are only comparing sim- ulated and observed flows for mid 2001 through 2004 when all datasets are available. However, it is important to spin up the model using realistic data inputs because we found that it takes between two and three years before the model results are independent of the starting conditions. 3.3.2.3 Climate Data Precipitation data was obtained from the NOAA gage at Hesperia, MI approximately 20 km NNW from the center of the Cedar Creek watershed (see locator map in Figure 3.1). This gage was chosen because lake effect precipitation is an important meteorological phenomenon in this area, and this gage lies at relatively the same distance from the Lake Michigan shoreline as the Cedar Creek watershed. NOAA data (shown in Figure 3.38) included hourly precipitation totals, as well as daily measurements of new snowfall and snow pack depth. Other climate data, including hourly temperature, relative humidity, wind speed, and incoming solar radiation (Figures 3.3C-F), were extracted from the Fremont, MI station of the Michigan Automated Weather Network (MAWN) (see Figure 3.1 for locations). This climate network is operated by the Michigan State University Exten- sion, the Michigan Agricultural Experiment Station and the Michigan Department of Agriculture. Since the MAWN data did not exist prior to 1996, from 1990-1995 we used the Julian-day average of the available data. 59 LAI mam are 00 co co Monthly 0 Weekly Average Windspeed (m/s) Precip. (mm) P 00:01 R.H. (%) 8 O i l I \ f.“- Weekly Average \s 1" ‘1 "31. . I ‘. h: i ‘ . ‘1... N O"! I: \A N": 0 Jan April I Figure 3.3: Time series data inputs to ILHM for 2003 and 2004: A) LAI of forest and agricultural land covers; B) monthly rain and snow along with daily average tempera- ture, a horizontal line at 0°C is included as a visual aid; C) weekly averaged windspeed and solar flux; D) inset of hourly windspeed and solar flux data from 6 / 21 / 04 through 6 / 24 / 04; E) weekly averaged relative humidity (RH) and temperature; F) inset of Jung (Sat Jan Aim: 20 3 I it. o'ct Jan 20 I hourly RH. and temperature from 6/ 21 / 04 through 6 / 24/04. 60 Weekly Average Solar Flux (kJ/m2) Weekly Average Temp (C) 3.3.3 Components of the Integrated Landscape Hydrology Model (ILHM) Figure 3.4 illustrates our conceptual model of the most important hydrologic processes in the Cedar Creek watershed, and diagrams the linkages between input datasets, ILHM modules, and model outputs. As mentioned earlier, this version of ILHM was developed by linking a novel landscape water balance model with a simple linear - delay unsaturated zone model and MODFLOW-2000 (Harbaugh et al., 2000), the most commonly used groundwater flow code. The landscape and near-surface portion of the ILHM combines several existing codes with a set of new modules, in order to speed development and to incorporate the full range of hydrologic processes. The canopy water balance model is based on equations published in Chen et al. (2005a). The surface hydrology model, including infiltration and runoff routing, are modified from the Distributed model for Runoff, Evapotranspiration, and Antecedent soil Moisture (DREAM) model by Manfreda et al. (2005). The snow pack is simulated using the UEB Snow Model by Tarboton and Luce (1996). Soil moisture accounting along with near-surface flows are handled by a set of codes we developed based on common unsaturated zone flow modeling methods. The ILHM suite calculates each term in the full water balance equation: AS=P-T—E—Pc+Tr—E:r—R (3.1) where AS is the changechange in surface soil moisture storage, P is watershed avail— able precipitation, T is transpiration, E is evaporation, R is precipitation excess runoff, and Pc is deep percolation beneath the root zone, Tr is lateral near-surface unsaturated flow called throughflow, and Er: is the exfiltration from each cell. For Equation 3.1 and the detailed equations presented in Appendix A, terms are in units of meters per unit time, unless specified otherwise. The landscape portion of our model sequentially calculates the water balance along the paths water takes as it is redis- tributed from precipitation to various subsurface 61 233 on» mo Ewuwflc xoa E9386 and 7:58 333280 in Sam—rm «$.53»... am>m COEMFE—QDm . .. _.Om £239.an— 62 and surface pathways. Incoming rainfall is first subjected to canopy interception, while snow is routed directly to the snow pack model. Next, canopy thoughfall and snowmelt are applied to the soil surface. These new inputs are then combined with any water stored in surface depressions and allowed to infiltrate into the soil. Any excess water at this point enters surface depression storage up to the available capacity. Infiltrated moisture is added to the existing surficial soil layer budget, where it can then percolate downward under the influence of hydraulic gradients. Any moisture within the first soil layer is then available for evaporation, along with any transpira- tion that may occur in any Of the biologically- active soil layers. Subsurface lateral throughflow is then calculated, which may cause moisture in down-gradient cells to exceed saturation. At this point, moisture in the lowest biologically-active soil layer may then percolate into the sediments beneath, where it becomes deep percolation. Remaining moisture in excess of saturation is exfiltrated back toward the surface where it also enters depression storage. Deep percolation is then delayed as a linear function of the thickness of the un- saturated zone, which is estimated based on a steady-state run of the regional MRW groundwater model. The delayed percolation then becomes recharge to the three-- dimensional transient groundwater flow model when it crosses the water table. Water stored in surface depressions is then subjected to direct evaporation. If depression stor— age capacity is exceeded, the excess water becomes surface runoff. Baseflow discharge from the groundwater model is then combined with the surface runoff and throughflow to produce the complete simulated stream hydrograph. The surface component of the landscape hydrology model is a 177x 153 grid at 106.3 m resolution. As shown in Figure 3.1, while the watersheds and groundwatersheds overlap for most of the modeling domain, some locations contribute only surface water or groundwater to Cedar Creek. To account for this, the landscape hydrology model is run for the entire domain while the unsaturated zone model only allows groundwater 63 Table 3.1: List of adjustable parameters in the UEB Snowmelt Model and their assumed values Parameter Units Value Snow Density kg m'3 200 Liquid Holding Capacity - 0.15 Thermally Active Soil Depth 111 1.0 Snow Thermal Conductance m hr‘1 0.2 recharge in active cells of the saturated groundwater model, and the stream routing module only includes areas within the surface-watershed of Cedar Creek. 3.3.3.1 Precipitation and Snowmelt Watershed available precipitation (P) is the sum of liquid rainfall and snowmelt. From late December through mid March, precipitation falls predominantly as snow in the Cedar Creek watershed. To model the storage and release of snow we used the UEB Snowmelt Model by Tarboton and Lace (1996), which is an explicit energy and water balance model designed to track three state variables: snow water equivalent, energy deficit (i.e., how much energy would be required to return the snow pack and soil layer to the 0 degree C reference condition), and the snow surface age. The model is computationally efficient because it assumes no temperature gradient within the snow pack and the layer of soil with which it interacts. For ease of integration the rest of the ILHM model suite, we ported the FORTRAN version of this code into MATLAB. The full UEB model requires air temperature, wind speed, relative humidity, and solar insolation. In addition it has several adjustable parameters, including the density of the snOWpack, the thermal conductance of the snow, liquid water holding capacity of the snowpack, and the depth of soil with which the snow interacts. Preliminary calibration to our snow depth data provided the parameter values given in Table 3.1. The ILHM suite only runs the UEB model if either the air temperature during a precipitation event is below freezing or the snow water equivalent of the snowpack is greater than 0.01 mm. Any water remaining in the snowpack at this point is then 64 applied to the surface as additional snowmelt. The watershed available precipitation (P) is defined, then, as the sum of all available snowmelt and liquid precipitation as determined by the UEB model. 3.3.3.2 Evaporation and Transpiration To calculate the E and T terms in the water balance equation, potential evaporation and transpiration were first calculated. All evaporation and transpiration potentials, (canopy, PEC; depression, PEd; soil, PBS; and transpiration, PT) are calculated using the modified Penman-Monteith equation (Monteith, 1965) presented by Chen et al. (2005a) and shown as Equation A.1 in Appendix A. Evapotransiration and transpiration rates vary according to land cover types through a stomatal conductance map which is assumed to be temporally invarient as discussed in Chen et al. (2005a) and with 8-day LAI scenes. Incoming rainfall is first subjected to interception up to the water holding capacity of the canopy, which is related to LAI in the cell. Evaporation (E) terms are cal- culated for the canopy, soil, and any surface depressions within each cell. The total transpiration in each cell T is calculated to be the sum of the root water uptake from each biologically active soil layer. The evaporation of moisture in depression storage is calculated after any infiltra- tion and exfiltration (described below) in a given time step have occurred. Total depression storage capacity is determined by land use, soil texture, and slope class; a tabular reference of storage capacities can be found in Manfreda et al. (2005). De- pression evaporation occurs at the potential rate until depression storage is depleted. Evaporation from depressions and directly from soil is allowed only from the propor- tion of soil that is exposed to solar radiation, thus assuming no soil evaporation from the portion shaded by canopy or covered with snow or ice. Direct soil evaporation is allowed only from the first soil layer, which is a reasonable 65 assumption in this relatively humid region. Soil evaporation occurs at the lesser of calculated potential rate or the soil exfiltration depth (discussed further in Appendix A following Equation A16). We are exploring alternative strategies for calculating evaporation using the model of Ritchie (1972). The total transpiration in each cell T is calculated to be the sum of the root water uptake from each biologically active soil layer. We assume that transpiration only occurs above a dormant threshold temper- ature, which was chosen to be 40 degrees F for this study. Stomatal conductance is also assumed to be constant for this case, although we plan to incorporate variations due to changes in temperature, carbon dioxide concentrations, and soil moisture as described in Chen et al. (2005a). 3.3.3.3 Infiltration, Percolation, Throughflow, and Exfiltration For this study we assume the biologically active soil can be described by two layers, with a total thickness calculated according to Equation A14 as the depth above which 90% of the root mass lies. The first soil layer, from which evaporation occurs and that controls infiltration capacity, is on the order of several centimeters thick. Infiltration capacity is calculated as the greater of either the soil-texture dependent saturated infiltration capacity, i sat; or the first layer moisture deficit from saturation. We chose the maximum infiltration rate to be (2 X i sat) which determines the choice of the first soil layer thickness. This formulation produces similar results to empirical infiltration rate descriptions, and has the advantage that it does not require storm event tracking or single storm event modeling. Here we assume that isat can take either its nominal soil-texture dependent value taken from literature values (see Appendix A), or i sat = 0 if the soil is frozen, which we only allow to occur in agricultural soils based on Schaetzl and Tomczak (2001). The soil is assumed to be frozen if its temperature is below —0.25°C, measured by the MAWN station in Fremont at a depth of 10 cm. 66 In addition, we assume that no infiltration occurs in cells classified as permanent water features. Although the water features in this watershed often cover only a small portion of each land-use cell, our assumption may nevertheless be fairly realistic. The true physical system process is more accurately described as rapid percolation to a shallow water table followed by equally rapid rise and subsequent relaxation of the water table. The effect is a temporary increase in groundwater discharge that generally does not modify what is more traditionally called stream baseflow. Our stream gage data indicate that the characteristic time of this response may be on the order of twice the surface runoff response, or perhaps 5—7 days. This does not allow for significant losses due to evapotranspiration, thus the combined increase in stream dis- charge due to percolation to the near-stream water table and direct overland flow would be nearly equal to that expected from an assumption of zero infiltration. Streamflows in our model would thus be expected to peak higher and return to baseflow levels more rapidly than observed. T hroughf low, defined here as lateral subsurface flow within the biologically active soil zone is calculated using a simplified Richards equation model. For a full devel- opment of the Richards equation, see (Hillel, 1980). For purposes of computational efficiency, and in order to assure that the subsurface redistribution of moisture occurs in only one dimension, we assume that flow only occurs parallel to the dip of the slope and thus cannot flow uphill on the ~100 meter scale of our model cells. For environments where these assumptions may be invalid, alternate two-dimensional for- mulations could be substituted for this ILHM module given adequate computational resources. The van Genuchten model was used to calculate all soil-moisture depen- dent properties (van Genuchten, 1980) with parameter values given in Appendix A. The downgradient cell is determined via the D8 flow direction function (ESRI, 2003). However, each cell can have more than one upgradient cell, thus throughflow is the summation of the shallow subsurface flow out of all upgradient cells. 67 3.3.3.4 Unsaturated Zone Delay Deep percolation beneath the biologically-active soil layer is delayed prior to becoming recharged as a linear function of the depth to the water table. The slope of this delay in units of days/meter, was determined from wells installed in the nearby Grand Traverse Bay Watershed, and was fixed at 2.5 for this study. It is important to note that this would not be the same as a solute transport time through the unsaturated zone. Delayed deep percolation then becomes recharge once it reaches the water table. The depth of the water table is assumed temporally invariant for the current version of the unsaturated zone delay module for ease of implementation. By fixing the depth of the water table for this purpose, we do not account for seasonal or trending differences in water travel times through the unsaturated zone. The seasonal differences in the depth of the water table are small (typically <1 meter) relative to average depths to water over most of the model domain. Locations very close to surface water features are an exception, but these areas comprise a small fraction of the total watershed area and thus will not significantly affect the dynamics of modeled recharge. 3.3.3.5 Groundwater Model The groundwater model for the expanded MRW was developed using a suite of MAT- LAB utilities that we developed to create input files from GIS layers for MODF LOW. A regular grid of 1798 x 1865 cells (106.3 meters on a side) was used so that each cell directly overlies 16 Digital Elevation Model (DEM) cells. The vertical domain of the saturated zone model was then subdivided into two layers with approximately equal saturated thicknesses based on the simulated water table in a single layer model. Au- tomated parameter estimation routines were applied to an early version of the Cedar Creek groundwater and soil balance model to estimate hydraulic conductivity values for aquifer sediments in geologic zones parameterized using a digital map created from Farrand and Bell (1982). 68 The Cedar Creek model is a single layer with dimensions of 184 x 162 and cell size of 100 meters. In this case, a single vertical layer provided an adequate description of flows through the Cedar Creek watershed, because it does not have significant vertical relief, extensive low-permeability subsurface layers, high-capacity pumping wells, or other features that tend to induce significant vertical head gradients. Groundwater discharge to streams is calculated with the Stream Flow Routing (SF R) package in MODFLOW that routes water via the kinematic wave equation (Prudic et al., 2004). 3.3.3.6 Runoff and Stream Routing Surface runoff is routed to the streams using an approach modified from that presented by Manfreda et al. (2005). In this version of the code, we assume that runoff cannot reinfiltrate once it is generated. It is routed overland and through streams according to the D8 flowdirection algorithm in ARC (ESRI, 2003) with runoff times given by the velocities in each cell along the flowpath. Runoff is assumed to travel overland at a velocity given by the Kerby time of concentration equation (K erby, 1959). Once the runoff enters the stream channel its velocity is calculated using Manning’s Equation. For this study, the hydraulic radius, r, for Manning’s Equation was determined as a function of discharge using low-flow channel geometry measurements in and around the Cedar Creek watershed along with geometries reported by the USGS for their stream gages. Wetted perimeter was assumed equal to 2 x depth + width, while area was simply depth X width. A power law fit to these data produced the empirical relationship for this watershed: r = 0.9046 - Q0283 (3.2) with a correlation coefficient R2 of 0.77 (see Figure 3.5). Q in the above equation is the measured stream discharge in m3/s. While streamfiow velocity can be dynamically calculated, for simplicity we have assumed a temporally constant vstream for each stream cell. To calculate Ustream, 69 E4. 0 (I) o .2 33' o e (I g 32’ 'C > I1, 10'4 10'2 100 1o2 Discharge (m3/s) Figure 3.5: Plot of measured discharge versus hydraulic radius in streams in the greater Muskegon River Watershed and adjacent watersheds. The hydraulic radius calculation assumes rectangular cross-section geometry. The power-law fit to the data produced a correlation coefficient, R2, of 0.77. we assume that each cell in the Cedar Creek watershed contributes a unit of runoff which is them routed to the gages using ARC’S fiowaccumulation function. To rescale the output to match a typical discharge event in a stream cell (Qij), we multiply the measured Q at the outlet by the ratio of the flowaccumulation value in each cell (i, j) to that of the outlet: f lowaccumulationi j ) (3 3) Qij = Qoutlet ( f lowaccumul ati onoutlet from each cell to the outlet. Here the weighting is the inverse of the velocity in seconds/meter. ARC then multiplies this “cost” by the distance traveled through each cell along the entire flowpath and outputs to the travel time to the outlet. The travel time grid is then used to transform the precipitation excess grid into a runoff hydrograph. 3.4 Results The results of the ILHM simulation for the Cedar Creek watershed are discussed in the context of the broader goals of the code, which are the prediction of temporal and 70 TU _ Table 3.2: List of calibrated parameters for the unsaturated and saturated ground- water models. Parameter Units Value Outwash conductivity m d'1 11.0 Till conductivity m d'1 4.4 Unsaturated zone delay d rn'1 2.5 spatial variations in recharge with very little direct calibration of model parameters using readily available remote-sensed and ground-based data sources. All parameters for this prediction were based on literature values (Tables 3.2 and A1—A3), except for the UEB model parameters (Table 3.1) as well as two hydraulic conductivity values and one unsaturated zone delay parameter that were calibrated using an early version of the model (Table 3.2). A plot of simulated versus observed heads (Figure 3.6) across this region shows a reasonable degree of agreement given the measurement uncertainty. Figure 3.6 shows no trending bias between simulated and observed heads, though a slight high-side bias is present at observed heads lower than approximately 200 meters. We would expect a higher degree of correlation between simulated and observed water levels in regions where pressure transducer data are available from wells, or if a parameter optimization were to be performed for this ILHM simulation. Detailed evaluation of modeled flows is hampered by the lack of long duration stream gages with stable channels in the basin. All observed stream discharges were collected via established methods, however the rating curves for the two stream gages are currently inadequate to account for temporal adjustments of the channel geometry after flood events. As is commonly the case, we have few high flow measurements, which limits the accuracy of our flows calculated from the rating curve during large floods. In addition, flows in the Lower Cedar gage from January through March appear to suffer from ice-induced over-pressurization not observed at the Upper Cedar gage. 71 230 . A 0’ ggzzo- . . g t ”210’ . £3: a” I .‘u 8 . 2 "j. 3.2.200 :, '. E 180 150 190 260 210 220 230 Observed Head (m) Figure 3.6: Uncalibrated simulated vs. observed groundwater heads plotted on top of a 1:1 line. Observations are from 96 residential wells installed in the Cedar Creek Watershed between 1990 and 2004 (see locations in Figure 3.1). Despite almost no calibration of parameters in the near-surface components of the ILHM code, the model provided a reasonable prediction of observed flows (Figure 3.7) for the two available gage sites in this 100 square kilometer watershed (Figure 3.1) during the fall and winter months. The ILHM also provided reasonable predictions of baseflow for this watershed system during the entire year. Because this prediction is based almost entirely on a set of widely avail- able meteorological inputs and GIS datasets combined with literature parameter values, the code appears suitable for directly simulating streamflows in ungaged basins. The close agreement between observed and simulated basefiow levels also suggests that the model is providing reasonable predictions of recharge, which provides base- flow in these streams during low flow periods. During May and October 2003, and January—March 2004, the ILHM- simulated total flows typically agreed with gaged values within 10% (when the lower Cedar gage was not affected by ice cover). Base- flow levels in the smaller upper Cedar catchment proved highly sensitive to hydraulic conductivity in the outwash sand / gravel zone (Figure 3.2A). The conductivity values presented in Table 3.2 should not be viewed as a fully calibrated parameter set, as optimization of the total stream-flow simulation is the subject of ongoing evaluation. 72 U'IO Discharge (m3/s) 0) A Discharge (m3/s) ll , if litiiaiutlk,‘ 0 May Sep Jan May Sep Jan May Sep Jan May Sep | 2001 | 2002 l 2003 | 2004 | Figure 3.7: A & C) Upper and B 81. D) Lower Cedar Creek stream discharges with simulated values shown in white and gage values in black calculated based on stage discharge relationships. Manually measured discharge values are shown as black cir- cles. 73 3 r: v ." fill. Y“ 14 ii ‘ Despite the limitations imposed by some parts the stream gage data, the values from April through December of 2003 are reasonable for quantitative flow compar- isons. During lower ET periods from April through early June and October through December, the ILHM-simulated total f luxes are approximately 20% higher than cal- culated from the flow gages installed at both the Upper and Lower Cedar gage sites. ILHM-simulated total fluxes during higher ET periods from June through October are much greater than observed, also likely due to an incomplete description of ET processes in this version of the ILHM code that is the subject of ongoing development. A scatter plot of simulated versus observed flows at the two gages on a log-log scale illustrates that most of the moderate to high flows in the system are reasonably described by the model (Figure 3.8). There is a larger degree of mismatch in the Up- per Cedar Creek site, as illustrated by a significant amount of scatter about the 1:1 line. Simulated peak discharge values are similar to those that have been measured, however the simulated discharge peaks are narrower than observed. This temporal offset at near-peak discharge, seen in Figure 3.7B, appears in Figure 3.8 as a tendency toward low simulated flows relative to observed values. These narrow simulated peaks are largely related to the simple nature of our stream routing package, and the unidi- rectional linkage between groundwater and surface water processes in this version of the ILHM code. The assumption of a temporally constant vstream results in average flow velocities lower than those that would be expected, with the effect that simulated peaks would be even sharper if vstream were a function of discharge. However, there are several known flow damping mechanisms, including f low through wetlands and bank exchange, which are not yet represented in this version of the code. A map of the simulated average annual recharge across the watershed implies that agricultural areas may have higher recharge than forested areas according to this simulation (Figure 3.9). Several inter—related factors combine to account for the simulated differences. Agricultural areas experience less canopy interception than 74 034’ 0 Upper Gage , E a Lower Gage . V o gmor ° 0 2...’ . . q ‘5 o .C D o -1 o E10 - 0 ° £9 0 E O o (9 O m - 10 - ‘ ‘ . 1 10‘2 10'1 10° Observed Discharge (m3/s) Figure 3.8: Uncalibrated simulated vs. observed flows for the Upper and Lower Cedar Creek gages plotted on top of a 1:1 line from approximately 20 manual discharge measurements taken at each of the upper and lower Cedar Creek gage locations. forested areas due to lower LAI values. Although this increases the infiltration into agricultural soils, a resulting increase in transpiration tends to narrow the difference in recharge between these land use types. Our model may also under-represent soil evaporation in agricultural soils because it does not incorporate solar heating of the shallow soil layer. This simulated difference is the subject of future evaluation across the much larger MRW where more flow data are available. The blocky nature of the simulated recharge in this map is mainly due to the large (1 km2) LAI cells. The effect of these coarse cells is to decrease forest LAI and increase agricultural LAI in regions with mixed land-uses. We plan to resolve this issue through downscaling the LAI information by assuming that the measured LAI value is a linear combination of the LAI of forest and agricultural land-uses represented by the much higher resolution IFMAP dataset. Thus unique “agricultural LAI” and “forested LAI” values can be approximated at the resolution of the IFMAP data constrained by the total measured LAI from MODIS. Areas with low hydraulic conductivity soils experience reduced recharge (Figure 3.9), such as portions of the upper watershed to the south side of the stream where loams and silty loams are common (compare with the soils map in Figure 3.2C). 75 Recharge (cm/yr) - 0 Figure 3.9: ILHM—simulated average annual recharge for the Cedar Creek watershed. Annual total precipitation averaged approximately 83 cm during the period of simu- lation. In contrast, recharge can be greatly enhanced in internally drained regions. In this simulation we deactivated the runoff mechanism in internally drained areas, thus potentially increasing infiltration and shallow subsurface flow. This may be very important in areas with moderate to low-conductivity soils, but the internally drained areas in this watershed tended to also be sandy so the effect is only localized. The Cedar Creek region experiences very little runoff from upland areas, which is consistent with the simulated map of precipitation excess (Figure 3.10). Nearly all the simulated precipitation excess in this watershed occurs in cells that are classified as “water” because there is no transpiration or percolation from those cells. The only cells with any significant precipitation excess that are not classified as water are in a region of lower conductivity sediments in the upper watershed. Despite a simple description of runoff processes, we do not expect significant runoff in most of the sediments across this watershed due to high infiltration capacities and saturated hydraulic conductivities. There is also very little simulated subsurface redistribution, 76 Precipitation Excess (cm I yr) -O-O.1 Figure 3.10: ILHM-simulated average annual precipitation excess runoff for the Cedar Creek watershed. Annual total precipitation averaged approximately 83 cm during the period of simulation. or throughflow, throughout most of the watershed due to the highly conductive soils and relatively gentle topography. This process is most active in areas with steep slopes or water tables very near the surface. The simulation provides evidence for a strong seasonality in recharge rates for the Cedar Creek watershed. The temporal variations in simulated deep percolation are shown in Figure 3.11. From September through March in the four illustrated years, the model predicts that approximately 70~80% of watershed available precipitation will percolate into the deep aquifer sediments where it eventually recharges groundwa- ter. In contrast, the simulations show virtually no deep percolation over the growing season from May through September for the same years, which is consistent with the statistical findings of Jayawickreme and Hyndman (2007). This simulation indicates that agricultural areas have more recharge in the fall months than forested areas, while the opposite occurs with higher relative forest recharge in the spring months. This is reasonable as forests tend to have less extensive frozen soils during snowmelt 77 I I I I I I I 100 150 A c a", -.= .-. %A100 '3 s E i 633?, 3- Q. L 3% 5°: 32 [l o .a... 30 I l I T j I I E 20 . - 2‘5 ‘° nil H ‘ L2 5' $5 0 “,2 U' W U — [:1 ll U... as .0 U” “U [j [j a) ' ' j f o_ . .20 1 l 1 1 l l l 1 Jan May Sep Jan May Sep Jan May Sep Jan May Sep Jan | 2001 | 2002 | 2003 | 2004 | Figure 3.11: A) Monthly average deep percolation (bars) in the Cedar Creek wa- tershed plotted with the ratio of percolation to total precipitation for each month (stars). B) Monthly difference in deep percolation between forested and agricultural land covers. periods Schaetzl and Tomczak (2001), and agricultural LAI often begins to decline earlier in the year than in forested areas. Although coniferous forests may transpire year round, they represent only a small percentage of the forested areas in this study region. Temporal variation in evaporation and transpiration are clearly the causes of most of the simulated variations in deep percolation because these are the primary loss pro- cesses. As Figure 3.12 illustrates, evaporation is generally a much smaller component of water loss in this watershed than transpiration, and this component is larger in forested land relative to agricultural land due to much higher forest canopy intercep- tion. Transpiration shows a stronger seasonal trend than evaporation, as it depends more strongly on LAI. Total agricultural transpiration is greater than that of forested 78 areas despite much greater potential transpiration in forested areas. Agricultural ar- eas experience less canopy interception than forests, and thus greater infiltration and higher average soil moisture. As a result, agricultural areas tend to transpire closer to their potential rate than forested areas. Unexpectedly, agricultural transpiration also rises in the spring more quickly than that of forested areas according to these model results due to the similarity in LAI values during early spring and higher stomatal conductance values for agricultural areas relative to forests. As the LAI of forests increases in the late spring, the transpiration in these areas becomes larger than that of agricultural areas, until they reach approximate equality in late June that continues through the rest of the summer. Also during the summer, soil moisture levels reach their lowest point and often approach the permanent wilting point. As deep percolation cannot occur until the field capacity of the soils is reached, most of the water that does infiltrate the soil is transpired. Thus deep percolation is almost non-existent during summer months according to these simulations. 3.5 Discussion and Conclusions We present the development and testing of a new suite of loosely coupled process-- based codes that we call the Integrated Landscape Hydrology Model (ILHM). This modeling framework has several advantages over existing coupled hydrology codes. It can simulate much larger domains than fully coupled process-based codes, with fewer data requirements. In addition, the ILHM accounts for the processes and mass balance in a more rigorous manner than semi-distributed codes, which tend to lump or oversimplify important watershed processes and use parameters that cannot be independently measured. The ILHM also facilitates model development via direct input of readily available GIS data, in contrast to the impractical level of manual data input required for large domains from some existing process-based models such 79 E (mm) T (mm) Difference, Forest-Agriculture (mm) E+T (mm) 0 3 - 7‘ Jan May Sep Jan May Sep Jan May Sep Jan May Sep Jan | 2001 | 2002 | 2003 l 2004 | Figure 3.12: Average monthly evaporation (E), transpiration (T), and evapotranspi- ration (E+T) plotted for both forested and agricultural landscape as stacked bars along with the difference (Forest-Agriculture) between the two land-use types plotted as circles. 80 as Mike-SHE. Finally, ILHM is well-suited for forecasting purposes because it allows forcing data and component process models to be interchangeable; thus a model developed and calibrated with current data can be rapidly converted to a forecast simulation by adding the appropriate component process code. This new modeling framework was designed to make development of models for large domains as simple as possible, while maintaining a rigorous fluid mass balance based on the primary processes that drive water movement over the landscape and through the subsurface. The approach is computationally efficient because it allows some processes to be simulated based on full numerical models while others can be described by simpler and thus faster water- and energy-balance approaches. Due to the loose—coupling framework, individual components can also be simulated at a variety of spatial and temporal scales appropriate to the individual processes. This framework also allows more rigorous simulation modules to be used in place of a simpler routine in cases where the additional computational burden provides necessary improvement in the model predictions. Alternatively, in cases where enough data exist to adequately describe a particular process, the data can be used in lieu of that process simulation module. As currently configured, ILHM is designed to simulate flows through regions with connected surface water and groundwater regimes such as Cedar Creek. This test watershed has a sub-humid and temperate climate, with flow through a glacio-fluvial aquifer, largely covered with deciduous forests, agricultural land, and small percentage of urban cover. Thus the code is expected to provide reasonable predictions for similar environments in the sub-humid Midwest. The general processes are the same in arid and montane regions; however alternate modules would likely provide more accurate simulations in such cases. In particular, some high-relief environments may require a full two—dimensional representation of overland flow, including depth-dependent velocities for sheet or rill flow. Areas with large proportions of urban land uses will 81 require additional modifications, especially when engineered storm water systems have a significant effect on the hydrograph shape after a storm event. In the Cedar Creek watershed, precipitation excess runoff routing and subsurface moisture redistribution are both largely inactive over most of the modeled domain and timeframe. As a result, the simulation results are similar even if these modules are not active for upland areas. Therefore, further testing in domains where these processes are responsible for a significantly larger percentage of the flow in a river system will be needed for these modules. The current unsaturated zone module is a very simple representation of hydrologic processes, thus we will explore the use of direct solution methods ranging from the Green-Ampt model through the full Richards equations. Additionally, using MODF LOW or any finite difference scheme has the disadvantage of requiring somewhat cumbersome rectangular grids that limit cell refinement at regional scales. However, ILHM can easily be altered to interface with a finite element code capable of representing and accounting for groundwater discharges to streams. The ILHM was tested in the Cedar Creek watershed because of the need for high- resolution flow simulations that provide the interface between land use change models and ecohydrology models in the near future. This first evaluation of the ILHM mod- eling framework demonstrated that these codes can reasonably predict groundwater recharge and streamf low through a 130 square kilometer watershed with very lit- tle calibration using readily available data. The simulation represented overall basin recharge accurately, but it appears to have slightly overestimated recharge in agricul- tural areas. This is likely due to an inadequate representation of soil evaporation that will be addressed in a future version of the ILHM. The simulated hydrograph peaks are too narrow and decline more rapidly than is observed because the current surface water/ groundwater linkage cannot represent bank storage and release processes, nor can the unidirectional coupling between surface water and groundwater fully repre— 82 sent near-stream processes at our chosen spatial scales. Nevertheless, the recharge and streamf low predictions provide reasonable descriptions of system behavior and will be further refined in future versions of the ILHM applied to much larger domains. Acknowledgements We would like to thank Mike Wiley and Paul Richards for supplying us with sur- face water flow measurements and stream levels from transducers, and Dushmantha Jayawickreme and Cheryl Kendall for their help with data processing and analy- sis. We are also grateful for the financial support of National Science Foundation grant (EAR-0233648), with supplemental support from both the Great Lakes Fish- eries Trust’s Muskegon River Initiative and Environmental Protection Agency grants (R830884). Any opinions, findings and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of the National Science Foundation or the Environmental Protection Agency. 83 CHAPTER 4 ILHM Development, Validation, and Simulation of Regional-Scale Stream Discharge 4.1 Introduction Climate and land use change are altering the terrestrial hydrologic cycle at all scales from hectare to continent . Over the last century, runoff from local to continental scales has increased in many areas of the world(Probst and Tardy, 1987; Labat et al., 2004), but the reasons for this are poorly understood (Gedney et al., 2006). Changes at the regional basin scale emerge from smaller-scale fluctuations within both me- teorological inputs and the catchments of higher-order tributaries that respond het- erogenously to underlying drivers. Those catchments themselves are collections of hectare-scale plots, each of which have a unique history of land use and management. Even at this scale, more localized alterations such as a re-planted hillslope or ripar- ian buffer strip on the order of several meters in width may significantly impact the quantity and quality of overland flow reaching a stream. Ultimately, predictive under- standing of the impacts of climate and land use change on the terrestrial hydrologic cycle requires models capable of explicitly simulating this entire range of scales. Historically, most hydrologic models have been developed for application across a 84 limited range of scales. Codes suited for regional to contential scales typically omit or simplify groundwater and include lumped parameters to explain much of the smaller- scale variability within model discretization units. The predictive capability of these models is limited outside of their calibrated set of conditions, as climate and land use change modify shallow groundwater systems in subtle and complex ways, and alter the values of lumped model parameters. Fully-coupled 3D numerical models are capable of predicting small catchment-scale fluxes but are generally either too computationally intensive, or too logistically cumbersome, to apply to large domains (Markstrom et al., 2008). In addition to the models discussed below, a number of models have appeared once or at most a few times in the literature (Manfreda et al., 2005; Quemer, 1997; Wigmosta et al., 1994; Chie'w et al., 1992). Originating out of a need to describe mass, energy, and momentum fluxes at the land surface withing general circulation models, land surface models have since evolved in complexity to describe much of the terrestrial hydrologic cycle. These models include, notably, BATS (Dickinson et al., 1993), LSM (Bonan, 1996), CLM (Yongjiu et al., 2001), LEAF (Walko et al., 2000), and IBIS (Foley et al., 1996). All share a common design goal: to accurately describe land surface fluxes within climate model cells on the order of 1Q x 1°, far too coarse for most hydrologic applications( Yu et al., 2006). Furthermore, the models either ignored or vastly simplified most hy- drologic processes including runoff generation, routing, and groundwater flow. Some incorporate an approach similar to TOPMODEL (Beven and K irkby, 1979), though at scales that preclude direct hydrograph simulation in all but the largest basins. To correct some of these shortcomings, VIC (Liang et al., 1994) was written primarily as a surface hydrologic model, though it, also lacks a comprehensive treatment of groundwater. A more recent efforts at coupling a detailed hydrologic models within a global climate model is detailed in Yu et al. (2006). Watershed hydrologists have developed a separate Class of models that describe the mass and energy fluxes within catchments at a variety of scales. These watershed models, including the Stanford model (Crawford and Linsley, 1966), SWAT (Neitsch et al., 2002), and HEC-HMS (USA CE-HEC, 1998), have evolved to suit the needs of soil conservationists and agricultural scientists studying water quality and quantity within catchments. These are semi-distributed, or lumped parameter models, whose parameter values are distributed within sub-aggregation units, sometimes referred to as hydrologic response units, or HRUS. Rainfall-runoff models, including HSPF (Johanson et al., 1984) and PRMS (Leavesley et al.), are similar towatershed models though are typically more finely discretized at the surface. These have arisen from the need to forecast floods within gaged basins. Both watershed and rainfall-runoff models all share a similar treatment of subsurface flows, which are routed through typically two (or more) subsurface pathways that are equated to near-surface throughflow and deeper saturated groundwater flow. These pathways are simple linear reservoirs, and must be calibrated for most applications. Agronomists and soil scientists required the capability to simulate the plot-scale, and evaluate the impacts local impacts such as tile drainage or tillage practices. They have developed models that simulate the full mass and energy balance at the surface. Examples of plot-scale models include HYDRUS-ID (Simunek et al., 2005), SHAW (Flerchinger, 2000), and RZWQM (Ahuja et al., 1999), which are fully-distributed, process-based (typically) 1D numerical models. These also include capabilities to assess nutrient fluxes through various storage reservoirs in the subsurface. Funda- mentally, however, these are limited for application to broader domains because of their construction as 1D models, and their lack of groundwater components. There are a number of models that include fully coupled surface and subsurface flow equations, with landscape process models of varying degrees of complexity. These include InHM (Vandeeraak and Loague, 2001), PIHM (Qu and Duffy, 2007), Hy- droGeoSphere (Therrien et al., 2004), MODHMS (HydroGeoLogic, 2002), GSSHA 86 (Downer and Ogden, 2003), and CLM-PARF LOW (Maxwell and Miller, 2005). Ap- plications of these models to small to mid-sized (1 - 100 km2) catchments demonstrate the value and efficiency of simultaneously solving fully-coupled systems of equations. However, upscaling these models can prohibitive, or require coarsening of the dis- cretization such that the value of the coupled system of equations approach is unclear. A more recent approach has been to couple 1D land surface models with full 3D subsurface models. This method scales efficiently to regional domains while still rep- resenting the dominant hydrologic processes at those scales . The recently-released GSFLOW is a coupling of the rainfall-runofl model PRMS with MODFLOW (Mark- strom et al., 2008), a 3D saturated groundwater model. MODFLOW has also been the subject of other integrated modeling efforts including HYDRUS-lD (560 et al., 2007), the VSF package for MODF LOW (Thoms et al., 2006), and HSPF (Said et al., 2005). These models incorporate land surface components of varying complexity and coupling schemes, though all are fully-coupled. MIKE-SHE is another fully integrated model, though its cost is considerable (Refsgaard and Storm, 1995). We have developed the Integrated Landscape Hydrology Model (ILHM) to simu- late fine-resolution hydrologic fluxes across regional to continental model domains. ILHM explicitly incorporates process modules that describe nearly the entire terres- trial hydrologic cycle, including groundwater. Its parameters are both physical and measurable, though generally available from literature. These attributes make it well suited for highly-calibrated applications within a single catchment or collection of catchments, modeling ungaged basins, and hindcast or forecast simulations. Here we present an application to a regional (7400 km2) catchment in Michigan, USA. This study demonstrates that ILHM, is capable of predicting hydrologic fluxes over a large domain at unprecedented resolution with minimal calibration. 87 4.2 The Integrated Landscape Hydrology Model ILHM is composed of four loosely coupled “domain” modules: 1) surface water rout- ing, 2) canopy and root zone moisture, 3) deep unsaturated zone, and 4) saturated groundwater (see Figure 4.1). The individual modules were selected based on the criteria that each includes: 1) fully—distributed physical representations of the most important hydrologic fluxes, 2) measurable parameters whose values are generally accessible through literature, 3) computationally efficient, highly scalable algorithms, from desktops to high performance computing (HPC) environments. Also, while not a formal criteria, existing open-source or published algorithms and solvers were chosen where available to reduce development time. Evaluating possible modules to include within the context of the study area that would be simulated resulted in the following four domain modules: 1) a simple, de—coupled kinematic wave surface routing model, 2) a comprehensive 1D / 2D water- balance canopy, root-zone, and wetland hydrology module, 3) a simple l-D unsatu- rated zone pulse delay module, and 4) MODFLOW, a widely-used 3D saturated zone groundwater model (Harbaugh et al., 2000). To facilitate rapid development while assuring platform independence, ILHM was written in the MATLAB computational environment (Mathworks, 2008). Also, a high-performance disk I/ O back end, HDF5 (HDF, 2008), was chosen for its speed, scalability, and portability. An alternate study area or set of questions to evaluate would have likely produced a different set of domain modules. Therefore, the code was written to incorporate alternate modules without modifying any existing code, and entirely new processes by modifying only a few lines. The choice of MATLAB greatly facilitated development due to its fully interactive development environment, however this imposes certain performance penalties. The structure of the code itself is highly modular (though not 88 explicitly object-oriented), allowing C/C++/FORTRAN/CUDA routines to replace computationally-intensive components. Furthermore, while the simulations presented in this study were all run in single-threaded desktop environments, the code may be readily extended and modified for HPC environments. The surface routing, canopy and root zone, and unsaturated zone modules are all “embarrasingly parallel”, and will scale to many nodes with little optimization. Figure 4.1 illustrates the components, dimensionality, and spatial discretization schemes employed by the four domain modules in ILHM. The dashed arrows repre- sent flux directions or dimensions that are partially supported. The root zone soil moisture module does allow for lateral cell to cell fluxes in the down SIOpe direc- tion. Additionally, vertical fluxes in the root zone are calculated if using the full 1D Richard’s Equation module, though that was not chosen for this study. Note that the root zone and saturated zone can consist of multiple layers, while the deep unsat- urated zone currently a single layer (though tracks multiple wetting fronts). Within ILHM, each domain can be independently discretized, though currently only struc- tured grids are supported. The root zone moisture module allows for a second level of discretization, allowing each cell to be subdivided into upland, wetland, and stream fractions. There are three general types of model input requirements: 1) static inputs, 2) dynamic inputs, and 3) observations data. Static inputs describe the model domain boundaries, thicknesses, layering, model stresses, and parameters. Dynamic inputs include climate and landscape data such as precipitation, solar radiation, relative hu- midity (or dew point), wind speed, air temperature, and leaf area index (LAI). Other dynamic inputs might be parameters whose values change throughout the simulation period, land use type, or anthropogenic stresses such as pumping wells or reservoirs. ILHM also accepts observation data for any quantity that it simulates, which may be used to implement sensitivity analysis or parameter estimation schemes, though 89 Hydrologic Flux Zone Dimensions SURFACE ROUTING A I noor ZONE 4...,3-5.) .71 DEEP UNSATURATED ZONE y SATURATED ZONE x 2 H Full Support 4..---.) Partial Support UPLANI WON] Figure 4.1: Graphic depiction of the spatial discretization scheme used in ILHM. these are not included at this time. 4.2.1 Model Process Flow and Coupling For this study, the canopy, root zone, and wetland hydrology modules run indepen- dently of the unsaturated, saturated, and surface routing modules, thus they are coupled serially rather than at run-time. Within each separate component, indepen- dent spatial and temporal discretization schemes may be employed. Currently, the surface processes run with hourly timesteps, while the unsaturated and saturated zone modules run on daily timesteps with weekly input stress periods. Each module writes outputs to disk, which are read by subsequent modules upon initialization. Figure 4.2 illustrates the flow of processes within ILHM. In general the processes follow the path of water, interacting with the canopy, infiltrating or running off, being redistributed in the shallow root zone or transpired, percolating into the deep 90 unsaturated zone, recharging the saturated groundwater, and finally discharging to surface water and being routed downstream. Each hour, ILHM reads dynamic inputs (climate, LAI, parameters) from disk, using a data buffer to minimize disk reads. The input file format separates the input data from how that data is distributed within the model. This allows a single input climate file, for instance, to be used for any number of parallel threads, or for smaller sub-domain models useful for testing and development. Thus the input controller performs three steps, it: 1) reads the data, 2) distributes those data to the model grid, or possibly just a subset of the grid within the current thread, and 3) stores the data in a central data controller that process modules may then request data from. The canopy components are run next. First, the model updates albedo values across the landscape and determines the shortwave radiation intercepted by the canopy, and upland or wetland portions of each cell. Next, solid and liquid precipitation are in- tercepted and stored as canopy water. Finally, canopy evaporation and transpiration are calculated. For the upland fraction of each cell, evaporation terms are then calculated for water stored in surface depressions and from the soil. Any canopy throughfall precipitation is first sent to the snowpack module. This module maintains both an explicit water and energy balance, accumulating snow in a single layer as well as calculating snow melt and sublimation. Snow melt or other liquid precipitation is then added to surface depression storage where it can either infiltrate into the soil or become overland flow runoff. Infiltration is then added to root zone soil moisture storage, transpiration is removed, and water is redistributed vertically and horizontally. Finally, if active, ex- filtrated water is passed to surface depressions where it may again runoff if depression storage capacity is exceeded. Potential evaporation from shallow and deep wetlands within the wetland fraction of each cell is calculated next. Then, for deep wetlands, the temperature is calculated 91 Surface Start Read Climate Input /LA Re ad / Controller Canopy Evolution l L_7/ Read Shortwave LU Radation l ‘ Evap. Intercept./ Storage l 1’ Snow ack §' Evap./ I p T . ; Upland ransp l '5 Controller Infiltration ICE, Evap. l l Soil Wetland Moisture Temp. l l ‘ DepressJ | Storage/ RUM“ ' ........ , Release , Runoff . l l Routing | Write '- - 4 Water Tab. Outputs Depth Evap./ Transp. Deep Unsat. Zone Perc Surface End l ,. _ - - _ - _ , / Sat. Zone I RUIIOII [I Gndwater T Controller Runoff Routing Write ‘ 7 Outputs / Figure 4.2: ILHM generalized process flow diagram, focusing on surface processes. Dashed outlines indicate optional items. 92 and an ice-pack accumulation and melt module is run. In addition to mass balance, the temperature and ice pack modules also maintain energy balance. The wetland water storage and release module calculates actual evaporation, transpiration, deep percolation, and runoff. Finally, the water table depth is updated. Following the physical process modules, outputs are written to disk for coupling with other domain modules or for comparison to observations and visualization. The model controller then loops through the remaining hours within the simulation. When the surface process modules are complete, the unsaturated zone, saturated zone, and runoff routing modules are then run in sequence. 4.2.2 Simulating Hydrologic Fluxes Hydrologic fluxes are calculated according to the equations in the Appendix A, how- ever note that for consistency, the symbols have been updated for this Chapter. There are numerous improvements and modifications to the code, which are described in detail in the sections below and Appendix B. 4.2.2.1 Canopy The canopy water balance is ASCAN = WPCP - WTHR - WEC (4-1) where ASC A N [m/s] is the change in water stored in the canopy SCAN [m], WPC p [m/s] is the total incoming precipitation, WT H R [m / s] is canOpy throughfall (WpCp— WT H R represents water intercepted by the canopy), and WEC [m/s] is canopy water evaporation. The canopy intercepts precipitation up to its storage capacity, which is a function of the Leaf Area Index (LAI) of each grid cell. Canopy evaporation and canopy transpiration, WTU [m/s], are calculated via a modified Penman-Monteith equation. ILHM treats the canopy as continuous across the upland, stream, and wetland fractions within each cell. 93 i=8 mo 98an 3&5 98:25 was 32v v53: Sm SEE .3 Ucaflaofio max—é 9&0?qu mo mason ”we BEE ’1.--------.- ”~----------------- 94 Sites , ionérm . . .. . . 2‘: mmnml 30m; «mzonmmmmnmq fivcs Total available energy for evaporation is a function of intercepted solar radiation, governed by the LAI, Stem Area Index (SAI), incident angle of solar radiation, relative proportions of beam and diffuse solar radiation, and albedo of the ground beneath the canopy. Resistance to vapor transport within the canopy depends on LAI. Resistance to momentum transport depends on input wind speed and canopy height (which varies as a function of LAI). Atmospheric vapor deficit depends on the input air temperature and vapor pressure. 4.2.2.2 Snowpack ILHM uses a modified version of the UEB snow accumulation and melt model by T arboton and Lace (1996). For coupling with ILHM, UEB was converted to MATLAB code from FORTRAN, and made fully distributed and vectorized. The snow model maintains the water balance AS5N0 = WTHR - WESN - WSNM- (41-?) Here, ASSNO [m/s] is the change in snowpack water storage SSNO [m], WESN [m/s] is evaporation (sublimation) from the snowpack, and W5 N M [m/s] is snow melt. SnOWpack sublimation is calculated using a mass-transfer approach, assuming neutral atmospheric conditions. Snow melt is the sum of liquid water in excess of the snowpack storage capacity, either from liquid precipitation or melt of the solid phase of the snOWpack. Currently, the snow model only runs in the upland portion of each cell. The wet- land and stream portions have an ice accumulation and melt routine that is roughly equivalent. This snow model does not explicitly track the temperature of the snow pack, but rather accumulates an internal energy deficit as a state variable. 4.2.2.3 Infiltration and Depression Storage The upland portion of each cell has a surface depression storage capacity that must be overcome prior to overland flow runoff. The capacity of this storage is a property of the land cover, soil texture, and average slope of each cell. The depression storage water balance is ASDEP = WTHR + l’VSNJlr + WEX — WINF - W EDP - W Roy, (43) where AS D Ep [m/s] is the change in surface depression storage S DE p [m], WIN 1: [m/s] is the rate of water infiltration into the soil, WE X [m/s] is the rate of ex- filtration from the soil due to near-surface interflow over saturation, WEDp [m/s] is depression evaporation, and WROU [m/s] is the overland flow runoff rate from uplands. Depression evaporation occurs at the potential rate, as calculated by the modified Penman-Monteith equation, using a resistance to vapor transport of 0, and a resistance to momentum transport determined by soil surface roughness. Soil infiltration depends on the average saturated hydraulic conductivity of the first soil layer, the extent of impermeable surfaces within each cell, and whether the first soil layer is frozen. The maximum infiltration rate is set equal to the average saturated conductivity, which underestimates total infiltration by not incorporating matric and ponding depth potentials. For this study, the soils are relatively coarse with high conductivity, thus direct overland flow runoff rarely occurs over much of the watershed. 4.2.2.4 Root Zone The overall root zone water balance in the upland portion of each cell is ASRZ = WINF — WEX - WESL — WTU — WINT — WDPU- (4-4) Here, ASRZ [m/s] is the change in root zone water storage S RZ [m], W ESL [m/s] is direct evaporation from the soil, IV INT [In/s] is near-surface interflow, or hori- 96 zontal redistribution between adjacent cells, and WDpU [m/s] is deep percolation to the unsaturated zone. Both interflow and over saturation due to variable vertical conductivity can generate WEX- For this study, W1 NT is a very minor term in the overall water balance, and was excluded for computational efficiency. The transpi- ration rate, WTU, is controlled by the total shortwave radiation intercepted by the canopy, canopy stomatal conductance, and soil moisture. Transpiration is partitioned among soil layers according the fraction of total root mass within each layer. Moisture is redistributed vertically between layers using either the 1-D Richard’s equation, or free-drainage at the rate of unsaturated hydraulic conductivity. For this study, the free-drainage method was chosen because of its numerical stability and computational performance. Water cannot enter the root zone from below, due to the one-way coupling between ILHM’s surface components and its deep unsaturated and saturated groundwater modules. Efforts are underway to couple these modules at run-time and incorporate capillary rise, saturated exfiltration, and other important fully-coupled processes. 4.2.2.5 Wetland Storage Shallow wetlands, deep wetlands, and streams, maintain the water balance ASWDorWS = WTHR + WICM — W DPW - Wall/,5 — Wrws - WRow, (4-5) where ASW D [m/s] and ASWS [m/s] are the changes in storage SW D [m] of deep wetlands, and SW3 [m] of shallow wetlands, WICM [m/s] is stored precipitation released by ice pack melt, WEW,S [m/s] is the rate of evaporation from surface storage in inundated wetlands, WTW,S [m/s] is the rate of transpiration from surface storage, and WROW [m/s] is the runoff or release rate. Inundated area within a watershed can vary dramatically seasonally, particularly in regions with shallow water tables. With a run-time coupling between the saturated groundwater and surface components of ILHM, direct calculation of inundated area is 97 possible. However, with the one-way coupling currently implemented in ILHM, a less direct means of calculating inundated area is required. ILHM currently implements a user-specifiable seasonal water table fluctuation function, which combined with input average depths of wetland complexes allows inundated area to vary. Runoff in wetlands is calculated differently depending on the inundation state. Inundated wetlands, i.e. those whose storage is positive or with a water table at or above the surface, release a portion of their storage each timestep. WROW for inundated wetlands is calculated as the product of a user-specified recession constant and the storage within each cell. For dry wetland cells, WROW is calculated as for upland runoff, if the rate of throughfall exceeds infiltration capacity in the wetland soils, runoff is generated. Wetland storage will be negative when the water table is beneath the surface. In this case, the minimum storage equals the depth of the water table, multiplied by the average porosity of the soil in between the surface and water table. When storage reaches this minimum, deep percolation no longer occurs. Deep percolation, WDpW from wetlands is calculated as the product of a specified wetland conductance and the storage head in each wetland cell. This head is the total head above the water table, including both sub-surface saturated water and ponded water at the surface. Potential transpiration is calculated with the Penman-Monteith equation. If the wetland is perched, transpiration WTW, S is removed from wetland storage provided water is available. Transpiration in connected wetlands is assumed to primarily occur within the phreatic zone, and thus is passed to the saturated groundwater module as WTW,P- The rate of evaporation is calculated via the Penman-Monteith equation for shallow inundated wetlands, and via the energy balance equation for deep wetlands whose temperature has been explicitly calculated. If water storage is positive, evaporation is removed from storage as WEW,S, otherwise this evaporative demand is passed to 98 'r' "1‘ ,_ .l.. ,> “J.- ‘r, 1.! r u q.- i... the runoff routing module where it will be removed (if available) from groundwater discharge to surface water or directly from the saturated groundwater module. If the wetland cell lies within an internally-drained zone (i.e does not drain to the host watershed), it is passed to the groundwater module as WEW, I [m/s], otherwise, evaporation WEW, E [m/s] is passed to the runoff routing module. Evaporation in wetlands with negative storage is calculated as soil evaporation with the Penman- Monteith equation. The difference between inundated and dry wetland evaporation is the surface resistance term, which is 0 for inundated wetlands and dependent on depth-to-water in dry wetlands. 4.2.2.6 Wetland Temperature and Ice The energy balance of a wetland can be expressed as DwPCpATW = QR + QA - QH - QE - QC, (4-6) where the terms on the left hand side are, in order, the depth of the wetland, density of the water, specific heat capacity of water, and change in temperature of the wetland. The energy flux terms are, in order, net radiative flux, advective flux, sensible heat flux, latent heat flux, and conductive heat flux. For computational efficiency, terms that are less important in shallow wetlands, including temperature change (i.e. heat storage), advective flux, and conductive flux may be ignored in ILHM. For deep wetlands, heat storage within the wetland becomes important. For those cells, a full multi-layer energy balance module that implements the 1D eddy diffusion model by (Hostetler and Bartlein, 1990). This module incorporates wind—generated diffusive eddies, as well as full convective mixing. A simple ice accumulation module has also been implemented that builds ice when the temperature of the lst wetland layer drops below zero. From there, the energy balance of the ice is solved via Equation 4.6, ignoring conductive fluxes to the water below. Because of this simplification, modeled ice thicknesses are greater observations, 99 however timing of melt is reasonably simulated which is important for calculating evaporation flux. The ice pack stores incoming precipitation separately from ice accumulated from below. This water is then released to the system proportionately as melt occurs. 4.2.2.7 Deep Unsaturated Zone The deep unsaturated zone module maintains the water balance ASUZ = WDP - WRCH- (4-7) Here, ASU Z [m/s] is the change in unsaturated zone storage SU Z [m]. Total deep percolation leaving the root zone and wetlands, WDp [m/s] is given by WDP = FUWDPU + (FW + F5) WDPW, (4-8) where FU, FW, and F S are the fraction of upland, wetland, and streams in each cell respectively. Individual wetting fronts representing the flux of water WD p travel vertically downward through the deep unsaturated zone at a rate controlled by the hy- draulic conductivity of the first saturated groundwater layer. The parameters for this were determined by fitting a 1-D Richard’s Equation model to water table fluctuation data collected in a nearby watershed. This simple routing greatly enhances the stability and speed of the model, but (1088 not allow vertical fluxes through gradients in matric potential, nor wetting front attenuation and broadening. An adequate description with a 1-D Richard’s Equation InoCiel was deemed too computationally extensive for a domain of this size, though it C()‘Jlld be easily implemented for smaller model areas. Methods of intermediate com- pleXity that may offer a more desirable trade off between speed and accuracy are under inveStigation, including transfer functions or wetting-front tracking as implemented in the UZF package in MODFLOW 2005 (Niswonger et al., 2006). After traveling through the unsaturated zone, water enters the saturated zone as grouIldwater recharge, WRC H [m / s]. The total travel time depends on rate of wetting 100 front movement and the thickness of the deep unsaturated zone. For this study, that thickness is assumed constant in time, and was set by extracting the average depth to the water table from a preliminary model run. Full run-time coupling of the surface, unsaturated, and saturated zone codes already under development will allow the thickness of this zone to vary. 4.2.2.8 Saturated Zone Saturated zone groundwater flow is calculated using MODFLOW-2000 (Harbaugh et al., 2000), and maintains the water balance A532 = WRCH — WSAT - WDIS - WTW,P - WEWJ — W WEL - WBDR, (4-9) Where ASSZ [m/s] is the change in saturated zone storage 33 Z [m], WDIS [m/s] is saturated discharge to water bodies, WW E L [m / s] are anthropogenic withdrawals, and WB DRY [m/s] is the sum of a variety boundary condition terms. ILHM prepares all input files for MODFLOW, and passes groundwater recharge and evapotranspirative fluxes to it after running the surface components. Groundwater discharge to drains, riVers, streams, or general head boundaries is then passed to a post-processing routine that routes these through the drainage network, and calculates remaining open-water eva«poration. A virtue of the one-way linkage to MODFLOW within ILHM is that any groundwa- ter Code may be used without modification of the source. Also, legacy models within thOSe codes may be used directly. This allows users to leverage the thousands of exist- ing groundwater models to rapidly build a more capable integrated hydrologic model Within ILHM. As already discussed, this one-way linkage does limit the ability to fully deScribe certain processes within ILHM. The importance of the full coupling can vary greatly depending on the degree of interaction between surface- and ground-water. ThuS even after completing the development of full-coupling, some applications (for Instance, parameter estimation) may still benefit from simpler one-way coupling. 101 4.2.2.9 Runoff Routing Total stream discharge generated in each cell is WRUN = As [FUWROU + (Fw + F3) WROW] + AG (Fw + F5) WDIs - ’EW,E- (4.10) Here, WRU N [m3/s] is total runoff from each cell , WUP [m3/s] is the up gradient discharge, A3 [m2] is the area of each cell in the surface components, and AG [m2] is the area of each cell in the groundwater module. Note that currently, ILHM does not allow loss along the flowpath. In relatively humid regions with shallow water tables, soil saturation increases downslope, and groundwater heads generally exceed heads in the stream. Thus streamflow routing can be decoupled from the surface and groundwater processes, and handled as a post-processing step. The total runoff generated in each cell is then routed along a flowpath generated by a D8 flow algorithm (O’Callaghan and Mark, 1984). The velocity along the flowpath depends on the kinematic wave velocity, Manning’s roughness of the cell and the average slope of the cell. Runoff is aggregated only at a user-specified set of points. Flow times between points are also specified. This simple routing approach is very fast, but tends to under-represent stream water velocities at high flow levels, thus flood peaks tend to be too narrow. In environments that produce more overland flow, enhanced flow routing algorithms could be implemented with little performance penalty, including Doo flow routing (T arboton, 1997), and stage-dependent channel velocity. 102 4.3 Model Domain and Input Data Prepa- ration 4.3.1 Study Area, Simulation Period, Model Domain, and Validation Data The Muskegon River Watershed (MRW) is located in central Lower Michigan (Figure 4.4), and covers approximately 7,400 km2. The Muskegon River ultimately drains into Lake Michigan, the second largest of the Laurentian Great Lakes. Its primary drainage, the Muskegon River, begins as an outflow from Higgins Lake, and flows through another large inland water body, Houghton Lake. It discharges to Lake Muskegon, which is connected to Lake Michigan via a narrow canal. There are several other sizable lakes along the Muskegon River, as well as a significant number of wet- lands, particularly in the upper portion of the watershed. According to other studies, the majority of the total flow in the watershed comes from groundwater discharge to streams (Jayawickreme and Hyndman, 2007). Several large control structures modify the watershed, including Hardy, Rogers, and Croton dams in the central portion of the watershed. The model boundary domain for this study has been extended well past the MRW to a total area of approximately 18,900 km2. The boundary was chosen to coincide with major hydrologic divides in neighboring basins. Integrated surface- and ground-water models must incorporate boundaries appropriate to both domains. In this region, particularly in low-relief areas, the watershed divide does a poor job of describing the “groundwatershed” of the Muskegon River. The model bottom is given by the contact between the surficial unconsolidated aquifer and the bedrock beneath, which is assumed to be impermeable. A significant portion of the bedrock contact is comprised of permeable aquifer materials, however for this study we have assumed that the fluxes between the surficial and bedrock aquifers are small relative to the surface to shallow- 103 ° A Model Boundary . e" , ”(W --i o USGS A/v NCDC (hrly/dly) e MAWN - Stream Discharge ~ Static Water Level I: Municipalities Figure 4.4: Map of the Muskegon River Watershed and surrounding model boundary. In part A, major municipalities along the drainage network are highlighted. Part B shows the locations of input gages and observation data locations. groundwater fluxes of interest here. This assumption is supported by two lines of evidence: 1) the freshwater-saline contact in bedrock aquifer materials beneath much of the MRW is relatively shallow, suggesting a small downward flux from the glacial sediments to the bedrock, and 2) the waters within the three primary bedrock aquifers are saline or brine (Westjohn and Weaver, 1998). A total of 28 years of hydrologic fluxes are simulated, from January I“, 1980, to December 31“, 2007. During this period, new input data sources became available, greatly increasing the quality of model input data (see Table 4.1). A new gage network came on line in 1996, the Michigan Automated Weather Network (MAWN) (MA WN, 2008), providing the first reliable measurements of total solar radiation in the region. Also beginning in 1996, NEXRAD precipitation data at 4 km, hourly resolution became available (Andresen, 2008). In 2000, Leaf Area Index (LAI) data became available from the Moderate Resolution Imaging Spectrometer (MODIS) on NASA’s TERRA and AQUA satellites (Knyazikhin et al., 1999). Prior to this, LAI had to be inferred from the Normalized Difference Vegetation Index (NDVI) products produced using data from instruments such as the Advanced Very High Resolution Radiometer (AVHRR). These heterogeneities of data sources and quality present some 104 of the greatest challenges to long-term change modeling. Below, (in Climate Data Preparation, Solar Radiation Preparation, and Precipitation Data Preparation) we present the extensive efforts taken to minimize systematic differences between datasets while still benefiting from newer, higher-quality data sources. Table 4.1: Input climate and LAI data sources. Start and end dates are generalized. Locations of gage data are shown in Figure 4.4. Source Type Units Start End Resolution Frequency AWOSl dew point F 1980 2006 15 stations hourly 2 m air temp. F 10 m wind speed knots sky conditions oktas pressure mbar MAWN2 relative humidity % 1996 2007 4 stations hourly 1.5 m air temp. C 3 m wind speed m/s total solar rad. kJ/m2 2 in soil temp. C 4 in soil temp. C NCDC3 dail precip. in-100 1980 2007 67 stations daily houry precip. in-100 1980 2007 10 stations hourly NEXRAD4 hourly precip. mm-10 1996 2007 4 km hourly GIMMS5 NDVI - 1981 2006 8 km 15 days AVHRR6 1989 2003 1 km 7 days MODIS7 LAI (v5) — 2000 2007 0.93 km 8 days 1(Andresen, 2007) 2(MAWN, 2008) 3(Center, 2001; NCDC Hourly, 2008) 4(Andresen, 2008; Fulton et al., 1998) 5(Tucker et al., 2005) 6(E'idenshinlc, 1992) 7(Knyazilrhin et al., 1999) In this study, two types of data will be used to validate the integrated model, groundwater heads and stream discharge. Static water levels were measured in nearly 50,000 wells during the 1980 - 2007 time period by drillers installing drinking water wells (see Figure 4.4). Though the quality of each individual measurement may be low, the dataset in aggregate presents an dataset for evaluating model-predicted groundwater heads. The USGS has maintained stream discharge gages at 6 locations, though no more than 5 were operating at any given time. These will be used to eval- 105 uate total simulated stream flow, along with errors in model fluxes. In addition to these, we have collected 105 stream discharge measurements during with a combina- tion of wading and bridge—board cross-section velocity methods. These measurements were collected at or near baseflow conditions in the system and will be compared to simulated stream baseflow. Table 4.2: Table of observations and validation data. Location abbreviations “MR” refers to “Muskegon River”. The “Area” column refers to the surface watershed area. Source/ Type Location Area (km?) Units Start End USGSl MR @ Newaygo 6,086 ft3/s 1980 1993 daily gage MR @ Croton 5,991 1995 2007 MR @ Evart 3,711 1980 2007 Little MR @ Oak Grove 894 1995 2007 Clam @ Vogel Center 629 1980 2007 Bear @ Muskegon 43 1980 2007 This study 42 locations 4.2 - 928 m3/s 2003/07/23—25 synoptic flow 28 locations 0.2 - 3,047 2004/ 08/ 21-24 35 locations 1.2 - 530 2004/10/09-10 MDEQ2 250,000 locations - ft 1980 2007 static water level 1(Usos, 2008) 2(MDEQ, 2008) 4.3.2 Static Inputs There are four primary static GIS inputs, shown in Figure 4.5. Many of the ILHM parameters depend on either soil texture, land cover type, or both. Here we use land cover from a 1998 update of the 1978 MIRIS aerial-photo classified land use map of the MRW region (Torbick et al., 2006). The two dominant land use classes in the MRW and surrounding region are forest and agricultural lands. During the 20th century, the extent of agriculture declined sharply, while forest, urban, and open/ grassland land covers moved into abandoned marginal agricultural lands. Three-dimensional soil textures to a depth of 3 m are taken from the digitized soil survey SSURGO data (USDA-NRCS, 2008). In the survey data, each map unit has a unique layering specified, based primarily on the soil horizons. To simplify 106 extracting the texture data, the texture values were mapped to four layers from 0-20, 20-50, 50-100, and 100-300 cm. To select the appropriate texture for each layer in the SSURGO map units, the majority texture was chosen. Soil hydraulic parameters were then mapped to these textures. Note, the layering within ILHM will not necessarily match the four layers extracted from the survey data as the ILHM total root zone thickness is variable, and the global number of layers adjustable. Thus on import, the soil parameters are mapped to ILHM layers based on the weighted average of the input layers within in each ILHM layer. The hydraulic conductivities of the deeper unsaturated and saturated zone sed- iments are assumed to be reasonably parameterized by the zones in the surficial, Quaternary geology map (Farrand and Bell, 1982). As illustrated in Figure 4.58, the dominant soil textures are sands, sandy loams, and loamy sands. These relatively coarse textured materials have high infiltration capacities and hydraulic conductivi- ties. As a result, most of the rainfall and snow melt within the watershed percolates into the soil, producing very little direct overland flow runoff. The important excep- tion to this is in flooded areas and low-lying uplands in close proximity to streams and lakes. These coarse-textured soils are produced from similarly coarse Quaternary sediments (Figure 4.5C). This map has three relatively distinct regions: 1) a low-relief area in the upper watershed comprised of alternating glacial outwash sandurs and finer-textured lacus- trine sediments from pro-glacial lakes, 2) an inter-lobate moraine region in the center of the watershed made up of chaotic assemblages of mostly coarse glacial till and subsequent incised alluvium, and 3) a lake shore region in the lower watershed of rel- atively high-energy near-shore deposits and coarse alluvium. The groundwater model consists of three layers, with conductivity parameters directly related to this map in the first layer, and a more generalized three-region conductivity map for the deeper two layers. A preliminary steady-state groundwater model was used to calibrated 107 hydraulic conductivity values, as discussed below. Surface flow routing directions and velocities are determined in part by the digital elevation map (DEM) (Gesch et al., 2002). The expanded model boundary encom- passes the highest elevation point in Lower Michigan, at roughly 500 meters. Lake Michigan defines base level of the watershed, at an average of approximately 176 meters during the simulation period. The highest relief portion of the watershed is located around Big Rapids, MI (see Figure 4.4A). Land Use Soil Texture Quaternary Geology Elevation (masl) - Urban - Deciduous [:1 Sands [:1 Outwash/Alluvium <22-EEE-> mAg - Coniferous - Till 325 425 El Shrub E] Wetland - Fines - Lacustrine A B ‘K' " Figure 4.5: Static GIS inputs for ILHM expanded model boundary Channel geometries were measured at each stream discharge measurement loca- tion. For measurements taken during low-flow conditions, discharge was correlated to stream width, average depth, and hydraulic radius with equations of the form Y = AX B , where Y is either width, depth, or radius, and X is low-flow discharge. Next, values from flow accumulation routine in ArcGIS (which counts the number of cells upstream of every point in the watershed) were correlated to low-flow discharge at each discharge measurement location, using an power-law equation. This equation was then applied to the flow accumulation map to create an average low-flow discharge map. This map, along with the channel geometry correlations, then calculated width, average depth, and hydraulic radius across the model domain. Lake and wetland depths are also required by the model for two purposes: 1) shallow 108 wetland depth information allows for a variable water table to affect inundation, and 2) deep wetlands and lakes require depths for temperature calculation. The Michigan Department of Natural Resources provides maps of the bathymetry of most of the large inland lakes in Michigan (MDNR, 2009). The Michigan Department of Environmental Quality has digitized some of these maps (MDEQ Bathymetry, 2008). For those that had not been digitized, the maps were imported into ArcGIS, georegistered, and digitized. Only wetlands larger than 0.25 km2 area were digitized. For smaller wetlands with bathymetry information, average depths were estimated from the maps. Large wetlands or lakes without bathymetry were assigned a single average depth across the model area. Smaller wetlands were assigned depths based on interpretation of the National Wetlands Inventory classification scheme (Table 4.3). Note here that negative average depths indicate wetlands that are inundated for only part of the year. Table 4.3: Wetland depths assigned based on NWI classifications. The depths are based on a sinusoidal water table fluctuation of ”2.5 m. NW1 Water Class Average Depth (m ) Typical Inundatz'on Saturated -1.5 Rarely Temporarily Flooded -1.0 Less than two weeks Seasonally Flooded -0.85 Dry by end of growing season Semi-permanently Flooded 0.25 Wet during growing season Artificially Flooded 1.0 Varies, controlled artificially Intermittently Exposed 1.1 Exposed during drought Permanently Flooded 2.0 Remains flooded in all years 4.3.3 Climate Data Preparation The first step in preparing climate input data from a heterogeneous set of data sources is to create a set of utilities to interpret the native format of those data, including all data and quality control flags, into a standard format for ILHM. As Table 4.1 illustrates, some data are delivered in English units only, while others included op- tions for output in SI units. ILHM uses the full SI standard set of units (meters, seconds, kilograms), so unit conversions were necessary in some cases. The AWOS 109 data included 10 In wind speed and dew point, where MAWN provides 2 m wind speed and relative humidity. While similar, these data were incompatible. 2 m wind speed values were approximated from the 10 m values using a power law wind profile (Peterson and Hennessey, 1978), and dew point was transformed to relative humidity using the equation from (Lawrence, 2005). The data were then subjected to a further level of quality control to identify three problems: 1) individual erroneous values, perhaps due to instrument error or subse- quent processing and data transfer, 2) erroneous periods within otherwise acceptable stations, 3) “bad” stations, whose measurements were consistent outliers. Individual bad measurements could be identified by specifying acceptable ranges of values; au- tomated outlier detection methods did not fare as well because of the relatively small number of stations. These ranges were determined by visual inspection of the time series, and by plotting the frequency histograms. Erroneous periods in data were identified primarily via visual inspection, and could only be done for the MAWN dataset because of the small number of stations and relatively short time period. Entire stations were omitted whose values consistently difl'ered from their neighbors past a specified threshold. After these steps were complete, there were hours within the model time period during which no stations reported a particular required input value. Any gaps in data, except for solar radiation and precipitation, were filled with the day and hour average of that data type from other years. Gaps were filled only for hours with no station data present; ILHM does not assume that any given station must have data defined for all time steps, as discussed in the next section. This method was chosen, as opposed to interpolating temporally between present values, because gaps where all stations were missing tended to be wider than a few hours, thus interpolation would not preserve diurnal cycles. 110 4.3.4 Spatial Distribution of Climate Data Data input to ILHM must specify three things: 1) the time of the measurement, 2) the value, and 3) how to distribute those data over the model domain. Spatially- distributed datasets, such as NEXRAD radar or satellite-sensed LAI inherently con- tain all three types of information. Gage data, however, must be interpolated to the ILHM domain. As each interpolation must be done every hour for over a dozen input data types, the scheme must be very fast, or be run prior to model execution and amortized over multiple model runs. Also, the scheme must handle stations dropping off- and coming on-line. An early version of ILHM applied to small sub-watershed of the MRW, published in Hyndman et al. (2007), accomplished this through a fast inverse-distance-weight method. Such a method, however, is fundamentally unsuitable for hourly interpo- lation at larger spatial scales. As weather systems move across the model domain, they cause abrupt changes in data values, for instance temperature and precipitation. If the gage spacing is greater than the distance across a weather front, as is often true for the MRW, any interpolation method that includes more than one station value and assumes stationarity will fail to accurately represent the true field. This is particularly true with precipitation data, as the intensity of precipitation is crit- ical to canopy interception and runofl' generation. A more appropriate method is to distribute a single gage data values within Theissen polygons. The precipitation data preparation section below describes how this method can be modified to better reproduce temporally-aggregated fields. 4.3.5 Solar Radiation Preparation Solar shortwave radiation drives the hydrologic cycle, and is a critical input to ILHM. However, for much of the model period only a single gage recorded solar radiation in the watershed, near Muskegon (see Figure 4.4A). Comparison of this gage with the 111 MAWN data beginning in 1996, shows that this gage does a poor job capturing total solar radiation, therefore it was not used in this study. To provide solar radiation prior to the availability of MAWN data, more fully represent its spatial variability, and to split total solar radiation into beam and diffuse components, a real-sky radiation model was used. This model requires air temperature, relative humidity, and Opaque cloud fraction. The only cloud fraction data available came as sky condition observations from airport weather stations. These observations, (clear, overcast, etc.), crudely represent cloud cover conditions in oktas (1 /8 increments). However, in 1996 the stations apparently changed their reporting, causing an abrupt shift in the distributions of cloud cover observations (Daz‘ et al., 2005). This shift had to be corrected to force consistency across the 28 model years. The distributions of cloud cover fractions from 1980 - 1996, and 1997—2007 were calculated. The distribution of data from 1980-1996 was shifted to match the later period by adjusting the mean and range of values within each quartile of the distribution. Cloud cover observations are used to calculate real-sky radiation from clear-sky modeled radiation. Even with the available observations, the model cannot represent the full range of values present in the gage data, even though the annual mean values were well simulated. To partially account for the monthly variability in cloud cover that is not fully captured by the sky condition observations, the gaged data were compared to simulated values, and a monthly-average weighting was calculated. For the period 1997—2007 individual monthly weightings were used, while for 1980-1996, an average monthly correction values from 1997—2007 were applied as no gage data were available for that period. After these corrections, the model captures much of the variability present in the gaged data, and provides additional spatial informa- tion not available from the more limited MAWN locations. Figure 4.6 plots annual model-averaged total solar radiation. The period of 1980-1996 includes only mod- 112 eled radiation, while 1997-2007 includes modeled and gaged values. Notice that the mean radiation is reasonably well simulated, but inter-annual variability is not fully captured prior to the availability of radiation data. .5 O) .5 -®®. _3 0| ‘9 _L 01 N W Annual Mean Solar Flux (W/mz) 61‘ 01 _L ‘45 (9“? 80 1984 1988 1992 1996 2000 2004 2008 Figure 4.6: Plot of annual model-averaged solar radiation, from 1980 through 2007. There are two primary periods, 1) only cloudiness observations are available for short- wave radiation modeling, 2) cloudiness observations along with limited solar radiation data are available. 4.3.6 Precipitation Preparation Within and around the model domain, many more stations stations report daily pre- cipitation than hourly (see Table 4.1). To incorporate the spatial variability present in the daily data into the hourly precipitation dataset, a four-step procedure was used. 1. For each model hour, hourly precipitation data are to interpolated to a grid using Theissen polygons (this is identical to a nearest-neighbor interpolation). This method was used rather than a scheme that uses multiple stations to predict unknown locations because of the issues with weather fronts mentioned previously. 2. Daily precipitation totals are not recorded at identical times at each station. Most of the stations recorded near midnight GMT, while others recorded near 113 noon GMT. Because of this, the daily data could not be directly interpolated to create a single daily precipitation grid. The daily data at each station were summed weekly, and interpolated to the grid with a inverse distance weighting scheme. 3. The hourly grids were summed weekly, a weight value for each grid cell was calculated as the ratio of the summed daily and summed weekly grids. 4. For each hour in the week, the weight grid was multiplied by the hourly precip- itation grid. The net effect of this procedure is that total precipitation values in each cell are fixed by inverse distance weighted interpolation of daily precipitation stations, while the hourly distribution is governed by the nearest hourly station. N EXRAD precipitation data is capable of providing relatively high resolution dis- tributed precipitation data at hourly intervals. For the MRW, a 4 km Level III calibrated, gage-adjusted dataset became available in 1996 (Sea, 1998). To determine whether NEXRAD data could be used in place of interpolated station values, monthly grids of each were compared. A time-series from 1996 - 2007 of the model-average monthly precipitation estimates are plotted in Figure 4.7. From visual inspection, there are three distinct periods of NEXRAD accuracy. From 1996 — 2000, model- average NEXRAD precipitation matches interpolated station values quite poorly. Beginning in April of 2000, growing-season NEXRAD values seem to provide rea- sonable monthly precipitation estimates. After March of 2003, NEXRAD and station values compare favorably in all months. This study assumes that if monthly NEXRAD precipitation compares favorably with interpolated station values, the NEXRAD dataset then likely provides more accurate hourly distributed estimates of precipitation over the watershed. Based on the comparison in Figure 4.7, no NEXRAD data prior to 2000 were used. Beginning 114 ® — Station ‘ — NEXRAD 4 N U! i G) N O I .5 .5 O 01 I l ... “(l ’1' "l*‘\\ : l 1 ii 1 L I l 1996 1997 1998 1999 2000 2001 2002 2003 2004 2005 2006 2007 2008 Monthly Precipitation (cm) OO‘I Figure 4.7: Plot of monthly model-averaged precipitation from daily interpolated NCDC precipitation gages, and N EXRAD radar from 1996 - 2008. The quality of the NEXRAD data can be differentiated into three periods, 1) all monthly observations are a poor fit to gage data, 2) NEXRAD provides a good match to months with liquid precipitation only, and 3) NEXRAD data and gage data match closely in almost all months. in April of 2000, NEXRAD data were used in place of station values for April - October. After April 2000, NEXRAD data were used exclusively. 4.3.7 Leaf Area Index (LAI) Preparation Leaf area index (LAI), a measure of the one-sided leaf area per unit area of horizontally-projected ground surface, drives seasonal vegetation change in ILHM. LAI controls vegetation height, canopy moisture and radiation interception, canopy water capacity, albedo, and canopy resistance to vapor and momentum transport. In 2000, a 1 km, 8-day LAI product (MOD15A2) was released, calculated from data captured by the Moderate Resolution Imaging Spectrometer (MODIS) instrument carried on-board NASA’s TERRA and AQUA satellites (Knyazz'khz'n et al., 1999). Prior to this, however, there are no readily-available distributed LAI products. To extend the LAI record as far backward as possible, a multi-step procedure for deriving LAI from normalized difference vegetation index (NDVI) values was devel- oped. NDVI is calculated using two spectral bands, the red (0.6 pm) and near infrared 115 (0.9 pm), and is calculated as NIB—RED . NDVI— NIR+RED' (4.11) Prior to MODIS, instances of the Advanced Very High Resolution Radiometer (AVHRR) instrument aboard several NOAA satellites have been in orbit since 1981. Two NDVI products derived from AVHRR data are readily available, the 15-day 8 km GIMMS global dataset, and the 7-day 1 km continental US dataset (see Table 4.1) (Eidenshz'nk, 1992; Tucker et al., 2005). The GIMMS dataset was specifically produced to correct certain errors present in the AVHRR data, including sensor and orbital drift (Tucker et al., 2005). Its coarse resolution, both temporally and spatially, are less than ideal for this application. The GIMMS grids were used to correct the 1 km N DVI. The mean value of the 1 km cells were calculated for each 8 km GIMMS grid cell. A weight grid was then calculated as the ratio of the 8 km values to the upscaled 1 km data. This weight grid was then multiplied by the 1 km data. This re-weighting allows the corrections present in the GIMMS dataset to be applied to the higher resolution 1 km, 7 day data. As discussed in Baret and Guyot (1991), the theoretical relationship between NDVI and LAI is an exponential of the form LAI = 0,4 - exp (CB - NDVI) (4.12) where C A and C B are constants that vary by biome type. Values of these constants for each land use class were Optimized by comparing, AHVRR-derived NDVI to MODIS LAI values. The 1 km NDVI dataset most closely matches the MODIS resolution, and four years of overlap between the two datasets were available for the following procedure: 1. Select all cells whose fraction of a given land use exceeds 90%. 116 2. Extract LAI and N DVI values for these cells. This requires re-interpolating one Of the datasets to match the other spatially. Here, the LAI data were resampled to 1 km resolution. 3. Temporally interpolate one Of the datasets to match the other for each model cell. Here, the LAI data were linearly interpolated tO match the dates Of the NDVI grids. 4. Optionally, randomly select a subset Of these points for computational efficiency. 150 points from each 7 day grid were chosen for this analysis. 5. Fit an equation Of the form shown in Equation 4.12 to the LAI and N DVI data, optimizing the values of C A and C B- Here, a bound optimization was conducted in MATLAB, using a mean-absolute residual objective function. 6. Calculate LAI from NDVI using Equation 4.12 and the constants optimized from step 5. Constrain the LAI values produced not to exceed the maximum and minimum values in the MODIS dataset. Optimized values of C A and CB are shown in Table 4.4, along with their R2 values. For all land use classes except urban, the empirical function provides reasonable estimates Of LAI from N DVI grids. Table 4.4: Table Of Optimized constants for the N DVI —» LAI mapping function 4.12, and coeflicient Of determination R2. Open water and barren are assumed to have LAI = 0 in ILHM. Land Use CA CB . R2 Urban 0.103 4.13 0.33 Agriculture 0.082 4.03 0.67 Open/Shrub 0.103 4.99 0.70 Deciduous 0.091 4.79 0.74 Coniferous 0.152 4.21 0.62 Wetland 0.096 4.39 0.67 Open Water - - - Barren - - - After calculating LAI grids from the 8 km and 1 km NDVI grids, the result is a heterogeneous resolution dataset from 1981 through 2007. To resample these grids to model resolution (425 m), the following downscaling procedure was applied: 1. Identify cells for each land use class that exceed 75% of that type, and extract LAI values for each Of these cells. 2. Interpolate these cells to a grid at model resolution using the nearest neighbor. This produces eight grids for each input LAI grid, one for each land use class, representing spatial distribution the LAI Of that land use at a point in time. 3. Add these eight grids together based on the fraction Of each land use class within each model cell. For example LAI = furbLAIurb'l'fagLAIag+f0penLAlopen'l'fco'nI/Alcmz+fdecLAIdeC'l-fuyet LAIwet where the sum of the land use fractions f is 1. 4. The result from Step 3 is a grid Of LAI values at model resolution. Now, correct the cell-by-cell values to match the coarse input dataset by resampling this grid to the original resolution (8 km, 1 km, or 0.92 km) and calculating the mean of the model grid values. Calculate a weight grid as the upsampled mean divided by the original LAI value. U! Resample this weight grid to model resolution, and multiply it by the model resolution LAI grid, resulting in a grid whose values match the original LAI data (or the NDVI-derived LAI values) at the coarse scale, but incorporate finer resolution information from land cover. The result Of these procedures is a 27 year continuous time series Of LAI values from two different instruments and three different sources. Consistent spatial biases remain between sources, primarily as a result of factors that control LAI other than NDVI. This can be easily visualized by calculating the 118 average annual maximum LAI from each dataset (not shown). In particular, areas Of mixed land use exhibit significant variations between the downscaled N DVI-derived and MODIS-derived LAI. However, the spatial bias varies consistently throughout the year, thus continuous spatial bias correction is needed. The following procedure was applied: 1. For the period of overlap between datasets (7 years for 8 km N DVI and MODIS, 4 years for 1 km NDVI and MODIS), temporally resample one Of the datasets linearly (they have already been downscaled to the same spatial resolution) to match the other. For this study, both datasets were resampled to the same 7 day intervals. 2. Calculate the ratio, LAI/NDVI for each grid. Then, calculate the Julian-day average ratio across all overlap years. The result is a set of ratio grids at 7 day intervals. 3. For the temporally-resampled NDVI-derived LAI grids, multiply each grid by the ratio grid for the corresponding Julian day. For instance, day 21 of 1992 is multiplied by the ratio grid for day 21 from step 2. Constrain the cell-by-cell maximum and minimum values to those Observed in the MODIS LAI grids. The result from this last procedure is that the long-term spatial averages of the N DVI—derived LAI grids match the MODIS LAI data. This was deemed essential for maintaining consistency within the model time period. Two more steps are then applied. The three data sources are composited, GIMMS data are used prior to 1989, following which the 1 km N DVI-derived LAI are used until MODIS data become available in 2000. Periods during which values were missing for over 1 month are filled with the Julian-day average from all model years. This includes the entire year Of 1980, and a several-month period in late 1989. Plotted in Figure 4.8 are the monthly model-average LAI values after all steps are complete. In 119 general, the three datasets are in good agreement, though peak LAI values derived from NDVI still exceed those from MODIS by 10—20% in some years. 3.3-y (l lll lllfll fill '1 i)” ll ll l2 E‘WUUUUUullltllllUlUl/UUUUUM £80 1982 1984 1986 1988-1990 1992 1994 1996 1998 2000 2002 2004 2006 2008 Figure 4.8: Plot of monthly model-averaged LAI from three sources: 1) 8 km GIMMS AVHRR data, converted from NDVI to LAI, 2) 1 km AVHRR data, uncorrected, converted from NDVI tO LAI, and 3) MODIS LAI data. 4.3.8 Model Parameters Model parameters in ILHM fall into four general classes: 1) physical constants with widely-agreed upon values (Table B4), 2) land use or soil texture dependent values with readily-available literature values that generally agree but vary from source to source, 3) descriptive physical parameters that are constant across the model, but are expected to be domain-independent, and 4) site-specific physical parameters whose values must be Specified by the modeler. Tables Of physical constants are provided in Appendix B, along with reference citations. Land use and soil texture dependent parameters are those that can generally be expected to be independent of the model domain location, thus literature values can reasonably constrain these parameters. Several sources Of soil texture parameters were used (Sexton and Rawls, 2006; Schaap et al., 1998; Stieglz'tz et al., 1997). Values of parameters in common among these sources were averaged. Some parameters varied little from source to source, including field capacity and total porosity that differed roughly 10% among sources. Saturated hydraulic conductivity and residual 120 porosity varied substantially, with individual sources differing more than 50%. All other parameters values came from a single source. Land use parameters are less commonly available, therefore a single source was used for these parameters. All parameter values are listed in Table 31. Domain-independent physical parameters relate to soil and canOpy processes, but are assumed to be independent Of land use or soil texture. Albedo and emissivity parameters fall into this category as well. These are listed in Tables 35 and B6. Table 8.7 lists those parameters whose values are domain-specific. These include the impermeable fraction of urban areas, wetland bottom hydraulic conductance, wet- land recession constant, deep unsaturated zone velocity, and saturated zone domain hydraulic conductivities. Urban areas were assumed to be 20% impermeable, a value that represents the paved portion of impermeability in residential urban areas (Ak- barz' et al., 2003). Other types Of imperviousness, like rooftops, quite often drain into the subsurface in residential areas. This value likely underestimates impermeability in dense urban, commercial, and industrial areas that make up a very small (and ungaged) portion Of the watershed. Future improvements will parameterize urban at Level 2 or 3 classification (Anderson et al., 1976)and incorporate this variability. Subsurface hydraulic conductivity values were preliminarily calibrated to a small subset of the static water level measurements using an early steady-state version Of the groundwater domain model. These are listed in Table B9. The velocity Of a wetting front in the deep unsaturated zone is closely related to these values, and was estimated using well hydrographs from the nearby Grand Traverse Bay Watershed (Table 8.8). A single wetland recession constant was derived from the Evart gage as the log— linear slope Of the recession portion of the hydrograph. There are two wetland bottom hydraulic conductance values, one for wetlands connected to the water table, another for disconnected wetlands. The connected wetland value was calculated assuming a 121 1.5 m (1‘1 horizontal hydraulic conductivity, with a vertical anisotropy Of 10 (i.e the ratio of horizontal to vertical conductivity), and a wetland bed sediment thickness Of 2 m. Conductivity values were similar to those calibrated in the groundwater model. The disconnected wetland conductance was arbitrarily chosen to be 8.6e—5 m d‘l. 4.4 Uncalibrated Model Results and Discus— sion ILHM calculates all near-surface hydrologic fluxes and storages with hourly timesteps for each grid cell, providing a wealth of information to validate the model’s predictions. There are two readily-available types Of data in many regional watersheds including the MRW, static water levels collected from drinking water wells, and stream discharge values measured by gaging stations. Each static water level measurement provides a snapshot Of the groundwater head at a single point in time and space. Stream discharge measurements integrate a much larger area and provide a means for bulk comparison of fluxes. Comparing model predictions to both types of data is essential for integrated surface- and ground-water modeling. Within the groundwater domain, hydraulic head is inversely sensitive to hydraulic conductivity and groundwater recharge Scan- lon et al. (2002). A model that reproduces measured heads will not necessarily ac- curately describe groundwater discharge. Similarly for the surface and near-surface domain, a model may reproduce stream discharge fluxes accurately while not accu- rately describing the partitioning Of the water budget between groundwater recharge, overland flow runoff, and near—surface throughflow. Here, uncalibrated model predictions Of daily stream discharge and hydraulic head are compared to data described in Table 4.2. With the exception Of a preliminary, steady-state manual calibration Of hydraulic conductivity values, the model has not 122 been calibrated to specifically match either fluxes or heads. Model calibration could substantially improve fits to measured data. Nevertheless, uncalibrated results are presented here to provide a more Open assessment of the capabilities of ILHM. Of the six USGS gages in the watershed with data records longer than a few years (Table 4.5), only three, Muskegon River at Evart, Clam River at Vogel Center, and Bear Creek near Muskegon, have continuous records for the entire model period, 1980 - 2007. Daily hydrograph comparisons for those three gages are plotted in Figure 4.9. Note that they are plotted in order Of descending catchment area. In all three catchments, ILHM describes the overall character of the hydrographs, particularly at the largest scale. Close inspection Of the hydrographs reveals that hydrograph recession is well- represented by ILHM. While ILHM has a site-specific recession constant for surface water bodies, the net slope of the hydrograph recession curve is the combined response of surface routing and evaporation. Particularly at the largest scales, rise times in ILHM are slower than in gage data and peak heights tend to be lower, especially dur- ing the spring. This suggests that a mechanism responsible for generating part Of the peak discharge response is not fully described. Because impermeable runoff is not a significant factor in much of the watershed, the likely reason for slower rise and lower peaks within the model is that seasonal inundation is not completely captured by the simple model described above. A run-time coupling Of the surface and groundwater domain modules will allow a much more detailed examination of this phenomenon. Stream baseflows, critical for ecosystems and water resources, are well-described by ILHM. Summer low flow values depend on both modeled groundwater discharge as well as calculated evaporation from the surface drainage network. Response tO pre- cipitation during the summer is an important indicator Of the sources of discharge. Greatly muted summer precipitation response relative to spring suggests that sea- sonal wetlands play an important role in generating peak flows. Notably, the flood 123 event of late summer 1986 is not well described in ILHM. This indicates that ground- water played an important role in the flOOd, as seasonal wetlands became inundated and contributed to a much larger runoff response. The one-way coupling currently implemented in ILHM cannot fully account for such an event. The broad base of runoff peaks are well described by this simulation, most Of which is groundwater dis- charge. This is another indication that baseflow, and therefore groundwater recharge, are accurately simulated. If T i I I I I T I T I 200 A Muskegon at Evart 100 a l l J . . . 'I'I " ' i I‘ ' ‘ 0\ oil .1 l . ' 1 l ' ll . . ‘ ‘ . l l . “.1 .t E I I I —I I I I 7 —I *I g, 40 B Clam at Vogel Center _ 2 _ .8 2° : D ' l l \ ."Vr',,\ . L» -- Ll. a\5h- _.,~ '- — ... £3: v x. E 0 ‘ (U 10 g ' ' ' 'Bear'near'Musk'egon' ' '—dage ‘ (0 ——Sim 5 —l ' _.__ _.._. _ -.. .5 .; h. —_ -. ...a _ a— —I—‘A~ ._ __ _ . ._ .‘l .— _ _‘,.‘.. .&-..A'~-< 1 80 1982 1984 1986 1988 1990 1992 1994 1996 19982000 2002 2004 2006 Figure 4.9: Simulated and Observed (gage) hydrographs for three USGS stream gages, A) Muskegon River at Evart, B) Clam River at Vogel Center, and C) Bear Creek near Muskegon. TO more quantitatively compare ILHM stream discharge predictions to gaged val- ues, Figure 4.10A plots simulated an Observed values using an unbiased residual cal- culated as (sim -— obs) / obs from the Evart gage. This figure demonstrates that ILHM predicts the bulk Of daily low-flow discharge to within +/-30%. This plot also shows the tendency for the model tO respond somewhat excessively tO summer precipitation events that tend to occur during low discharge periods. Similarly, the higher discharge periods that most Often occur during the spring are not fully captured. 124 As shown in Table 4.2, three synoptic stream discharge datasets were collected at or near annual low-flow levels across the MRW. These data are plotted in Figure 4.1OB, along with simulated model baseflow. Note that stream discharge varies across three orders of magnitude in this figure. In general, there is no significant trend in the data. Data collected during August of 2004 were within the central portion of the watershed. ILHM over-predicts baseflow during this period, suggesting either that the model over-predicts recharge or that the simulated water table relaxes tOO slowly during the growing season. Note a portion of the October 2004 and July 2003 data are under-predicted. This is likely due to small precipitation events shortly before the data were collected that continued to contribute some surface-originated discharge. In Figure 4.100, simulated hydraulic heads are compared to Observed static water levels. A best-fit line (not shown) deviates from the 1:1 line by less than 5%, indicating no significant head-dependent trend. The region between 230-260 m, where ILHM tends to underestimate head, is due to a sequence Of large impoundments along the main channel Of the Muskgeon River (Croton, Hardy, and Rogers dams) that are not included in the groundwater model. Other discrepancies between simulated and observed heads, such as a tendency toward over-prediction at the highest head levels, are likely due to incompletely calibrated hydraulic conductivity values. In general , ILHM simulates hydraulic heads within +/- 5 m, with little bias (Figure 4.10D). Here, the median annual head residual Sim—obs is shown as a heavy black line on top of a shaded region that outlines +/- 1 standard deviation in annual residual from the median. ILHM slightly under-predicts heads during the 19808, and slightly over-predicts heads during for the remainder of the model period. Overall, the mean residual bias is less than 0.20 m. A common quantitative estimate of model performance is the Nash-Sutcliffe effi- Unbiased Residual (%) Simulated Head (m) 300 A mi" E 2 . °°. a ‘8 U) (U m 13 2 .‘2 D 100 9' ' - g ’ o 100 200 300 ‘0 Observed Discharge (male) 400 -.c A E 350 a 1?: 300 '7) as 250 U .» — 3 200 ’f I 260 250 300 350 400 Observed Head (m) AJul2003 ,, .. B °Aug 2004 30 ° nOct2004A 00 D A030“ :0 3‘ at: 0° 5 1 A 6° 43 “GA ° 1 0'3 10'2 10'1 100 1o1 Observed Discharge (ma/s) -+/— std D — Median —1G 19 80 1985 1990 1995 2000 2005 Figure 4.10: Plots of simulated and observed streamflows and groundwater head, Observation data are described in Table 4.2. A) Observed daily discharge at the Evart, MI gage vs. unbiased residual (sim — obs) /obs, B) observed discharge from 105 points across the MRW during low-flow conditions vs. simulated baseflow , C) observed groundwater head vs. simulated head, and D) annual median and +/— 1 standard deviation of the head residuals sz'm — obs. All horizontal and diagonal lines are reference lines, not linear fits. 126 n. aw . '1‘. 9'- gut . ciency E Nash and Sutclz'fie (1970), _ 221(01— a)? — 2 232:1 (0i — 0) where O,- are Observed values (Of a time series) and P,- are predicted or simulated (4.13) values. Essentially, E is a measure Of the deviation of predictions from Observed values normalized by the variance of the Observed dataset. Here, a log-transform Of the Observations and predictions are used in order to more evenly weight the efficiency metric toward all flow levels (K rouse et al., 2005). Values Of E > 0 indicate that the predictions describe more Of the variance in the Observations than the mean Of the observations alone. Annual E values Of log-transformed daily stream discharge are plotted in Figure 4.11 for each of the six USGS gages. The three gages with the largest catchments, Evart, Croton, and Newaygo, have efliciencies in the range -0.2 to 0.6, with most years near 0.4. The Clam and Little Muskegon River gages, particularly from 1994 - 2007, generally have negative efficiency values. This is because of an apparent overestimation Of groundwater recharge within these basins. The gage on Bear Creek, with the smallest catchment, has positive efficiencies for much Of the model period, which is notable because random errors in model parameters and climate inputs can be much more significant at this scale. Plotted in Figure 4.12A are water-year (Oct-Sep) total discharge flux error values (Sim — obs) normalized by the area Of each catchment, and expressed in cm. This flux error is directly comparable to total precipitation, which averages approximately 85 cm across the watershed. Thus, an 8.5 cm flux error in stream discharge is approx- imately 10% error of the total water budget. For most of the gages, excepting the early period in Bear Creek, variability in flux error is less than 10 cm annually, and most Of the gages fall within 10 cm most years. In particular, the larger catchments Of MR at Evart and Croton are simulated more accurately than the smaller catchments, suggesting that random errors either in input fluxes or model parameters account for 127 i ’6 0.5 * /% L g 3' . f l' i . ' 1 «:2. ° « ‘ ' l V WU m... ms... t *? —Clam g —1 _ —Litt|e MR / . 5 —MR Evart l [41—1 5 —MFi Croton I —MR Newaygo ’2 1985 1990 1995 2000 2005 Water Year Figure 4.11: Plot Of water-year (Oct - Sep) Nash-Sutcliffe efficiency Of the log- transformed simulated and gage streamflow values for all six USGS gages. “MR” refers to “Muskegon River”. much Of the error in smaller gage basins. During the early simulation period, Bear Creek (Figure 4.9) responded to precip- itation events with rapid hydrograph rise. This has, over time, lessened. This may be to the requirement that new development be accompanied by retention ponds or infiltration basins tO reduce peak runoff response. As more Of these were installed, the system responds more slowly to precipitation, and with less intensity. This kind Of time-varying land use eflect can be difficult to incorporate in a regional scale model, and is a source of error across the basin. Average monthly flux errors at each gage (Figure 4128) exhibit the same general behaviors discussed for Figure 4.9; spring flows are under-predicted, and summer flows slightly over-predicted. For all gages, except the Bear Creek gage which is likely influenced by anthropogenic water withdrawals near Muskegon, the general shape of monthly flux errors are similar. This implies that the processes generating monthly error variability are relatively uniform across the model. Indeed, they correspond quite closely if Oflset with a scalar constant at each gage. The gage-dependent bias is likely related tO spatially distributed model parameter values, though no parameter sensitivity or calibration was conducted here. 128 ANOD 000 1 .l Annual Error (cm) 0 -10. -20 -30 i A n A . 1985 1990 1995 2000 2005 Water Year 2 r a B —L O I —L I —— Little MR _2_ —MR Evart 1 —MR Croton _3 . g.— MR Neway90_ Average Monthly Error (cm) JEMAMJJASOND Calendar Months Figure 4.12: Plots of stream discharge error sim - obs for all six USGS gages normal- ized by the area of each gage catchment. Plotted in A) are water-year error values, and B) are average monthly error values across the 27 water years in the simulation. 129 Over the nearly three decade model period, the influence Of time-variant inputs on model error is reduced. Table 4.5 lists daily efficiency values for the entire period, along with total flux errors at each gage. Of the gages active during the entire period, the largest flux error Observed is 7.9 cm, under 10% Of the water budget. The larger catchment areas, including the MR @ Evart and Croton, had flux errors of 1.7 and 3.1%, respectively. This level Of error compares favorably to other methods of estimating regional water budgets (Scanlon et al., 2002; Arnold and Allen, 1996; Winter, 1981). The total daily efficiencies are similar to values reported elsewhere for uncali- brated distributed rainfall-runoff models. The first Distributed Model Intercompari- son Project evaluated 12 different hydrologic models for 8 catchments ranging between 65 and 2484 km2 (Smith et al., 2004). Model efficiencies varied among models and basins, but in general, larger basins had better efficiencies. Average uncalibrated model performance for large basins was generally less than 0.5, and in some cases negative. Calibration improved all model efficiencies, generally increasing them by 0.1 to 0.2 (Reed et al., 2004). Table 4.5: Total Nash-Sutclifl'e efficiency and average annual flux error values for each of the 6 USGS gages. Efficiency values are calculated using log-transformed simulated and gaged values. Average annual flux error is sim — obs divided by the catchment area Of each gage. The final column is average annual error normalized by average annual precipitation for that catchment. Gage Area Avg. Precip. E Avg. Error Avg. Error (kmz) (cm) log-transformed (cm) (70 precip) MR @ Newaygo 6,086 81.1 0.43 -4.43 -5.46 MR @ Croton 5,991 85.4 0.43 1.24 1.45 MR @ Evart 3,711 79.0 0.45 -4.16 -5.27 Little MR 894 89.2 -0.32 12.01 13.47 Clam 629 82.4 0.11 3.69 4.48 Bear 43 87.6 0.25 -8.12 -9.26 130 5 4.5 Conclusions The Integrated Landscape Hydrology Model (ILHM) is capable Of fine-resolution, large-domain simulations. It is a fully distributed, process-based integrated hydrologic model designed to be computationally eflicient and highly flexible. While there are numerous watershed and integrated hydrologic modeling codes, few are capable of fully-distributed simulations over large domains. Originally detailed in (Hyndman et al., 2007), the code has been significantly improved. It now incorporates a full range of hydrologic processes covering upland and wetland domains. Currently, the coupling between surface and deeper subsurface domain modules is one-way, though full run-time coupling is under development. ILHM is structured to handle both gaged and fully-distributed input data sources at a variety Of resolutions. Heterogeneous sources can be included to provide con- tinuous records over many decades. Extensive efforts were undertaken to provide a consistent, quality-controlled climate and landscape input dataset. Coarse cloud cover Observations along with a sparse set Of total solar radiation data were modeled to provide a temporally- and spatially-variable real-ground solar radiation beam and diffuse components. Methods for compositing gage and radar precipitation data are detailed, as well as a discussion of proper spatial distribution of gaged datasets. Leaf area index data that describe phenological development within each model cell came from three separate sources, including two sets of data that provided only vegetative indexes. These index values were transformed to LAI with locally-Optimized empir- ical functions. Following this, a land use-based downscaling method was detailed to provide higher resolution LAI estimates within each coarse data cell. Comparisons Of uncalibrated ILHM predictions of stream discharge and ground- water heads tO Observations indicate that the code simulates most of the important hydrologic processes, and accurately describes fluxes both regionally at locally. Gen- eral hydrograph shapes are well-represented, with tendencies toward under-prediction 131 of peaks in the spring and over-prediction Of peaks during the summer months. Dur- ing all periods, baseflow is reasonably estimated, with regionally-dependent over- prediction present at some gages during low—flow periods. Simulated hydraulic heads display no general biases, with modest residuals at this spatial scale. Total model flux errors are small, even for this uncalibrated simulation. Based on these results, the current ILHM process modules are well suited for sim- ulating unconsolidated, loosely-coupled surface— and ground-water systems. Other environments will dictate modification or development of new process modules, which the code is structured to allow relatively easily. Moreover, multi-decadal simulation periods can be run on modest desktop machines, though the surface and deep unsatu- rated zone modules are highly parallel, and the code is written to allow subdivision to many nodes on a cluster. The need to calibrate very few domain-specific parameters allows ILHM to be readily and rapidly applied to ungaged and remote basins, with reasonable confidence in predictions. Acknowledgments Funding for this work has been provided by the Great Lakes Fisheries Trust Muskegon River Initiative, the National Science Foundation through its Water Cycles program (EAR-0233648), and the Michigan State University Center for Water Sciences. We Would like to thank former students including Cheryl Kendall, Jefl Eggleston, Nick- laus Welty, and Dushmantha Jayawickreme for their assistance in collecting stream discharge data and preparing input datasets, including the digitization of inland lake maps. 132 CHAPTER 5 Influence Of Land Use, Climate, and Soils on Recharge and Evapotranspiration Across a Regional Watershed 5.1 Introduction Climate and land use change are altering the timing, distribution, and quantity Of hy- drologic fluxes worldwide. During the last century, air temperatures (NCDC Temp, 2008), precipitation (Center, 2001), and stream discharge (USGS, 2008) have in- creased in the Muskegon River Watershed (MRW), a regional catchment in central lower Michigan, USA. Over the same period, land use shifted dramatically; forests QXpanded to cover much Of their pre-settlement range, urban areas expanded consid- erably, and the extent of agricultural declined by half Piajnowski et al. (2007). Predicting the heterogeneous local and regional responses to these changes is critical to managing water resources for both human uses and natural ecosystems. A relative abundance Of shallow groundwater in Michigan is a valuable resource for community and agricultural water supplies, but also provides the bulk Of annual stream discharge in this baseflow—dominated system . Numerous segments of stream channels within 133 the MRW provide habitats for economically-important cold-water spawning fishes due to the robust supply Of groundwater discharge during summer months (O’Neal, 1997). One of the primary consumptive water uses in many basins, irrigated agriculture, is not currently widespread in the MRW (Morenz et al., 2005). Spring snowmelt provides a reliable source Of soil moisture storage for the early growing season, complemented by relatively wet months later in the growing season. Furthermore, much Of the surface soil textures across the MRW are sands or loamy sands with little plant-available water capacity. These marginal soils are farmed in areas of low relief, but are not well suited to crop systems such as seed corn that generate demand for irrigation in this region (Morenz et al., 2005; USDA-NRCS, 2008). During the 21St century, temperatures across the region are forecast to increase between 3 and 4 QC under the A1B emissions scenario, while annual precipitation may increase 10% (Christensen et al., 2007). Higher temperatures would hasten the arrival Of the growing season, and extend the period Of moisture extraction from soil storage. Despite a forecast annual increase in precipitation, climate model ensembles show no change in precipitation during the growing season. Increased demand for moisture with a static supply may lead to drier soils throughout the region and threaten the viability of some non-irrigated farming practices. Recent demand for biofuel crop production may provide economic incentive for marginal lands to shift crop systems, or perhaps be put pack into production (NASS, 2007). Expanded production Of corn and soybeans into sandier soils may increase the demand for irrigation to compensate for poor soil moisture storage capacity. Biofuel extraction and processing facilities are also major water-users, consuming 40 L of water per L of corn ethanol produced (Pimentel and Patzek, 2005). Aside from shifting the consumptive demand for water resources, land use and cli— mate change will alter the water budget Of the region. Because of coarse-textured soils 134 an associated high infiltration capacities, most Of precipitation is either evapotran- spired (ET), or enters the groundwater as recharge and is subsequently discharged to surface water systems. Impervious surfaces in urban areas Of the MRW increase runoff, and alter the balance between ET and groundwater recharge in upland ar- eas. The succession of abandoned agriculture to shrub, and ultimately forests, also modifies the water budget through increased canopy interception, decreased soil evap- oration, and other changes. Effective management Of existing water resources and the sustainability Of land use practices requires a quantitative understanding of the interactions among climate, land use, and terrestrial hydrology. These interactions are facilitated and regulated by various aspects Of the landscape, including soil textures, relief, and surface wa- ter bodies. The complexity of interactions among these factors require more than Simply a conceptual or theoretical understanding of the impacts. Furthermore, the heterogeneity of landscape, land use, and climate demands a predictive capacity that addresses all scales, from hectare tO basin. The Integrated Landscape Hydrology Model (ILHM) is a fully-distributed hydro- logic model capable Of fine-resolution simulations over large-domain. It consists of four domain modules for surface water, canopy and root zone, deep unsaturated zone, and saturated groundwater (Hyndman et al., 2007). Within those domains, sets Of individual physical process models simulate hydrologic fluxes such as evaporation, transpiration, infiltration, groundwater recharge, overland flow runoff, and ground- Water flow. ILHM is written to run on desktop computers, but be rapidly ported to high performance clusters if necessary. Chapter 4 described the improvement of ILHM since the publication of Hyndman et al. (2007) and the development Of new Capabilities including more sophisticated wetland process modules. Also detailed in Chapter 4 are a set Of methods to develop spatially and temporally consistent input datasets to drive a fine-resolution, regional—scale model. 135 L'_) I Uncalibrated predictions of streamflow and groundwater head for the MRW com- pared favorably tO several different sources Of observation data. Over the simulation period 1980 - 2007, the model simulated total stream discharge to within 5% Of pre- cipitation for the largest basins, and between 5 and 15% for smaller catchments. N ash-Sutcliffe efficiencies, a measure Of the ability Of the model simulate Observation dataset relative to its variance, were generally around 0.45 for the larger basins, and between -0.3 and 0.25 for smaller watersheds. Nearly 50,000 measurements Of static water table levels, collected by drinking water well drillers, were used to determine the accuracy Of hydraulic head predictions within the model. The mean head residual was less than 20 cm. This Chapter examines the spatial and temporal variability of key hydrologic fluxes including groundwater recharge and evapotranspiration. Precipitation and temper- ature vary significantly over the model domain, due largely to the effects Of nearby Lake Michigan. Imprinted on these regional scale trends are land use and soil texture variability at the hectare scale. ILHM is capable Of resolving variability across these scales, which is necessary to quantify the influences Of climate, land use, and soil textures on the terrestrial hydrologic cycle. 5.2 Methods and Study Site The Muskegon River Watershed (MRW) is within the Lake Michigan/ Huron basin, part of the Laurentian Great Lakes (Figure 5.1). The MRW covers approximately 7,400 km2, and encompasses several large lakes (Houghton, Higgins, Mitchell, Cadil- lac, Muskegon) and wetland complexes. The main stem Muskegon River is a fifth order stream, with a single fourth order tributary, the Little Muskegon River. Stream flow within the MRW is dominated by groundwater, which contributes approximately 85 percent of mean annual discharge (see Chapter 4). Its major urban centers are 136 Muskegon (pop. ~40,000), Big Rapids (”11,000), and Cadillac (710,000). The MRW has relatively little relief, with elevations between 175 and 500 meters. During the late 19th century, the MRW underwent a rapid period of deforestation leaving a virtually denuded landscape by the late 18903. By the first decade of the 1900s, less than 20% of the watershed was forested, while nearly 50% was managed agriculture (Piajnowski et al., 2007). Subsequent abandonment of marginal agricul- tural lands led to significant reforestation with only moderate urbanization during the remainder of the 20th century. -90° -86° -82° —78° Figure 5.1: Map of the Muskegon River Watershed and surrounding model domain boundary, with an inset showing the location Of the model domain within lower Michigan, USA. Filled symbols on the lower map are climate gage locations for data listed in Table 5.1. The current land use within the MRW is predominantly forest (58%), with a sig- hificant fraction of agriculture (17%), some shrub or Open lands (9.5%) and urban 137 (7.4%), and a significant portion of Open water or wetlands (8.4%) (Torbick et al. (2006), see Figure 5.2). Soil textures from SSURGO data (USDA-NRCS, 2008) are predominantly sands (61%) and loamy sands (19%), with considerable areas of sandy loams (7.2%), loams (6.2%) and finer textures (6.8%). Outside the MRW, soil textures are somewhat finer, particularly in the southern model domain (see Figure 5.2). Land Use - Major Land Use - Minor Soil Texture Average Precip. (cm)- [:1 Agriculture B Urban E] Sands <- E! El-> - Deciduous E] Shrub B Loams 74 79 84 89 94 - Coniferous - Wetland - Fines . l D Figure 5.2: Maps of select ILHM input datasets. Sub-figures A and B show distri- butions of major and minor land use classes, C illustrates the dominant soil textural classes, and D plots the average annual precipitation across the model domain. 5.2.1 Input Climate and Landscape Data The climate Of the MRW is temperate and humid, classified according to the KOppen- Geiger scheme as Dfa (snow climate, fully humid precipitation, and warm summer temperatures) ( K ottek et al., 2006). Figure 5.3 plots monthly average temperatures, precipitation, solar radiation, and leaf area index (LAI) for the 1980 - 2007 simula— tion period. Average daily temperatures are typically below freezing from December- March, which typically results in a persistent seasonal snowpack that stores approx— imately 10 cm Of water. That snowpack usually melts between late March and mid April, when daily maximum temperatures rise above freezing (Figure 5.3A). Winter is typically the driest season, and summer the wettest, with significant potential for very wet months during the spring and autumn (see Figure 5.3B). Frontal 138 systems moving over Lake Michigan generate lake effect precipitation, particularly in late fall and early winter. Lake effect enhances average annual precipitation in a band near the lake by 10-20 cm annually (see Figure 5.2). Monthly average shortwave solar radiation varies by as much as 20% from year to year during the spring and fall, with a smaller annual range during the winter months (Figure 5.3C). Watershed-average leaf area index values are greatest in forested areas (peak growing season N5), and are lowest in agricultural areas (peak growing season N2, and vary typically less than 10% year-tO-year (Figure 5.3D). Solar Rad. (W/m2) Precip. (cm) Air Temp. ('C) g a " ¥--iorest ',...:T I? A ——0 en - ~L ufban.’ ‘.\ _ —a / — - S 2.5 9 ’7'". ...... >x..\ /‘/ \'\ “J FMAMIJ JASOND Calendar Month Figure 5.3: Monthly values Of select climate and leaf area index (LAI) inputs. The heavy solid line is the median monthly value during the simulation period, which is surrounding by a shaded range indicating the 5% and 95% values. Plotted in A) monthly temperatures, B) monthly precipitation, C) monthly average solar radiation, and D) monthly average LAI for four land use classes. 139 Ea" ml pk- C" Significant effort was invested tO create a composite climate dataset from heteroge- neous sources while maintaining spatial and temporal consistency (see Table 5.1 for data sources). Quality assurance, gap filling, and data consistency procedures that were used are detailed in Chapter 4 Of this thesis. TO spatially distribute the data, Theissen polygons were constructed for each unique combination of on-line gages. Weekly precipitation, aggregated from daily data at over 60 stations, was used to re- weight hourly precipitation values available at only 9 stations. Distributed datasets, including LAI and NEXRAD precipitation data, were bilinearly downsampled tO the model resolution. Table 5.1: Input climate and LAI data sources with the start and end dates of each source. Locations of gage data are shown in Figure 5.1. Source Type Start End AWOSl gaged climate 1980 2006 MAWN2 gaged climate, solar radiation 1996 2007 N CDC3 gaged precipitation 1980 2007 N EXRAD4 radar precipitation 1996 2007 GIMMS5 8-km NDVI 1981 2006 AVHRR6 1-km NDVI 1989 2003 MODIS7 LAI (v5) 2000 2007 1(Andresen, 2007) 2(MAWN, 2008) 3(Center, 2001; NCDC Hourly, 2008) 4(Andresen, 2008) 5(Tucker et al., 2005) 6(Eidenshink, 1992) 7(Knyazikhin et al., 1999) 5.2.2 The Integrated Landscape Hydrology Model ILHM is loosely-coupled set of fully-distributed physical process models, grouped into four domains: 1) surface water storage and routing, 2) canopy and root zone, 3) deep unsaturated zone, and 4) saturated zone groundwater. Its internal parameters are virtually all physically-measurable, and generally accessible through literature, requiring no site-specific calibration for reasonable flux estimates. This makes it well suited both for application in ungaged basins and for calibrated modeling within highly instrumented catchments. 140 Individual domain modules within ILHM may have different spatial and temporal discretization schemes. Currently, all domains use structured grids. For this study, the surface routing, canOpy and root zone, and deep unsaturated zone domains share the same hourly timesteps, with square grid cells 425 meters on a side. The saturated zone module is discretized at 1/4 the grid size Of the surface domain, or approximately 106 meters, and is run with daily time steps and weekly input stress periods. The root zone is divided into 7 vertical layers, varying in thickness between 5 and 70 cm. Deep lakes have 10 vertical layers, between 30 cm and tens Of meters thick. The deep unsaturated zone is a single layer that tracks multiple discrete wetting fronts. Three vertical layers are used for the MRW saturated groundwater model, ranging between 3 meters and several hundred meters thick. Neighboring cells do not necessarily have identical vertical layer structures. Each cell in ILHM is composed of three “landunits”, uplands, wetlands, and streams. The fraction Of each cell covered by the landunits are Fup, Fweta and Fstr- Within each landunit, hydrologic fluxes and storage volumes are independently maintained, and the total flux in a cell (for shared flux types, such as deep percolation) is calculated as Q = QupFup‘l'QwetFwet +QstrFstr- Only the root-zone domain module is subdivided in this fashion. For this study, an 18,900 km2 model domain of was established around the MRW. The boundaries Of this domain extend outward to significant hydrologic divides, pri- marily the main channels of adjacent similar-order watersheds. The base Of the shal- low unconfined aquifer, a contact with bedrock, occurs at depths between 15 and 400 m. This contact was interpolated from thousands Of Oil and gas wells. Within the model, the root zone domain extends to 2.5 m depth or shallower if the water table (taken as the average annual low level from preliminary modeling) is closer tO the surface. Below this, the deep unsaturated zone domain extends to the water table that defines the upper limit Of the saturated groundwater domain. 141 5.2.3 Climate Teleconnection Indices Numerous studies have identified correlations between low-frequency climate oscilla- tion, or teleconnection, patterns and precipitation or temperature in the Great Lakes region (Rodionov and Assel, 2000; Rajagopalan et al., 2000; Assel, 1992; Leathers et al., 1991). Some of the most influential such patterns in the Northern Hemisphere are the North Atlantic Oscillation (NAO, Visbeck et al. (2001); Walker and Bliss (1932)), Pacific/ North American (PNA, Wallace and Gutzler (1981)), Arctic Oscilla- tion (A0, Thompson and Wallace (1998)), and El Nino Southern Oscillation (ENSO, Trenberth (1997)). These patterns are summarized numerically via indices measuring the departure from normal Of geopotential heights within certain geographic zones (NA 0, A0, or PNA), or sea surface temperatures (ENSO). These index values are typically normal- ized such that monthly average values are small, typically within +/- 2 for indices used here. For this study, monthly index values were obtained from the NOAA Cli- mate Prediction Center. Note, the ENSO index was calculated from SST anomalies in the so-called 3.4 region, NAO and PNA using 50 hPa height anomalies, and A0 using 100 hPa anomalies. TO understand the influence Of the climate indices on fluxes in this region, cor- relations between monthly, seasonal, and annual index values and precipitation and temperature were evaluated. The observed correlations are discussed below, all of which are consistent with other studies. This study examines correlations between those indices and the resultant landscape hydrologic fluxes, particularly evapotranspi- ration and groundwater recharge. Understanding these linkages will help tO further our understanding Of the processes that influence the variability Of terrestrial hydro- logic fluxes. 142 5.2.4 Linear Regression Modeling Multiple linear regression modeling is applied to extract additional information from the complex temporal and spatial variability present in the ILHM outputs. Here, annual cell-by-cell input variables are correlated to flux estimates from ILHM to produce models capable Of estimating the spatial distribution Of those annual fluxes. A nested linear modeling approach was used along with variables including precipitation and temperature during both the growing and non-growing seasons, field capacity Of the soils, fraction Of the landunit covered by each land use class, and the annual average value of the PNA and NAO indices. Linear model inputs are annual or seasonal cell-by-cell averages, except for annual mean climate index values that are spatially uniform. Due to memory limitations in the 32-bit version Of the statistical package R (R Development Core Team, 2005) that we used, 30,000 randomly selected ILHM cells were used for the regression analyses. To determine the quantitative improvement Of each successive variable addition, the root-mean square residual and coeflicient Of determination, R2, between the an- nual linear model and the ILHM estimates were calculated. To compare the relative influence of each variable on the linear model, the model coefficient was multiplied by the range of values expressed by each variable. In this fashion, the maximum and minimum contributions to the linear model for each variable can be directly assessed. 5.3 Model Results and Discussion Cell-by-cell monthly fluxes Of water budget components were saved from ILHM for the 1980 - 2007 model period. The discussion below examines the spatial and temporal variability of those fluxes, to better understand the expected influences of climate, soil texture, and land use on terrestrial hydrology. In addition to direct spatial and temporal aggregation Of the model results, multiple regression analysis was used to 143 extract information from the complex ILHM outputs. Model parameters, as discussed in Chapter 4, were largely derived from literature values. With the exception Of preliminary calibration Of hydraulic conductivity values with a steady state groundwater flow simulation, the model parameters have not been adjusted to better match measured values. This uncalibrated simulation displayed no significant net basin-wide bias is stream discharge, although individual catchments had annual errors ranging from -4.3% to 13.6% Of precipitation over the 28 year model period. A wide range Of model parameters could be estimated to reduce this bias and mini- mize the residuals between simulated and Observed values, however here we choose to present uncalibrated results for several reasons. Residuals from a process-based model inform the modeler Of improvements that can be made to the sources or preparation Of input data, process modules, and physical parameterization of a system. This study seeks to better understand the interactions among hydrologic processes and climate and landscape inputs. Calibration in place of properly simulating physical processes can mis-represent the role Of individual processes within the overall system. Finally, by presenting uncalibrated results and explicitly acknowledging the model bias, the reader can make a more honest assessment Of ILHM’s capabilities tO predict flows in ungaged basins as well as for application to other problems. The results below focus primarily on the near-surface processes that drive spatial and temporal variability Of hydrologic fluxes in the MRW. Because of this deep per- colation (DP), or the flux Of water below the root zone, is presented rather than groundwater recharge to the water table. In the coarse-textured, high-conductivity soils of the model domain, the near-surface water budget is dominated by DP and evapotranspiration (ET). Note also that all annual fluxes are calculated for 26 com- plete water years, from October - September (e.g., Oct 1980 - Sep 1981 is the 1981 water year). 144 5.3.1 Spatial and Temporal Variability Of DP and ET Figure 5.4 shows the average cell-by-cell annual ET and DP along with the coefficient of variation Of the annual fluxes. The coefficient Of variation q, is a measure of the inherent variability Of a dataset, and is calculated as cv = arm/i, where on; is the standard deviation Of the dataset, and Li: its mean. Note that, due to lateral subsurface flows, the sum Of ET and DP may exceed precipitation. ET (Figure 5.4A) displays no significant regional trends, with the exception of high values in the upper part Of the MRW due tO a concentration Of wetlands. Large lakes and concentrated urban areas have lower ET, as does the broad floodplain Of the Muskegon River main channel. Areas Of higher ET at the 10—50 km scale (see Figure 5.1 for a scale) tend to have finer-textured soils (see Figure 5.2). Patterns in mean ET due to land use are more subtle, and relative relationships between classes depend on underlying soil textures. By contrast, deep percolation (Figure 5.48) has a strong regional trend. From the southwest - northeast, DP decreasing by 50% within cells Of similar land use and soil type. The underlying factors driving this regional trend are discussed in greater depth below, but the broad spatial pattern Of DP is indicative Of lake—effect precipitation during the fall and winter. ET responds much more directly to the more spatially- uniform growing season precipitation. Certain wetland classes also produce significant DP, including the complex in the northeast corner Of the watershed. Coarse-textured soils and sediments, and concentrated urban areas also lead to higher DP. In general, annual DP is more variable than ET (Figure 5.4C&D); model average c1, Of DP is 32.4%, and ET is 14.0%. The DP cu map exhibits significantly more vari- ability at the 5-10 km scale than ET as well. ET cu exhibits some geometric patterns that are artifacts of the nearest-neighbor spatial interpolation Of solar radiation and temperature. Regions with high cv Of DP tend to have low cv Of ET. In uplands, deep percolation 145 Mean Modeled ET (cm) Mean Modeled DP (cm) ----’ <----> 1530456075 1525354555 “ Coefficient of Variation V 7 Coefficient of Variation Annual ET (%) Annual DP (%) <----> <----> 11 131619 22 22 27 32 37 42 Figure 5.4: Maps of total cell mean annual modeled A)evapotranspiration (ET), B) mean annual modeled deep percolation (DP), C) coefficient of variation (ratio of standard deviation to the mean) of annual ET, D) coefficient Of variation, on of annual DP. Note the color scales differ for each sub-figure. 146 occurs only when root zone soil moisture exceeds field capacity; whether and to what degree this occurs is controlled largely by spring and autumn rains, which vary significantly inter-annually. Much Of ET is determined by the storage Of moisture in the soil root zone, nominally the difl'erence between the field capacity and wilting point of a soil. Most years, soils are near or at field capacity at the start of the growing season, leaving summer precipitation as the remaining factor governing inter-annual variability for ET. Regions with low cu Of ET tend to have finer textured soils, and thus greater soil moisture storage, reducing dependence on summer precipitation. These same regions also require more soil moisture to exceed field capacity, decreasing the likelihood Of significant DP occurring during any given year. TO further examine the temporal variability Of ET and DP fluxes across the MRW, Figure 5.5 plots annual watershed-averaged fluxes. In Figure 5.5A, total watershed ET, storage and runoff are plotted. ET dominates the watershed water balance, which is typically 60-70% Of precipitation (Figure 5.5C). The annual watershed runoff time series is quite steady, compared to the variability in other water budget components. Steady annual watershed runoff is facilitated by storage Of water within the water- shed, primarily in the saturated groundwater system. Storage is the most variable component Of the MRW water budget, and is capable Of storing and releasing as much as 25% of annual precipitation in any given year. Discussing longer-term trends within a time span of just a few decades can be difficult, particularly when the model input data sources change during the simulation period (see Chapter 4 Of this thesis). Nonetheless, simulated watershed ET seems to be declining steadily, while runoff is increasing. Some Of this apparent trend may be driven by the years 2004 - 2007, which are the four highest runoff years in the simulation period. As a fraction Of the total water budget, runoff seems tO be increasing as well. The water budget of upland areas is plotted in Figure 5.5B & D. The only two 147 Annual Flux (cm) % of Precipitation Watershed Average: Upland Root-Zone Average: -——-ET ----Runoff —ET —Overland Flow ———-Storage —Storage —DP 75 a . . . . . . . A ,-. 60' B ’ ‘—\ . 'n‘ I ‘Tl "_‘~ \t’.\ 1’“ so " ~ g” V - 45 \ ¢’\\ /\ 1’ \""\,~—— 1“], 30 25,” - ~ ~ [l I‘\ III I“ A Ill| A 1" ‘., l\‘ 1"v\ 0 v "| I; -' \ J' ‘U’ I \ ’l ‘1’ -25 100 - r- C ., 75 ' 1‘. r“ I“ l 5’-‘ ’z‘ I /~’ I“ ‘4' ‘4' . I ‘ ‘. ‘iv/a_. 50* .. ’ I; I‘ I". ,‘ ’1‘“ 15- G 'I \‘V'I “ ”‘\" \\\I 1‘ I V ‘ T l 5 A AA A A A A 1", v \ \.I U V 'V \I VV V \ —25..‘«.. . _,5 n..n.-n. 82 85 88 91 94 97 00 0306 82 85 88 91 94 97000306 Water Year Water Year Figure 5.5: Plots of annual watershed-averaged fluxes for the MRW. On the left, watershed-average ET, storage, and runoff for the MRW are plotted, while on the right upland ET, overland flow, root-zone storage, and DP within upland landunits are shown. Parts A & B plot annual fluxes in cm. Those quantities, divided by annual precipitation, are plotted in parts C & D. Note, for storage, negative values indicate release of water from storages in the model. 148 dominant components at the watershed scale are ET and DP, runoff is a small fraction except in urban land uses. Storage within the root zone is an important component of the annual variability in fluxes, but is minimal long-term. Within uplands, ET remains remains the largest component of the water budget, though deep percolation exceeds ET in some years. Over the entire model period, ET makes up 53% Of the water budget, and DP 43%. The 2004, 2006, and 2007 water years experienced the greatest DP during the 28 year model period, as a result Of heavy spring and fall rains. In absolute terms, DP is annually more variable than ET, though the Opposite is true as a fraction Of precipitation. Strong 2-3 year cycles are evident in both ET and DP that may be indicative Of the influence Of broader climatic patterns. Table 5.2 summarizes a nested linear regression model of annual watershed-averaged upland ET and DP with annual means Of the four climate teleconnection indices PNA, NA 0, ENSO, and A0. Four models were tested, each adding one additional index. The ordering was chosen based on single regressions for each dataset, with highly correlated indices included earlier in the nested linear models. The final five years (water years 2003-2007) were excluded from this analysis because Of behavior that is significantly diflerent from the remainder of the model period. A model including just NAO + PNA explains between 15% and 24% of the vari— ability in average annual upland DP and ET, respectively. The addition Of ENSO has little effect on either model, while adding A0 seems to significantly increase overall model fit. The significance Of this, however, is suspect as the total number Of data points (22) may be over-parameterized by a four-component linear model. Note, since the range Of index values is approximately equal for all indexes, the values of coefficients from the models indicate the relative influence Of each parameter. The positive phase Of NAO (and the closely-related A0) brings warmer than aver- age annual temperatures tO the region, along with slightly increased average precipi- 149 tation. In addition to the continental effects induced by modifications to mean storm tracks, discussed in Harrell et al. (2003), increased average temperatures may pro- mote decreased ice cover and increased lake effect precipitation. According to these ILHM simulations, both annual upland DP and ET tend to be higher when NAO is in its positive phase. A positive phase Of PNA corresponds to reduced upland ET, but slightly increased upland DP. As Observed in Rodionov and Assel (2000); Leathers et al. (1991), winter and spring temperatures tend to be cooler during a positive PNA cycle, while sum- mer and fall precipitation is reduced. These lower winter temperatures prolong the seasonal snowpack. which leads to a more significant pulse Of meltwater infiltration in the spring. Reduced summer and fall precipitation limits ET significantly during a period when surface moisture storage is typically at or near annual minimums. Table 5.2: Tables Of multiple regression models for annual watershed-average A) Up- land ET, and B) Upland DP. Columns include the root-mean square residual between the multiple regression model and ILHM, along with the correlation coefficient R2, and the model parameters. Note, the final five years Of the simulation were excluded from this analysis because DP seems to behave exceptionally during this period com- pared the rest Of the simulation. A (Upland ET) Coefficients Model RMSR (an) R2 Intercept NAO PNA ENSO A0 NAO 4.44 0.20 46.3 6.88 + PNA 4.33 0.24 46.8 5.81 -293 +ENSO 4.31 0.25 46.7 5.68 -310 0.64 +A0 3.47 0.51 45.5 17.5 -4.87 1.07 -11.8 B (Upland DP) Coefiicients Model RMSR R2 Intercept NAO PNA ENSO A0 NAO 5.20 0.08 33.9 4.79 + PNA 5.00 0.15 33.3 6.34 4.24 +ENSO 5.00 0.15 33.2 6.32 4.22 0.08 +A0 4.90 0.18 32.8 10.9 3.53 0.24 -452 There is disagreement over which Of the A0 and NAO patterns better represent fun- damental atmospheric modes (Ambaum et al., 2001; Thompson and Wallace, 1998). In any case, the two indices are highly cross-correlated in this study (Table 5.3) and 150 others (Ambaum et al., 2001). NAO is virtually uncorrelated with either PNA or ENSO, whereas A0 is somewhat correlated with both. This strong cross-correlation between A0 and NAO partially explains the behavior Of the fourth model for Up- land ET in Table 5.2A. As expected, including both NAO and A0 alters the model coefficients of NAO because Of interaction between the two variables. Including A0 also increases the model R2 for upland ET dramatically, while having little effect for upland DP. Based on Tables 5.2 and 5.3, only PNA and NAO seem to Offer unique correlative information about upland ET and DP. Table 5.3: Cross-correlation values between monthly climate indices. PNA NAO ENSO A0 PNA | 1 |-0.003| 0.257 L-0.264] NAO |-0.003| 1 | 0.001 1 0.608J ENSO L0.257 I 0.001 | 1 j-O.125j AO [-0.264 1 0.608 L—0.125 ] 1 | Average monthly storage values provides an informative view of the interactions among model domain processes (Figure 5.6). From January through March, water is stored both in the snowpack and root-zone soil. At the onset of Of the growing season, a large pulse Of deep percolation (Figure 5.6B) balances extraction of water by ET. As this pulse of deep percolation and, eventually, groundwater recharge drains from the system, water is removed from groundwater storage until August. At this point, streams are near annual low flows, and ET has depleted most Of the soil moisture storage capacity. Fall precipitation at the end of the growing season then begins adding water to root-zone storage and deep percolation. On a monthly basis, upland ET and DP vary inversely (Figure 5.68). Most DP oc- curs following spring snowmelt in late March and early April, while a second, smaller peak occurs in October and November after leaf senescence. Deep percolation (and therefore groundwater recharge) is small during the winter months because the snow- pack stores precipitation, and during the growing season when transpiration and soil 151 evaporation consume most of the available moisture. Watershed runoff (Figure 5.6A) is the sum of groundwater discharge and overland flow runoff, minus evaporation from the drainage network. Direct overland flow runofl' accounts for a small fraction Of total runoff from the MRW. While still weighted towards higher runoff during the spring, storage in unsaturated and saturated zones dampen the sharp DP peaks significantly. Watershed Average: Upland Root-Zone Average: ----ET ----Runoff -—ET —Overland Flow ---~Storage ——Storage —DP 1 . ---. E 0 A ”I \ 3 ll \\ X I, \\ 2 I \ LL 5’ I, \‘ 2‘ [I \\ \ 5 1“\ \I” \\ C ~~~~~ 7 ‘ __________ /*:"'\ g >\—-—’::“‘ ’I \\_- Q) 0 \\ I, c» x ,I ‘3 \ ,I g \\ I, —2 < \ ,I __5 . . . . 7": 1 1 a . _4 1 a . . n . 1 . . n JFMAMJJASOND JFMAMJJASOND Calendar Month Calendar Month Figure 5.6: Average monthly fluxes Of A) watershed ET, runoff (overland + ground- water discharge), and storage, and B) upland ET, overland flow, root-zone storage, and DP within upland landunits. Evapotranspiration is the sum six transpiration and evaporation components, as shown in Figure 5.7, and Tables 5.4 & 5.5. Upland and wetland transpiration are the dominant components, followed by canopy evaporation, soil evaporation, snow- pack sublimation, and finally open water evaporation. Canopy evaporation during the winter months is due in large part to evergreen coniferous interception Of snowfall. Snowpack sublimation peaks in February because the cold, dry winter air has a high vapor pressure deficit. Open water evaporation during the late summer and autumn occurs primarily in large lakes that accumulated heat during the summer. This sim- 152 ulation may under-represent open water evaporation in these lakes, as modeled lake temperatures tend to be lower than Observed. -Open Water Evap. -Soil Evap. -Sublimation “Wetland Transp. ElCanopy Evap. -Upland Transp. 1o . . . . _ . . . . Monthly Flux (cm) OJFMAMJJASOND Calendar Month Figure 5.7: Filled area plot of watershed-average monthly evaporation and transpira- tion components for the MRW. Shown are upland transpiration, wetland transpira- tion, soil evaporation, canopy evaporation, snow and ice sublimation, and Open water evaporation. Listed in Tables 5.4 & 5.5 are watershed-average annual mean, standard deviation, and coefficient of variation, cu, of hydrologic flux components. Each table consists of three groups of rows. The first is annual mean precipitation, the second are the components of the water budget for the watershed (Table 5.4) and for upland lan- dunits (Table 5.5). The third group includes notable fluxes within ET and runoff components in the second group, ranked in order Of mean annual flux. Standard deviation, here, is calculated as the standard deviation of annual watershed—average or landunit—average fluxes. As in Figure 5.5, storage is the most variable component of the water budget apart from precipitation (Table 5.4), over the model period, total change in storage is minimal (0.3 cm). Most of this storage occurs within the saturated groundwater domain, and is controlled largely by the variability in annual deep percolation (Table 153 Table 5.4: Watershed-average annual mean, standard deviation, and coefficient of variation, cmof hydrologic fluxes. Flux Type Mean (cm) Std (cm) cv % Precip 84.6 12.4 14.6 Total ET 55.1 5.26 9.5 Total Runoff 29.3 4.64 15.9 Total Storage 0.26 8.60 3308 Baseflow 25.6 4.15 16.2 Upland T 23.2 2.36 10.2 Wetland T 14.1 1.66 11.8 Canopy E 8.67 1.48 17.1 Soil E 3.71 0.45 12.1 Snow and Ice Sub. 3.06 1.92 62.7 Wetland E 2.37 0.45 19.0 Wetland Run 2.73 0.57 20.9 Baseflow 0.96 0.18 18.8 Table 5.5: Watershed-averaged mean, standard deviation, and coeflicient Of variation, on, Of upland water budget components. Flux Type Mean (cm) Std (cm) on % Precip 84.6 12.4 14.6 Upland ET 46.3 5.10 11.0 Upland DP 37.6 9.95 26.5 Upland Run 1.15 0.21 18.6 Upland T 29.0 2.95 10.2 CanOpy E 8.8 1.15 17.1 Soil E 4.71 0.57 12.0 Snow Sub 3.84 2.40 62.5 154 5.5). With the exception Of snow and ice sublimation, the av values for all ET and runoff sub-components are between 10 and 20%. 5.3.2 Influence of Climate, Soils and Land use on Upland Hydrologic Fluxes Lake Michigan plays a major role in the hydrology of the MRW and the surrounding region. TO examine the effect Of proximity to Lake Michigan on upland water budgets, the distance Of each cell in the model to Lake Michigan along the prevailing wind direction was calculated. Wind rose diagrams from the MAWN gages (not shown) demonstrate that the prevailing winds are roughly from the SSW at an azimuth Of 60°. Figure 5.8 illustrates the lake effect on average annual precipitation, ET, and DP. Shown are monthly upland fluxes within 25 km distance bands downwind from the lake. Here, all cells containing greater than 75% deciduous forests, and sand as the dominant soil texture both horizontally and vertically are used. There are two distinct lake effect precipitation modes (Figure 5.8A), warm-lake and cool-lake Andresen and Winkler (2009). The warm-lake effect, when the lake is warmer than the air above, occurs most strongly in late September through early December. During these times, Lake Michigan significantly enhances precipitation near the lake shore, to approximately 150 km in the downwind direction. When the lake is cooler than the air above, which occurs from late May through early September, near-lake precipitation to 50 km is reduced as the cool water inhibits the formation of convective thunderstorm cells. The cool-lake eflect is visible in Upland ET (Figure 5.88) during the summer months. ET is roughly uniform from September - April. Then, the combined effects Of decreased precipitation, and slightly cooler temperatures near the lake reduce ET by roughly 10% within 50 km of the lake shore. Monthly Flux (cm) 012345678910 A Precipitation D N 53 C O A EM —A 8M F J B Upland ET Calendar Month L11§>§LL>CDOZU c Upland DP Calendar Month c.11§>Zr_L>CDOZU 0 _25 50 75 100125150175 200 225 Distance from Lake Mlchigan (km) along Prevailing Wind Direction Figure 5.8: Filled value plots of average monthly A) precipitation, B) upland evap- otranspiration (ET), and C) upland deep percolation (DP) vs. distance from Lake Michigan along the prevailing wind direction. All values are averages Of only cells within the model domain containing greater than 75% deciduous forests, with sand as the dominant texture. Values are averaged over the entire model period. 156 Upland DP exhibits both cool-lake and warm-lake effects (Figure 5.8C). From Octo- ber - December, and again in March - May, deep percolation is strongly biased toward the lake shore. Near-lake DP may be 50% greater than inland values during those months because Of the warm-lake effect. January and February do not exhibit strong trends. From June - September, near-lake DP is slightly reduced by the cool-lake effect, though not as significantly as ET during the same period. Although not specifically examined here, the role that Lake Michigan plays in the water budget of the MRW may be changing. Recent declines in lake ice cover (Assel et al., 2003) may extend the period Of the warm lake effect later into the winter. This may, in turn, increase winter precipitation and lead tO greater DP (Burnett et al., 2003). The cool-lake effect may decline as summer water surface temperatures warm more rapidly than the air above, as has been Observed over the last few decades at Lake Michigan tide gages (Austin and Colman, 2007). Because of the strong influence of proximity to the lake on upland fluxes, Figures 5.9 & 5.10 are calculated within the distance bands 75—150 km from the lake shore. That band was chosen tO provide relatively uniform ET and DP conditions according to Figure 5.8, while still providing enough variability in soil texture and land use classes to examine the influence Of those two variables. However, most Of the sandy loam and loam soils are in agricultural areas as these soils are more productive. Figure 5.9 plots upland ET and DP as a function Of the majority horizontal and vertical soil texture in cells with greater than 75% agricultural land use within the 75-150 km downwind distance band. TO reduce the impact Of vertical discontinuities in soil texture, cells with mismatched majority horizontal and vertical textures were not included. Finer-textured soils influence hydrologic fluxes in two ways: 1) reduced hydraulic conductivity slows the loss Of water out of the root zone as deep percolation, and 2) higher field capacity stores more water after the spring. This increased storage capacity is clearly evident in Figure 5.9.4. Finer-textured soils evaporate and transpire 157 much more water than coarser sands. Note that the fines class in Figure 5.9A includes soils with a wide range Of textures including silty loams, muck, and sandy-clay loams, so the behavior is less uniform than the other three texture classes. 'Upland ET -Agriculture 1 .5 h A, .5.5 ON °/o of W.Y. Precip N-hOJCD . Upland DIP -'Ag'riculture - B —— Sands — Sandy Loam , — Loam .5 o CD O) °/o of W.Y. Precip A N 0 1 1 1 1 1 1 1 1 1 1 J F M A M J J A S O N D Calendar Month Figure 5.9: Plots Of monthly average upland ET and DP across generalized soil texture classes as a percent Of water-year (Oct - Sep) precipitation. Values are from cells within the model domain containing greater than 75% agricultural land use between 75-150 km from the lake shore along the prevailing wind direction. Soil texture influences both the quantity and timing Of deep percolation, as shown in Figure 5.98. Less water is lost during spring snowmelt in finer textured soils. For textures except the fines class, DP during the early growing season is somewhat greater as silt and clay content increases because the reduced hydraulic conductivity in these soils slows transport of moisture through the root zone profile. As the growing season progresses, greater utilization Of water by plants again reduces DP in finer textures relative to sands. Regardless of soil texture, by the end Of the summer soil moisture is near the wilting point most years (not shown here). At the end Of the growing season, a greater amount Of water is required to restore fine-textured soils to their field capacity. This has the eflect Of dampening the autumn DP peak seen in the sandy soils. Table 5.6 presents the annual average upland ET and DP fluxes plotted in Figure 5.9. Compared to sandy-textured soils, finer soil textures can increase upland ET by as much as 30%, and decrease upland DP by almost 40%. Note here than the fines tex- ture class has the fewest members, though the differences shown are still statistically significant. This very strong dependence on soil texture highlights the importance of using higher-resolution (preferably vertically-variable) soil texture information avail- able from the digitized SSURGO data. Many studies rely on the Older, coarser, but more manageable STATSGO data ( USDA-SCS, 1993; Wang and Melesse, 2006) that lacks much of the detailed spatial variability present in SSURGO. Table 5.6: Average annual upland ET and DP by soil texture in agricultural areas between 75—150 km from the lake shore along the prevailing wind direction. Texture Upland ET (cm) Upland DP ( cm) Sands 47.6 37.3 Sandy Loam 52.8 32.4 Loam 58.5 27.6 Fines 62.0 21.9 range 14.4 15.4 TO examine the influence of land use on the terrestrial water budget, upland ET and DP are plotted in Figure 5.10 are as a function Of the dominant land use class within sandy soils (sub-figures A and B) and loam or sandy loam soils (sub-figures C and D). For each cell to be included in this analysis, at least 75% Of it had to be covered by a single land use class. Sandy soil textures dominate the model domain, and tend to be more vertically homogeneous than loam soils, so the the interrelationships among land use classes are perhaps more robust for the sands. As in Figure 5.9, only model 159 cells within the 75-150 km distance band are included. Land use affects the ILHM—simulated fluxes through a variety Of mechanisms, most prominently through canopy interception Of moisture and radiation, and through transpiration differences. Note, there are several factors not yet included in ILHM that may influence the land use relationships shown here, including sub-canopy mi- croclimate effects, soil temperature variation, ground litter, and differences among land uses in plant soil water extraction efficiency. Furthermore, the influence that a particular land use class has on terrestrial hydrology may vary depending on the texture Of soil beneath it (compare Figures 5.10 A & B to C & D). Within both soil textures, coniferous forests tend to evapotranspire more than other land use classes, with the exception Of sandy agriculture. This is due primarily to interception of precipitation during periods when plants in other land covers are dormant and their leaves senescent. Urban areas, in both sands and loams had the least ET and DP, due to several factors. ET in urban areas may be under-predicted due to the difficulty of retrieving urban LAI from remote sensing platforms. Even so, LAI in urban areas tends to be lower than shrub or forested land uses, as space is appropriated for structures and transportation. Furthermore, runoff in urban areas is significant due to impermeable surfaces, which for this simulation were assumed to cover 20% the urban fraction Of each cell. Next tO urban, shrub land uses have the lowest ET. This is perhaps because the shrub land use class includes rangeland and other managed vegetation that maintain persistently low LAI throughout the year (see Figure 5.3), while the resistance to transpiration in woody shrubs in particular is much higher than in agriculture crops. Low runoff in shrub land uses leads to higher DP than other classes. Agriculture and deciduous forests, the two dominant land use classes in the watershed, tend to fall in between other classes, with the exception Of DP in sandy-soil agriculture. Though a less important factor than soil texture, understanding the impacts of 160 Upland ET — Sands .A' .5.5 ON % of W.Y. Precip ON-hOim .5.5 ON % of W.Y. Precip ONAODCD A Calendar Month JFMAMJJASOND Upland DP - Sands Upland DP — Loams urban -—- agriculture -— shrub —— deciduous ‘ coniferous 011111411111 JFMAMJJASOND Calendar Month Figure 5.10: Plots Of monthly average upland ET and DP across land use classes as a percent Of water-year precipitation. Values are from cells within the model domain located between 75-150 km from the lake shore along the prevailing wind direction with either sands (A & B) or loams (C & D) as the dominant soil texture class both vertically and horizontally. 161 land use on hydrologic fluxes is critical because relative changes in fluxes are more important both ecologically and for water resources than the absolute values Of fluxes. Listed in Table 5.7 are annual average upland ET and DP fluxes across land use classes within sandy soils in the 75—150 km down wind distance band. At the extremes, land use in sand soils produces a range Of 7.8 cm for upland ET and 14.0 cm for upland DP. From Table 5.6, soil textures induced a range within agricultural areas Of 14.4 cm for both upland ET and DP. For the most part, the influence Of soils is important only spatially, while land use has both spatial and temporal implications. Table 5.7: Average annual upland ET and DP by land use in sandy soils 75-150 km from the lake shore along the prevailing wind direction. Land Use Upland ET (cm) Upland DP (cm) urban 40.3 31.2 agriculture 47.6 37.3 shrub 39.8 45.2 deciduous 44.7 39.2 coniferous 45.8 36.6 range 7.8 14.0 5.3.3 Linear Regression Analysis Of ILHM Outputs Traditionally, one of the few tools available for mapping the spatial and temporal variability of groundwater recharge and ET over regional domains has been linear regression modeling. The only available estimates Of spatially-variable groundwater recharge in the MRW were generated by correlating baseflow-separated recharge es- timates from stream gages with land cover, geographic, topographic, and climatic variables within catchments (MDEQ, 2009; Holtschlag, 1996). Two problems with these estimates are that 1) they rely on baseflow separation to determine groundwater recharge, a method that can be problematic in catchments with significant wetland area, and 2) variables that co-vary at scales smaller than the resolution Of catchment-based methods will inevitably confuse the linear models. An example of this second case is agricultural lands and soil textures. Within the MRW 162 and surrounding model domain, agricultural land use is highly positively correlated with field capacity Of the soil (correlation coefficient Of 0.46). Linear modeling of ILHM-calculated ET and DP suffers neither of those two prob- lems. ILHM calculates ET and groundwater recharge (or DP) directly, and at the scale of the ILHM model cells in this study, roughly 425 meters (cell area Of 18 hectares), there are sufficient instances Of all variables in all combinations within cells that the co—variation problem is greatly reduced. Therefore, linear regression models derived from the ILHM outputs can potentially provide broadly- and readily-applicable annual estimates Of ET and DP. Furthermore, these models provide additional information about the relative influence of each variable without needing to painstakingly isolate each variable as in the previous section. We refer tO this analysis as process-inferred statistical modeling. Here five nested linear models are calculated separately for upland ET and DP. The first includes only seasonal precipitation, the second adds seasonal temperatures. Cell-average soil field capacity values are added next, followed by within-landunit land use fractions, and finally climate index values. A sixth model, labeled “P” to indicate “parsimonious” has been calculated omitting seasonal temperatures and climate cycles. Note that the climate variables vary both in space and time (except the indices, which are annually single-valued), while the landscape parameters only vary spatially. As shown in Tables 5.8 and 5.9, each additional set Of variables Offered an im- provement in model fit, ultimately resulting in the linear models explaining 44 and 62% of the annual variability in upland ET and DP, respectively. The root-mean square residuals for DP are somewhat higher, due to the inherently more variable nature Of this quantity. The influence Of each variable relative to the mean field (the model intercept) is given in part C Of Tables 5.8 and 5.9, and is calculated as the product of the range Of the parameter and its model coefficient. The most influential 163 variable for DP by a significant margin is non-growing season precipitation, followed by growing season precipitation, field capacity Of soils, coniferous forests, and urban land use. The least influential are the climate indices NAO and PNA, and air tem— peratures. Upland ET is most influenced by soil field capacity, followed by growing season precipitation, forest, and agricultural land use. Here, non-growing season air temperature, urban land use, and the two climate indices were the least influential. The primary limitations of a linear model is that the system may respond non- linearly to input variables, and threshold responses are difficult to describe. Such responses are evident in the ILHM outputs, particularly for soil texture and pre— cipitation. Additionally, the variables may interact with each other, which a linear model cannot accurately represent. Ideally, the residuals between the linear model and ILHM would be randomly distributed. Non-linearities and interactions among parameters will produce non-random distributions. Figure 5.11 maps water-year average upland ET and DP predicted by the fifth linear model, along with the cell-by-cell root—mean square annual residuals. Parts A and B capture most Of the spatial patterns evident in the upland ET and DP maps (Figure 5.4 is somewhat comparable, but includes wetland ET and DP components as well). The residual maps (Figure 5.11C&D) display a combination Of both random and non-random residuals. The non-random portion of the residuals are associated with both finer soil textures and agricultural land use, suggesting that these variables either interact, or behave non-linearly, or both (likely). There is also a region Of high residual in the southern portion of the model domain similar to the high 0., region present in Figure 5.4C that may be an artifacts of a single outlier gage or climate station. Overall, however, the linear model does a reasonable job Of reproducing the spatial and temporal variability in the ILHM outputs. The mean RMSR of annual upland ET and DP are 7.7 cm and 8.8 cm respectively. 164 Sfimeem...aewm team 8.8 3.: Sea 23. .3. Sam :3 on: a..- a was- 1...; 88 8.: as 2.2 was 8.8 Re- 3a was 3.: as..- a. £8 8.: a; as. a... £8 ems- cm. on: 3.3 8.8. a. 83 was- fie 83 8mm 82- a New. Se ma. 2.2 was- a a: 8mm 5.: . SE 9;. 1.5.3 .88 as. aw as be 5.52.. we... 825. men. 88:. as: .83 853.8 88883. O 1 1 1 1 1 1 1 8\EO Do Do 80 80 3.25 an. a... 8.. 8.. 8.. 8.. 8.. 8.8 83 Se 5.3 $8 sees. 38 82 .50 .86 .18 a... as be. m3... 8.. 821. men. eases. emcee $888.9. m wmdm mos. .md e.g.... .m...V 3.. mod mm... pm... .ed 2% 0. Sam- 8.. .mdm ops. med w..©. was 3.. mo..- mod o..o .md 3.3..- we... no... u. emdm mums. Ema and. «we 3.. mo..- ow. mod hm... Snow- med .w.w v mod .m..- we. mod hm... m2..- mm... we... «1. mm..- em. o..o mm... worm- m..o mm... m 8... nmd 5.5 ...o .ad. N SE 82 .189 .38 .115. a... as be. 52.. 8.. 821. men. ease. am .55. 82.... as: $35.). :o.mm8wom1o.8..:.2 < 688.6an 8m. 883.0 mm: .88. :8 was 5.880 Soc :8 m. bk 53.8388 18382.88 .88 8053.888 808% wageméo: .88 888cm w:.Bo..w 86. mQZ 1m .88 m6 1.“ $1020. 60% .< 8.3 88.. ESSEmoO .258 a... .88 586888. some .0 $88 2.”. .o 8.68.. 2: mm 88.5.8 1.8.35 US$5.88 638 88.88. 23 Ema :38 2: 3 3.8.8 85888. no.8 .0 8.83.8 2.”. @505 AU £83888. 3.58 838. some .8. 836» .0 89:8 2.. 83 Am. 8 .338. £3.58 .88.. 2.. was 3.58.. .22.. 8.... 88153 AQWEEV $8.368 8.8.8 888180.. .88 .mucmsamg 828.880 6826580 .3383... .3308 2: 9.33% A14. Hm. .839. .8. £0.58 828888163588 .88.. .3887. ”w...“ 658... 165 .. m. w w m N. m m a o. . m 9.8... 3.8. .m....1 .m.w1 $.31 mi..- and? made modm cc... n. ma... 8&1 Eda- mm..- mm.w1 mw.m.1 8.21 .mdm- a»... 3.1 mm... Exam 3.3 m. madm- 8.21 3.x- 8.....- coda- and? was on..- 3.5 oadm 3.3 V o..©m1 3... mm... mm... m..om eqw m. ...w mm..- 3.8 amdm E. w flame mmdm 9.....- . $5.. 6.9. .50 gm. gem. $.\ .95 D... .262... WUB mUZnN $60. 2.3:. .253. .80. 282...... 8.5882. 0 1 1 1 1 1 1 1 8 \ So 05 Do 8. 85 2.25 ww. .2. co. co. co. co. oo. .53 we... on... .o..w 3.3 5:8. SE 22 $5 .85 .28 .3. .3. E .82.. 8.. 82. men. $3585. 5252 5.588... m 3.81 3.5.1 Sam. 38.1 $.31 m..- cm... mvd we... .3. med A. m... 8.1 SAN. .5.5- an..- mw.m.1 86.1 m..- h... .o..- mm... 3... mwdm mod av... m. 3.81 8.5.1 3&1 0.3.1 8.81 m...1 a... E... mm... .5... 3.3 mm... .21.... V 3...- om. .561 mm... 3... co.w mm... 3... m mm. we..- .3. .3. Sum ow... .m... m wwd mvd 35.1 m3. hm... . DR 9025 MOB WUZNN mum $5. 6.9. .250 SQ .12.... .3.\ .95 8.2:. am. it.» Ems... .25... 225.). 8285.52.28.83. 4 .w.m 2.5... 8 mm 2252. 885.00 .< .82. 80.. 2822.88 .258 2.. .85 8.58.82. :28 .o 52.2 2.. .o .5.5... 2.. mm 55.2.5.3 5.5882. .55 .o 252...... 53.28 2.. .o 228 5.5.82. 2.. 952... .0 £588.22. .258 .82... .55 5. 25.? .o 82.8 2.. 5.8 .m 8 .532. .2258 82... 2.. .28 2.52 3...... 2.. 883.5 226.8. 88.8 8881.5. .55 £822.38 82.2255 £822.88 £58.58. .258 2.. 9.2%... A... .n.Q .529. 5. 2258 80.823.128.38 82... .5287. ”a...” 2pm... 166 Mean Linear Model Upland DP (cm) <----> 20 27 34 41 48 I Mean Linear Model Upland ET (cm) ---> 32 39 46 53 60 RMSR RMSR Annual ET (cm) Annual DP (cm) <----> <----> 4681012 5791113 Figure 5.11: Maps of linear regression modeled annual average A) Upland ET and B) Upland deep percolation (DP) for 1980 - 2007. Parts C) and D) show mean root mean square residuals (RMSR) between the ILHM annual fluxes and the linear regression modeled values of ET and DP respectively. 167 5.4 Conclusions This study primarily examines the spatio—temporal variability of groundwater recharge (or deep percolation) and watershed evapotranspiration which dominate the water budget of uplands in regions with high-conductivity soils. The region surround- ing the Muskegon River Watershed was simulated using the Integrated Landscape Hydrology Model (ILHM), a novel hydrologic modeling suite. Twenty-eight years of hourly hydrologic fluxes were simulated at 425 m resolution over a 19,000 km2 area. In the MRW, groundwater recharge is more variable than ET both spatially and temporally. Over the period 1980 - 2007, cell-by-cell recharge (or DP) varied ap- proximately 32% annually, with ET varying somewhat less at 14%. In the sandy soils of the MRW, virtually all precipitation in excess of soil field capacity becomes groundwater recharge. In contrast, ET variability is buffered by root-zone storage of moisture to capacity following spring snowmelt, and the subsequent summer rains. Regionally, recharge exhibited a strong decreasing trend away from the lake shore. Two climate teleconnection indices, the Pacific/North American (PNA) and the North Atlantic Oscillation (NAO) were observed to correlate with both ET and recharge. In the MRW, a positive PNA phase indicates cooler winter and spring temperatures, and reduced summer and fall precipitation. This leads to increased an- nual groundwater recharge in a positive PNA phase as the snowpack persists longer in the season, and decreased annual ET due to late-summer moisture deficits. Posi- tive NAO index values are associated with warmer average annual temperatures, and slightly increased annual precipitation. As a result, both annual ET and recharge correlate positively with NA 0. The proximity of Lake Michigan is critical to the hydrologic behavior of the MRW. Because of lake effect precipitation during the fall and winter, groundwater recharge near the lake shore is nearly 50% greater than inland values. If the lake effect is considered, soil texture plays the next most significant role in controlling the spatial 168 variability of recharge and ET. ET and DP are oppositely affected by finer-textured soils; ET increases by up to 30% and recharge decreases by as much as 39% relative to sands. Land use plays an important role as well, particularly because unlike soil texture, land use changes over decadal timescales. Variability between land use classes in sandy soils is almost 20% of upland ET and 45% of recharge. Spatially-variable linear models of annual recharge and ET were calculated from the ILHM outputs. These process-inferred statistical models incorporate seasonal precipitation and temperature, soil field capacity, land cover, and climate index values. They explained approximately 40% and 60% of the spatial and temporal variability of ET and recharge, respectively. Critically, the resolution and domain size of this simulation allowed for strongly co-variate parameters like soil texture and agricultural land cover to be resolved independently. The three most important factors in the linear models were soil field capacity, seasonal precipitation, and coniferous land cover. Recharge and ET exhibited considerable variability at all scales, bounded only by the extent of the model domain and the resolution of the grid cells. The presence of significant spatial and temporal variability even at the smallest model scales high- lights the importance of explicitly modeling the fully—distributed, non-linear processes that drive landscape hydrologic fluxes. Future efforts will explicitly examine the im- pact of model resolution (both through parameterization and numerical solutions) on predicted fluxes to quantify importance of finer resolution simulations for regional and larger scale modeling. Acknowledgments Funding for this work has been provided by the Great Lakes Fisheries Trust Muskegon River Initiative, the National Science Foundation through its Water Cycles program (EAR-0233648), and the Michigan State University Center for Water Sciences. 169 APPENDIX A ILHM Model Development A.1 Evaporation and Transpiration The basis for our potential evaporation and transpiration calculations is the modified Penman-Monteith equation (Montez'th, 1965)presented by Chen et al. (2005a): 85—8 AF + pop 7of ,. . A, (A+7(1+;§s)) Variables appearing in Equation A.1 and others that are not explicitly defined in the PE, PT = (A.1) text are explained in Table A.1. In Equation A.1 , F is the net radiation flux (W 111—2), which is the product of total solar radiation measured at the MAWN gage multiplied by the albedo of each cell. Here we assume that albedo varies seasonally from a leaf-off “brown albedo”, ab, condition to a peak growing season “green albedo”, ag value. When LAI = 0, albedo equals ab and increases linearly to ag until canopy closure is complete, which we assume occurs at LAI = 3. Thus albedo ag LAI 2 3 a. = (A2) £341 (0.9-ab) +ab LAI<3 3) Values of ab and 0.9 are provided in Table A.2. 9 (kg m— is the density of moist air Patm as = , A.3 p (1%?“ + BUT) ( ) given by the ideal gas law: 170 The barometric pressure, Patm (Pa) is calculated according to the barometric formula (Berbemn-Santos et al., 1997): —1W 2 Patm = P0 exp <—R £7)? ) (AA) 9 where z is the elevation relative to mean sea level (m) given by the DEM. 68 (Pa) is calculated from the Goff-Gratch equation (Gofir and Gratch, 1946): log10 e, = —7.90298 (Tst/T — 1) + 5.02808 log10 (Tst/T) . (A.5) — 1.3816 x 10‘7(1011-344(1—T/Tst) — 1) + 8.1328 x 10—3 (10—3-4914ngst/T-U — +log10 est) A (Pa K— 1) is the slope of the saturated vapor pressure- temperature curve, calculated as the numerical derivative of the Goff-Gratch equation; 6, the product of es and measured fractional relative humidity, is the ambient water vapor pressure; Av is the latent heat of vaporization for water (J kg-l) (Harrison, 1963): A, = 103 (2500.5 — 2.359T) (A.6) and '7' is the psychrornetric coefficient (Pa K‘l) (Brunt, 1952) '7 = 1-610ppatm/Av (A.7) The canopy resistance to vapor transport, rci, (m s‘l) is calculated as 1 - = A. TC? gsLAI ( 8) where gs is the stomatal conductance (s m‘l). Values for maximum stomatal con- ductance were taken from Schulze et al. (1994). Unlike Chen et al. (2005a), the aerodynamic resistance, Tm“ (m s‘l) is calculated based on canopy properties and height-adjusted gaged wind speed Allen et al. (1998) h — h h v ‘_ h _1 Ta,- = loge (Jig—9) ~10ge (JET—0) - (Hum) (A.9) m 'U 171 $83 .3 a $22.85 :83 838 as is? 388 Eases mam: eases: as 95388.53. 63: Ease as has... 8%: s a seam“ £8: a a 5:3 v mica Xm.m TmmE Bag “353 a8 EBUEmoo QUEBEC 53832 9Q as m - mica—9.89 mom a m mid - E20580 33 SEE Ewmmaémémmmmvfia EoEEEm 3.6 H 36 - 33250 :mESMéQV 0‘ - m - 53.3 ummbo Ewflmgéwémmmmwfik b39354 ox. - m E voommpfia mo Emma: 8582:ng o8: - m E 35:55 @3938 Ho Emma: “.582:ng 3:: m mmmJoH mam masuwuoaaou .EBQéBBm 2: am no “so m madam vm 83.3383 pfioméamum “Mum o 82 TKme H. be m6 m0 pawn uEqum 3585438200 m imam 5135 n he 5w mo 3.3980 mew H302 9% N. mwwowd mlm 8 A20 mad “—25 nosmuflmuuw EnoSm£>EU om o mmod mx he. Ev mo max: .832 3 o mmmJoH am 62x 5% :38 as 3:985 in moamummmm o& o mvaov Hlvrlmx n Bass. 533 mo S3250 mew am e 3.5% Hivrlwx a he b6 mo 3338 $0 at vasom 33> .333 :osEmoQ 699mm Ema waSSBEo 5 com: mawemfifimm ”#4. £nt 172 Because windspeed was only measured at one height, the effective measurement height is adjusted for canopy height. Here we assume that the height to which windspeed is adjusted, hmw is given by h7nliw = max (hmO, thc) (A.10) where [1C is the canopy height assumed constant for a given land cover (Table A2), and f0 is a factor to move the adjusted wind height some distance above the canopy. The zero displacement height ho is assumed to be 2/3hc (Allen et al., 1998). lm is the roughness length for momentum transport (m) taken as (Allen et al., 1998) (same as the zero displacement), and lv is the roughness length for vapor and heat transport (m) assumed to be (Allen et al., 1998). um is the measured wind speed (m s‘l) adjusted for measurement height according to the wind profile power law assumption ’Uwa = Uw(%n1;:)g) in) (All) (Elliot et al., 1986). where nu, is the raw measured wind speed (m 3’1). In order to calculate evaporation from leaf surfaces EC, surface depressions (Ed), and the soil (E3), Equation A.1 is used but the two conductance terms are modified. 7‘0, is set to 0 for evaporation from leaf surfaces and surface depressions, and Ta,- is calculated using a crop height of 2 cm for surface depression evaporation. To calculate surface soil layer evaporation, TC,- is replaced by rs given by Choudhury and Montez'th (1988): T16 = .12 7's where [6 is the depth from the surface to the top of the evaporative layer of water (m), here assumed to be half the depth of the top soil layer, and (I) is the total porosity. The total transpiration in each cell T (m) is the sum of the root water uptake from each biologically active soil layer, T1 and T2. Following Manfreda et al. (2005), we assume in this version of the ILHM code that the actual transpiration is calculated 173 383 .3 3 £235 68: .3 s :33: 388 $353 as 33: a $8: .3 s 33%? v 88d :53 2 mod and e3 3 - $6628 nosoééoemo 3‘ m 83 33 we was was 3. one - 882 555 § m 83 83. 2.0 3d 3 mg no - 832 520 as m 85 3o 23 35 $2 $3 $3 - Sam 80m n m 8d 8d 2 93 ca 3 wd E Ewsm 3250 3 H od 3. 3 9m 3 9: 3: TE .0. 853300 338% mm 850m Sam H335 $5335 awoken Edam w< 9213 EEO :03“:on Basra 23H. $5-984 mhmumEmBQ $3-923 ”max 2an 174 from the potential value by linearly interpolating between 0 at the permanent wilting point and the potential rate at 75% of saturation according to: 1i Zr 12' where S- is the soil moisture m of the it" soil layer, (IL is the permanent wiltin ‘ z 33 8 4 T2“ 2 (Si > (I)_33) - min [1, 35,“, ((1)210—1] -PT- (A.13) point of the soil (Table A3) and Q is the thickness of the ith soil layer (m). The term S,- > (I)_33 is a logical statement that returns a value of “1” if true and “0” if false. Porosity values, (Dz- are taken as a function of soil type as given by Table A.3. Biologically active soil thickness is calculated as the depth above which 90% of the root mass lies using the asymptotic equation (Gale and Grigal, 1987): _ d y — 1 — B (A.14) where y is the cumulative root fraction at depth d. = 0.9 (cm) for this study; [J’ is a land cover-dependent parameter (Table A.1). We use Equation A.14 to solve for d with a fixed cumulative root fraction 3; = 0.9, 1 = Z 11- =1ogfl(0.9). (A.15) i Total evaporation E (m) is the sum of canopy evaporation, EC, soil evaporation E3, and surface depression evaporation Ed. Soil evaporation is calculated according to Chen et al. (2005a) as E3 = min (PES, d3) - (1 — fc) (A.16) where d3 is the soil-controlled exfiltration depth (m) calculated by d3 = 36 -At"1/2 (A.17) and At is the model timestep length (3), while Se is the soil desorptivity (m 8-1/2), calculated as in Entekhabi and Eagleson (1989): = 8®1ksat¢b ”2 S(m/2+2) Se 3(1+3m)(1+4m) 0 (M8) 175 where ksat is the saturated hydraulic conductivity of the first soil layer (m s "1), m is the pore size distribution index assumed to be a function of soil texture (Table A.3), and SD = 61/ ((1)1l1) is the fractional saturation of the first soil layer. As in Manfreda et al. (2005), the closed canopy fraction ( fc) is defined by the empirical relationship (Eagles-0n, 1982) fc = 1 — e‘fl'L’“ (A.19) where ,u is a constant for a given land cover type given by Table A1 Calculating EC requires a full canopy water balance model. The canopy water balance is calculated using: where ASC is canopy water storage (In). Incoming rainfall is first subjected to in— terception up to the water holding capacity of the canopy (m) given by (Dickinson et al., 1991): Smax = 1 x 10—4LAI (A21) where LAI is the leaf area index (mZ/m2). The available interception capacity of the canopy is then given by 55 = scmax — 351. (A22) Additionally, we modify the model of Chen et al. (2005a) to allow some water to penetrate the canopy at all times based on the assumption that the canopy is not completely closed. Interception at time t is then Int = min (65, P - fc). (A23) Canopy evaporation, EC is then calculated as in Manfreda et al. (2005) with EC 2 min [(Sc/Scmax)2/3 - PEG, Sc] . (A24) Surface depression evaporation, Ed occurs only when water is stored in surface depressions. At each time step, any water stored in surface depressions Sd from the 176 previous timestep is added to throughfall from the canopy, or snowmelt from the UEB model, such that precipitation excess runoff, Re is given by R8 = min (0,1D + 35-1 — Inf — Int) (A25) where infiltration, Inf, is calculated as discussed below. Sd is then calculated as Sd = min (Sdmax, Re + Exl) (A26) where Exl is the exfiltration out of the first soil layer and the depression storage capacity. Sdmax is assumed to be constant for a given combination of slope, land cover, and soil type (see table in Manfreda et al. (2005)). Depression evaporation, Ed is then given by Ed 2 min [8d, (1 — FC) - PEd] . (A27) A2 Infiltration, Percolation, Throughfiow, and Exfiltration The next three terms of the water balance (Equation 3.1), percolation, Pc; throughf low, Tr; and exfiltration Ex are calculated within the soil water balance model. First, the outputs of the canopy model, snowmelt model, and depres- sion storage model are used to calculate infiltration into the surface soil layer Inf = min (P — Int + 53*, ism“) . (A28) The infiltration capacity, is max is a function of the moisture content of the surface soil layer and is calculated according to Zsmax = max ((1)111 ‘ ASE—1:33.901) (A29) where ll is the thickness of the first soil layer as defined previously, 31 is the moisture stored within the first soil layer, and isat is the saturated infiltration capacity, which can vary with time due to the influence of impermeable frozen soils. 177 H 32 Sad 33 was $3 83 - aoHsEEmHQ mam 28H E H $2 NH; 82 $3 $3 HE; E 252m wsznsm as m wsH 8H H3 9: mg mg - Meagan 83850 as, z m. ah 9m ON m.» v.3 93 HIE 588553 53380 :9, 6 m 2.9 HH.o 2d 23 So So - Eon wag? 8o 2'9 a mg mg was 2. 26 Ed - 3688 25 male m aHH a.” 3 a: EH was. Tm Ea: x $688 nosesafi :5 H 83. Rod 23. H2; £3 83 - 3380 $33 Hsfiem .6 H wand 82. Home 32 E3 53 - .Hzmeom H309 e H mHH Em aH m3 EH mam. Him :6on £38358 523% was moi—om €32 EH34 Smog ham Ewoq 355m vcmm >595 95m 3:5 aoEEwmQ Bantam $20 83me :om moshmmoa mom ”mas @558 178 Infiltration is applied to the first soil layer, which can then percolate into the second layer. First, the soil moisture storage at the end of each timestep in the first layer is calculated as 31+ = 55* + Inf — Es — T1 + Tn — P1 + 153:2 — E331 (A.30) where P1 is the percolation of water from the first soil layer to the second. Note that calculating T1 requires knowing S 1. To avoid having to solve the coupled equations, T1 is calculated using an intermediate value of S?“ = Si-1+is—Es. Then SI = Si‘ —T1, and Si+ = S 1 + T71 — P1 + E1132 — Eccl. Given S], percolation into the second layer, P1 is given by P1 = max (31 - SlmaXa PDREAM) (A31) where 51mg,X = (1)111 and PD R E A M is the percolation calculated according to Man- freda et al. (2005) given by: 0 St“) S —33lz' PDREAM = 1_ 1/(1—7) . . At-k ('y—l) S- 7 . . 52 max (SI “ I: imax + (Sir—31;) :l ) St“) > (I)-33lz (A32) where 7 = (2 + 3m)/m. PD RE AM effectively allows percolation only when soil moisture exceeds the field capacity given by (bili. The second-layer soil moisture at the end of the timestep, (52), is calculated simi- larly to EquationA.30: 55+ = 53* + P1 — P2 — T2 + Trg — E32 (A.33) where T2 is calculated from the values of 83—1 at the previous timestep. PC = P2 = PDREAM (A34) . . . ,t_ , _ is calculated from an intermediate value of 51 = i. f 1 + P1 — T2- 179 Throughflow out of a cell is calculated as 2 AZ Alb 63 down _ 03' . : l- , A . . __ —_ . ’— A. Tron” 7' :L‘ (l/kfli + 1/k0i,d0wn) [Ax + (A6 )z' < A]; )i ( 35) where the effective unsaturated hydraulic conductivity for subsurface flow in layer 2' is taken as the harmonic average of the unsaturated hydraulic conductivity in the cell, km- (m 3.1) and the down slope value kflwown; A2: is the model cell resolution (m); Az / A1: is the vertical gradient in the down-slope direction; (Aw/A6),- is the slope in meters of the moisture retention curve in layer 2', and _ Si(t)/Si max " 6r 9 _ q) _ 0T (A.36) where GT is the residual volumetric moisture content assumed to be a soil-texture dependent property (see Table A.3). The assumption that flow only occurs parallel to the dip of the slope requires that Trout 2 0. Tr,- is then calculated as TT‘i = T7'out,'up - T7'out,z'- (A37) Finally, Ecrgis calculated as the soil water in excess of saturation given by E172 = max (0, 5'2 + Trg — 82 max). (A.38) This is then applied to the first layer prior to calculating E121 = max (0, 31 + TT1+ EIL‘Q — 31 max)- (A39) A.3 Runoff Routing Water exfiltrated from layer one is then applied to the surface depression model, thus R (m) is simply R. = Re — S’d (A.40) which is then routed to the streams using an approach modified from that presented by Manfreda et al. (2005). Once generated, runoff cannot infiltrate and is instead 180 routed overland and through streams according to the D8 flowdirection algorithm in ARC (ESRI, 2003) with runoff times given by the velocities in each cell along the flowpath. Runoff is assumed to travel overland at a velocity given by the Kerby time of concentration equation (Kerby, 1959) I. ll _ TL 0.467 Leell = 86.735 (if) (A41) where tee); is the time required to completely traverse a model cell (5), lcell is the length of the model cell (m), n is the dimensionless Manning’s Roughness coefficient (values from McCuen (2004)) and s is the fractional slope of the cell in the downslope direction. The velocity (111 3’1) is then 1 ll "land 2 teen. (A42) C8 Once the runoff enters the stream channel its velocity is calculated using Manning’s Equation (McCuen, 2004) I vstream = gT2/331/2 (A43) where 7* is the hydraulic radius (m) given by the ratio of the stream cross-sectional area to the wetted perimeter. 181 APPENDIX B ILHM Regional Upscaling This appendix provides additional equations from those in Appendix A and Chapter 4. In some instances, variable names and symbols have been changed from Appendix A to be more consistent. B.1 Surface Domain B. 1.1 Shortwave Radiation Incoming total shortwave solar radiation is provided to the model as two constituents: beam and diffuse radiation, Q33 and QSD [W/mz]. For this study, the the beam and diffuse fractions were obtained from real ground solar radiation modeling (see Solar Radiation Preparation, below), though one could use the method of Chen et al. (2005a). The canOpy then intercepts a portion of each both constituents, transmitting to the ground beam radiation QGB = Q33 X TCBand diffuse radiation QGD = Q S D x TCthere TC B is the canopy transmissivity for beam radiation, and TC D for diffuse. These transmissivities are defined as Fa = q Bu]. TCB COS (63) and ( ) PG 4 = B2 TCD COS (913V ( ) 182 where 63 is the zenith angle (degrees from vertical) of beam radiation, and 6D is the effective diffuse radiation angle, defined as OD = 0.537 + 0.025L(Chen et al., 2005a), and ngs the gap fraction of the canopy, defined as FG = exp [—0.5§2 (LAI + SAI)] (8.3) where Q is the canopy clumping index, a measure of the distribution of leaves within the canopy. Q values near 0.5 indicate that leaves are distributed in rows, and values closer to 1 that the leaves are randomly distributed within the canopy. LAI is the leaf area index [m2/m2] of the canopy representing the leaf area per unit surface area, and SA] is the stem area index. The ground then reflects back to the canopy a portion of the transmitted shortwave radiation, given by QGC = (FW + F5) X (OWBQGB + GWDQGD) + FU (aUBQGB + aUBQGD) (34) where FW S is the fraction of each cell covered by streams and wetlands, and FU the fraction by uplands; these terms are related by 1 = FW + F S + FU. The a terms refer to shortwave albedos (weighted average of visible and near IR) for wetlands and uplands to beam and diffuse radiation. It is assumed that the canopy then reflects back to the ground a negligeable portion of this radiation, thus the infinite series of reflected radiation terms is truncated after the first reflection (ground - canopy). The net radiation absorbed by the canopy in each cell is then given by QCN = (1 - GOD) >< (Q31) - QGD + QGC) + (1 - 003) X (QSB - QGB)- (B-5) Here, the ground-reflected shortwave radiation QGC is assumed to be diffuse radia- tion. The ground absorbs QGN = (QGB + QGD) - QGC- (B-G) 183 From this equation, the net shortwave radiation in the upland portion of each cell is QU N = QGN FU, while for the wetland and stream portions QW N = QGN (FW + F5)- For this study, the beam and diffuse radiation terms for the canopy are assumed to be equal, and are given as a linear mixing of two end members: the brown and green albedo of a given vegetative canopy type as in Walko and Tremback (2005). %l (0031' — aC,bI‘) + (1035,. LAI < 3 with aC,gr and ac br are the albedos of the canopy when fully green and fully brown, respectively. Frozen water stored on the canopy can have a siginficant effect on the albedo, which is described by the following two equations: 003 = aCB (1 - F05) + GSNFCS, and (8.8) 000 = 000 (1 - Foss/2) + O‘SN (PCS/2)- In Equation 8.8, 0 SN is the albedo of new snow, and FCSis the fraction of canopy covered by snow, calculated as F05 = SC / SC M where SC is the total water stored on the canopy, which is assumed to be frozen when the air temperature is below freezing, and SC M is the maximum canopy water storage. The factor of 2 in the denominator of diffuse albedo calculation arises from the fact that snow only covers the top half of each leaf, and only a small portion of the total stem area. Since beam radiation primarily sees the tops of leaves and stems, no modification is necessary. Net upland albedo is calculated according to the net soil albedo and the beam and diffuse components of snow albedo. Bare ground upland albedo is calculated via the linear combination of dry and wet soil albedos, asaSO = O‘SO,dry (1 — F 3 AT,1) + o' SQwet F S AT,1- (150 is the net soil albedo, aSO,dry is the dry soil albedo, 080,2uet is the wet soil albedo, and F S AT, 1 is the soil saturation fraction in the first soil layer. The soil saturated fraction in layer 2' is defined as F 5 A71,- = S 3’, / S Swami where S 3,,- is the soil water storage in a given layer [m], and SS.max,z' is the maximum water 184 storage. SS,max,z' = (195,, 03,2" where (119,,- is the total soil porosity in a given layer and D5,, is the thickness of that soil layer [m]. Diffuse and beam snow albedo on the upland portion of a cell, oz 3 N D and a S N B: are a function of the age of the snow at the air interface. They are calculated as in the UEB snow model Tarboton and Luce (1996), which itself is derived partially from the land surface model BATS Yang et al. (1997). This formulation includes a modification of the albedo to account for shallow snow and partial transmission of the visible and near IR components of shortwave radiation. The net upland albedo is calculated as 0.903 DSN = 0 “U8 = OSNB DSN > 0 (B 9) ar301) DSN = 0 (IUD = aSND DSN > 0 Here, D5 N is the depth of snow in the upland portion of each cell. In this version of the model, snow cover in wetland and stream portions of cells is ignored, thus wetland albedo depends only on the albedo of water and any ice cover. The diffuse albedo of open water is (1W0 D, and the beam albedo is 0W0 B = 0W0 D / (cos(63)1'7 + 0.15)(Bonan, 1996). Ice albedo is assumed constant, aw I D = aw I B- The net wetland albedo is then aW013 DICE = 0 O'WB = aWIB DICE > 0 (B 10) OWCD DICE = 0 awn = aWID DICE > 0 185 cam: .3 3 seems 388 :3 a :25: 62: .228me 38m $3525 as 33: Va 62: re e :33? m. a. m m m m N a 85% em. 3o :2 was as :2 :2 was 205.68 03 8d 85 Be So as a; 83 mean 3 a; and Ed 2d 2: 84 23 3233 OS Se 86 3o 2:. as :3 83 as; 3. a; one «.3 86 0.2 9.2 83 2.820% E :5 03 mg 8d we 3 Sad wages e: mg one 23 mg as 3 Sad 233%? S: a; Q3 and OS 3 3 E3 5%: xexéagi Ia Ii TE? Tee... sees E59 In ..E...3 22.: E whoeoawnma Suwanee om: ESQ ”mm 2nt 186 Table 8.2: Depression capacity, S D E P,max [m x 103]. Values are from Manfreda et al. (2005), assuming no soil-texture dependence, and taking values for clay soil textures. Slope [‘70] Land Use 0 - 0.5 0.5 - 5.0 5.0 - 10.0 > 10.0 urban 1.80 1.37 0.94 0.51 a riculture 1.50 1.08 0.66 0.23 s rubland 2.00 1.46 0.93 0.39 deciduous 2.50 1.88 1.25 0.63 water 0 0 0 0 wetland 0 0 0 0 barren 1.00 0.71 0.42 0.12 coniferous 2.50 1.88 1.25 0.63 B.1.2 Snowpack Accumulation and Melt In Appendix A, the snow depth was calculated at only a single point, and averaged across the watershed. Here, the code was modified to be fully-explicit, and snow depths are calculated for each cell in the model. B.1.3 Infiltration and Root-zone Moisture Vertical Redistri- bution Water falling through the canopy may infiltrate, according to WINE = min (VVTHRFIMa WINEmaX) a (311) where WIN F‘max is the maximum infiltration capacity, defined as KSAT T1 2 TSF , WINFm’iax = 7 (B.12) O-I'KSAT T1mfinw was: 85.38% .mumuoafiaq uaowaoaow mnzaxgfiom ”0m 2an 188 Table 8.4: Table listing physical constants Symbol Variable Description Units Value pair density of dry air @ 10 °C kg m-3 1.247 Pice density of ice kg III—3 920 p.,,.,.,.,, density of mineral soil matrix kg m-3 2650 phgo density of water kg m‘3 1000 Hair gas constant dry air J kg‘lK—1 287.05 Rhgo gas constant water vapor J kg_1K_1 461.5 g gravitational acceleration @ msl and 44.5°N m 8‘ 9.80665 Cpm-T heat capacity of moist air J kg-1K_1 1013 Cm“a heat capacity of ice J kg‘lK“1 2114 Emmi heat capacity of mineral soil matrix J kg‘1K_1 754.7 'ng heat capacity of organic soil matrix J kg’lK”1 1923 Cpfigo heat capacity of water J kg‘lK'“1 4181.3 um, kinematic viscOsity of air m23-1x105 1.5 H f latent heat of fusion J kg"1x10’5 3.34 H s latent heat of sublimation J kg—1 x10“5 28.34 mm, molecular weight dry air kg mol’1 x102 2.8964 whgo molecular weight water kg mol"l x 102 1.8015 0 Stefan Boltzmann constant J s"lm_2K‘1x108 5.6704 K (Lair thermal conductance air W m-1K_1 0.025 K q,ice thermal conductance ice W m’lK’1 2.2 K (I’mml thermal conductance mineral soil matrix W m"1K"1 2.9 qug thermal conductance organic soil matrix W m‘1K_1 0.25 K (#120 thermal conductance water W m-1K_1 0.57 kv K von Karman’s constant - 0.41 Table B.5: Table listing albedo and emissivity parameters, and sources Symbol Variable Description Units Value Source (ISA/ms albedo snow vis - 0.95 1 a S Nmir albedo snow nir - 0.72 1 0W0 D albedo water diff - 0.06 2 0W0 B albedo water beam - 0.06 2 QWI’dTy albedo ice dry - 0.50 3 avg/[Met albedo ice wet - 0.25 4 e S N emissivity snow - 0.98 5,6 (30 emissivity soil - 0.96 5 6 1,20 emissivity water - 0.96 4 Eice emissivity ice - 0.97 4 DC, 5 N albedo extinction depth, snow m 0.10 6 1(Aoki et al., 2003) 2(Hostetler and Bartlein, 1990) 3( Yongji'u et al., 2001) 4(Patterson and Hamblin, 1988) 5(Bonan, 1996) 6(Tarb0ton and Lace, 1996) 189 Table 8.6: Table of domain-indepedent physical parameters and sources Symbol Variable Description Units Value Source h R,wet roughness height, wetland m 0.10 - hR,so roughness height, soil In 0.01 - T501: soil freezing temperature 0C -0.25 - TC D canopy dormant temperature °C 1.0 1 500 specific canopy capacity mm m‘2 (leaf) 0.15 2 77W wetland light attenuation coefficient m‘1 0.70 3 P0 wetland neutral Prandtl number - 1.0 4 a wind-offset power - 0.14 5 lg, wind-offset factor - 2.0 - 1(Chen et al., 2005a) 2Average of Dickinson et al. (1991) and Yongjiu et al. (2001) 3Modified from Bonan (1996) 4(Hostetler and Bartlein, 1990) 5(Allen et al., 1998) Table B.7: Table of domain-dependent physical parameters, values are explained in Chapter 4 Symbol Variable Description Units Value F150 soil frozen impermeable factor - 0.90 [U R B specific urban impermeability 0.20 CW ET D conductivity, disconnected wetlands m m”1(head) hr—1x106 3.60 C ETC conductivity, connected wetlands m m'1(head) hr"1><103 3.125 WT water table amplitude m 1.25 FRW ET recession constant, wetlands m m-1(head) hr—1x103 7.0 190 Movement of water vertically in the shallow root zone can be calculated according to one of two methods: a 1—D Richard’s equation (Hillel, 1980), or via a free-drainage model similar to that in Appendix A. For this study, the free-drainage model was chosen because of its stability and computational efficiency. For this model, the movement of water from one soil layer to the next is calculated as: WS’i = min (max (55,1; — SSFC,iv 0) ,ng) . (B.13) Here, S 59,;is the water stored in layer i [m], S S FCJ is the field capacity of the soil layer [m], and K9,,- is the unsaturated hydraulic conductivity in that layer. K9,,- is calculated via a modified form of the Maulem-vanGenuchten equations(Maulem, 1976; van Genuchten, 1980; Schaap and Leij, 2000) M 2 K9 = K06’ [1 — (1 — (WM) [ . (3.14) In Equation B.14, K0 is a soil-textural dependent value that is less than K S AT: which provides a more accurate description of unsaturated conductivity according to Schaap and Leij (2000).6’ = (6 - OR)/(63 — 03) where 0 is the volumetric soil moisture [m3/m3], HR is the residual volumetric soil moisture, and 03 = (133 is the total porosity in a given soil layer. Also, in Equation B.14, M = 1/N, where N is the van Genuchten empirical parameter. Flux out of the bottom soil layer into the deep unsaturated zone is terms deep percolation, WDpU. B. 1.4 Canopy Transpiration Upland root-zone transpiration has been slightly modified from Appendix A to ac- count for the relative fraction of root mass within each soil layer. The new formulation is min(1,4/3 - FSAT) ° WTp ° FSR,i SS,i > SSW,i WTU,I = , (315) 0 35¢ S SSW,2' 191 where WTP is potential canopy transpiration [m/s], calculated as in Appendix A, FSR,i is the fraction of root mass in the root zone within layer 2', and SSW is the wilting capacity of the layer [m]. The fraction of the root mass within each layer is calculated by assuming that root mass declines exponentially with depth, according to the following equation, modified from Gale and Grigal (1987) FSR’z’ = [3(DS,i—1+DS,i)/2/100 _ [8(DS,i+DS,i+1)/2/100, (BIC) dealing appropriately with the first and last soil layers. It can be necessary to nor- malize F313,,- by the sum 221:5; F513,,- for cases when soils are thin relative to the root depth. The parameter [3 is a land-cover dependent parameter, and D 5,, referes to the cell-center depth of layer i [m]. In wetlands, transpiration occurs at a rate equal to the potential, WpT, unless lim- ited by available moisture. For wetlands connected to the water table, transpiration is removed directly from the phreatic zone at the potential rate as WTW, p = WT p [m / s]. Transpiration is removed surface storage in disconnected wetlands up to the limit of available moisture, WTP 3W5 > SWS,min 0 5W5 S Sws,min B. 1.5 Inundated Area Extent Inundated areas are defined as wetland landunits containing ponded water at the surface. This can occur in two situations: 1) wetland sediments have low conductance, and are disconnected from the water table, and 2) the water table in connected wetlands is above the bed elevation. The first situation is handled simply by tracking the storage in disconnected wetlands. The second situation requires some capability to describe the transient depth of the water table. Provided a fully-coupled surface 192 and groundwater system, this becomes straightforward. As this capability is not yet available in ILHM, another approach is required. Generally speaking, the water table depth varies annually roughly sinusoidally, peaking in the spring and early summer in response to snowmelt and spring infil- tration, and reaching an annual minimum sometime near October. Locations with thicker unsaturated zones exhibit smaller amplitudes of variability and phase shift relative to shallow unsaturated zone areas. For shallow unsaturated zone areas, the water table is assumed to vary sinuoisdally as: 27r (DD — 0190)) (8.18) D = A ° WT WT Sln ( 366 where DWT[m[ is the depth to the water table, AWT [m] is the amplitude of water table fluctuation in areas of shallow water table, DD [day] is the decimal day of the year, and DDO [day] is an offset to control the phase shift of the water table. For this study, DDO = 0, AWT = —1.25 which results in a peak water table elevation (a minimum depth) in April, and a maximum depth in October. To then calculate which wetlands are inundated as a function of water table depth, each wetland is assigned an average depth to water (Table 4.3). The depth to water in a particular wetland is then was = DWA,z' — DWT — Dwss, (319) with depth to water in wetland cell 2' equal to DWW,z‘ [m], DW A [m] is the average depth to water in the wetland cell i, and DW S is the thickness of water stored in that wetland cell, calculated as: WWS,z' WWS,i Z 0 DWS’, = . (8.20) WWS,i(I)1,i WWS,i < 0 Note that WW S, the water stored in wetlands, may be negative, indicating water stored in the soil within a non-inundated wetland cell. In Equation B20, (FL,- is the 193 soil porosity in the wetland landunit in the first layer of cell i. If the depth to water DWW < 0, the wetland cell is inundated. B.1.6 Wetland Evaporation, Lake Temperature, Ice Pack Wetland evaporation is calculated directly from the lake temperature, air temperature and aerodynamic resistance of each wetland cell. First, the sensible heat flux is calculated as pm aiGC.air QH = ’ ' B21) Tai (Tair — TW) ( where Q, [Vi/m2] is the sensible heat flux, pmm'r [kg/m3] is the moist air density, Cnair is the heat capacity of air, rm- [5 / m] is the aerodynamic resistance (calculated via Equation A9), Tair [OC] is the input air temperature, and TW [OC] is the wetland surface layer temperature. The moist air density is calculated via Equation A.3. After calculating sensible heat flux, latent heat flux Q L [W/m2] is calculated by QL = AevPEVV (13.22) with Aev the latent heat of evaporation, and PEW [m/s] the potential evaporation rate of wetlands. The latent heat of evaporation is As es where A6 is the gradient of the vapor-pressure temperature curve, es is the saturated vapor pressure at the water temperature TW, and Rth is the gas constant of water vapon Lake temperature is calculated by solving the energy balance Equation 4.6 using the method of Hostetler and Bartlein (1990). Ice pack thickness is calculated in this process module as discussed in Chapter 4. 194 Table 38: Lookup table of unsaturated zone wetting front velocities. Ksat [m/dl VWF [d/ml 1.7 4.18 2.5 3.13 5.0 2.47 10 1.81 20 1.62 B.2 Deep Unsaturated Zone Unsaturated zone routing has been modified from Appendix A, adding variable wetting-front velocity. Wetting fronts are assumed to move with constant velocity through the deep unsaturated zone. All moisture leaving the root-zone is accumu- lated for the length of the groundwater domain stress period, currently 7 days. This pulse of water then moves through the unsaturated zone at a fixed rate, depending on the calibrated saturated conductivity of the deep sediments. No data from monitoring wells in the MRW were available for this study, so wells instrumented in the nearby Grand Traverse Bay Watershed were used to calibrate a 1-D Richard’s equation model. For this model, a saturated conductivity of 5 m d ’1 was assumed. van Genucthen equation parameters a and N were calibrated such that the observed mean wetting front velocity VW F [cl/m] (determined by the mean arrival time of a pulse of water traveling through a simulated 10 m homegeneous column) matched that observed in wells across the GTBW, roughly 2.5 d m-l. These values are a = 1.76 and N = 3.2. These values were then help constant as the saturated hydraulic conducitivity, K sat, was varied to build a lookup table spanning the range of values used in the MRW saturated groundwater domain module. Table 38 lists the values obtained through this procedure. 195 Table 8.9: Table of prelimmarily—calibrated saturated hydraulic conductivities for deep sediments Sediment Type Conductivity [m/d] Water 100 Peat and Muck 5 Dune Sand 50 Lacustrine Fines 5 Lacustrine Sand and Gravel 20 Outwash and Alluvium 15 Ice-Contact Outwash 20 Morainal Tills 8 Layers 2 and 3 2.5 B.3 Saturated Zone Discharge from each “drain” type cell in the MODF LOW model is recorded at ev- ery timestep. After completing the MODFLOW run, these discharge values are then summed by sub-basin (there are approximately 40 sub-basins in the MRW). Evapora- tion demand not satisfied by surface storage is then passed to this step. This demand is summed at the same sub-basin level as discharge, and subtracted from discharge as as detailed in Equation 4.10. 196 BIBLIOGRAPHY Ahuja, L., J. Hanson, K. Rojas, and M. Shaffer (Eds) (1999), The Root Zone Water Quality Model, Water Resources Publications, CO USA. Akbari, H., L. Rose, and H. Taha (2003), Analyzing the land cover of an urban environment using high-resolution orthophotos, Landscape and Urban Planning, 63,1—14. Allen, R., L. Pereira, D. Raes, and M. Smith (1998), Crop evapotranspiration - guidelines for computing crop water requirements - FAO irrigation and drainage paper 56, Tech. rep., Food and Agriculture Organization (FAO) of the United Nations, Rome. Ambaum, M., B. Hoskins, and D. Stephenson (2001), Arctic Oscillation or North Atlantic Oscillation?, Journal of Climate, 14 (16), 3495—3507. Anderson, J ., E. Hardy, J. Roach, and R. Witmer (1976), A land use and land cover classification system for user with remote sensor data, Tech. rep., US Geological Survey Professional Paper 964. Andresen, J ., and J. Winkler (2009), Michigan Geography and Geology, chap. Chapter 19: Weather and Climate, Pearson Custom Publishing, Boston, MA. Andresen, J. S. (2004), Level-III NEXRAD data, Archived from the National Weather Service. Andresen, J. S. (2007), Automated Weather Observing System (AW OS) climate data, Archived from the National Weather Service. Andresen, J. S. (2008), Level-iii nexrad data, Archived from the National Weather Service. Aoki, T., A. Hachikubo, and M. Hori (2003), Effects of snow physical parameters on shortwave broadband albedos, Journal of Geophysical Research, 108(D19), 4616. Arnold, J ., and P. Allen (1996), Estimating hydrologic budgets for three illinois watersheds, Journal of Hydrology, 176(1-4), 57—77. Arnold, J., P. Allen, and G. Bernhardt (1993), A comprehensive surface- groundwater flow model, Journal of Hydrology, 142, 47—69. 197 Assel, R. (1992), Great Lakes winter-weather 700-hPa PNA teleconnections, Monthly Weather Review, 120( 9), 2156w2163. Assel, R., K. Cronk, and D. Norton (2003), Recent trends in Laurentian Great Lakes ice cover, Climatic Change, 57, 185—204. Austin, J ., and S. Colman (2007), Lake Superior summer water temperatures are increasing more rapidly than regional air temperatures: A positive ice- albedo feedback, Geophysical Research Letters, 34, L06,604. Avnir, D., O. Biham, D. Lidar, and O. Malacai (1998), Is the geometry of nature fractal?, Science, 279, 39—40. Bach, M., and T. Meigen (1999), Do’s and dont’s in Fourier analysis of steady- state potentials, Documenta Ophthalmologica, 99, 69—82. Baret, F., and G. Guyot (1991), Potentials and limits of vegetation indices for LAI and APAR assessment, Remote Sensing of Environment, 35, 161—173. Berberan-Santos, M., E. Bodunov, and L. Pogliani (1997), On the barometric formula, American Journal of Physics, 65(5), 404—412. Beven, K., and M. Kirkby (1979), A physically-based variable contributing area model of basin hydrology, Hydrology Science Bulletin, 24 (1), 43—69. Blackman, R., and J. Tukey (1958), The measurement of power spectra, from the point of view of communications engineering, New York: Dover Publications. Bogena, H., R. Kunkel, T. Schobel, H. Schrey, and F. Wendland (2005), Dis- tributed modeling of groundwater recharge at the macroscale, Ecological Mod- elling, 187, 15—26. Bonan, G. (1996), A Land Surface Model (LSM version 1.0) for ecological, hydrological, and atmospheric studies: Technical description and user’s guide, Tech. Rep. NCAR/ TN—41 7+5 TR, NCAR Technical Note. Bouraoui, F., G. Vachaud, L. Li, H. Le Trent, and T. Chen (1999), Evaluation ' of the impact of climate changes on water storage and groundwater recharge at the watershed scale, Climate Dynamics, 15(2), 153—161. Boutt, D., D. W. Hyndman, B. Pijanowski, and D. Long (2001), Identifying potential land use-derived solute sources to stream baseflow using ground water models and GIS, Ground Water, 39, 24—34. Bras, R., and I. Rodriguez-Iturbe (1985), Random Functions and Hydrology, New York: Dover Publications. 198 Brunt, D. (1952), Physical and Dynamical Meteorology, Cambridge University Press: Cambridge. Burnett, A., M. Kirby, H. Mullins, and W. Patterson (2003), Increasing Great Lake-effect snowfall during the twentieth century: A regional response to global warming?, Journal of Climate, 16', 3535—3542. Cengel, Y., and M. Boles (2001), Thermodynamics: An Engineering Approach, McGraW-Hill. Center, N. C. D. (2001), NESDIS, NOAA, US. Department of Commerce, U.S. daily surface data (DS—3200 and DS-3210), Tech. rep., National Climatic Data Center, Asheville, NC. Chen, J ., X. Chen, W. Ju, and X. Geng (2005a), Distributed hydrological model for mapping evapotranspiration using remote sensing inputs, Journal of Hydrol- ogy, 305, 15—39. Chen, J ., C. Menges, W. Ju, and X. Geng (2005b), Global mapping of foliage clumping index using multi-angular satellite data, Remote Sensing of Environ- ment, 97, 447—457. Chiew, F., T. McMahon, and I. O’Neill (1992), Estimating groundwater recharge using an integrated surface and groundwater modelling approach, Journal of Hydrology, 131, 151—186. Choudhury, B., and J. Monteith (1988), A four-layer model for the heat budget of homogeneous land surfaces, Quarterly Journal of the Royal Meteorological Society, 114, 373—398. Christensen, J., et al. (2007), Climate Change 2007: The Physical Sciences Basis, Contribution of Working Group I to the Fourth Assessment Report of the Intergovernmental Panel on Climate Change, chap. Regional Climate Pro- jections, Cambridge University Press: Cambridge, United Kingdom and New York, NY, USA. CODATA (2002), CODATA internationally recommended values for the funda- mental physical constants. Cooley, J ., and J. Tukey (1965), An algorithm for the machine calculation of complex fourier series, Mathetmatics of Computation, 19, 297—301. Coulibaly, P., and D. Burn (2004), Wavelet analysis of variability in annual Canadian streamflows, Water Resources Research, 40(3), W03,105. 199 Crawford, N ., and R. Linsley (1966), Digital simulation in hydrology: Stanford Watershed Model IV, Stanford University Department of Civil Engineering, Tech Report 39, Palo Alto, CA. Dahlstedt, K., and H. Jensen (2004), Fluctuation spectrum and size scaling of river flow and level, Physica A, 348, 596—610. Dai, A., K. T.R., B. Sun, and K. Trenberth (2005), Recent trends in cloudi- ness over the United States: A tale of monitoring inadequacies, Bulletin of the American Meteorological Society, 87, 597—606. Dickinson, R., A. Henderson-Sellers, C. Rosenzweig, and P. Sellers ( 1991), Evap- otranspiration models with canopy resistance for use in climate models, a re- view, Agriculture and Forest Meteorology, 54, 373—388. Dickinson, R., A. Henderson-Sellers, and P. Kennedy (1993), Biosphere- Atmosphere Transfer Scheme (BATS) version 1e as coupled to the NCAR Community Climate Model, Tech. rep., NCAR Technical Note NCAR/TN- 387+STR, Boulder, CO. Downer, C., and F. Ogden (2003), GSSHA user’s manual: Gridded surface subsurface hydrologic analysis, version 1.43 for WMS 6.1, Tech. rep., US. Army Corps of Engineering, Engineering Reseaerch and Development Center. Eagleson, P. (1982), Climate, soil and vegetation, 5: A derived distribution of storm surface runoff, Water Resources Research, 18(2), 325—340. Eidenshink, J. (1992), The 1990 conterminous U.S. AVHRR data set, Pho- togrammetric Engineering and Remote Sensing, 58, 809—813. Elliot, D., C. Holladay, W. Barchet, F. HP, and W. Sandusky (1986), Wind energy resources atlas of the United States, Tech. rep., Pacific Northwest Lab- oratory: Richland, WA. Entekhabi, D., and P. Eagleson (1989), Land surface hydrology parameteri- zation for the atmospheric general circulation models including subgrid scale spatial variability, Journal of Climate, 2, 816—831. ESRI (2003), ArcInfo workstation, Tech. rep., Environmental Research Systems Institute: Redlands, CA. Farrand, W., and D. Bell (1982), Quaternary geology of southern Michigan, The University of Michigan. Fleming, 8., A. Lavenue, A. Aly, and A. Adams (2002), Practical applications of spectral analysis to hydrologic time series, Hydrological Processes, 16, 565—574. 200 Flerchinger, G. (2000), The Simultaneous Heat and Water (SHAW) Model: Technical Documentation, Northwest Watershed Research Center, USDA Agri- cultural Research Service. Foley, J ., I. Prentice, N. Ramankutty, S. Lewis, D. Pollard, S. Sitch, and A. Hax- eltine (1996), An integrated biosphere model of land surface processes, terres- trial carbon balance, and vegetation dynamics, Global Biogeochemical Cycles, 10(4), 603*628. Frigo, M ., and S. Johnson (1998), F F TW: An adaptive software architecture for the FFT, in Proceedings of the International Conference on Acoustics, Speech, and Signal Processing, vol. 3, pp. 1381—1384. Fulton, R., J. Breidenbach, D.-J. Seo, and D. Miller (1998), WSR-88D rainfall algorithm, Weather Forecasting, 13, 377—395. Gale, M., and D. Grigal (1987), Vertical root distributions of northern tree species in relation to successional status, Canadian Journal of Forest Research, 17, 829—834. Gaucherel, C. (2002), Use of wavelet transform for temporal characterization of remote watersheds, Journal of Hydrology, 269, 101—121. Gedney, N ., P. Cox, R. Betts, O. Boucher, C. Hunting‘ford, and P. Stott (2006), Detection of a direct carbon dioxide effect in continental river runoff records, Nature, 439(7078), 835—838. Gesch, D., M. Oimoen, S. Greenlee, C. Nelson, M. Steuck, and D. Tyler (2002), The National Elevation Dataset, Journal of the American Society for Pho- togrammetry and Remote Sensing, 68(1), 5—11. Goff, J ., and S. Gratch (1946), Low-pressure properties of water from -160 to 212f, Transactions of the American Society of Heating and Ventilating Engi- neers, 52, 95—122. Halliday, D., R. Resnick, and J. Walker (2004), Fundamentals of Physics, Wiley. Hameed, S. (1984), Fourier analysis of Nile flood levels, Geophysical Research Letters, 1(9), 843—845. Harbaugh, A., E. Banta, M. Hill, and M. McDonald (2000), MODFLOW-2000, the us. geological survey modular ground-water model- user guide to modu- larization concepts and the ground water flow process, Tech. rep., Open-File Report 00-92, USGS: Reston, VA. 201 Harris, F. (1978), On the use of windows for harmonic analysis with the Discrete Fourier Transform, Proceedings of the IEEE, 66, 67—67. Harrison, L. (1963), Humidity and moisture, Vol 3, chap. Fundamental concepts and definitions relating to humidity, Reinhold Publishing Co., NY. HDF (2008), HDF5 User’s Guide, The HDF Group. Healy, R., and P. Cook (2002), Using groundwater levels to estimate recharge, Hydrogeology Journal, 10(1), 91. Hillel, D. (1980), Introduction to Soil Physics, Academic Press Inc.: San Diego, CA. Holtschlag, D. (1996), A generalized estimate of ground-water—recharge rates in the lower peninsula of Michigan, Tech. rep., USGS Water-Supply Paper 2437. Hostetler, S., and P. Bartlein (1990), Simulation of lake evaporation with appli- cation to modeling lake level variations of Harney-Malheur lake, Oregon, Water Resources Research, 26 (10), 2603—2612. Hurrell, J ., Y. Kushnir, M. Visbeck, and G. Ottersen (2003), The North Atlantic Oscillation: Climate Significance and Environmental Impact, chap. An overview of the North Atlantic Oscillation, pp. 1—35, AGU Geophysical Monograph Series 134. HydroGeoLogic (2002), MODHMS—MODFLOW-based hydrologic modeling system: Documentation and user’s guide, Tech. rep., HydroGeoLogic Inc., Hern- don, VA. Hyndman, D. W., A. D. Kendall, and N. R. H. Welty (2007), Subsurface Hydrol- ogy: Data Integration for Properties and Processes, chap. Evaluating Temporal and Spatial Variations in Recharge and Streamflow Using the Integrated Land- scape Hydrology Model (ILHM), pp. 121—141, AGU Geophysical Monograph Series 171. IPCC (2001), Contribution of Working Group I to the Third Assessment Report of the Intergovernmental Panel on Climate Change, chap. IPCC 2001, Climate Change 2001: The Scientific Basis, Cambridge University Press: Cambridge. Jackson, R., J. Canadell, and J. Ehleringer (1996), A global analysis of root distributions for terrestrial biomes, Oecologia, 108, 389—411. Jayawickreme, D., and D. W. Hyndman (2007), Evaluating the influence of land cover on seasonal water budgets using Next Generation Radar (N EXRAD) rainfall and streamflow data, Water Resources Research, 43(2), W02,408. 202 Johanson, R., J. Kittle, and A. Donigian (1984), Hydrological Simulation Program-FORTRAN (HSPF): User’s Manual for Release 8.0, EPA-600/3—84- 066 US. Environmental Research Laboratory, Athens, CA. J ukic, D., and V. Denic—Jukic (2004), A frequency domain approach to ground- water recharge estimation in karst, Journal of Hydrology, 289, 95—110. Kendall, A. D., and D. W. Hyndman (2007), Subsurface Hydrology: Data Inte- gration for Properties and Processes, chap. Examining Watershed Processes Us- ing Spectral Analysis Methods Including the Scaled-Windowed Fourier Trans- form, pp. 183-200, AGU Geophysical Monograph Series 171. Kerby, W. (1959), Time of concentration for overland flow, Civil Engineering, 29(3), 60. Kirchner, J. (2005), Aliasing in 1/falpha noise spectra: Origins, consequences, and remedies, Physical Review, 71 (6), 066,110(16). Kirchner, J., X. Feng, and C. Neal (2000), Fractal stream chemistry and its implications for contaminant transport in catchments, Nature, 4 03, 524—526. Knyazikhin, Y., et al. (1999), MODIS Leaf Area Index (LAI) and Fraction of Photosynthetically Active Radiation Absorbed by vegetation (F PAR) product (MOD15) algorithm, Tech. rep., Theoretical Basis Document, NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA. Kottek, M., J. Grieser, C. Beck, B. Rudolf, and F. Rubel (2006), World map of the Koppen-Geiger climate classification updated, Meteorologische Zeitschrift, 15(3), 259—263. Krause, P., D. Boyle, and F. Base (2005), Comparison of different efficiency creteria for hydrological model assessment, Advances in Geosciences, 5, 89—87. Labat, D., Y. Godderis, J. L. Probst, and J. L. Guyot (2004), Evidence for global runoff increase related to climate warming, Advances in Water Resources, 27, 631—642. Lawrence, M. (2005), The relationiship between relative humidity and the dew- point temperature in moist air: A simple conversion and applications, Bulletin of the American Meteorological Society, 86(2), 225-233. Leathers, D., B. Yarnal, and M. Palecki (1991), The Pacific/ North American teleconnection pattern and the United States climate. part i: Regional temper- ature and precipitation associations, Journal of Climate, 4(5), 517—528. 203 Leavesley, G., R. Lichty, B. Troutman, and L. Saindon (), Precipitation-Runofl' Modeling System—User’s manual, US Geological Survey Water—Resources In- vestigations Report 83-4238. Lee, J ., and K. Lee (2000), Use of hydrologic time series data for identifica- tion of recharge mechanism in a fractured bedrock aquifer system, Journal of Hydrology, 229, 190—201. Liang, X., D. Lettenmaier, E. Wood, and S. Burges (1994), A simple hydrologi- cally based model of land surface water and energy fluxes for general circulation models, Journal of Geophysical Research, 99(D7), 14,415—14,428. Manfreda, S., M. Fiorentino, and V. Iacobellis (2005), DREAM: a distributed model for runoff, evapotranspiration, and antecedent soil moisture simulation, Advances in Geosciences, 2, 31—39. Mann, M. (2004), On smoothing potentially non-stationary climate time series, Geophysical Research Letters, 31, L07,214. Markstrom, S. L., R. G. Niswonger, R. S. Regan, and D. E. Prudic (2008), GSFLOW—Coupled Ground- Water and Surface- Water Flow Model Based on the Integration of the Precipitation-Runoff Modeling System (PRMS) and the Modular Ground- Water Flow Model (M 0DFLOW-2005), chap. D1, USGS Tech- niques and Methods 6. Mathworks (2008), MATLAB Computational Environment, The Mathworks, Inc. Maulem, Y. (1976), A new model predicting the hydraulic conductivity of un- saturated porous media, Water Resources Research, 12, 513—522. MAWN (2008), Hourly climate data, http://www.agweather.geo.msu.edu/mawn. Maxwell, R., and N. Miller (2005), Development of a coupled land surface and groundwater model, Journal of Hydrometeorology, 6, 233—247. McCuen, R. (2004), Hydrologic Analysis and Design, Pearson Prentice Hall: Upper Saddle River New Jersey. MDEQ (2008), WelLogic, statewide ground water database. MDEQ (2009), Michigan estimated groundwater recharge. MDEQ Bathymetry (2008), Inland lake bathymetry. 204 MDNR (2001), Integrated forest monitoring assessment and prescription (IF MAP), Tech. rep., Michigan Department of Natural Resources: Lansing, MI. MDNR (2009), Inland lake maps. Monteith, J. (1965), Evaporation and environment, Syrnposia for the Society of Experimental Biology, XIX, 205—234. Morenz, M., R. van Til, and C. Luukkonen (2005), Water use for irrigation in Michigan, 2001, USGS Fact Sheet FS 2005-3079. Myeni, R., et al. (2002), Global products of vegetation leaf area and fraction absorbed PAR from year one of MODIS data, Remote Sensing of Environment, 83,214—231. Naff, R., and A. Gutjahr (1983), Estimation of groundwater recharge parame- ters by time series analyses, Water Resources Research, 19(6), 1531—1546. Nash, J., and J. Sutcliffe (1970), River flow forecasting through conceptual models. part i: A discussion of principles, Journal of Hydrology, 10, 282—290. NASS (2007), Prospective plantings, Tech. rep., National Agriculture Statistical Service, USDA. NCDC Hourly (2008), Hourly precipitation gauge data, http: / / www.ncdc.noaa. gov / . NCDC Temp (2008), Daily temperature data, http://www.ncdc.noaa.gov/. Neitsch, S., J. Arnold, J. Kiniry, R. Srinivasan, and J. Williams (2002), Soil and Water Assessment Tool User’s Manual Version 2000, GSWRL Report 02- 02 BRC Report 02-06 Texas Water Resources Institute TR—192, College Station, TX. Niswonger, R., D. Prudic, and R. Regan (2006), Documentation of the Unsaturated-Zone Flow (UZFl) package for modeling unsaturated flow between the land surface and the water table with MODFLOW-2005, in U. S. Geological Techniques and Methods Book 6, Chapter A19, p. 62. O’Callaghan, J., and D. Mark (1984), The extraction of drainage networks from digital elevation data, Computer Vision, Graphics, and Image Processing, 28, 3239344. O’Neal, R. (1997), Muskegon River Watershed assessment, Tech. rep., Michigan Department of Natural Resources, Fisheries Special Report 19. 205 Panday, S., and P. Huyakorn (2004), A fully coupled physically-based spatially- distributed model for evaluating surface/subsurface flow, Advances in Water Resources, 27(4), 361—382. Patterson, J ., and P. Hamblin (1988), Thermal simulation of a lake with winter ice cover, Limnology and Oceanography, 33(3), 323-338. Pelletier, J ., and D. Turcotte (1997), Long-range persistence in climatological and hydrological time series: Analysis, modeling and application to drought hazard assessment, Journal of Hydrology, 203, 198—208. Percival, D., and A. Walden (1993), Spectral Analysis for Physical Applications, Cambridge: Cambridge University Press. Peterson, E., and J. Hennessey, Jr. (1978), On the use of power laws for esti- mates of wind power potential, Journal of Applied Meteorology, 17, 390—394. Piajnowski, B., D. Ray, A. Kendall, J. Duckles, and D. Hyndman (2007), Using backcast land-use changes and groundwater travel-time models to generate land- use legacy maps for watershed management, Ecology and Society, 12(2), 25 [online]. Pimentel, D., and T. Patzek (2005), Ethanol production using corn, switchgrass, and wood; biodiesel production using soybean and sunflower, Natural Resources Research, 14(1), 65—76. Press, W., S. Teukolsky, W. Vetterling, and B. F lannery (1992), Numerical Recipes in FORTRAN 77, Cambridge: Cambridge University Press. Probst, J ., and Y. Tardy (1987), Long range stream-flow and world continental runoff fluctuations since the beginning of this century, Journal of Hydrology, 94, 289—311. Prudic, D., L. Konikow, and E. Banta (2004), A new Stream-Flow Routing (SFRl) package to simulate stream-aquifer interaction with MODFLOW-2000, Tech. rep., Open-File Report 2004-1042, USGS: Reston, VA. Qu, Y., and C. Duffy (2007), A semidiscrete finite volume formulation for mul- tiprocess watershed simulation, Water Resources Research, 43, W08,419. Querner, E. (1997), Description and application of the combined surface and groundwater flow model MOGROW, Journal of Hydrology, 192(1-4), 158-188. R Development Core Team (2005), R: A language and environment for statisti- cal cornputing, R Foundation for Statistical Computing, Vienna, Austria, ISBN 3-900051-07-0. 206 Rajagopalan, B., E. Cook, U. Lall, and B. Ray (2000), Spatiotemporal vari- ability of ENSO and SST teleconnections to summer drought over the United States during the twentieth century, Journal of Climate, 13, 4244—4255. Reed, S., V. Koren, M. Smith, Z. Zhang, F. Moreda, D.-J. Seo, and D. Partici- pants (2004), Overall distributed model intercomparison project results, Journal of Hydrology, 298(1-4), 27-60. Refsgaard, J., and B. Storm (1995), Computer models of watershed hydrology, chap. MIKE SHE, pp. 809—846, Water Resources Publication Co. USA. Ritchie, J. (1972), Model for predicting evaporation from a row crop with in- complete cover, Water Resources Research, 8(5), 1204—1213. Rodionov, S., and R. Assel (2000), Atmospheric teleconnection patterns and the severity of winter in the Laurentian Great Lakes basin, Atmosphere-Ocean, 38(4), 601—635. Said, A., D. Stevens, and G. Sehlke (2005), Estimating water budget in a re- gional aquifer using HSPF -MODF LOW integrated model, Journal of the Amer- ican Water Resources Association, 41(1), 55—66. Saxton, K., and W. Rawls (2006), Soil water characteristic estimates by texture and organic matter for hydrologic solution, Soil Science Society of America Journal, 70, 1569—1578. Scanlon, B., R. Healy, and P. Cook (2002), Choosing appropriate techniques for quantifying groundwater recharge, Hydrogeology Journal, 10(1), 18—39. Schaap, M., and F. Leij (2000), Improved prediction of unsaturated hydraulic conductivity with the Maulem-van Genuchten model, Soil Science Society of America Journal, 64, 843-851. Schaap, M., F. Liej, and M. van Genuchten (1998), Neural network analysis for hierarchical prediction of soil water retention and saturated hydraulic conduc- tivity, Soil Science Society of America Journal, 62, 847—855. Schaetzl, R., and D. Tomczak (2001), Wintertime temperatures in the fine- textured soils of the Saginaw Valley, Michigan, The Great Lakes Geographer, 8(2), 87—99. Schulze, E., F. Kelliher, C. K6rner, J. Lloyd, and R. Leuning (1994), Relation- ships among maximum stomatal conductance, ecosystem surface conductance, carbon assimilation rate, and plant nitrogen nutrition: a global ecology scaling exercise, Annual Reviews of Ecological Systems, 25, 629—660. 207 Seo, D. (1998), Real-time estimation of rainfall fields using radar rainfall and rain gage data, Journal of Hydrology, 208, 37—52. Seo, H., J. Simunek, and E. P. Poeter (2007), Documentation of the HYDRUS package for MODFLOW-2000, the U.S. Geological Survey modular ground- water model, Tech. Rep. GWMI 2007-01, Int. Ground Water Modeling Center. Simunek, J ., M. van Genuchten, and M. Sejna (2005), The HYDRUS-1D Soft- ware Package for Simulating the One-Dimensional Movement of Water, Heat, and Multiple Solutes in Variably-Saturated Media, Department of Environmen- tal Sciences, University of California Riverside. Smith, M., D.-J. Seo, V. Koren, S. Reed, Z. Zhang, Q. Duan, F. Moreda, and S. Cong (2004), The distributed model intercomparison project (DMIP): motivation and experiment design, Journal of Hydrology, 298(1-4), 4—26. Smith, S. (1997), The Scientist and Engineer’s Guide to Digital Signal Proce- sesing, California Technical Publications. Sophocleous, M., and S. Perkins (2000), Methodology and application of com- bined watershed and ground-water models in Kansas, Journal of Hydrology, 236, 185—201. Stieglitz, M., D. Rind, J. Famiglietti, and C. Rosenzweig (1997), An efficient approach to modeling the topographic control of surface hydrology for regional and global climate modeling, Journal of Climate, 10, 118—137. Tang, Z., B. Engel, B. Pijanowski, and K. Lim (2005), Forecasting land use change and its environmental impact at a watershed scale, Journal of Environ- mental Management, 76, 35-45. Tarboton, D. (1997), A new method for the determination of flow directions and upslope areas in grid digital elevation models, Water Resources Research, 33(2), 309—319. Tarboton, D., and C. Luce (1996), Utah Energy Balance snow accumulation and melt model (UEB), computer model technical description and user’s guide, Tech. rep., Utah Water Research Laboratory and USDA Forest Service Inter- mountain Research Station. Tessier, Y., S. Lovejoy, P. Hubert, and D. Schertzer (1996), Multifractal analysis and modeling of rainfall and river flows and scaling, causal transfer functions, Journal of Geophysical Research, 101 (D21), 26,427—26,440. 208 Therrien, R., R. McLaren, E. Sudicky, and S. Panday (2004), HydroGeoSphere, A Three-Dimensional Numerical Model Describing Fully-Integrated Subsurface and Surface Flow and Solute Transport, Groundwater Simulation Group, Uni- versity of Waterloo, Ont. Canada. Thompson, D., and J. Wallace (1998), The Arctic Oscillation signature in the wintertime geopotential height and temperature fields, Geophysical Research Letters, 25(9), 1297—1300. Thorns, R., R. Johnson, and R. Healy (2006), User’s guide to the Variably Saturated Flow (VSF) process for MODFLOW, Tech. rep., U.S. Geological Survey Techniques and Methods 6-A18. Torbick, N., J. Qi, G. Roloff, and R. Stevenson (2006), Investigating impacts of land-use land cover change on wetlands in the Muskegon River Watershed, Michigan, USA, Wetlands, 26(4), 1103—1113. Torrence, C., and G. Compo (1997), A practical guide to wavelet analysis, Bulletin of the American Meteorological Society, 79, 61—78. Trenberth, K. (1997), The definitoin of El Nino, Bulletin of the American Me- teorological Society, 78(12), 2771—2777. Tucker, C., J. Pinzon, M. Brown, D. Slayback, E. Pak, R. Mahoney, E. Vermote, and N. El Saleous (2005), An extended AVHRR 8-km N DVI data set compat- ible with MODIS and SPOT vegetation NDVI data, International Journal of Remote Sensing, 26{ 20), 4485—5598. USACE-HEC (1998), HEC-HMS Hydrologic Modeling System user’s manual, U.S. Army Corps of Engineers Hydrologic Engineering Center, Davis, CA. USDA-NRCS (2008), Soil survey geographic (SSURGO) database. USDA—SCS (1993), State Soil Geographic ( STA TSGO) Database, Miscellaneous Publication No. 1492, U.S. Government Printing Office, Washington, DC. USGS (2005), National Water Information System (NWISWeb), data avialable on the World Wide Web. USGS (2008), National Water Information System (NWISWeb), data avialable on the World Wide Web. van Genuchten, M. (1980), A closed-form equation for predicting the hydraulic conductivity of unsaturated soils, Soil Science Society of America Journal, 44, 892—898. 209 Vandeeraak, J. (1999), Numerical simulation of flow and chemical transport in integrated surface-subsurface hydrologic systems, Ph.D. thesis, Dept. of Earth Sciences, University of Waterloo: Ontario, Canada. Vandeeraak, J ., and K. Loague (2001), Hydrologic-response simulations for the R—5 catchment with a comprehensive physics-based model, Water Resources Research, 37 (4), 999—1013. Visbeck, M., J. Hurrell, L. Polvani, and H. Cullen (2001), The North Atlantic Oscillation, past present and future, Proceedings of the National Academy of Sciences, 98(23), 12,876—12,877. Vorosmarty, C., P. Green, J. Salisbury, and R. Lammers (2000), Global water resources: vulnerability from climate change and population growth, Science, 298, 284—288. Voss, C., and A. Provost (2002), A model for saturated-unsaturated variable- density ground-water flow with solute or energy transport, Water-Resources Investigations Report 02-4231, U.S. Geological Survey. Walker, G., and E. Bliss (1932), World weather V, Memorandum of the Royal Meteorological Society, 4, 53—84. Walko, R., and C. Tremback (2005), Modifications for the transition from LEAF-2 to LEAF-3, Tech. rep., ATMET Technical Note. Walko, R., et al. (2000), Coupled atmosphere-biophysics—hydrology models for environmental modeling, Journal of Applied Meteorology, 39, 931—944. Wallace, J ., and D. Gutzler (1981), Teleconnections in the geopotential height field during the northern hemisphere winter, Monthly Weather Review, 109(4), 784-812. Wang, X., and A. M. Melesse (2006), Effects of STATSGO and SSURGO as inputs on SWAT model’s snowmelt simulation, Journal of the American Water Resources Association, 42(5), 1217—1236. Westjohn, D., and T. Weaver (1998), Hydrogeologic framework of the Michigan Basin regional aquifer system, Tech. rep., US Geological Survey Professional Paper 1418. Wigmosta, M., L. Vail, and D. Lettenmaier (1994), A distributed hydrology- vegetation model for complex terrain, Water Resources Research, 30(6), 1665— 1679. 210 Wiley, M., and P. Richards (2006), Unpublished streamflow data, cedar Creek. Winter, T. (1981), Uncertainties in estimating the water balance of lakes 1, Journal of the American Water Resources Research, 17(1), 82—115. Yang, Z., and D. Han (2006), Derivation of unit hydrograph using a transfer function approach, Water Resources Research, 42, W01,501. Yang, Z.-L., R. Dickinson, A. Robock, and K. Vinnikov (1997), Validation of the snow sub-model of the Biosphere-Atmosphere Transfer Scheme with Russian snow cover and meteorological observational data, Journal of Climate, 10, 353— 373. Yeh, G., and G. Huang (2003), A numerical model to simulate water flow in watershed systems of 1-d stream-river network, 2-d overland regime, and 3-d subsurface media (WASH 123D: Version 1.5), Technical report, Dept. of Civil and Environmental Engineering University of Central Florida, Orlando, Florida. Yongjiu, D., Z. Xubin, and R. Dickinson (2001), Common Land Model (CLM): Technical Documentation and User’s Guide, National Center for Atmospheric Research. Yu, Z., D. Pollard, and L. Cheng (2006), On continental-scale hydrologic simu- lations with a coupled hydrologic model, Journal of Hydrology, 331 (1-2), 110— I24. Zhang, Y., and K. Schilling (2004), Temporal scaling of hydraulic head and river base flow and its implications for groundwater recharge, Water Resources Research, 40(3), W03,504. 211 IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII 4Hll(HIJILIIJILLILIIINH(HlflllllfllllHUHIH 62 5986