LAND-USE LEGACY EFFECT: COMBINING SPATIAL AND TEMPORAL DRIVERS IN STATISTICAL AND MECHANISTIC MODELS OF LAKE WATER CHEMISTRY By Sherry L. Martin A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of DOCTOR OF PHILOSOPHY Fisheries and Wildlife 2010 ABSTRACT LAND-USE LEGACY EFFECT: COMBINING SPATIAL AND TEMPORAL DRIVERS IN STATISTICAL AND MECHANISTIC MODELS OF LAKE WATER CHEMISTRY By Sherry L. Martin Lake Classification: A classification system is often used to reduce the number of different ecosystem types that governmental agencies are charged with monitoring and managing. We compare the ability of several different hydrogeomorphic (HGM) classifications to group lakes for water chemistry/clarity. We ask three questions: (1) Which approach to lake classification (regionalization, landscape position, lake-specific, or some combination) is most successful at classifying lakes for similar water chemistry/clarity? (2) Which HGM features are most strongly related to the lake classes? and, (3) Can a single classification successfully classify lakes for all of the water chemistry/clarity variables examined? We use classification and regression tree (CART) analysis of HGM features to classify six water chemistry/clarity variables from 151 minimally disturbed lakes in Michigan USA. We developed two CART models for each water chemistry/clarity variable: HGM characteristics alone and HGM characteristics combined with regionalizations and landscape position. The combined CART models had the highest strength of evidence (ωi range 0.92-1.00) and maximized within class homogeneity (ICC range 36-66%) for all water chemistry/clarity variables except water color and chlorophyll a. The most successful single classification in our study was on average 20% less successful in classifying other water chemistry/clarity variables. Thus, our results show that no single classification maximizes success for all lake variables examined. Therefore, we suggest that the most successful classification is (1) specific to one response variable, and (2) capable of incorporating information at multiple spatial scales and from a variety of different sources (regionalization and local HGM variables). Land use legacies: The recognition of legacy effects from historical land use/land cover (LULC) is a conceptual advance that has clarified the relationship between LULC and ecosystem responses. Legacy effects can be defined as effects which perpetuate beyond an expected or perceived endpoint in time. The goal of our research was to investigate LULC legacy effects on lake water chemistry. Water chemistry and five time steps of LULC data were collected from 35 lakes in the Huron River Watershed, Michigan. We took both a correlational and mechanistic approach to represent how temporal changes in LULC influence lake water chemistry. We used principal components of LULC over time to build hierarchical regression models linking to water chemistry. We also created a mechanistic groundwater flow model to estimate spatiallyexplicit groundwater travel times. The groundwater travel time was used to create a legacy LULC map for subsequent regression modeling. Our correlative models show that some water chemistry characteristics show a stronger link to legacy LULC than others and may be explained by the solubility and reactivity of the chemical. Our mechanistic models offer insights about how groundwater interacts with LULC change to create legacy effects and show how naturally occurring conservative tracers can provide a basis for comparison against nutrient relationships to the landscape. By categorizing the chemistry variables by their key characteristics of solubility and reactivity, we are better equipped to explore other mechanisms that are important for the physical transport and biogeochemical transformations of these chemicals. DEDICATION From FCG to PhD iv ACKNOWLEDGEMENTS Thanks to the following personnel for their contributions to the development of the landscape databases: Michael Belligan, Patricia Sorrano, Mary Bremigan, Kendra Cheruvelil, Jim Breck, Howard Wandell, Gary Weissman, Stephen Bowman, Tyler Rosa, Remy Brim, Cassie Meier, Dave Meyers and Sarah Wills. I wish to thank Erica Mize, Jason Karl, Gary Roloff, Matt Baker, and John Van Sickle for invigorating discussions about statistical design; Lidia Szabo-Kraft, Arthur Cooper, and Anthony Kendall for assistance with GIS; Dave Weed, Lauren Bailey, Bryan Burroughs, Jordan Burroughs, Emi Fergus, Jon Hansen, Geoff Horst, Justin Miller, Emily Norton, Elizabeth St James, and Dan Wieferich for assistance with field collections. I thank my committee – Dan Hayes, Brian Maurer, Jan Stevenson, and Dave Hyndman – for their guidance during my PhD experience. Your challenges and support helped to create my success. Dr. Horton H. Hobbs III mentored me during my time at Wittenberg University. He has been a great supporter all along the way. I am very grateful to have been a member of the Order of the Claw and the Limnuts. Some of my fondest memories are from your classes and field trips! Thank you for starting me down this path. v TABLE OF CONTENTS LIST OF TABLES…………………………………………………………………...…….….. viii LIST OF FIGURES…………………………………………………………………………....... ix CHAPTER 1 COMPARING HYDROGEOMORPHIC APPROACHES TO LAKE CLASSIFICATION: ISSUES OF SPATIAL SCALE AND PRACTICALITY Introduction……………………………………………………………………….……… 2 Methods………………………………………………………………………….……….. 6 Hydrogeomorphic characteristics…………………………………………….….. 7 Classification frameworks………………………………………………….……. 8 Comparing classifications…………………………………………………...….. 10 Results………………………………………………………………………………...… 11 Comparing classifications……………………………………………………..... 11 Relationships between hydrogeomorphic features and lake classes………...….. 13 Evaluation of CART tree stability…………………………………………...…. 14 A single classification for all lake water characteristics……………………...… 15 Discussion…………………………………………………………………………...….. 16 Comparing classifications…………………………………………………...….. 17 Relationships between hydrogeomorphic features and lake classes………...….. 18 A single classification for all lake water characteristics…………………...…… 21 Applications to ecosystem management…………………………………...…… 22 Appendix……………………………………………...……………………………...…. 34 References…………………………………………………………………………...….. 40 CHAPTER 2 THE LAND-USE LEGACY EFFECT: ADDING TEMPORAL CONTEXT TO UNDERSTANDING LAKE WATER CHEMISTRY Introduction…………………………………………………………...………………… 48 Methods…………………………………………………………………………...…….. 50 Study Area………………………………………………………………...……. 50 Description of data……………………………………………………...………. 51 Statistical analyses…………………………………………………………...…. 53 Results………………………………………………………………………...………… 54 Changes in land use/cover…………………………………………………….... 54 Comparing regression models……………………………………………...…… 56 Discussion……………………………………………………………...……………….. 57 Management Implications…………………………………………………...….. 60 Appendix……………………………………………………………………………..…. 69 References…………………………………………………………………………...….. 83 vi CHAPTER 3 THE LAND-USE LEGACY EFFECT: A MECHANISTIC INVESTIGATION INTO HOW GROUNDWATER DELIVERY, WETLAND PROCESSING, AND RIPARIAN DYNAMICS AFFECT LAKE WATER CHEMISTRY Introduction………………………………………………………………………...…… 89 Methods………………………………………………………………………………..... 91 Study Area……………………………………………………………...………. 91 Water chemistry data…………………………………………………………….92 Land use/cover data………………………………………………………...…... 93 Groundwater travel time…………………………………………………...…… 93 Creating the legacy map………………………………………………...………. 96 Statistical analyses…………………………………………………………...…. 96 Results……………………………………………………………………………...…… 99 Groundwater travel time………………………………………………...……… 99 Characterizing LULC…………………………………………………………… 99 Relationships with lake water chemistry…………………………………..... 100 Discussion……………………………………………………………………………. 102 Wetland processing………………………………………………………….… 103 Riparian zones……………………………………………………………….… 105 Conclusions……………………………………………………………….…… 107 References…………………………………………………………………….……….. 120 vii LIST OF TABLES Table 1.1. Summary of lake water chemistry/clarity and hydrogeomorphic characteristics for 151 minimally disturbed lakes in Michigan USA. Variables are arranged into broad categories. Abbreviations are listed in parentheses. PCU, platinum cobalt units. CA:LK, catchment area to lake area ratio. CUCA:LOCA, cumulative catchment area to local catchment area ratio…………………………………………………………………….. 23 Table 1.2. Summary of model selection statistics for candidate classifications per lake water characteristic. Classification type is indicated as regionalization, landscape position, or CART. Individual classification names are indicated within each classification type. Intra-class correlation coefficients (ICC) are presented as percent of total variance that is among the classes. A high ICC indicates high within class homogeneity and thus, high classification success. The number of classes per classification model (K) is presented. ∆AICC is the difference between the AICC for each model and the minimum AICC for each lake water characteristic. The ∆AICC will equal 0 for the best model per lake water characteristic. The Akaike weights (ωi) sum to 1 for each lake water characteristic and is interpreted as the likelihood that a given model is the best model relative to others included in the analysis. n/a, not applicable……………………………………………. 25 Table 1.3. Summary of the intraclass correlation coefficients (ICCs) for CART models across lake water characteristics. CART models are listed by the original dependant variable and type of classification. Mean ICC across lake water characteristics is computed. Rank of mean ICC is listed. See Table 2 for acronyms……………………………………… 27 Table 1.4. Pearson product-moment correlation matrix for lake water characteristics………… 28 Table 2.1. Eigenvectors and proportional variance explained from principal components analysis of land use/cover types. Extreme values for each PC are indicated in bold………….… 62 Table 2.2. Regression statistics for best models. Dependant variables are listed under broad categories. Details of regression models include: time step(s) represented in final model, number of regression parameters fit (p), Akaike Information Criteria (AIC), Akaike weight (ωi), and coefficient of determination (R2)…………………….………………. 63 Table 2.3. Characteristics of the study lakes including water chemistry variables as categorized into a priori groups and lake area. Minimum, maximum, mean, and coefficient of variation are shown for each variable……………………………………………………70 Table 2.4. Regression coefficients from best models as reported in Table 2.2……………….....71 Table 3.1. Size and water chemistry characteristics of study lakes, including minimum, maximum, mean, and coefficient of variation. Phosphorus species include total phosphorus (TP), and soluble reactive phosphorus (SRP). Nitrogen species include total nitrogen (TN), nitrate (NO3), and ammonia (NH4). Reactive ions include calcium (Ca), viii silica (Si), and sulfate (SO4). Conservative ions include chloride (Cl), potassium (K), magnesium (Mg), and sodium (Na)…………………………….……………………. 108 ix LIST OF FIGURES Figure 1.1. Map of the upper and lower peninsula of Michigan. Lakes included in the analysis are shown as solid dots. Boundaries are shown for each regionalization: A) 8-digit USGS hydrologic units (HUC), B) ecological drainage unit (EDU), and C) hydrological landscape region (HLR)………………………………………………………………… 29 Figure1.2. Results from the CART analysis of local hydrogeomorphic features (HGM) in classifying lake water characteristics: (A) alkalinity, (B) water color, (C) Secchi, (D) TN, (E) TP, and (F) Chl a. Each split is labeled with the splitting variable (see Table 1 for abbreviations), and proportional reduction in error (PRE). Branches are labeled with splitting value. Terminal nodes (rectangles) represent lake classes and are labeled with an alphabetical class, class mean, and number of lakes per class (in parentheses)………. 30 Figure 1.3. Results from the CART analysis combining regionalization, landscape position and local hydrogeomorphic features (HGM+) in classifying lake water characteristics: (A) alkalinity, (B) water color, (C) Secchi, (D) TN, (E) TP, and (F) Chl a. Each split is labeled with the splitting variable (see Table 1 for abbreviations), and proportional reduction in error (PRE). Branches are labeled with splitting value. Branches split using categorical variables (e.g., HUC, EDU) are detailed in appendices. Terminal nodes (rectangles) represent lake classes and are labeled with an alphabetical class, class mean, and number of lake per class (in parentheses)……………………..…………………… 32 Figure 1.4. Details of HUC and EDU categorical splits in the alkalinity HGM+ model (refer to Figure 3): A) first node splitting further into class A and towards classes B and C, and B) second node splitting further to classes B or C. Grey areas indicate HUCs which split to the left. Striped areas indicate HUCs which split to the right……………………..……. 35 Figure 1.5. Details of HUC categorical split in the water color HGM+ model (refer to Figure 3). Black areas indicate HUCs which were not represented at the node. Grey areas indicate HUCs which split to the left into class B. Striped areas indicate HUCs which split to the right towards classes C and D……………………………………………………………36 Figure 1.6. Details of HUC categorical splits in the Secchi HGM+ model (refer to Figure 3): A) left split from WRT splitting further to classes B and C, and B) right split from WRT splitting further to classes D and E. Black areas indicate HUCs which were not represented at the node. Grey areas indicate HUCs which split to the left into classes B or D. Striped areas indicate HUCs which split to the right into classes C or E…………….37 Figure 1.7. Details of HUC categorical split in the TN HGM+ model (refer to Figure 3). Grey areas indicate HUCs which split to the left towards classes A, D, E, and F. Striped areas indicate HUCs which split to the right towards classes B and C……………………..….38 x Figure 1.8. Details of HUC categorical split in the TP HGM+ model (refer to Figure 3). Grey areas indicate HUCs which split to the left into class B. Striped areas indicate HUCs which split to the right into class C………………………………………………………39 Figure 2.1. Map showing the Huron River Watershed within the state of Michigan, including detailed hydrography features along with outlines for the cities of Ann Arbor and Ypsilanti. For interpretation of the references to color in this and all other figures, the reader is referred to the electronic version of this thesis (or dissertation)………...……. 64 Figure 2.2. Average proportion of each land use/cover type for the study lakes calculated for each time step in the study. Solid fine line represents agriculture. Solid heavy line represents urban. Dashed fine line represents forest. Dashed heavy line represents open. Dotted fine line represents water. Dotted heavy line represents wetland………………. 65 Figure 2.3. Change in land use/cover over time represented by principal component 1 (PC1), principal component 2 (PC2), and principal component 3 (PC3). Each dot represents the average score for that time step………………………………………...…………….. 66 Figure 2.4. Akaike weights (ωi) from chrono-sequence of regression models. (A) TP, solid circles with heavy solid line; and, SRP, open circles with dashed line. (B) TN, solid circles with heavy solid line; NO3, open circles with dashed line; and, NH4, X inside circle with dotted line. (C) Ca, open circle with heavy solid line; Si, open triangle with dashed line; and, SO4, open square with dotted line. (D) Cl, open with circle heavy solid line; K, open triangle with dashed line; Mg, open square with dotted line; and, Na, open diamond with fine solid line………………………………………...………………….. 67 Figure 2.5. Coefficient of determination (R2) from chrono-sequence of regression models. (A) TP, solid circles with heavy solid line; and, SRP, open circles with dashed line. (B) TN, solid circles with heavy solid line; NO3, open circles with dashed line; and, NH4, X inside circle with dotted line. (C) Ca, open circle with heavy solid line; Si, open triangle with dashed line; and, SO4, open square with dotted line. (D) Cl, open with circle heavy solid line; K, open triangle with dashed line; Mg, open square with dotted line; and, Na, open diamond with fine solid line……………………………………………………………. 68 Figure 2.6. Comparisons between regression results using 2 PCs and 3 PCs for calcium, chloride, potassium, magnesium, sodium, ammonium, nitrate, silica, sulfate, soluble reactive phosphorus (SRP), total nitrogen (TN), and total phosphorus (TP). Akaike weights (yaxis) were computed across regression results for each water chemistry variable from both 2 PCs and 3 PCs to total ten regression models………………………………...…..72 Figure 2.7. Change in LULC over time as represented in multivariate space. Each line represents the LULC trajectory of a lake in the study. Each of the first three PC’s are shown……..78 Figure 2.8. Regression diagnostics for best models as reported in Table 2.2. Observed values (yaxis) are plotted against values predicted from the regression equation (x-axis)………..81 xi Figure 3.1. Map showing the Huron River Watershed within the state of Michigan, including detailed hydrography features along with outlines for the cities of Ann Arbor and Ypsilanti……………………………………………………………………….………. 109 Figure 3.2. Map showing land use/cover (LULC) from 1938, 1955, 1968, 1978, and 1996 within the Huron River Watershed…………………………………...………………………. 110 Figure 3.3. Surficial geology zones and optimized hydraulic conductivity (K, m/d) values used in the final model. Glacial Outwash = glacial outwash sand and gravel and postglacial alluvium. Ice Contact = ice contact outwash sand and gravel. Percent of each type within the Huron River groundwater sourceshed (HRGW) upstream of Ypsilanti, and the mean percent for the study lake groundwater sourcesheds are shown. The HRGW is dominated by glacial outwash. The study lakes had similar geology as the HRGW, in general…. 111 Figure 3.4. Observed versus predicted water levels (m) recorded for 15,581 wells used to estimate values for hydraulic conductivity (K). The line of best fit and coefficient of determination are shown………………………………………………...…………….. 112 Figure 3.5. Steps to create a land use legacy map. Step 1: The simulated groundwater travel times in years were calculated for each model location down gradient to a discharge point (e.g. river) using Darcy’s Law, as detailed in Methods. Step 2: Groundwater travel time year was reclassified to bracket times with available LULC data (in this case, 1996, 1978, 1968, 1955, 1938, or Pre-1900). Step 3: LULC type for each cell was assigned using time categories from step 2. Step 4: LULC from 1996, showing areas that differ with the legacy LULC with black cross-hatching…………...………………………………….. 113 Figure 3.6. Groundwater travel times reclassified to represent the 5 time steps of land cover data. Time classes span from midpoint to midpoint around a specified year. Percent of each groundwater travel time category within the Huron River groundwater sourceshed (HRGW), and the mean percent for the study lake groundwater sourcesheds are shown. The HRGW is dominated by groundwater needing less than 10 years to travel from it’s source to a surface water delivery pathway. The study lakes had similar groundwater travel times as the HRGW, in general…………………………………………..…….. 114 Figure 3.7. Land use/cover within the model area showing A) 1996 and B) legacy land use/cover. Percent of each LULC type within the Huron River groundwater sourceshed (HRGW), and the mean percent for the study lake groundwater sourcesheds are shown. Maximums are shown in bold. The HRGW is dominated by urban LULC type in 1996, but agriculture dominates in the legacy map. The study lakes had similar LULC representation as the HRGW, in general…………………………………………….… 115 Figure 3.8. Areas within the Huron River groundwater sourceshed (HRGW) with differences between 1996 LULC and legacy LULC. The model boundary is indicated with the bold black line. Study lake groundwater sourcesheds are indicated with dark grey lines and light grey shading. Red areas indicate where legacy land cover differs from 1996 land xii cover (10.5% of the area within the HRGW, 8.4% of the area within the study lake groundwater sourcesheds)…………………………………………………………….. 116 Figure 3.9. Coefficients of determination (R2) from regression models of total phosphorus (TP), soluble reactive phosphorus (SRP), total nitrogen (TN), nitrate (NO3), ammonia (NH4), silica (Si), calcium (Ca), sulfate (SO4), chloride (Cl), potassium (K), magnesium (Mg), and sodium (Na) using correlative legacy models from Chapter 2, LULC characterized by the current time step (current LULC), and LULC characterized by the legacy map. Asterisks indicate model chosen by AIC weights……………………………...……… 117 Figure 3.10. Coefficients of determination (R2) from regression models of total phosphorus (TP), soluble reactive phosphorus (SRP), total nitrogen (TN), nitrate (NO3), ammonia (NH4), silica (Si), calcium (Ca), sulfate (SO4), chloride (Cl), potassium (K), magnesium (Mg), and sodium (Na) using LULC characterized by the GW legacy models measuring wetlands in three ways: all wetlands in the watershed (GW legacy), area of wetland connected to rivers, and length of wetland connected to rivers. Wetland area includes the entire area of the wetland polygons which connect to rivers. Wetland length differs from wetland area by only accounting for the area of the connected wetland polygons within a 50m buffer of rivers…………………………………………………………………… 118 Figure 3.11. Coefficients of determination (R2) from regression models of total phosphorus (TP), soluble reactive phosphorus (SRP), total nitrogen (TN), nitrate (NO3), ammonia (NH4), silica (Si), calcium (Ca), sulfate (SO4), chloride (Cl), potassium (K), magnesium (Mg), and sodium (Na) using correlative legacy models from chapter 2, LULC characterized by the legacy map, LULC within the 50m riparian zone buffer only, and legacy LULC combined with LULC from the 50m riparian zone buffer. Asterisks indicate model chosen by AIC weights. The asterisk for nitrate indicates the correlative legacy model………………………………………………………...………………………… 119 xiii CHAPTER 1: COMPARING HYDROGEOMORPHIC APPROACHES TO LAKE CLASSIFICATION: ISSUES OF SPATIAL SCALE AND PRACTICALITY 1 INTRODUCTION Ecosystem structure and function are controlled in large part by the hydrology, geology, land cover, and climate characteristics of that ecosystem. These natural landscape features have been used to identify similarities in temperate (Host et al. 1996) and tropical (Mora and Iverson 2002) forests, rangelands (Kunst et al. 2005), streams (Frissell et al. 1986) and their riparian zones (Vidon and Hill 2004), wetlands (Brinson 1993), lakes (Winter 1977, Riera et al. 2000), and coral reefs (Rodgers 2005). Brinson (1993) outlined an approach to wetland classification based on hydrologic and geomorphic features such as precipitation, groundwater flow, and landscape position. This hydrogeomorphic (HGM) approach was intended to provide a flexible classification framework based on knowledge of how HGM factors drive ecosystem structure and function. The accumulation of HGM data and the advancement of analytical techniques capable of handling complex datasets have extended the capacity of the HGM approach, making it possible to include more characteristics in the ecological classification of more ecosystem types, beyond wetland classification (Host et al. 1996). Researchers have been classifying lakes since the early 1900’s; in fact, lake classification was the major focus of the International Congress for Limnology in 1956 (Moss et al. 1994), and interest in this topic has continued to the present. Although a wide variety of approaches to classification has been adopted, one feature common to many has been to group lakes based on the statistical similarity of specified water chemistry variables. Such a classification approach has been applied to lake water chemistry data in Canada (Pitblado et al. 1980, Zimmerman et al. 1983), Northeastern U.S. (Hunsaker et al. 1986, Newton and Driscoll 1990, Young and Stoddard 1996, Momen and Zehr 1998, and Jenerette et al. 2002) and Sweden (Hakanson 1996, Hakanson 2005). However, this site-specific approach can group only those lakes for which water 2 chemistry data are available; therefore, it has limited application for other regions or where data are lacking. To move beyond this limitation, we apply an HGM approach to lake classification which groups systems according to relationships between important ecosystem characteristics (e.g., lake water chemistry) and HGM features. With technological advancements in remote sensing and coordinated data collection strategies, large stores of HGM data are readily available in geographic information systems (GIS) for many areas of the world (Johnson and Gage 1997). Because an HGM-based classification can be applied to previously unsampled lakes without logistically challenging and expensive field collections, its usefulness can be extended to encompass many more lakes than can be physically sampled (Brinson 1993, Young and Stoddard 1996), thereby allowing inferences and predictions to be made for individual lakes across broad geographic regions. Regional land classifications (e.g., regionalizations) take advantage of the wealth of HGM data to group large geographic regions based on the similarity of physiographic, climatic and terrestrial features (Omernik 1987, Bailey et al. 1994, Albert 1995). Interestingly, unlike these land classifications, few classifications of aquatic systems have taken an HGM approach (but see Hakanson 1996, Hershey et al. 1999, Higgins et al. 1998, Wolock et al. 2004). One hydrologically driven example of regionalization is the USGS hydrologic units (HUC; Seaber et al. 1987). HUCs are delineated using topographical boundaries specific to a surface drainage area and have been used as management units by many agencies. More recent examples of HGM-based hydrological regionalizations include hydrologic landscape regions (HLR, Winter 2001, Wolock et al. 2004) and ecological drainage units (EDU, Higgins et al. 2005). The concept of hydrologic landscapes provides an aquatic analog to the land-based regionalization 3 approaches, delineating land areas with similar HGM-drivers of surface and ground water movement and storage, specifically land-surface form, geologic texture, and climatic setting (Winter 2001, Wolock et al. 2004). Alternatively, Higgins et al. (2005) delineated EDUs by combining HUC watersheds with similar climate and landscape features. Although these regionalizations are conceptually appealing, the few studies testing such regionally-based lake classifications report that some critical lake water characteristics, such as productivity, are not always similar among lakes within these aquatic and terrestrial regions (Jenerette et al. 2002, Cheruvelil et al 2008). The concept of lake landscape position describes the local hydrologic landscape of a system. Lake landscape position quantifies the hydrologic connectivity and spatial arrangement of various aquatic systems to infer similarity in ground and surface water hydrology (Kratz et al. 1997, Soranno et al. 1999, Riera et al. 2000, Martin and Soranno 2006). Several metrics of landscape position have been derived measuring various combinations of local hydrologic connectivity to other ecosystems (e.g., streams, lakes, wetlands). Each of these metrics of landscape position has shown significant relationships with important lake ecosystem characteristics, such as acid neutralizing capacity, dissolved organic carbon and nitrogen to phosphorus ratio (Kratz et al. 1997, Martin and Soranno 2006). However, many other measures of water clarity and productivity have not shown significant relationships with landscape position metrics (Riera et al. 2000, Quinlan et al. 2003, Martin and Soranno 2006). Although many of the above ecological classification schemes have demonstrated some success in classifying lakes, they do so with little regard to other important HGM features. For instance, although some regionalizations successfully group lakes with similar water quality, mechanisms that act though local scale variables, such as lake morphometry, are not 4 incorporated into such regionalizations and are likely important for lake classification success (Cheruvelil et al. 2008). Indeed, others have emphasized the importance of including both regional and local scale variables simultaneously in analyses of stream and lake characteristics (Seelbach et al. 1997, Goransson et al. 2004, Stendera and Johnson 2006) and call for an approach that combines regional and local features (Pyne et al. 2007, Cheruvelil et al. 2008). A combined approach, however, has inherent technical demands; it must be able to incorporate both continuous and categorical data, and account for local scale variation concurrently with regional scale phenomena. To date, the majority of statistical techniques that have been employed for classification development have used traditional linear models, such as principal components analysis and clustering (Emmons et al. 1999, Bryan 2006). Classifications created with these linear methods are limited statistically when including categorical variables, such as regionalization or landscape position. By including such spatially-explicit categories into a larger classification framework, additional variation in water characteristics may be captured that local HGM data alone may miss. In addition, although linear approaches have been found to accurately represent some ecological relationships, these approaches may not effectively represent non-linear relationships and may mask the true character of the data by forcing it to conform to a linear arrangement (De’Ath and Fabricius 2000, Robertson et al. 2006, Soranno et al. 2008). To date, few classification efforts have taken advantage of advances in statistical methods that alleviate these shortcomings (but see Zimmerman et al. 1983, Emmons et al. 1999, Olden and Jackson 2002, Robertson et al. 2006). The goal of our study is to build an ecological classification for lake water characteristics that is built using variation in HGM features over multiple spatial scales. More specifically, we incorporate the phenomena captured by regional summaries of HGM features (i.e., 5 regionalizations) with local HGM features that are intrinsic to each lake (e.g., lake morphometry). We strive to create a classification approach that: (a) maximizes within class homogeneity and between class heterogeneity for lake water characteristics, (b) is developed from natural landscape features that are temporally stable on the scale of decades to centuries, (c) minimizes the confounding effects of non-natural landscape features (e.g., human disturbances), and (d) provides an example of a broadly applicable approach for other ecosystems. We ask three questions: (1) Which approach to lake classification is most successful at grouping lakes with similar water characteristics (regionalization, landscape position, lake-specific, or some combination)? (2) Which HGM features are most strongly related to the lake classes? and (3) Can a single lake classification successfully group lakes for all of the water chemistry characteristics examined? METHODS Our dataset includes 151 minimally disturbed lakes in Michigan, U.S.A. that are greater than 20 hectares in area. We define minimally disturbed lakes as those with no dam or water control structure and less than 25% human land use/cover (i.e., agriculture and urban) in the cumulative catchment (detailed below). Our study lakes had only an average of 8% human land use/cover in the cumulative catchment and were surrounded mostly by forest (mean 80% forested land use/cover). We chose to limit our dataset to these lakes in order to reduce the confounding effects of human disturbances and maximize our ability to detect relationships with HGM characteristics (D’arcy and Carignan 1997, Stoddard et al. 2006). We obtained data on lake water characteristics during the time period of 1975 through 1982 from the U.S. EPA Storet database. The Michigan Department of Environmental Quality 6 sampled the epilimnion of each lake during summer stratification (July, August, and September) for a wide range of limnological variables: alkalinity, water color, Secchi disk depth, total nitrogen (TN), total phosphorus (TP), and chlorophyll a (Chl a). The study lakes vary widely in all variables (Table 1.1). Hydrogeomorphic characteristics We created a digital HGM database for our study lakes from geospatial data including bedrock geology, surficial geology, land use/cover, and climate/hydrology. These landscape features were summarized for each lake using a 500 m buffer. Bedrock geology data were obtained from the Geologic Survey Division of the Michigan Department of Environmental Quality. Sedimentary clastic is the dominant bedrock type in Michigan and for many of our study lakes, but our dataset also includes lakes dominated by other bedrock types (Table 1.1). Surficial geology data were provided by the Michigan Natural Features Inventory and Michigan Department of Natural Resources. Outwash and moraine surficial geology types dominate the study lakes (Table 1.1). Land use/cover data were obtained from the Michigan Resource Information Service (MIRIS 2000) based on aerial photo interpretation of photos taken between 1978 and 1985. Average annual precipitation for the period 1971-2000 was obtained from the Spatial Climate Analysis Service (www.ocs.oregonstate.edu). Average annual runoff for the period 1951-1980 and mean base-flow index (BFI) were obtained from USGS (http://water.usgs.gov). Base-flow index provides an estimate of groundwater input relative to surface water input. Cumulative lake catchments (CUCA) were delineated to include the catchment area associated with all lakes and streams draining into the lake using 1:100,000 resolution stream hydrography data, digital elevation models (30 m resolution) and topographic 7 maps. Using the above data, we also delineated local catchments (LOCA) as the portion of the cumulative catchment downstream from any upstream lake greater than 0.2 km2. Wetland land cover was summarized using the 500 m buffer as well as a 100 m buffer (to represent a lake’s riparian zone), and for the local and cumulative lake catchments. Lake area, perimeter, shape, mean depth, and maximum depth were gathered from bathymetric maps. Mean depth was calculated by taking the average depth of approximately 100 points evenly spaced across each bathymetric map (Omernik and Kinney 1983). Lake shape was calculated as the ratio of shoreline perimeter to the circumference of a circle of the same area (Wetzel and Likens 2000). Water residence time (WRT) was estimated as: [(lake area*mean depth) ÷ (cumulative catchment area*runoff)]. Area, shape and slope were measured for both local and cumulative catchment. Classification frameworks We included three regionalization frameworks: 1) USGS 8-digit hydrologic units (HUC, Seaber et al. 1987), 2) ecological drainage units (EDU, Higgins et al. 2005), and hydrologic landscape regions (HLR, Winter 2001). The location of each study lake within each region determined the class membership. Our study lakes were located within 19 HUCs, 6 EDUs, and 5 HLRs (Figure 1.1). We included three metrics of landscape position that can be easily measured from existing data using GIS (described in brief here, see Martin and Soranno 2006 for more detailed descriptions). Lake hydrology (LH; n=7 categories) is a general measure of lake surface hydrologic connections, incorporating both connections to streams and lakes. Lake network number (LNN; n=5 categories) measures the degree of surface connectivity to other lakes. Lake 8 network complexity (LNC; n=4 categories) is a measure of the complexity of connections to other lakes (e.g., dendritic or linear chain). We created two new classifications for each lake water characteristic using: 1) local HGM features for each study lake (HGM), and 2) local HGM features combined with regionalization and landscape position categories (HGM+). These classifications were created using classification and regression tree (CART) analysis. We chose to use CART models because they: 1) maximize class homogeneity, 2) do not penalize for including many independent variables, 3) handle high-order interactions among variables, and 4) accommodate both continuous and categorical data (De’ath and Fabricius 2000). All CART models were built using the recursive partitioning algorithm “rpart” in the R software system (R Development Core Team, http://www.R-project.org). CART trees were grown using 10-fold cross-validation and subsequently pruned using the 1-SE rule (Breiman et al 1984, Venables and Ripley 1999). Terminal nodes (i.e., lake classes) were required to have a minimum of five observations (i.e., lakes). The proportional reduction in error (PRE) for each split was summed to produce an overall PRE for each tree. Output detailing splitting decisions from each CART tree was reviewed to assess tree stability and correlations among independent variables. Independent variables maximizing class homogeneity and PRE were always selected as the primary splitter. The top five independent variables for a primary split, measured by class homogeneity, were considered as competitor splits. The top five independent variables grouping lakes into classes similar to the primary split, measured by percent similarity, were considered as surrogate splits (R Development Core Team, http://www.R-project.org). We assessed tree stability using information about competitor and surrogate splits, in combination. A split was considered 1) stable if there were no competitor 9 splits within 3% reduction in error from the primary splitter, 2) somewhat stable/unstable if there were competitor splits within 3% reduction in error from the primary splitter but these competitors were also surrogates, or 3) unstable if there were competitor splits within 3% reduction in error from the primary splitter but these competitors were not surrogates. Therefore, given small changes in input data 1) stable trees are not likely to change in tree structure or class membership, 2) somewhat stable/unstable trees may split on different independent variables yet yield similar class membership, and 3) unstable trees would likely yield different tree structure and class membership. Comparing classifications We compared the success of the three regionalization frameworks, the three landscape position metrics, and our two CART models for classifying each of the six lake water characteristics included in this study. The HGM and HGM+ CART models were parameterized for each lake water characteristic independently, and thus, we could determine the success of a classification developed from one variable for another. For example, an HGM model and a HGM+ model were built specifically for alkalinity. Classification success of these two alkalinityspecific models was then assessed for of each of the other lake water characteristics (e.g., water color, Secchi, TN, TP, and Chl a). However, because individual lakes are assigned to a category within the regionalization and landscape position classification systems, without respect to individual lake water characteristics, comparison across lake water characteristics was not appropriate. Two model selection statistics were used to compare among the candidate classifications for each lake water characteristic. First, we took an information-theoretic approach for multi- 10 model comparison, using the corrected Akaike information criteria (AICC) for small sample sizes (Burnham and Anderson 2002), computed in SAS (version 9.1) using PROC MIXED (SAS Institute Inc.). We compared the relative support for each classification using Akaike weights (ωi). These weights sum to equal 1 and are interpreted as the probability that a model is the best model relative to others included in the analysis (Johnson and Omland 2004). Secondly, we used the intra-class correlation coefficient (ICC) to compare the ability of each of the classifications to maximize class homogeneity for each dependent variable (Donner and Koval 1980, Cheruvelil et al. 2008). We calculated the ICC from the error terms of a one-way ANOVA with random effects: Yij = γ00 + rij + u0j where, Yij = observation of dependent variable for lake i in lake group j, γ00 = grand mean of the dependent variable, rij = random error term for lake i in lake group j, where rij ~ N (0, σ 2 ) and σ 2 represents the within-group error in the dependent variable, u 0 j = random error term for lake group j, where u0 j ~ N (0,τ 00 ) and τ 00 represents the among-group error in the dependent variable. The ICC is the amount of the total variance that is among groups: ∧ ∧ ∧ 2 ICC = τ 00/ ( τ 00 + σ ) A successful classification has a high ICC, meaning that a large amount of the variation is among the groups created from the classification, maximizing class homogeneity. All variables used in linear techniques (i.e., ANOVA) were transformed to meet normality assumptions. RESULTS 11 Comparing classifications Classification success, indicated by class homogeneity as represented by the ICC, ranged from 0% to 66% across all classification approaches and all lake water characteristics (Table 1.2). Across lake water characteristics, alkalinity was classified most successfully (ICC mean 42%, range 14% - 66%), followed by water color (ICC mean 21%, range 4% - 54%). Secchi disk depth and measures of lake productivity were classified least successfully (ICC mean <15%) with 6 classification failures (ICC = 0%). For each lake water characteristic evaluated, we observed that one model received an Akaike weight greater than 0.9, and all other models received very low weights, less than 0.1 (Table 1.2). Thus, only one classification was supported by the data for each lake water characteristic, with supported models differing among lake water characteristics. No regionalization or metric of landscape position alone was supported as a suitable classification of our data, as AICC values were substantially higher than most CART models. Among all classification approaches and all lake water characteristics, CART models had the highest strength of evidence (ωi range 0.92 - 1.00) and were the most successful at maximizing within class homogeneity (ICC range 36% - 66%). HGM+ models had more AICC support for a majority of lake water characteristics. Only two lake water characteristics (water color and Chl a) had a higher weight of evidence for HGM models. The HGM+ model for Chl a did not differ from the HGM model and, therefore, was not included in comparisons of model fit (indicated by “n/a” in Table 1.2). Class homogeneity (ICC) was not always maximized by the most parsimonious model, as indicated by AICC (Table 1.2). According to weight of evidence, the best classification for water 12 color, Secchi, and TP had ICC’s 12%, 4%, and 8%, lower, respectively, than the maximum ICC for that lake water characteristic. Relationships between hydrogeomorphic features and lake classes HGM models divided the study lakes into between 2 and 4 lake classes, capturing between 16% and 53% of the variation among lakes (Figure 1.2). Measures of lake morphometry were the most frequent classifiers across HGM models (4 of 6 models). Mean depth, in particular, was the most important feature driving HGM models of lake productivity, with water residence time and maximum depth also included in some models. Various measures of catchment morphometry were important in classifying alkalinity and water color. The proportion of the local catchment in wetlands was the most important classifier for water color and Chl a. Geology and climate variables were present in only one model each (Chl a and alkalinity, respectively). HGM+ models divided the study lakes into between 3 and 6 lake classes, capturing between 30% and 60% of the variation among lakes (Figure 1.3). All HGM+ models (except Chl a) explained more variation than HGM models: alkalinity by 10%, water color by 12%, Secchi by 11 %, TN by 29%, and TP by 8%. All HGM+ models (except Chl a) included the regionalization framework HUC as an important classifier (Figures 1.2 and 1.3). No landscape position metrics were included as important classifiers in any of these models. As with HGM models, measures of lake morphometry were frequently important classifiers across HGM+ models (4 of 6 models), followed by catchment morphometry and wetlands (2 of 6 models each). The proportion of clastic bedrock geology type was present in the HGM+ model for TN. Climate was not included as an important splitting variable in any of the HGM+ models. 13 The tree structure for most lake water characteristics was similar when comparing among HGM models and HGM+ models (Figures 1.2 and 1.3). Most notably, the HGM+ model of Chl a did not include any regionalizations and was, therefore, identical to the HGM model. HGM and HGM+ models for Secchi and TP shared initial structure and classification variables, differing only by the addition of HUC as the last classifier. Alternatively, the HGM features driving the two CART classifications of water color were from the same broad categories; however, different variables represented these categories. More specifically, for HGM models of water color, catchment morphometry was represented by both the ratio of cumulative catchment area to local catchment area (CUCA:LOCA) and the ratio of cumulative catchment area to lake area (CA:LK), whereas in HGM+ models, catchment morphometry was represented by one variable (cumulative catchment shape, CUCA shape). Despite these similarities, there were also some striking differences between HGM models and HGM+ models. For example, all HGM features included in the HGM model for alkalinity were completely replaced by regionalizations in the HGM+ model, with the 1st split (HUC) explaining 50% of the variation (Figures 1.2 and 1.3). In another example, the tree structure for the HGM+ model for TN is quite different than the HGM model, although lake morphometry is still important for classifying TN in both models. Catchment morphometry and bedrock geology are additional HGM features included in the HGM+ model of TN, increasing the number of classes created from 2 in the HGM model to 6 in the HGM+ model. Evaluation of CART tree stability Evaluating competitor and surrogate splitting options available in the detailed output from CART analyses can give a sense of the stability of a classification model. Some 14 classification splits can be labeled as unstable (multiple competitors with none being surrogates, further detail in Methods), thereby indicating a sensitivity of the resulting classification to the particular study lakes used to build the classification. For example, HUC is the most important classifier in the HGM+ model for TN. However, two other variables explain approximately 3% less variation than HUC and do not serve as surrogates (data not shown). Using our methods, this split can be labeled “unstable” and highly dependant on the input data. In contrast, we label the 2nd split in the TN HGM+ model as “somewhat stable” because mean depth explained only slightly less variation (approx. 2%) than maximum depth, the primary splitter at this node. In this case, however, mean depth acts as a surrogate for maximum depth since the majority (94%) of lakes would follow the same splitting path under either scenario. Therefore, this split is likely to be less dependent upon the specific dataset used in the analysis and can be considered stable. Over all 32 splits created in the CART models, 28% of splits were stable, 31% were somewhat stable/unstable, and 41% were unstable. A single classification for all lake water characteristics Classification success of CART models was compared across lake water characteristics to determine the model most successful at classifying lakes for all lake water characteristics. The mean ICC across water characteristics for the Secchi-HGM model ranked highest at 38% (Table 1.3). Class homogeneity varied slightly for most lake water characteristics when classified by the Secchi-HGM model in comparison to the variable-specific CART models. Homogeneity increased for water color (1%), Secchi (4%), and Chl a (2%), and decreased only a moderate amount for TN (6%) and TP (8%). However, the class homogeneity decreased much more for 15 alkalinity (50%). Simple correlations indicate that Secchi is significantly correlated to all lake water characteristics, except alkalinity (Table 1.4). Although most lake water characteristics were significantly correlated with one another (Table 1.4), classification success varied widely when a single classification was used for all lake water characteristics (Table 1.3). For example, alkalinity was significantly correlated with water color, TP, and Chl a (Table 1.4). However, neither of the alkalinity CART models successfully classified any other lake water characteristic (mean ICC excluding alkalinity: HGM 6%, HGM+ 5%; Table 1.3). In another example, water color had the highest correlation with all lake water characteristics (except TP, Table 1.4), yet water color CART models ranked low for overall classification success (HGM rank 8, HGM+ rank 5; Table 1.3). Therefore, correlations among lake water characteristics did not predict classification success across lake water characteristics. DISCUSSION There are three main conclusions that follow from our research questions. First, although some lake water characteristics were classified most successfully by local HGM features alone, most lake water characteristics were best classified when models included both lake-specific information and one or more regionalization. Second, lake and catchment morphometry plays a dominant role in structuring the classifications for the lake water characteristics. Third, using a single classification severely erodes the classification success for most water characteristics. Overall, because it is important for management agencies to balance the logistics and the effectiveness of classification, we suggest that the most successful classification system is (1) specific to one response variable, and (2) capable of incorporating information at multiple spatial scales and from a variety of different sources (regionalization and local HGM variables). We 16 found that CART models effectively modeled the complex interrelationships among our explanatory variables and are thus a useful tool for the classification of ecosystems, particularly those with multiple data types and/or non-linear relationships Comparing classifications Our results agree with previous studies showing that most lake water characteristics are not similar within many areas delineated by regionalizations (Jenerette et al. 2002, Cheruvelil et al. 2008). Regionalizations alone had poor classification success for all lake water characteristics, except alkalinity, which was classified well by HUCs and EDUs. Moreover, our alkalinity HGM model split on some characteristics that are more indicative of regional-scaled processes, such as climate. In another study of Michigan lakes, the authors found strong relationships between landscape position and alkalinity (Martin and Soranno 2006). However, our results show that this relationship is comparatively weak in contrast to larger-scale regionalizations. For example, landscape position (specifically, LH), explained 21% less variation than HUC in the 1st split of the alkalinity HGM+ model. These results may indicate that as the spatial scale of the classifying feature grows (LH50%) is needed before detecting significant relationships with water chemistry characteristics (Prepas et al. 2001). However, other studies report much lower thresholds (6-25%) beyond which wetland presence becomes important (Dillon et al. 1991, D’arcy and Carignan 1997, Canham et al. 2004). We also compared wetland cover measured over four spatial scales (cumulative catchment, local catchment, 500m buffer, and 100m buffer) and found that the proportion of wetlands in the local lake catchment was the only scale represented in any of the final CART classifications (water color and Chl a). However, analysis tree stability shows that other spatial scales act as competitor and/or surrogate splits for all lake water characteristics (except alkalinity). In some cases, there were only small losses in explanatory power when choosing other spatial scales. For example, cumulative catchment wetland cover explained 2% and 0.5% 20 less than local catchment wetland cover in the HGM models for water color and Chl a, respectively. Moreover, cumulative catchment wetland cover had 96% and 98% classification similarity with local catchment wetland cover. Although our results show that wetland presence at the local catchment scale is the strongest classifier for our study lakes, our results support other studies finding little difference in explanatory power between wetland cover measured at different spatial scales (Gergel et al. 1999, Strayer et al. 2003, Canham et al. 2004). Therefore, additional investigations are needed to more fully understand the scale and magnitude of wetland presence important for lake ecosystem dynamics. A single classification for all lake water characteristics The ultimate lake classification would successfully classify lakes for all water characteristics. Our results show that such a classification likely does not exist as this approach severely erodes the success of lake classification for most water characteristics. We found that the most successful single classification (Secchi HGM CART) for the lake water characteristics that we analyzed was on average 20% less successful in classifying other water characteristics and as much as 50% less successful in classifying alkalinity. However, when compared to regionalizations, the Secchi HGM CART was on average 16% more successful in classifying characteristics and 29% better when alkalinity was not included in the analysis. Thus, our results demonstrate that no single classification scheme maximizes success for all lake water characteristics because each classification depends on a different suite of local and regional HGM variables. However, comparisons such as ours should help guide the application of different approaches to lake classification and allow for management agencies to make choices between logistical practicality and ecological robustness. 21 Applications to ecosystem management Our use of detailed CART ouput allows for the evaluation of tree stability and splitting decisions. Alternative splitting decisions can be compared to allow for an assessment of practicality to enter into the classification process. For example, our results show that susceptibility of lakes to acidification may be adequately captured at larger spatial scales, as lakes within regions had similar alkalinity. On the other hand, management of nitrogen inputs to lakes should benefit from a more complex classification combining regionalizations and local HGM features. This increase in classification success, however, comes at the cost of increasing the number of lake classes. Management agencies can use this information to evaluate the tradeoffs involved in choosing different models. Most importantly, our approach to lake classification combines the strengths of a regionalization approach and a local HGM approach with analytical advances in multivariate statistics. Our approach can fulfill the needs of management agencies for an ecologically-based classification system which will allow for robust trend detection through time by reducing variation in natural HGM features within each class. Therefore, any resulting trend may be attributed to other factors such as changes in land use/cover or resource use (e.g., public access). The classification of other ecosystem types should also benefit from taking a multi-scale HGM approach by building upon foundational relationships between ecosystem function and hydrogeomorphic setting which can be measured over different spatial scales. 22 Table 1.1. Summary of lake water chemistry/clarity and hydrogeomorphic characteristics for 151 minimally disturbed lakes in Michigan USA. Variables are arranged into broad categories. Abbreviations are listed in parentheses. PCU, platinum cobalt units. CA:LK, catchment area to lake area ratio. CUCA:LOCA, cumulative catchment area to local catchment area ratio. units Min. Max. Mean SD mg/L PCU m µg/L µg/L µg/L 1 1 0.8 88 1.0 0.2 206 99 7.9 1044 32.0 29.0 72 14 3.4 461 11.7 4.3 55 16 1.3 192 6.0 4.4 % % % % 0.0 0.0 0.0 0.0 100.0 100.0 100.0 100.0 12.5 48.3 20.5 18.7 32.3 48.4 38.5 37.1 % % % % % % % 0.0 0.0 0.0 0.0 0.0 0.0 0.0 100.0 37.5 100.0 100.0 100.0 100.0 36.4 3.2 0.9 17.6 9.8 27.7 36.0 0.6 16.9 4.9 32.8 27.4 41.2 42.5 3.9 km2 unitless m m year 0.20 1.1 1.2 3.0 1.2* 70.38 6.3 21.8 58.5 31.6 2.84 1.9 4.9 14.4 2.5 9.34 0.7 3.3 9.2 4.2 Area km2 Shape unitless Slope % Cumulative catchment morphometry (CUCA) 0.2 1.1 0.6 1759.3 3.0 4.9 52.0 1.7 2.4 182.0 0.3 1.0 km2 unitless % 0.2 0.0 0.6 1948.3 2.6 5.7 86.0 0.6 2.8 276.7 0.4 1.1 Water chemistry/clarity Alkalinity Water color Secchi Total nitrogen (TN) Total phosphorus (TP) Chlorophyll a (Chl a) Bedrock geology Carbonate (B-Carb) Clastic (B-Clast) Hardrock (B-Hard) Iron (B-Iron) Surficial geology Bedrock (S-Bed) Dune (S-Dune) Glacial till (S-Till) Lacustrine (S-Lacu) Moraine (S-Mora) Outwash (S-Outw) Peat and muck (S-PeMu) Lake morphometry Lake Area (LK) Shape Mean depth Maximum depth (Max. Depth) Water residence time (WRT) Local catchment morphometry (LOCA) Area Shape Slope 23 Table 1.1 (cont’). CUCA:LK CUCA:LOCA Climate/hydrology Precipitation Runoff Baseflow index (BFI) Wetlands CUCA LOCA 500m buffer 100m buffer ratio ratio 0.9 1.0 2655.1 14.4 63.5 1.6 247.1 2.2 cm/year cm/year % 72.7 20.3 55 90.8 50.8 89 81.7 35.8 71 4.1 5.5 9 % % % % 0 0 0 0 25 23 40 53 5 4 7 8 5 4 7 11 24 Table 1.2. Summary of model selection statistics for candidate classifications per lake water characteristic. Classification type is indicated as regionalization, landscape position, or CART. Individual classification names are indicated within each classification type. Intra-class correlation coefficients (ICC) are presented as percent of total variance that is among the classes. A high ICC indicates high within class homogeneity and thus, high classification success. The number of classes per classification model (K) is presented. ∆AICC is the difference between the AICC for each model and the minimum AICC for each lake water characteristic. The ∆AICC will equal 0 for the best model per lake water characteristic. The Akaike weights (ωi) sum to 1 for each lake water characteristic and is interpreted as the likelihood that a given model is the best model relative to others included in the analysis. n/a, not applicable. Lake water characteristic Alkalinity Type Regionalization Landscape Pos. CART Water color Regionalization Landscape Pos. CART Secchi Regionalization Landscape Pos. CART TN Regionalization Landscape Pos. Name HUC EDU HLR LH LNN LNC HGM HGM+ HUC EDU HLR LH LNN LNC HGM HGM+ HUC EDU HLR LH LNN LNC HGM HGM+ HUC EDU HLR LH LNN ICC 63 50 14 34 32 22 56 66 17 14 15 9 4 9 46 54 0 4 0 5 2 10 54 50 13 5 0 2 0 25 K 19 6 5 7 5 4 4 3 19 6 5 7 5 4 4 4 19 6 5 7 5 4 3 5 19 6 5 7 5 AICC 1545 1586 1626 1598 1609 1617 1542 1515 404 396 403 402 409 401 361 366 522 523 524 521 524 520 474 454 182 187 185 187 185 ∆AICC 30 71 111 83 94 102 27 0 43 35 42 41 48 40 0 5 68 69 70 67 70 66 20 0 57 62 60 62 60 ωi 0.00 0.00 0.00 0.00 0.00 0.00 0.00 1.00 0.00 0.00 0.00 0.00 0.00 0.00 0.92 0.08 0.00 0.00 0.00 0.00 0.00 0.00 0.00 1.00 0.00 0.00 0.00 0.00 0.00 Table 1.2 (cont’). CART TP Regionalization Landscape Pos. CART Chlorophyll a Regionalization Landscape Pos. CART LNC HGM HGM+ HUC EDU HLR LH LNN LNC HGM HGM+ HUC EDU HLR LH LNN LNC HGM HGM+ 1 21 47 13 3 2 7 0 8 47 39 13 9 7 1 0 1 36 n/a 26 4 2 6 19 6 5 7 5 4 2 3 19 6 5 7 5 4 4 n/a 187 172 125 291 294 295 291 293 291 250 245 382 378 382 386 384 386 355 n/a 62 47 0 46 49 50 46 48 46 5 0 27 23 27 31 29 31 0 n/a 0.00 0.00 1.00 0.00 0.00 0.00 0.00 0.00 0.00 0.08 0.92 0.00 0.00 0.00 0.00 0.00 0.00 1.00 n/a Table 1.3. Summary of the intraclass correlation coefficients (ICCs) for CART models across lake water characteristics. CART models are listed by the original dependant variable and type of classification. Mean ICC across lake water characteristics is computed. Rank of mean ICC is listed. See Table 2 for acronyms. CART model Alkalinity Alkalinity-HGM 56 Alkalinity-HGM+ 66 Water color-HGM 1 Water color-HGM+ 22 Secchi-HGM 16 Secchi-HGM+ 25 TN-HGM 5 TN-HGM+ 22 TP-HGM 13 TP-HGM+ 17 Chl a-HGM 6 Lake water characteristic Water color Secchi TN 7 5 12 9 0 5 46 22 12 54 50 7 47 54 41 39 50 31 14 19 21 38 36 47 31 29 24 28 23 18 31 19 15 27 Mean TP 5 5 12 9 31 29 26 40 47 39 20 Chl a 1 8 11 11 38 27 12 12 23 24 36 Rank 14 16 17 26 38 34 16 33 28 25 21 11 10 8 5 1 2 9 3 4 6 7 Table 1.4. Pearson product-moment correlation matrix for lake water characteristics. *** p<0.001, ** p<0.01, * p<0.05, NS, not significant (p>0.05). Lake water characteristic Alkalinity Secchi Water color TN TP Chl a Alkalinity … 0.10NS -0.22** -0.02NS -0.17* -0.21** Lake water characteristic Secchi Water color TN … -0.66*** -0.45*** -0.50*** -0.46*** … 0.48*** 0.53*** 0.52*** 28 … 0.59*** 0.38*** TP … 0.39*** Chl a … Figure 1.1. Map of the upper and lower peninsula of Michigan. Lakes included in the analysis are shown as solid dots. Boundaries are shown for each regionalization: A) 8-digit USGS hydrologic units (HUC), B) ecological drainage unit (EDU), and C) hydrological landscape region (HLR). A. HUC B. EDU C. HLR 29 Figure1.2. Results from the CART analysis of local hydrogeomorphic features (HGM) in classifying lake water characteristics: (A) alkalinity, (B) water color, (C) Secchi, (D) TN, (E) TP, and (F) Chl a. Each split is labeled with the splitting variable (see Table 1 for abbreviations), and proportional reduction in error (PRE). Branches are labeled with splitting value. Terminal nodes (rectangles) represent lake classes and are labeled with an alphabetical class, class mean, and number of lakes per class (in parentheses). A) Alkalinity (overall PRE 0.53) B) Water color (overall PRE 0.38) CUCA Area (0.36) LOCA wetlands (0.17) < 28.6 km2 > 28.6 km2 < 4.4% > 4.4% Runoff (0.11) A: 127 (39) A: 8.5 (89) CUCA:LOCA ( 0.16) > 14 cm/yr LOCA Area (0.06) < 1.0 < 14 cm/yr B: 86 (32) < 4.6 km2 CUCA:LK ratio (0.05) C) Secchi (overall PRE 0.33) Mean depth (0.24) < 1.1 yr D: 21.4 (26) D) TN (overall PRE 0.16) Mean depth (0.16) > 7.8 m > 3.6 m < 3.6 m A: 392 (83) B: 546 (68) A: 5.2 (18) > 1.1 yr B: 2.8 (70) > 8.1 C: 8.6 (20) C: 23 (42) D: 58 (38) WRT(0.09) B: 37.9 (16) < 8.1 > 4.6 km2 < 7.8 m > 1.0 C: 3.6 (63) 30 Figure 1.2 (cont’). E) TP (overall PRE 0.22) Mean depth (0.22) F) Chl a (overall PRE 0.30) LOCA wetlands (0.11) < 3.2% > 5.1 m A: 2.8 (74) < 5.1 m > 3.2% Max. depth (0.10) > 15.4 m < 15.4 m B: 2.6 (22) S-Till (0.09) < 41% A: 7.7 (50) C: 6.0 (46) B: 13.7 (101) 31 > 41% D:12.0 (9) Figure 1.3. Results from the CART analysis combining regionalization, landscape position and local hydrogeomorphic features (HGM+) in classifying lake water characteristics: (A) alkalinity, (B) water color, (C) Secchi, (D) TN, (E) TP, and (F) Chl a. Each split is labeled with the splitting variable (see Table 1 for abbreviations), and proportional reduction in error (PRE). Branches are labeled with splitting value. Branches split using categorical variables (e.g., HUC, EDU) are detailed in appendices. Terminal nodes (rectangles) represent lake classes and are labeled with an alphabetical class, class mean, and number of lake per class (in parentheses). A) Alkalinity (overall PRE 0.59) B) Water color (overall PRE 0.49) HUC (0.50) LOCA Wetlands (0.17) < 4.4% > 4.4% A: 8.5 (89) HUC (0.16) A: 133 (43) EDU (0.09) B: 14.0 (39) CUCA Shape (0.16) > 0.384 B: 35 (77) C: 22.4 (15) C: 78 (34) C) Secchi (overall PRE 0.45) < 0.384 D: 56.6 (8) D) TN (overall PRE 0.46) HUC (0.17) Mean Depth (0.24) < 7.8 m > 7.8 m Max. Depth (0.09) Max. Depth (0.08) WRT (0.09) A: 5.2 (18) > 16.6 m < A: 298 (29) > 10.4 m < < 1.1 yr > 1.1 yr HUC (0.07) HUC (0.05) C: 3.6 (20) C: 681 (23) CUCA:LOCA (0.06) B: 481 (23) < 1.0 > 1.0 D: 413 (56) B-Clast (0.06) E: 4.1 (27) > 0.4% B: 2.4 (50) E: 405 (8) D: 3.2 (36) 32 < 0.4% F: 657 (12) Figure 1.3 (cont’). E) TP (overall PRE 0.30) Mean depth (0.22) F) Chl a (overall PRE 0.30) LOCA wetlands (0.11) < 3.2% < 5.1 m > 3.2% A: 2.8 (74) > 5.1 m Max. depth (0.10) > 15.4 m A: 7.7 (50) HUC (0.08) < 15.4 m B: 2.6 (22) S-Till (0.09) < 41% B: 11.0 (38) C: 6.0 (46) C: 15.3 (63) 33 > 41% D:12.0 (9) APPENDIX 34 Appendix 1. Additional figures for Chapter 1 Figure 1.4. Details of HUC and EDU categorical splits in the alkalinity HGM+ model (refer to Figure 3): A) first node splitting further into class A and towards classes B and C, and B) second node splitting further to classes B or C. Grey areas indicate HUCs which split to the left. Striped areas indicate HUCs which split to the right. A. B. 35 Figure 1.5. Details of HUC categorical split in the water color HGM+ model (refer to Figure 3). Black areas indicate HUCs which were not represented at the node. Grey areas indicate HUCs which split to the left into class B. Striped areas indicate HUCs which split to the right towards classes C and D. 36 Figure 1.6. Details of HUC categorical splits in the Secchi HGM+ model (refer to Figure 3): A) left split from WRT splitting further to classes B and C, and B) right split from WRT splitting further to classes D and E. Black areas indicate HUCs which were not represented at the node. Grey areas indicate HUCs which split to the left into classes B or D. Striped areas indicate HUCs which split to the right into classes C or E. A. B. 37 Figure 1.7. Details of HUC categorical split in the TN HGM+ model (refer to Figure 3). Grey areas indicate HUCs which split to the left towards classes A, D, E, and F. Striped areas indicate HUCs which split to the right towards classes B and C. 38 Figure 1.8. Details of HUC categorical split in the TP HGM+ model (refer to Figure 3). Grey areas indicate HUCs which split to the left into class B. Striped areas indicate HUCs which split to the right into class C. 39 REFERENCES 40 REFERENCES Albert, D.A., 1995. Regional Landscape Ecosystems of Michigan, Minnesota, and Wisconsin: A working Map and Classification, USDA Forest Service North Central Forest Experiment Station General Technical Report NC-178. Bailey, R.G., P.E. Avers, T. King, and W.H. McNab. 1994. Ecoregions and subregions of the United States (map) (supplementary table of map unit descriptions compiled and edited by McNab, W.H. and Bailey, R.G.): Washington, D.C., U.S. Department of Agriculture Forest Service, scale 1:7,500,000. Breiman, L., J.H. Friedman, R.A. Olshen, and C.J. Stone. 1984. Classification and regression trees. New York: Chapman and Hall. Brett, M.T., and M.M. Benjamin. 2008. A review and reassessment of lake phosphorus retention and the nutrient loading concept. Freshwater Biology 53:194-211. Brinson, M.M. 1993. A hydrogeomorphic classification for wetlands. Technical Report WRPDE- 4, U.S. Army Corps of Engineers, Waterways Experiment Station, Wetlands Research Program, Washington, DC, USA. Bryan, B.A. 2006. Synergistic techniques for better understanding and classifying the environmental structure of landscapes. Environmental Management 37(1): 126–140. Burnham, K.P., and D.R. Anderson, 2002. Model Selection and Multimodel Inference: A Practical-Theoretic Approach, 2nd ed. Springer-Verlag. ISBN 0-387-95364-7. Canham C.D., M.L. Pace, M.J. Papaik, A.G.B. Primack, K.M. Roy, R.J. Maranger, R.P. Curran, and D.M. Spada. 2004. A spatially explicit watershed-scale analysis of dissolved organic carbon in Adirondack lakes. Ecological Applications 14(3): 839–854. Cheruvelil, K.S., P.A. Soranno, M.T. Bremigan, T. Wagner, and S.L. Martin. 2008. Grouping lakes for water quality assessment and monitoring: the roles of regionalization and spatial scale. Environmental Management 41: 425-440. D’Arcy, P., and R. Carignan. 1997. Influence of catchment topography on water chemistry in southeastern Québec Shield lakes. Canadian Journal of Fisheries and Aquatic Sciences 54: 2215–2227. De’ath, G., and K.E. Fabricius. 2000. Classification and regression trees: a powerful yet simple technique for ecological data analysis. Ecology 81: 3178-3192. Detenbeck, N.E., C.A. Johnston, and G.J. Niemi. 1993. Wetland effects on lake water quality in the Minneapolis/St. Paul metropolitan area. Landscape Ecology 8(1): 39-61. 41 Dillon, P.J., L.A Molot, and W.A Scheider. 1991. Phosphorus and nitrogen export from forested stream catchments in Central Ontario. Journal of Environmental Quality 20: 857–864. Donner, A., and J.J. Koval. 1980. The estimation of intraclass correlation in the analysis of family data. Biometrics 36: 19-25. Emmons, E.E, M.J. Jennings, and C. Edwards. 1999. An alternative classification method for northern Wisconsin lakes. Canadian Journal of Fisheries and Aquatic Sciences 56(4): 661-669. Fee, E. J. 1979. A relation between lake morphometry and primary productivity and its use in interpreting whole-lake eutrophication experiments. Limnology and Oceanography 24: 401–416. Frissell, C.A., W.J. Liss, C.E. Warreb, and M.D. Hurley. 1986. A hierarchical framework for stream habitat classification: Viewing streams in a watershed context. Environmental Management 10(2):199-214. Gergel, S.E., M.G. Turner, and T.K. Kratz. 1999. Dissolved organic carbon as an indicator of the scale of watershed influence on lakes and rivers. Ecological Applications 9: 1377–1390. Goransson, E., R.K. Johnson, and A.Wilander. 2004. Representativity of a mid-lake surface water chemistry sample. Environmental Monitoring and Assessment 95: 221–238. Griffith, G.E., A.J. Kinney and J.M. Omernik. 1987. Interpreting patterns of lake alkalinity in the Upper Midwest Region USA. Lake and Reservoir Management 3(1): 329 - 336. Hakanson, L. 1996. Predicting important lake habitat variables from maps using modern modelling tools. Canadian Journal of Fisheries and Aquatic Sciences 53(Suppl. 1): 364– 382. Hakanson, L. 2005. The importance of lake morphometry and catchment characteristics in limnology – ranking based on statistical analyses. Hydrobiologia 541: 117-137. Halsey, L.A., D.H. Vitt, and D.O. Trew. 1997. Influence of peatlands on the acidity of lakes in northeastern Alberta, Canada. Water, Air, and Soil Pollution 96: 17-38. Hershey, A.E., G.M. Gettel, M.E. McDonald, M.C. Miller, H. Mooers, W.J. O'Brien, J. Pastor, C. Richards, and J.A. Schuldt. 1999. A geomorphic-trophic model for landscape control of arctic lake food webs. BioScience 49(11): 887-897. Higgins, J.V., M. Lammert, M. Bryer, M. DePhilip, and D. Grossman. 1998. Freshwater conservation in the Great Lakes basin: Development and application of an aquatic community classification framework. Report to The George Gund Foundation, Cleveland, Ohio (Grant Number 95-65). 42 Higgins, J.V., M.T. Bryer, M.L. Khoury, and T.W. Fitzhugh. 2005. A freshwater classification approach for biodiversity conservation planning. Conservation Biology 19: 432-445. Host G.E., P.L. Polzer, D.J. Mladenoff, M.A. White, and T.R. Crow. 1996. A quantitative approach to developing regional ecosystem classifications. Ecological Applications 6(2): 608-618. Hunsaker, C.T., J.L. Malanchuk, R.J. Olson, S.W. Christensen, and R.S. Turner. 1986. Adirondack headwater lake chemistry relationships with watershed characteristics. Water, Air, and Soil Pollution 31(1-2): 79-88. Jenerette, G.D., J. Lee, D.W. Waller, and R.E. Carlson. 2002. Multivariate analysis of the ecoregion delineation for aquatic systems. Environmental Management 29(1): 67–75. Johnson, L.B. and S.H. Gage. 1997. Landscape approaches to the analysis of aquatic ecosystems. Freshwater Biology 37: 113-132. Johnson, J.B., and K. S. Omland. 2004. Model selection in ecology and evolution. Trends in Ecology and Evolution 19(2): 101-108. Kratz, T.K., K.E. Webster, C.J. Bowser, J.J. Magnuson, and B.J. Benson. 1997. The influence of landscape position on lakes in northern Wisconsin. Freshwater Biology 37: 209–217. Kunst C., E. Monti, H. Perez, and J. Godoy. 2005. Assessment of the rangelands of southwestern Santiago del Estero, Argentina, for grazing management and research. Journal of Environmental Management 80: 248–265. Martin, S.L., and P.A. Soranno. 2006. Lake landscape position: Relationships to hydrologic connectivity and landscape features. Limnology and Oceanography 51(2): 801-814. MIRIS (Michigan Resource Information System). 2000. [Internet] Michigan Department of Natural Resources [accessed 2001 May]. Available from http://www.ciesin.org/datasets/mirislcover/miriscov.html Momen, B., and J.P Zehr. 1998. Watershed classification by discriminant analyses of lakewaterchemistry and terrestrial characteristics. Ecological Applications 8(2): 497-507. Mora, F. and L. Iverson. 2002. A spatially constrained ecological classification: rationale, methodology and implementation. Pant Ecology 158:153-169. Moss B., P.Johnes, and G. Phillips. 1994. August Thienemann and Loch Lomond — an approach to the design of a system for monitoring the state of north-temperate standing waters. Hydrobiologia 290: 1-12. Newton, R.M., and C.T. Driscoll. 1990. Classification of ALSC lakes. Pages 2–70 to 2–91 in J. P. Baker, S.A. Gherini, S.W. Christensen, C.T. Driscoll, J. Gallagher, R.K. Munson, 43 R.M. Newton, K.H. Reckow, and C.L. Schofield, editors. Adirondack Lakes survey: an interpretive analysis of fish communities and water chemistry, 1984–87. Adirondack Lakes Survey Corporation, Ray Brook, New York, USA. Olden, J.D., and D.A. Jackson. 2002. A comparison of statistical approaches for modeling fish species distributions. Freshwater Biology 47: 1976-1995. Omernik, J.M., and A.J. Kinney. 1983. An improved technique for estimating mean depth of lakes. Water Research. 17: 1603–1607. Omernik, J.M. 1987. Ecoregions of the conterminous United States. Annals of the Association of American Geographers 77, 118-125. Pitblado, J.R., W. Keller, and N.I. Conroy. 1980. A classification and description of some northeastern Ontario lakes influenced by acid precipitation. Journal of Great Lakes Research 6(3): 247-257. Prepas E.E., D. Planas, J.J. Gibson, D.H. Vitt, T.D. Prowse, W.P. Dinsmore, L.A. Halsey, P.M. McEachern, S. Paquet, G.J. Scrimgeour, W.M. Tonn, C.A. Paszkowski, and K. Wolfstein. 2001. Landscape variables influencing nutrients and phytoplankton communities in Boreal Plain lakes of northern Alberta: A comparison of wetland- and upland-dominated catchments. Canadian Journal of Fisheries and Aquatic Sciences 58: 1286–1299. Pyne M.I., R.B. Rader, and W.F. Christensen. 2007. Predicting local biological characteristics in streams: a comparison of landscape classifications. Freshwater Biology 52: 1302–1321. Quinlan R., A.M. Paterson, R.I. Hall, P.J. Dillon, A.N. Wilkinson, B.F. Cumming, M.S.V. Douglas, and J.P. Smol. 2003. A landscape approach to examining spatial patterns of limnological variables and long-term environmental change in a southern Canadian lake district. Freshwater Biology 48: 1676–1697. R Development Core Team. 2008. R: a language and environment for statistical computing. R Foundation for statistical computing, Vienna, Austria. Rasmussen, J.B., L. Godbout, and M. Schallenberg. 1989. The humic content of lake water and its relationship to watershed and lake morphometry. Limnology and Oceanography 34: 1336–1343. Riera, J.L., J.J. Magnuson, T.K. Kratz, and K.E. Webster. 2000. A geomorphic template for the analysis of lake districts applied to the Northern Highland Lake District, Wisconsin, USA. Freshwater Biology 43: 301–318. Robertson, D.M., D.A. Saad, and D.M. Heisey. 2006. A regional classification scheme for estimating reference water quality in streams using land-use-adjusted spatial regressiontree analysis. Environmental Management 37(2): 209–229. 44 Rodgers, K. 2005. Evaluation of nearshore coral reef condition and identification of indicators in the main Hawaiian islands. PhD Dissertation University of Hawai‘i. Seaber, P.R., F.P. Kapinos, and G.L. Knapp. 1987. Hydrologic Unit Map, USGS Water-Supply Paper 2294. Seelbach, P.W., M. Wiley, J.C. Kotanchik, and M.E. Baker. 1997. A landscape-based ecological classification system for river valley segments in Lower Michigan (MI-VSEC version 1.0). Michigan Department of Natural Resources Fisheries Division, Fisheries Research Report 2036. Soranno P.A, K.E. Webster, J.L. Riera, T.K. Kratz, J.S. Baron, P.A. Bukaveckas, G.W. Kling, D.S. White, N.Caine, R.C. Lathrop, and P.R. Leavitt. 1999. Spatial variation among lakes within landscapes: Ecological organization along lake chains. Ecosystems 2: 395–410. Soranno, P.A, K.S. Cheruvelil, R.J. Stevenson, S.L. Rollins, S.W. Holden, S. Heaton, and E. Torng. 2008. A framework for developing ecosystem-specific nutrient criteria: integrating biological thresholds with predictive modeling. Limnology and Oceanography 53(2): 773-787. Stendera, S. and R.K. Johnson. 2006. Multiscale drivers of water chemistry of boreal lakes and streams. Environmental Management 38:760–770. Stoddard J.L., D.P. Larsen, C.P. Hawkins, R.K. Johnson, and R.H. Norris. 2006. Setting expectations for the ecological condition of streams: The concept of reference condition. Ecological Applications 16(4): 1267–1276. Strayer, D.L., R.E. Beighley, L.C. Thompson, S.Brooks, C.Nilsson, G.Pinay, and R.J. Naiman. 2003. Effects of land cover on stream ecosystems: roles of empirical models and scaling issues. Ecosystems 6: 407–423. Venables W.N., and B.D. Ripley. 1999. Modern applied statistics with S-PLUS. Third edition. New York: Springer. Vidon, P.G.F., and A.R. Hill 2004. Landscape controls on the hydrology of stream riparian zones. Journal of Hydrology 292:210-228. Vollenweider, R.A., 1968. The Scientific Basis of Lake Eutrophication, with Particular Reference to Phosphorus and Nitrogen as Eutrophication Factors. Technical Report DAS/DSI/68.27, OECD, Paris, 159 p. Webster, K.E., P.A. Soranno., K.S. Cheruvelil, M.T. Bremigan, J.A. Downing, P.D. Vaux, T.R. Asplund, L.C. Bacon, and J. Connor. 2008. An empirical evaluation of the nutrient-color paradigm for lakes. Limnology and Oceanography 53(3): 1137–1148. 45 Wetzel R.G. and G.E. Likens. 2000. Limnological analyses, 3rd ed. Springer-Verlag. Winter, T.C. 1977. Classification of the hydrologic settings of lakes in the north-central United States. Water Resources Research 13: 753-767. Winter, T.C. 2001. The concept of hydrologic landscapes. Journal of the American Water Resources Association 37: 335-349. Wolock, D.M., G.M. Hornberger, K.J. Beven, and W.G. Campbell. 1989. The relationship of catchment topography and soil hydraulic characteristics to lake alkalinity in the northeastern United States. Water Resources Research 25: 829–837. Wolock, D.M., T.C. Winter, and G. McMahon. 2004. Delineation and evaluation of hydrologiclandscape regions in the United States using geographic information system tools and multivariate statistical analyses. Environmental Management 34, S71-S88 Young, T.C., and J.L. Stoddard. 1996. The Temporally Integrated Monitoring of Ecosystems (TIME) project design 1. Classification of northeast lakes using a combination of geographic, hydrogeochemical, and multivariate techniques. Water Resources Research 32(8): 2517–2528. Zimmerman, A.P., K.M. Noble, M.A. Gates, and J.E. Paloheimo. 1983. Physicochemical typologies of south-central Ontario lakes. Canadian Journal of Fisheries and Aquatic Sciences 40(10): 1788-1803. 46 CHAPTER 2: THE LAND-USE LEGACY EFFECT - ADDING TEMPORAL CONTEXT TO UNDERSTANDING LAKE WATER CHEMISTRY 47 INTRODUCTION The conversion of land to serve human purposes can alter lake ecosystems. In particular, many human activities drive lake eutrophication through direct (e.g. sewage and industrial wastes) and indirect (e.g. agricultural run-off) inputs of excess nutrients (Dillon and Kirchner 1975, Omernik 1976, Holtan et al. 1988). Although many studies have related land use/cover to lake and stream characteristics, much variation remains to be explained (Allan 2004, Declerck et al. 2006). For example, agriculture is a known source of nitrogen and phosphorus to aquatic systems. However, large differences exist in nutrient export to streams among watersheds dominated by agriculture (Arbuckle and Downing 2001, Vanni et al 2001, Knoll et al. 2003) as well as forest (Findlay et al. 2001, Brett et al. 2005). In contrast, some studies report finding no relationship between riverine nutrient concentrations and watersheds with contrasting land cover (Pellerin et al. 2004, Burcher and Benfield 2006). Such inconsistencies indicate that the scientific understanding of ecosystem response to land use/cover (LULC) is still developing (Allan 2004). The recognition of legacy effects from historical LULC is an important conceptual advance that has clarified the relationship between LULC and ecosystem responses in terrestrial (Foster et al. 1998, Foster et al. 2003,Chauvat et al. 2007), and aquatic systems (Harding et al. 1998, McTammany et al. 2007). Legacy effects can be defined as effects which perpetuate beyond an expected or perceived endpoint in time. In other words, an ecosystem response to LULC disturbance depends on contemporary LULC cover as well as effects from prior LULC that are still propogating through the system. Harding et al. (1998) provide a prominent example of legacy effects in aquatic systems. They found that historic LULC had a stronger relationship with stream macroinvertebrate and fish community diversity than current LULC and called for further research investigating land use legacy effects. 48 The value of adding a temporal component to analyses of LULC impacts has been discussed in many studies (Allan 2004, Pellerin et al. 2004, Opperman et al. 2005, Goodman et al. 2006, Burcher and Benfield 2006, Davies and Jackson 2006, Foy and Lennox 2006). Specifically, Goodman et al. (2006) report no differences in stream macroinvertebrate community diversity between forested and clear-cut sites, which they attributed to legacy effects from previous agricultural land use. In another example, Foy and Lennox (2006) indicated that a land-use legacy effect may be the mechanism driving increasing riverine phosphorus export despite temporally stable LULC and reductions in point source phosphorus input. Legacies have been suggested to play roles in other small scale (Reed-Anderson et al. 2000, Tomer and Burkhart 2003, Schilling and Spooner 2006) as well as global nutrient budgets (Bennett et al. 2004). Despite the many studies that suggest the importance of legacy effects, few studies have been designed explicitly to investigate these temporal effects of LULC changes on aquatic systems. Moreover, studies that have addressed legacies have thus far focused on biological endpoints, such as measures of community diversity. Because aquatic biological communities integrate environmental stresses over their lifetime and respond indirectly to some aspects of LULC change, measuring the effects of land use legacies on water chemistry may provide a more direct mechanistic link to aquatic ecosystem condition. There are at least three potential mechanisms for land-use legacy effects in lakes: 1) import through surface water runoff, 2) delayed import through groundwater pathways, and 3) internal recycling from lake sediments. Knowledge of biogeochemical cycles can help differentiate among these mechanisms. Differences in biogeochemical activity and therefore cycling periods will likely produce different signals in the data for dissolved versus particulate 49 substances and for reactive versus conservative ions. Phosphorus tends to adsorb to soil, being delivered primarily as particulate matter though surface run-off (Wetzel 1983). In contrast, nitrogen is highly mobile, readily dissolved, and delivered through both surface and groundwater pathways (Wetzel 1983). Therefore, differences between nitrogen and phosphorus legacy signals may indicate differences between surface and groundwater delivery. Furthermore, because lake systems retain a majority of water column nutrients and export a majority of conservative ions (Larsen et al. 1981, Kling et al. 2000), comparing the legacy signals of conservative ions to those of nutrients can help differentiate between external delivery and internal recycling. The goal of our research was to investigate temporal LULC effects (i.e. legacy effects) on lake water chemistry to increase understanding of how changes in LULC over time influence lake ecosystem responses, such as eutrophication. We used multiple linear regression, where the response variable was one of a suite of water chemistry variables and the predictor variables were time-specific LULC represented by principal components. We then compared model fit and model predictive ability between models using LULC from a single time, and legacy models, which use LULC from multiple time periods. By comparing the legacy signals from a wide range of limnological variables, we can better understand the processes that differentiate the legacy signal among these chemicals. METHODS Study area This study was conducted within the Huron River Watershed (HRW) in Michigan (Figure 2.1). This watershed is approximately 2,359 km2 in size and ranges from 390m to 173m in elevation. The HRW contains numerous rivers, streams and lakes supporting a diverse 50 assemblage of aquatic species (Hay-Chmielewski et al. 1995). Due to close proximity to major metropolitan areas, such as Detroit and Ann Arbor, and the development of interstate highways, many people live, work, and enjoy the recreational opportunities the area provides. The watershed has undergone extensive LULC change over the past century, shifting from an agrarian to suburban society (Hay-Chmielewski et al. 1995, Rutledge and Lepczyk 2002). However, numerous natural areas have been maintained in various publicly owned parks. These characteristics make the HRW a prime area to investigate legacy effects as the watersheds cover a range of LULC trajectories. Description of data Water samples were collected from 35 lakes in the HRW during 2008 spring mixing. Lakes were chosen based on their accessibility, including both private and public lakes, and ranged in size from 0.05 to 2.6 km2. All water samples were taken from 1m below the surface at the deepest portion of the lake. Samples were analyzed for a range of limnological attributes. Cation (calcium, magnesium, potassium, and sodium) and anion (chloride, nitrate, and sulfate) concentrations were determined using membrane-suppression ion chromatography (Wetzel and Likens 2000). Silica concentrations were determined using the molybdate colorimetric method (Wetzel and Likens 2000). Total nitrogen concentrations were determined using the 2nd derivative of the absorbance curve at 224 nm following persulfate digestion. Ammonia concentrations were determined following the indophenol-blue method. Soluble reactive phosphorus (SRP) and total phosphorus (TP) concentrations were determined spectrophotometrically following persulfate digestion method. 51 Water chemistry variables were organized into a priori groups based largely on biological reactivity and solubility (Appendix 1). Phosphorus species were the most highly reactive variables in our study (Wetzel 1983). Nitrogen species are also highly reactive, but less so than phosphorus due to the relative abundance of nitrogen versus phosphorus (Wetzel 1983). We group Ca, Si, and SO4 together under the heading of “reactive variables” due to their role as micronutrients in the metabolisms of macrophytes, diatoms, and bacteria, respectively (Wetzel 1983). Conservative ions (i.e. Cl, K, Na, and Mg) play a much smaller role in biological metabolism and therefore are influenced more by external supply than internal conversions (Wetzel 1983). Land use/cover data (30m resolution) were available from five time steps (1938, 1955, 1968, 1978, and 1996), classified into six categories (urban, agriculture, open, forest, water, and wetland) based a modified version of the Anderson et al. (1976) LULC classification scheme. All data were compiled into a multi-temporal GIS database for analysis (full details in Rutledge 2001). Briefly, digitized land use/cover data for 1978 were provided by the Michigan Department of Natural Resources Michigan Resource Information System (MIRIS) and served as a base for digitizing all other time steps. Aerial photos from each of the other time steps were scanned (150 dpi) to create digital images. These images were then registered and rectified to the 1978 data using the placement of county roads. Land use/cover polygons were then digitized. This technique increased the consistency of polygon locations. Even though LULC data from more recent times were available for the study area, the source images and LULC classification techniques made the datasets incomparable to our historical data. The proportion of each land use/cover type during each time step was calculated for the surface drainage area (watershed) of each study lake. 52 Statistical analyses The relationship between time-specific land use/cover and lake water chemistry was analyzed using a combination of principal components analysis and multiple regression. Principal components (PC) analysis on the covariance matrix of the six LULC types was used to reduce the dimensionality of the data, reduce collinearity among the LULC variables, and summarize the LULC composition for lake watersheds for each time step. Land use/land cover composition for the watershed of each lake (n = 35) entered the PC analysis separately for each time step, totaling 175 observations entering into the PC analysis (35 lakes by 5 time steps). The six LULC types were combined via PC analysis to represent each lake watershed for each time period, preserving changes in watershed composition through each of the five time steps as well as differences among watersheds. This approach to PC analysis of LULC shows differences in watershed LULC between lakes and between times, highlighting different LULC trajectories. We built regression models to test the effects of land-use legacies for each water chemistry variable following a chronological hierarchy of LULC following the general form: yi = β0 + β1,t PC1t + β2,t PC2t + β3,t PC3t + εi where, yi is the value of water chemistry variable i, β1,t is the coefficient for PC axis 1 in time step t (1996 to 1938), and PC1t is the value of PC axis 1 at timestep t (1996 to 1938). Model building continued until all time steps were included in a single regression model (max 15 independent variables) as shown below: Model 1: yi = β0 + β1,1PC11996 + β2,1PC21996 + β3,1PC31996 + εi 53 Model 2: yi = β0 + β1,1PC11996 + β2,1PC21996 + β3,1PC31996 + β1,2PC11978 + β2,2PC21978 + β3,2PC31978 + εi … Model 5: yi = β0 + β1,1PC11996 + β2,1PC21996 + β3,1PC31996 + β1,2PC11978 + β2,2PC21978 + β3,2PC31978 + β1,3PC11968 + β2,3PC21968 + β3,3PC31968 + β1,4PC11955 + β2,4PC21955 + β3,4PC31955 + β1,5PC11938 + β2,5PC21938 + β3,5PC31938 + εi By building the regression models in this way, our approach recognizes the contribution from all time periods and makes the explicit assumption that current LULC is a stronger driver of lake water chemistry than older LULC. It was not our intent to interpret the effect of individual LULC types over time, as represented by individual regression coefficients, but rather we focused on the aggregate effect over time (see Appendix 2). We used Akaike Information 2 Criteria (AIC), AIC weights (ωi) and the coefficient of determination (R ) (following Burnham and Anderson 2002) to compare the performance among the five models for each water chemistry variable. We also used AIC weights to compare models built from two PC’ to models built from 3 PC’s (Appendix 3). All analyses were computed using 3 PC’s in SAS using PROC PRINCOM and PROC REG (SAS). RESULTS Changes in land use/cover 54 Averaged across all lake watersheds, the largest changes in LULC occurred in the agriculture and urban categories (Fig. 2.2). Agriculture predominated in 1938, comprising an average of 40% of the study watersheds. Over time, agriculture steadily declined in area to represent only about 15% of the watershed areas by 1996. As agriculture declined, urban LULC steadily increased from less than 5% in 1938 to over 20% by 1996. Wetland LULC type also declined steadily from 12% to 8% over the study period. Both forest and open LULC types increased initially, reaching 23% and 21% of the watershed areas in 1968 and 1978, respectively. However, these LULC types then decreased for the remainder of the study period, comprising 20% and 19%, respectively, of the watershed areas by 1996 (Fig. 2.2). The first three PC’s explained approximately 89% of the variation in the dataset (Table 2.1). The first PC represented almost half of the variation in the dataset and was weighted most heavily by agriculture and urban LULC types (Table 2.1). The pattern of urban expansion and agricultural contraction, seen using average watershed LULC (Fig 2.2), can also be seen using a multivariate representation of the LULC data (Fig. 2.3). Agriculture dominated the landscape in 1938, shifting steadily to urban dominance between 1955 and 1968 (Fig 2.3). The second PC was weighted most heavily on urban and forest (Table 2.1), and showed a slight shift between forest and urban LULC between 1938 and 1968 (Fig. 2.3), as the pace of urban growth exceeded that of forest (Fig. 2.2). A more dramatic shift in forest and urban LULC types occurred between 1968 and 1996, as the growth of urban increased as forest waned. The third PC was weighted most heavily on open and forest LULC types (Table 2.1). The pattern of LULC change represented with the third PC is more complex than what was seen in the first two PC’s. Open areas were slightly more prevalent than forest in 1938, switching to forest prevalence for 1955 and 1968. Forest and open areas were roughly equal for the last two time steps of the study. 55 Plotting the change in LULC of each lake watershed over time showed that our dataset included lakes that vary in initial conditions, rate of change, and final composition (Appendix 4). The individual trajectories of LULC change represented by the first PC were nearly parallel, showing a change from agriculture to urban dominance. However, some of the lakes changed more rapidly than others, some changed more than others, and some changed very little. The second PC showed a narrower range of initial conditions then the first PC, but also showed largely parallel trajectories for the majority of the lakes with few divergences. The trajectory of the third PC was more erratic than the other two PCs, representing more lake to lake variation in the timing of shifts in open and forest LULC. Comparing regression models The regression models using only current LULC provided the best model fit only for phosphorus species (Fig. 2.4). The explanatory power, as measured by R2, for these models ranged from 19% to 64%, but had low Akaike weights (mean 8%) for all other water chemistry variables. Our results showed that legacy models provided a better model fit for all other lake water characteristics. Including additional information about historical LULC in regression models always improved explanatory power (Fig. 2.5), but did not always improve model fit enough to account for the added complexity (Fig. 2.4). The explanatory power of the simplest model for TP and SRP (the 1996-only model) was only 25% and 22%, increasing to 45% and 35%, respectively, when all of the time steps were included (Fig. 2.5). Among the five regression models evaluated, the only model supported by Akaike weights for TP or SRP was the 1996-only model (Fig. 2.4). 56 All other lake water characteristics were best modeled by some combination of current and past LULC data (Fig. 2.4). All species of nitrogen were best modeled by including historical LULC information through the 1978 time step (Fig. 2.4), increasing the explanatory power from 18-27% for the current-LULC-only regression approach to 38-43% for the model including 1978 (Fig. 2.5). Calcium, silica, and sulfate were also best modeled by a 1978 legacy model (Fig. 2.4), accounting for 47-62% of the variation in the data (Fig. 2.5). Other legacy models for these variables also received some support by AIC weights. The 1968 legacy models for calcium and silica were supported by AIC weights of 0.27 and 0.25, respectively. Conservative ions showed the longest legacy effect and the highest explanatory power of all water chemistry characteristics in our study. The best model for potassium included historical LULC information through the 1955 time step (Fig. 2.4), accounting for 93% of the variation in the data (Fig. 2.5). Explanatory power increased the most (approximately 30%) when 1978 LULC was added to the regression model, but continued to increase through the addition of 1955 (Fig. 2.5). Chloride, magnesium, and sodium were best modeled by the full legacy model that included LULC information for all time steps back to 1938 (Fig. 2.4). Although the 1996-only model for sodium was also supported with an AIC weight of 0.23, the full legacy model received the majority of support, with an AIC weight of 0.59. The 1938-legacy models of chloride, magnesium, and sodium improved explanatory power by 21%, 35%, and 20%, respectively, over the 1996-only model (Fig. 2.5). These models accounted for 82-85% of the variation in the data (Table 2.2). All regression models show a roughly linear relationship (Appendix 5). DISCUSSION 57 Among key ecological principles essential in sustainable land planning, Dale et al. (2000) describe the relevance of ‘time’ in terms of ecological legacies, time-lag feedbacks, and future constraints. Our study sought to quantify this temporal component of ecosystem response to land conversion (i.e. land use legacy effect). The method we used shows both the presence and timescale of legacy effects in models relating LULC to most lake water characteristics. Furthermore, by comparing lake water characteristics with vastly different biogeochemical cycles, our study gives mechanistic insight into legacies. We report finding the first evidence of land use legacy effects in lakes. The traditional approach to studying ecosystem responses to LULC compares LULC from a single time period to a measure of ecosystem response from a similar time period. A legacy approach simply modifies the traditional approach by adding historical LULC to models. We found that taking a legacy approach improved model fits for total and dissolved nitrogen, reactive ions, and conservative ions. However, the added complexity inherent to legacy models did not substantially improve explanatory power for total and soluble reactive phosphorus. These results seem to support a conclusion that phosphorus responds most strongly to the most recent LULC. Another plausible explanation is that lake phosphorus concentrations are governed more by internal recycling from lake sediments rather than delivery pathways (Stauffer 1987, Soranno et al. 1996). The low explanatory power for our models of phosphorus also indicates that other mechanisms, such as internal recycling, should be considered in future analyses of phosphorus. Results from other studies linking watershed landscape context to phosphorus also show that external landscape factors fail to predict lake phosphorus well (Riera et al. 2000, Martin and Soranno 2006, Martin et al. In Review). 58 In contrast to phosphorus, legacy effects were more apparent for nitrogen species. Nitrogen is also a highly reactive nutrient; however, nitrogen is readily held in multiple biologically available forms and remains dissolved much more easily than phosphorus. Previous studies provide evidence that nitrogen species are delivered primarily as dissolved forms through subsurface hydrologic pathways (Jordan et al. 1997, Vanni et al. 2001), whereas phosphorus binds to soils and is delivered primarily through surface runoff (Jordan et al. 1997, Johnes and Heathwaite 1997). Differences in the rate of delivery between these two routes, surface delivery being fast and subsurface delivery being slow, provide another potential mechanistic explanation for differences in legacy timescales for nitrogen and phosphorus. Carpenter and Turner (2000) discuss how interactions between the proverbial tortoises and hares of ecological processes work to structure ecosystem responses. Our results provide empirical evidence supporting the importance of both slow (soil saturation) and fast (surface runoff) processes in structuring lake ecosystem response to external nitrogen loads. By grouping our water chemistry variables into a priori groups based on biological reactivity, we are able to show that the level of biological reactivity can also be used to understand differences in legacy timescales. As major nutrients for aquatic primary producers, phosphorus and nitrogen have high levels of reactivity. Calcium (Alstad et al. 1999), silica (Martin-Jezequel et al. 2000) and sulfate (Holmer and Storkholm 2001) are also important in biological processes, but less than nitrogen (Wetzel 1983). On the other hand, Cl, K, Mg, and Na are grouped as conservative ions given their low relative reactivity in aquatic ecosystems (Wetzel 1983). Our results indicate that legacy timescales are negatively related to biological reactivity: highly reactive elements, such as nutrients have relatively short legacy timescales and conservative ions have the longest legacy timescales. From these results, we can estimate that the 59 signal from nitrogen and reactive ions are still apparent from historical LULC occurring 30 years ago. In comparison, conservative ion concentrations, which are not as strongly linked to biological metabolism, are still responding to historical LULC from 70 years ago. It is possible that if LULC data were available prior to 1938, many of the conservative ions would have shown an even longer legacy. The few biogeochemical studies of legacy effects have found similar legacy timescales for nutrients (Harding et al. 1998, Chauvat et al. 2007) and conservative ions (Boutt et al. 2001). We also see an increase in explanatory power as reactivity decreases, with forms of phosphorus explained poorly (R2 22 and 25%), nitrogen explained fairly well (R2 37-43%), reactive ions explained well (R2 47-62%), and conservative ions explained extremely well (R2>80%). Kling et al. (2000) found that biological reactivity was positively related to the coefficient of variation calculated from multiyear sampling on ten lakes, with higher reactivity showing higher interannual variability. This high variability may also explain differences in explanatory power of the legacy models, in that high interannual variability of a variable in a lake will likely reduce the predictability (hence model power) for that variable. Management Implications A recent renewal of interest in eutrophication (Conley et al. 2009, Smith and Schindler 2009, Wagner and Adrian 2009) and calls for nutrient reductions (Carpenter 2008) echo the calls of the past (Edmonson et al. 1956, Edmonson 1970, Schindler 1974). Many management strategies, such as changes in land use policies and reductions in point source inputs, have been put into place to reduce anthropogenic effects on lakes (Schindler 2006). Setting realistic expectations is important part of the management cycle and failure to recognize legacy impacts 60 may lead to unachievable expectations. While the occurrence of legacy effects may be more widely known or accepted, the quantification of legacy timescales is important for setting realistic goals for the restoration of ecosystem services (Chauvat et al. 2007). Our study shows a nitrogen legacy signal from more than a decade prior. On the other hand, our study did not show a legacy effect for phosphorus. Therefore, management actions aimed at reducing nitrogen impacts may take longer to show improvements than actions aimed at reducing phosphorus impacts. The long legacy that conservative ions exhibited may also be indicative of the legacy timescale of other soluble substances, such as persistent toxic chemicals In conclusion, relationships between LULC and lake water characteristics go beyond that which is currently being measured. It is essential to include the temporal context of LULC in modeling important ecosystem dynamics. By adding an analysis of conservative ions to the suite of limnological variables, we were able to show differences in inferred legacy mechanisms that would have otherwise been missed. Furthermore, ecosystem responses to future changes in land cover and/or climate are likely linked to historic landscape context. 61 Table 2.1. Eigenvectors and proportional variance explained from principal components analysis of land use/cover types. Extreme values for each PC are indicated in bold. PC1 PC2 PC3 PC4 PC5 Agriculture -0.8663 -0.1345 -0.0426 -0.0509 -0.2457 Urban 0.3302 -0.7274 -0.0216 -0.4412 0.0075 Forest 0.1567 0.5363 0.5599 -0.4196 -0.1779 Water 0.2239 -0.1552 0.2985 0.7826 -0.2399 Open 0.2416 0.3533 -0.7711 0.0099 -0.2355 Wetland -0.0862 0.1276 -0.0231 0.1191 0.8916 Eigenvalue % variance explained 0.043 0.024 0.013 0.007 0.003 48% 27% 14% 7% 4% 62 Table 2.2. Regression statistics for best models. Dependant variables are listed under broad categories. Details of regression models include: time step(s) represented in final model, number of regression parameters fit (p), Akaike Information Criteria (AIC), Akaike weight (ωi), and coefficient of determination (R2). Regression model p ωi R2 TP 96 3 0.872 0.255 SRP 96 3 0.930 0.216 TN 96+78 6 0.639 0.376 NO3 96+78 6 0.788 0.400 NH4 96+78 6 0.722 0.435 Ca 96+78 6 0.494 0.496 Si 96+78 6 0.611 0.471 SO4 96+78 6 0.672 0.616 Cl 96+78+68+55+38 15 0.926 0.850 K 96+78+68+55 12 0.870 0.933 Mg 96+78+68+55+38 15 0.968 0.847 Na 96+78+68+55+38 15 0.593 0.815 Dependant Phosphorus Nitrogen Reactive Ions Conservative Ions 63 Figure 2.1. Map showing the Huron River Watershed within the state of Michigan, including detailed hydrography features along with outlines for the cities of Ann Arbor and Ypsilanti. For interpretation of the references to color in this and all other figures, the reader is referred to the electronic version of this dissertation. Ann Arbor Ypsilanti 64 Figure 2.2. Average proportion of each land use/cover type for the study lakes calculated for each time step in the study. Solid fine line represents agriculture. Solid heavy line represents urban. Dashed fine line represents forest. Dashed heavy line represents open. Dotted fine line represents water. Dotted heavy line represents wetland. 0.5 Urban Agriculture 0.4 Open Forest Proportion Wetland 0.3 0.2 0.1 0 1938 1955 1968 65 1978 1996 Figure 2.3. Change in land use/cover over time represented by principal component 1 (PC1), principal component 2 (PC2), and principal component 3 (PC3). Each dot represents the average score for that time step. 0.15 PC1 PC2 PC3 0.1 0.05 0 1938 1955 1968 -0.05 -0.1 -0.15 -0.2 66 1978 1996 Figure 2.4. Akaike weights (ωi) from chrono-sequence of regression models. (A) TP, solid circles with heavy solid line; and, SRP, open circles with dashed line. (B) TN, solid circles with heavy solid line; NO3, open circles with dashed line; and, NH4, X inside circle with dotted line. (C) Ca, open circle with heavy solid line; Si, open triangle with dashed line; and, SO4, open square with dotted line. (D) Cl, open with circle heavy solid line; K, open triangle with dashed line; Mg, open square with dotted line; and, Na, open diamond with fine solid line. 1 0.8 0.6 ωi 0.4 0.2 0 1 0.8 0.6 ωi 0.4 0.2 0 96 96 78 96 78 68 96 78 68 55 96 96 78 68 55 38 67 96 78 96 78 68 96 78 68 55 96 78 68 55 38 2 Figure 2.5. Coefficient of determination (R ) from chrono-sequence of regression models. (A) TP, solid circles with heavy solid line; and, SRP, open circles with dashed line. (B) TN, solid circles with heavy solid line; NO3, open circles with dashed line; and, NH4, X inside circle with dotted line. (C) Ca, open circle with heavy solid line; Si, open triangle with dashed line; and, SO4, open square with dotted line. (D) Cl, open with circle heavy solid line; K, open triangle with dashed line; Mg, open square with dotted line; and, Na, open diamond with fine solid line. 1 A. Phosphorus B. Nitrogen C. Reactive Ions D. Conservative Ions 0.8 0.6 R2 0.4 0.2 0 1 0.8 0.6 R2 0.4 0.2 0 96 96 78 96 78 68 96 78 68 55 96 96 78 68 55 38 68 96 78 96 78 68 96 78 68 55 96 78 68 55 38 APPENDIX 69 Appendix 2. Additional information for Chapter 2. Table 2.3. Characteristics of the study lakes including water chemistry variables as categorized into a priori groups and lake area. Minimum, maximum, mean, and coefficient of variation are shown for each variable. Lake characteristic 2 Lake area (km ) Phosphorus TP (µg/L) SRP (µg/L) Nitrogen TN (mg/L) NO3 (mg/L) NH4 (µg/L) Reactive Ions Ca (mg/L) Si (mg/L) SO4 (mg/L) Conservative Ions Cl (mg/L) K (mg/L) Mg (mg/L) Na (mg/L) Minimum Maximum Mean CV 0.06 2.60 0.65 0.93 5.4 0.33 45.0 3.45 24.1 1.15 0.37 0.56 0.53 1.85 0.94 0.30 0.00 1.03 0.20 1.20 3 143 36 0.96 16 0.00 92 5.61 52 1.53 0.38 0.88 1.9 88.4 26.1 0.83 2 0.58 3.9 0.9 242 7.53 24.7 77.3 66 2.09 15.4 24.0 0.90 0.57 0.26 0.88 70 Table 2.4. Regression coefficients from best models as reported in Table 2.2. Ca Intercept PC1_96 PC2_96 PC3_96 PC1_78 PC2_78 PC3_78 PC1_68 PC2_68 PC3_68 PC1_55 PC2_55 PC3_55 PC1_38 PC2_38 PC3_38 Cl K Mg Na 49 -42 -145 335 -23 147 -353 -44 -705 -377 491 1362 5.6 338 -657 298 -681 590 -603 436 -746 255 -471 1.4 -13 -11 14 13 11 -8.5 1.0 -3.3 -2.6 -4.1 1.8 -0.1 17 4.7 -4.2 24 7.2 -23 -11 -11 27 -13 -20 -27 -23 19 33 37 -15 -228 -129 352 451 34 -60 -241 107 -232 245 -267 169 -284 70 -200 NH4 NO3 37 0.352 -416 -6.4 -175 -2.2 -731 0.5 503 5.5 159 1.7 751 -0.2 71 Si SO4 SRP TN TP 1.2 -13 -14 1.5 11 20 -6.6 22 135 -61 527 -212 73 -536 1.3 -1.0 1.1 0.4 1.1 -8.3 -3.6 -0.6 7.2 3.1 1.2 27 -27 -2.4 -12 Figure 2.6. Comparisons between regression results using 2 PCs and 3 PCs for calcium, chloride, potassium, magnesium, sodium, ammonium, nitrate, silica, sulfate, soluble reactive phosphorus (SRP), total nitrogen (TN), and total phosphorus (TP). Akaike weights (y-axis) were computed across regression results for each water chemistry variable from both 2 PCs and 3 PCs to total ten regression models. AIC weights 100% Calcium 80% 60% 40% 20% 0% 2 PC's 3 PC's AIC weights 100% Chloride 80% 60% 40% 20% 0% 2 PC's 3 PC's 72 Figure 2.6 (cont’). AIC weights 100% Potassium 80% 60% 40% 20% 0% 2 PC's 3 PC's AIC weights 100% Magnesium 80% 60% 40% 20% 0% 2 PC's 3 PC's 73 Figure 2.6 (cont’). AIC weights 100% Sodium 80% 60% 40% 20% 0% 2 PC's 3 PC's AIC weights 100% Ammonium 80% 60% 40% 20% 0% 2 PC's 3 PC's 74 Figure 2.6 (cont’). AIC weights 100% Nitrate 80% 60% 40% 20% 0% 2 PC's 3 PC's AIC weights 100% Silica 80% 60% 40% 20% 0% 2 PC's 3 PC's 75 Figure 2.6 (cont’). AIC weights 100% Sulfate 80% 60% 40% 20% 0% 2 PC's 3 PC's AIC weights 100% SRP 80% 60% 40% 20% 0% 2 PC's 3 PC's 76 Figure 2.6 (cont’). AIC weights 100% TN 80% 60% 40% 20% 0% 2 PC's 3 PC's AIC weights 100% TP 80% 60% 40% 20% 0% 2 PC's 3 PC's 77 Figure 2.7. Change in LULC over time as represented in multivariate space. Each line represents the LULC trajectory of a lake in the study. Each of the first three PC’s are shown. 0.5 0.4 0.3 0.2 PC1 0.1 0 -0.1 -0.2 -0.3 -0.4 -0.5 -0.6 1938 1955 1968 78 1978 1996 Figure 2.7 (cont’). 0.6 0.4 PC2 0.2 0 -0.2 -0.4 -0.6 -0.8 1938 1955 1968 79 1978 1996 Figure 2.7 (cont’). 0.4 0.3 0.2 0.1 PC3 0 -0.1 -0.2 -0.3 -0.4 -0.5 1938 1955 1968 80 1978 1996 Figure 2.8. Regression diagnostics for best models as reported in Table 2.2. Observed values (yaxis) are plotted against values predicted from the regression equation (x-axis). 81 Figure 2.8 (cont’). 82 REFERENCES 83 REFERENCES Allan, J.D. 2004. Landscapes and riverscapes: the influence of land use on stream ecosystems. Annu. Rev. Ecol. Evol. Syst. 35: 257–84. Alstad, N.E.W., L. Skardal, and D.O. Hessen. 1999. The effect of calcium concentration on the calcification of Daphnia magna. Limnology and Oceanography 44, 2011–2017. Anderson, J. R., E. H. Harvey, J. T. Roach, and R. E. Whitman.1976. A land use sensor and land cover classification system for use with remote sensor data geological survey. Professional Paper 964. US Government Printing Office. Arbuckle, K.E., and J.A. Downing. 2001. The influence of watershed land use on lake N: P in a predominantly agricultural landscape. Limnology and Oceanography., 46(4): 970–975. Bennett, E.M., S.R. Carpenter, and M.K. Clayton. 2004. Soil phosphorus variability: scaledependence in an urbanizing agricultural landscape. Landscape Ecology 20: 389-400. Boutt, D.F., D.W. Hyndman, B.C. Pijanowski, and D.T Long. 2001. Identifying potential land use-derived solute sources to stream baseflow using ground water models and GIS. Groundwater 39(1): 24-34. Brett, M.T., G.B. Arhonditsis, S.E. Mueller, D.M. Hartley, J.D. Frodge, and D.E. Funke. 2005. Non-point-source impacts on stream nutrient concentrations along a forest to urban gradient. Environmental Management 35(3): 330-342. Burcher, C.L., and E.F. Benfield. 2006. Physical and biological responses of streams to urbanization of historically agricultural watersheds. Journal of the North American Benthological Society 25(2): 356–369. Burnham, K.P., and D.R. Anderson, 2002. Model Selection and Multimodel Inference: A Practical-Theoretic Approach, 2nd ed. Springer-Verlag. ISBN 0-387-95364-7. Carpenter, S. R. 2008. Phosphorus control is critical to mitigating eutrophication. Proceedings of the National Academy of Sciences 105(32):11039-11040. Carpenter, S.R. and M.G. Turner. 2000. Hares and Tortoises: Interactions of Fast and Slow Variables in Ecosystems. Ecosystems 3(6): 495-497. Conley, D.J., G.E. Likens, D.C. Buso, L. Saccone, S.W. Bailey, and C.E. Johnson. 2008. Deforestation causes increased dissolved silicate losses in the Hubbard Brook Experimental Forest. Global Change Biology 14: 2548-2554. 84 Dale, V. H., S. Brown, R. A. Haeuber, N. T. Hobbs, N. Huntly, R. J. Naiman, W. E. Riebsame, M. G. Turner, And T. J. Valone. 2000. Ecological principles and guidelines for managing the use of land (ESA Report). Ecological Applications 10(3): 639–670. Davies, S.P., and S.K. Jackson. 2006. The biological condition gradient: a descriptive model for interpreting change in aquatic ecosystems. Ecological Applications 16(4): 1251–1266. Declerck, S., T. De Bie, D. Ercken, H. Hampel, S. Schrijvers, J. Van Wichelen, V. Gillard, R. Mandiki, B. Losson, D. Bauwens et al., 2006. Ecological characteristics of small farmland ponds: Associations with land use practices at multiple spatial scales. Biological Conservation 131(4): 523-532. Dillon, P. J., and W.B. Kirchner. 1975. The effects of geology and land use on the export of phosphorus from watersheds. Water Research 9: 135–148. Edmonson, W.T., G.C. Anderson, and D.R. Peterson. 1956. Artificial eutrophication of Lake Washington. Limnology and Oceanography 1(1): 47-53. Edmonson, W.T. 1970. Phosphorus, nitrogen, and algae in Lake Washington after diversion of sewage. Science 169(3946): 690-691. Findlay, S., J.M. Quinn, C.W. Hickey, G. Burrell, M. Downes. 2001. Effects of land use and riparian flowpath on delivery of dissolved organic carbon to streams. Limnology and Oceanography 46(2): 345-355. Foster, D. R., G. Motzkin, and B. Slater. 1998. Land-Use History as Long-Term Broad-Scale Disturbance: Regional Forest Dynamics in Central New England. Ecosystems 1(1): 96-119. Foster, D., F. Swanson, J. Aber, I. Burke, N. Brokaw, D. Tilman, and A. Knapp. 2003. The importance of land-use legacies to ecology and conservation. Bioscience 53(1): 77-88. Foy, R. H., and S. D. Lennox. 2006. Evidence for a delayed response of riverine phosphorus exports from increasing agricultural catchment pressures in the Lough Neagh catchment. Limnology and Oceanography 51(1, part 2): 655–663. Goodman K.J., A.E. Hershey, and K. Fortino. 2006. The effect of forest type on benthic macroinvertebrate structure and ecological function in a pine plantation in the North Carolina Piedmont. Hydrobiologia 559: 305-318. Harding, J.S., E.F. Benfield, P.V. Bolstad, G.S. Helfman, and E.B.D. Jones III. 1998. Stream biodiversity: The ghost of land use past. Proceedings of the National Academy of Sciences of the United State of America 95(25): 14843-14847. Hay-Chmielewski, E. M., P. W. Seelbach, G. E. Whelan, and D. B. Jester, Jr. 1995. Huron River assessment. Fisheries Special Report No 16. Michigan Department of Natural Resources. 85 Holmer, M. and P. Storkholm. 2001. Sulphate reduction and sulphur cycling in lake sediments: a review. Freshwater Biology 46: 431-451. Holtan, H., L. Kamp-Nielsen, and A.0. Stuanes. 1988. Phosphorus in soil, water and sediment: an overview. Hydrobiologia 170: 19-34. Johnes P.J. and A.L. Heathwaite 1997. Modelling the impact of land use change on water quality in agricultural catchments. Hydrological Processes 11: 269-286. Jordan, T.E., D.F. Whigham, K.H. Hofmockel, and M.A. Pittek. 1997. Nutrient and sediment removal by a restored wetland receiving agricultural runoff. Journal of Environmental Quality 32: 1534-1547. Kling, G.W., G.W. Kipphut, M.M. Miller, and J. O’Brien. 2000. Integration of lakes and streams in a landscape perspective: the importance of material processing on spatial patterns and temporal coherence. Freshwater Biology 43: 477-497. Knoll, L.B., M.J. Vanni, and W.H. Renwick. 2003. Phytoplankton primary production and photosynthetic parameters in reservoirs along a gradient of watershed land use. Limnology and Oceanography, 48(2): 608–617. Larsen, D.P., D.W. Schults, and K.W. Malueg. 1981. Summer internal phosphorus supplies in Shapawa Lake, Minnesota. Limnology and Oceanography 26(4): 740-753. Martin-Jezequel, V. M. Hildebrand, and M.A. Brzezinski. 2000.Silicon metabolism in diatoms: Implications for growth. Journal of Phycology 36: 821-840. Omernik, J.M., 1976. The influence of land use on stream nutrient levels. USEPA Ecological Research Series EPA-60013-76-014. US Environmental Protection Agency, Corvallis, OR. Opperman, J.J., K.A. Lohse, C. Brooks C, N.M. Kelly, and A.M. Merenlender. 2005. Influence of land use on fine sediment in salmonid spawning gravels within the Russian River Basin, California. Canadian Journal of Fisheries and Aquatic Sciences 62(12): 2740-2751. Pellerin, B.A., W.M. Wollheim, C. S. Hopkinson, W.H. McDowell, M.R. Williams, C.J. Vörösmarty, and M.L. Daley. 2004. Role of wetlands and developed land use on dissolved organic nitrogen concentrations and DON/TDN in northeastern U.S. rivers and streams. Limnology and Oceanography 49(4): 910–918. Reed-Andersen T., S.R. Carpenter, R.C. Lathrop. 2000. Phosphorus flow in a watershed-lake ecosystem. Ecosystems 3(6): 561-573. Rutledge, D.T. 2001. Changes in land cover and wildlife habitats in two watersheds in the lower peninsula of Michigan. Ph.D. dissertation, Michigan State University, East Lansing, MI. 86 Rutledge, D.T., and C.A. Lepczyk. 2002. Landscape change: patterns, effects, and implications for adaptive management of wildlife resources. In: Integrating Landscape Ecology into Natural Resource Management. J. Liu, and W.W. Taylor Ed. New York, NY: Cambridge University Press. p 312-333. Schindler, D.W. 1974. Eutrophication and Recovery in Experimental Lakes: Implications for Lake Management. Science 184: 897-898. Schindler, D.W. 2006. Recent advances in the understanding and management of eutrophication. Limnology and Oceanography 51(1, part 2):356-363. Schilling, K.E., and J. Spooner. 2006. Effects of watershed-scale land use change on stream nitrate concentrations. Journal of Environmental Quality 35:2132–2145. Smith, V.J. and D.W. Schidler. 2009. Eutrophication science: where do we go from here? Trends in Ecology and Evolution. 24(4): 201-207. Soranno, P.A., S.L. Hubler, S.R. Carpenter, and R.C. Lathrop. 1996. Phosphorus loads to surface waters: a simple model to account for spatial pattern of land use. Ecological Applications 6(3): 865-878. Stauffer, R.E. 1987. Vertical nutrient transport and its effects on epilimnetic phosphorus in four calcareous lakes. Hydrobiologia 154: 87-102. Tomer, M.D., and M.R. Burkart. 2003. Long-term effects of nitrogen fertilizer use on ground water nitrate in two small watersheds. Journal Environmental Quality 32:2158–2171. Vanni M.J., W.H. Renwick, J.L. Headworth, J.D. Auch, and M.H. Schaus. 2001. Dissolved and particulate nutrient flux from three adjacent agricultural watersheds: A five-year study. Biogeochemistry 54(1): 85-114. Wagner, C. and R. Adrian. 2009. Exploring lake ecosystems: hierarchy responses to long-term change? Global Change Biology 15: 1104-1115. Wetzel, R.G. 1983. Limnology, 2nd ed. Harcourt College 87 CHAPTER 3: THE LAND-USE LEGACY EFFECT - A MECHANISTIC INVESTIGATION OF HOW GROUNDWATER DELIVERY, WETLAND PROCESSING, AND RIPARIAN DYNAMICS AFFECT LAKE WATER CHEMISTRY 88 INTRODUCTION There have been numerous studies linking land use/land cover (LULC) to ecosystem responses. The typical approach involves correlating current LULC within a specified area to current ecosystem condition(s). The main difference among these studies is the way in which the contributing zone has been delineated, including fixed-distance (Strayer et al. 2003, King et al. 2005) and flow distance buffers around a sampling station (Brazner et al. 2007), fixed-distance buffers around a stream (Hunsaker & Levine 1995, Van Sickle et al. 2004, Floyd et al. 2009) with varying distance upstream of a station (Sponseller et al. 2001, Frimpong et al. 2005), and, finally, different variations on “watersheds” (Soranno et al. 1996, Hollister et al. 2008). Regardless of spatial scale, these approaches assume that materials originating from different cover types are delivered to the ecosystem over the same time period that the LULC is measured. This is likely a reasonable assumption for surface delivery. Burcher (2009) estimated that the maximum travel time for overland flow from precipitation was on the order of days. However, groundwater delivery often takes longer than a decade and can exceed a century (Pint et al. 2003, Pijanowski et al. 2007). Therefore, it is likely that much of the water delivered through groundwater pathways is not representative of the current LULC. In systems where groundwater provides the dominant source of water, this temporal mismatch can obscure relationships between LULC and ecosystem responses. While there have been several studies that highlight the importance of representing groundwater geochemistry (Wayland et al. 2003) and transport delays in models of aquatic ecosystem response (Baker et al. 2006, Fraterrigo & Downing 2008, Kelly et al. 2008), few have linked the dynamics of changing LULC to a mechanistic understanding of flow paths and travel times (Boutt et al. 2001, Wayland et al. 2002). Pijanowski et al. (2007) coupled a land 89 transformation model with a groundwater travel time model to produce a temporally adjusted mosaic of LULC that they term the land-use legacy map. They found large discrepancies between current LULC and legacy LULC in their study watershed. For example, there was only a 22% agreement in urban LULC type between current and legacy maps. Additionally, 11% of their study area that was currently in a human dominated LULC type (i.e., urban and agriculture) was assigned to forested cover type in the legacy map. Using a legacy map of LULC, as described above, should improve models of ecosystem response to LULC. In addition to groundwater flow paths, wetland processing and riparian zone dynamics modify how current LULC impacts aquatic ecosystem condition. Peterjohn and Correll (1984) showed that riparian forests can effectively capture surface and subsurface nutrient fluxes from adjacent agricultural land to stream systems. In recent reviews of riparian buffer characteristics, Zhang et al. (2009) and Yuan et al. (2009) agree that buffer width and composition are important in determining the flux of non-point pollutants. Other studies have shown that stream and lake water chemistry characteristics are related to wetland extent and proximity (Osborne & Wiley 1988, Johnston et al. 1990, Detenbeck et al. 1993). However, the appropriate spatial scale for wetland measurement, whole watershed proportion versus proportion within a specified distance, remains an active area of inquiry (Gergel et al. 1999, Xenopoulos et al. 2003, Moreno-Mateos et al. 2008). The goal of this study is to combine the knowledge of groundwater delivery, wetland and riparian processing in mechanistic and statistical models to better understand the relationship between LULC and lake water chemistry characteristics. We created a mechanistic groundwater flow model to estimate spatially-explicit groundwater travel times for a watershed to surface water bodies. We used this information to create a legacy LULC map, which accounts for delays 90 in materials delivered via groundwater pathways from past LULC. Three approaches for incorporating wetland processing were compared: 1) wetland extent, 2) area of wetland connected to rivers, and 3) length of wetland connected to rivers. Riparian processes are incorporated by adding measures of riparian LULC to the mechanistic groundwater legacy model. We use regression models to relate lake water chemistry variables to various LULC representations. We applied these modeling approaches to twelve water chemistry variables, ranging from nutrients to conservative ions, to better understand the roles that biological reactivity and solubility play in connections between LULC and aquatic ecosystem function. We expect that: 1) chemicals of high solubility (e.g. SRP, NO3) will have a stronger relationship with groundwater flow paths than chemicals of low solubility (e.g. TP), which will be more dependent upon surficial transport processes, such as retention in riparian zones, than groundwater flow paths; 2) chemicals of low biological reactivity (e.g. conservative ions) will have a stronger link to groundwater flow paths than chemicals of high biological reactivity (e.g. nitrogen and phosphorus), which will be more dependent on wetland processing than groundwater flow paths. METHODS Study Area 2 This study was conducted within the 2,359 km Huron River Watershed (HRW) in Michigan (Figure 3.1), which ranges from 390m to 173m in elevation. The HRW contains numerous rivers, streams and lakes supporting a diverse assemblage of aquatic species (Hay91 Chmielewski et al. 1995). Due to close proximity to major metropolitan areas, such as Detroit and Ann Arbor, and the development of interstate highways, many people live, work, and enjoy the recreational opportunities of the area. The watershed has undergone extensive LULC change over the past century, shifting from an agrarian to suburban society (Hay-Chmielewski et al. 1995, Rutledge and Lepczyk 2002). However, numerous natural areas have been maintained in publicly owned parks. Water chemistry data Water samples were collected from 35 lakes in the HRW during 2008 spring mixing. Lakes were chosen based on their accessibility, including both private and public lakes ranging 2 in size from 0.05 to 2.6 km . All water samples were taken from 1m below the surface at the deepest portion of the lake and analyzed for major ions and nutrients. Cation (calcium, magnesium, potassium, and sodium) and anion (chloride, nitrate, and sulfate) concentrations were determined using membrane-suppression ion chromatography (Wetzel and Likens 2000). Silica concentrations were determined using the molybdate colorimetric method (Wetzel and Likens 2000). Total nitrogen concentrations were determined using the 2nd derivative of the absorbance curve at 224 nm following persulfate digestion. Ammonia concentrations were determined following the indophenol-blue method (Wetzel and Likens 2000). Soluble reactive phosphorus (SRP) and total phosphorus (TP) concentrations were determined spectrophotometrically following persulfate digestion. Water chemistry variables were organized into a priori groups based largely on biological reactivity and solubility (Appendix 3.1). Phosphorus species were the most highly reactive variables in our study; nitrogen species are also highly reactive, but less so than 92 phosphorus due to the relative abundance of nitrogen versus phosphorus (Wetzel 1983). We grouped Ca, Si, and SO4 together as “reactive variables” due to their role as micronutrients in the metabolism of macrophytes, diatoms, and bacteria, respectively (Wetzel 1983). Conservative ions (i.e. Cl, K, Na, and Mg) play a much smaller role in biological metabolism and are thus influenced more by external supply than internal conversions (Wetzel 1983). Land use/cover data Land use/cover data (30m resolution) from five available time steps (1938, 1955, 1968, 1978, and 1995) were classified into six categories (urban, agriculture, open, forest, water, and wetland) based a modified version of the Anderson (1976) LULC classification scheme (Figure 3.2). All data were compiled into a multi-temporal GIS database for analysis (full details in Rutledge 2001). Briefly, digitized land use/cover data for 1978 from the Michigan Resource Information System (MIRIS) served as a base for digitizing all other time steps. Aerial photos from each of the other time steps were scanned (150 dpi) to create digital images that were then registered and rectified to the 1978 data using the placement of county roads. Land use/cover polygons were then digitized. Groundwater travel time calculation and modeling We modeled groundwater travel times following procedures similar to Boutt et al. (2001) and Pijanowski et al. (2007). This approach is based on Darcy’s Law of groundwater flow: Q = -KAi where, Q is discharge in m3/day, K is hydraulic conductivity in m/day, A is the cross sectional area in m2, and i is the hydraulic gradient. Groundwater flow velocity (v, m/day) is calculated as: 93 v= Q An where n is porosity. Travel time is then calculated by dividing the flow length by flow velocity. We used the Groundwater Modeling System preprocessor (BYU 1994) to facilitate the discretization of geospatial data for input to MODFLOW-2000 (Harbaugh et al. 2000) to solve for the groundwater flow equations. We used a one-layer model, with 180,000 cells with rd approximately 109 x 109 m dimensions. We used elevations from USGS 1/3 arc second National Elevation Dataset (NED, resolution approx. 10m) to represent surface topography. We also used the NED surface elevations to model stream coverage using procedures and tools in ArcInfo 8.3 (ESRI, Inc): FILL SINKS, FLOW DIRECTION, and FLOW ACCUMULATION. The location of stream cells were defined according to a threshold of 15,000 cells using the STREAM DEFINITION tool, which produced stream cell locations comparable to known hydrography features. Stream cells were combined with cells along lake edges to create a complete representation of hydrography features. These features were modeled as drains in the groundwater model. Bedrock elevations were provided by MSU Department of Geography (D. Lusch personal communication). Recharge was estimated using measured flows from the Ypsilanti USGS gauge on the Huron River (#4174800) for 1974-1994. Low flow values for each year in this period of record were averaged and divided by the drainage area, providing the recharge estimate of 0.00138 m/d (49 cm/yr). This is roughly 55% of the mean annual precipitation measured in the near-by city of Ann Arbor. Static water levels from 15,581 wells recorded in the Michigan Department of Environmental Quality Statewide Groundwater Database, called Wellogic, were used to 94 interpolate water table elevations across an initial model boundary (HRW plus 1km buffer). We filtered this dataset to only include wells with the following characteristics: 1) physically located within the model boundary, 2) reported surface elevation within 2m of NED surface elevation at that location, and 3) casing deeper than the reported static water level. This groundwater surface was then used to delineate a groundwater sourceshed for all river cells upstream of Ypsilanti, as all of the study lakes were located upstream of Ypsilanti. We use the term sourceshed in parallel to how watershed is used for surface water drainage. Therefore, the groundwater sourceshed represents the area contributing groundwater to a particular point (e.g. study lake). The modeled groundwater sourceshed, which we will refer to as the Huron River Groundwater Sourceshed (HRGW), was then used as the model boundary for the MODFLOW-2000 groundwater model. The interpolated static water levels were also used as starting heads in the model. A digital version of the Farrand and Bell (1982) map detailing the location of surficial geology types was obtained from the Michigan Geographic Data Library (http://www.mcgi.state.mi.us/mgdl). The area within the HRGW has some variation in surficial geology types, but is dominated by glacial outwash and end moraine deposits (Figure 3.3). Hydraulic conductivity (K) values for each surficial geology type were optimized using PEST 10.0 (Doherty 2004) (Figure 3.3). Simulated groundwater elevations were in reasonable agreement with observations from the Wellogic database (Figure 3.4). We found that wells with high residual values were surrounded by other wells with groundwater elevations similar to the 2 modeled values, had low influence on the R when removed, and therefore were kept in the welldataset despite being clear outliers. Groundwater flow velocities and travel times were calculated following Darcy’s Law, as described above, from the simulated groundwater 95 elevations and optimized hydraulic conductivity values, using FLOWDIRECTION and FLOWLENGTH tools in ArcGIS (Pijanowski et al. 2007). Creating the legacy map The legacy map was created by combining simulated groundwater travel times with interpreted historical LULC categories for each model cell. Groundwater travel times were calculated using Darcy’s Law, as described above. The calculated travel times (Figure 3.5 – step 1) were grouped around available LULC data (Figure 3.5 – step 2), with 1996 being the most recent land cover map available in this area. The mid-point between the LULC time steps was used to define category thresholds. For example, there is an eighteen year gap between the 1996 and 1978 LULC maps, thus we assigned LULC from the 1996 time step to all model cells with a groundwater travel time of zero (1996) to nine years earlier (1987). These reclassified travel times (Figure 3.5 – step 2) were then combined with the LULC maps developed from air photo analysis for particular years to produce the legacy LULC map (Figure 3.5 – step 3). Therefore, the legacy map is a spatially explicit representation of the LULC, corresponding to groundwater delivery to each study lake. Finally, we compared the legacy LULC to the 1996 LULC, highlighting areas where the two differ (Figure3. 5 – step 4). Statistical analyses The relationships between LULC and lake water chemistry variables were analyzed using multiple linear regression, implemented in PROC REG (SAS). Each water chemistry variable was regressed against the proportional cover of the six LULC classes within the study lake groundwater sourcesheds, following the equation: 96 yi = β0 + βurbxurb + βagxag + βopenxopen + βforxfor + βwaterxwater + βwetxwet + εi where, yi is the value of water chemistry i; β0 is the constant (y-axis intercept); βp is the regression coefficient for xp; xp is the proportional cover of LULC p; and εi is the error associated with the model of yi. The proportion of each LULC type within the sourcesheds was determined for the 1996 time step (current LULC) and the legacy map (GW legacy). We modified the general regression equation (above) to incorporate wetland and riparian processes. Wetland processes were represented by replacing the single regression parameter describing wetlands (xwet) with two terms for wetlands: 1) proportional area of wetlands that border drain cells (xconwet), and 2) proportional area of wetlands that are unconnected to drain cells (xunconwet), following the form: yi = β0+βurb xurb+ βag xag+βopen xopen+βfor xfor+βwater xwater+βconwet xconwet+ βunconwet xunconwet+ εi Connected and unconnected wetlands were summarized within two spatial extents: 1) proportion within the groundwater sourceshed that represents processing by wetland area, and 2) proportion within a 30m river buffer that represents processing along the wetland length. Riparian processes were incorporated by adding parameters to the general regression equation for the proportional cover of each of the six LULC classes as represented in 1996 within a 50m river buffer. Therefore, the regression equations modeling overland flow has six LULC types measured over two extents, equaling twelve regression parameters, as follows: yi = β0 + LEGACY[βurbxurb + … + βwetxwet] + RIPARIAN[βurbxurb + …+ βwetxwet] + εi 97 We used Akaike Information Criteria (AIC), AIC weights (ωi) and the coefficient of 2 determination (R ) (following Burnham and Anderson 2002) to compare the parsimony and explanatory power of each model. Specifically, to understand if groundwater pathways are the dominant mechanism for legacy effects, we compared results from: 1) current LULC models, 2) GW legacy models, and 3) a correlative estimation of legacy land cover (as described fully in Chapter 2). Briefly, the correlative legacy models use regression to relate LULC, as represented through principal components analysis, to water chemistry variables. These models were built sequentially, starting with the most recent LULC then adding the next most recent LULC and so on until each of the five time steps of LULC were included in the model. For the current study, we limited our comparison to only the correlative models with highest support as demonstrated by AIC weights. The correlative legacy models represented legacy effects through a correlational “black box”, in that no specific mechanisms were specified. In comparison, GW legacy models specified groundwater pathways as the dominant mechanism for legacy effects. The current LULC models represented the traditional approach linking LULC to ecosystem responses and hypothesizes that historical LULC has no consequence for our study lakes. Comparisons were also made with models measuring wetland and riparian processes. Hypotheses inherent to the models were: wetland presence throughout the entire lake watershed was important to lake water chemistry (wetland extent tested); the entire area of riparian wetlands was important to lake water chemistry (wetland area tested); only the area of riparian wetland in immediate contact with the river was important to lake water chemistry (wetland length tested); the LULC within only the riparian zone was important to lake water chemistry (riparian LULC only tested); there is an additive effect of riparian LULC to the watershed LULC (legacy + riparian tested) . 98 RESULTS Groundwater travel time Modeled groundwater travel times in the HRGW ranged from <1 yr to >400 years, encompassing each of the six time classes of available LULC data (Figure 3.6). Just over half of the HRGW has a groundwater travel time between 0-9 years, following the surface drainage network in the watershed. The longest travel times occurred in the central and northeastern portions of the HRGW. In general, the groundwater travel times in the study lake groundwater sourcesheds match those of the larger HRGW, but have more area within the 0-9 yr travel time category than the overall watershed (Figure 3.6). Characterizing LULC Urban and agriculture land cover dominated the HRGW in 1996, with open and forested land cover types comprising most of the remaining area (Figure 3.7A). The legacy LULC map shows agriculture as the dominant land cover, followed by urban, forest and open cover types (Figure 3.7B). The largest differences between the two LULC maps were in human dominated cover types, as agriculture increased by approximately 5% and urban decreased by approximately 4% of the model area, when legacy LULC was used. The 1996 LULC within the studied groundwater sourcesheds also show an urban dominance (Figure 3.7A), however forest was the second most common LULC type in the groundwater sourcesheds, and there was 12% less agriculture cover than in the HRGW as a whole. Legacy LULC within the studied groundwater sourcesheds had almost equal 99 representation of forest, urban, and open LULC types, and agriculture was still underrepresented relative to the HRGW (Figure 3.7). Overall differences in LULC composition between 1996 and the legacy map within the HRGW and the groundwater sourcesheds were relatively small (Figure 3.8), with only 10.5% of the area of the HRGW where legacy LULC differed from current LULC. This is mainly due to the high hydraulic conductivities of the sediments, coupled to the relatively small size of the watershed. Differences in the studied sourcesheds were even less, at 8.5%. Relationships with lake water chemistry Current LULC and GW legacy models had similar explanatory power (Figure 3.9). 2 Across all water chemistry variables, the R values from these models differed by a maximum of 4% (TN, NH4, and SO4) but were within 1.4% of each other on average. In general, models showed a pattern of high correlation with lower-reactivity soluble ions, such as Cl, K, Mg, and Na, and decreasing correlation for more biologically reactive chemicals such as nitrogen and phosphorus Correlative legacy models received most support from AIC weights for 9 out of 12 of the water chemistry variables (Figure 3.9). Correlative legacy models for TP, NO3, and Si explained 8-11% more variation than current LULC or GW legacy models. However, there were larger differences between the models for all other water chemistry variables. The correlative legacy model explained on average 21% more variation than the GW legacy or current LULC models for all other water chemistry variables, ranging between 19% and 28% more for TN and K, respectively. Models using current LULC or legacy LULC received more AIC support for NH4, 100 Ca, and SO4, and explained 11%, 21%, and 7% more variation, respectively, than the correlative legacy models. Overall, these results provide evidence that adding a mechanistic representation of land cover through delayed groundwater delivery was not sufficient to capture legacy effects for most water chemistry variables in our study lakes. Adding different measures of wetlands (all wetland, wetland area, and wetland length) to the GW legacy model changed the explanatory power of most models very little, changing R 2 less than 3% in most cases (Figure 3.10). However, the explanatory power increased by 14% and 11% for TP and SRP, respectively, when wetland connections were specified in the regression. Interestingly, the method for representing wetland processing differed between TP and SRP; adding the length of connected wetland performed better for the TP model, whereas the area of connected wetland performed better for SRP. Overall, our comparison of different ways to represent wetlands shows that there was little difference for most water chemistry variables. However, forms of phosphorus did show a stronger relationship to LULC when wetland connections were included in the regression. Because wetland length provided similar or slightly better model fits for most water chemistry variables, we only show results for GW legacy models using this term in our analyses of riparian zones. However, wetland area was used in the GW legacy models for SRP. In comparison to adding wetland processing to the GW legacy model, adding factors representing riparian processes increased R2 values for every water chemistry variable (Figure 2 3.11). Models combining GW legacy with riparian processes improved R values by 15% on average and up to 26% (Si) and 43% (SRP) from models of GW legacy alone. Many of these models received more AIC support than correlative legacy models (6 out of 12), with the largest 101 2 increase in R for SRP (46% over correlative legacy). However, some nutrients (TP, TN, and NO3) and conservative ions (K and Mg) were still better described by correlative legacy models than GW legacy combined with riparian processes. Models including riparian LULC alone have little support from AIC weights and had R 2 values similar to or lower than GW legacy for most water chemistry variables. Only SRP and Si show more explanatory power with riparian only models versus GW legacy models. Differences between GW legacy and riparian only models were greater for conservative ions than nutrients. DISCUSSION Legacy effects have been shown to be important in terrestrial (Foster et al. 1998, Foster et al. 2003, Chauvat et al. 2007), stream (Harding et al. 1998, McTammany et al. 2007), and lake (Martin et al. Chapter 2) ecosystems. These studies have all used a correlational approach, relying solely on statistical relations to detect legacy effects. This study moves beyond a simple correlational approach by combining temporal and spatial changes in LULC with a mechanistic model of groundwater flow paths to create a process-informed representation of legacy LULC. We show how this legacy map can be used to connect changes in LULC to important ecosystem characteristics, such as water chemistry in lakes. We also incorporate other mechanisms known to be important, such as wetland and riparian processes, for a more complete analysis of the relationships between LULCs and lake ecosystem dynamics. This mechanistic approach is generalizable across ecosystem types and will likely increase the accuracy and predictability of other models linking LULC with ecosystem responses. This study shows how delays in groundwater travel time can be successfully coupled with changes in LULC as a mechanism for land use legacies. Unfortunately, our study area was 102 not the ideal location for groundwater legacies to play a major role. Rapid changes in LULC occurred in the HRW between 1938 and 1968 (Martin et al. Ch2), a time period represented by only 17% of the groundwater sourceshed. LULC largely stabilized during the time represented by the majority of the groundwater, showing only small declines in agriculture (2%) and increases in urban (5%) area. Therefore, it is not surprising that our results show little difference in explanatory power between legacy LULC and current LULC. Legacy effects through groundwater delivery would be more apparent in 1) areas where groundwater travel times are long enough to “reach back further in time”, or 2) areas where LULC is currently undergoing rapid conversions within areas of shorter delivery times. Wetland processing Wetlands are a particularly important land cover type in biogeochemical cycles, greatly altering the chemistry of water from input to output (Howard-Williams 1985). However, there is no consensus about the most appropriate method for representing the impact of wetland processing within the landscape. Most studies have measured wetland extent as the proportional cover within an area, usually a watershed. Whigham et al. (1988) offer an early meta-analysis showing the variable nature of nitrogen and phosphorus removal along a gradient of increasing wetland extent. Since then, other methods for measuring wetlands have been used, such as proportional cover within a specified distance of the target ecosystem (Johnston et al. 1990, Weller et al. 1996, King et al. 2005). We add to this research by offering a comparison of three measures of wetland influence: watershed extent, area of river connected wetland, and length of river connected wetland. 103 Our results agree with previous studies that show a stronger relationship between phosphorus and riparian wetlands in comparison to wetlands located further from receiving bodies (Johnston et al. 1990, Weller et al. 1996). Furthermore, our analyses provide evidence that wetland area and wetland length are important in different ways for nutrient cycling. Models including the area of connected wetland had higher explanatory power for dissolved phosphorus (SRP) whereas models including the length of connected wetland had higher explanatory power for total phosphorus (TP, Figure 3.10). Given that dissolved forms move more freely through groundwater pathways, it makes sense that the entire area of a connected wetland will influence SRP. On the other hand, overland flow has been found to be the primary transport mechanism for particulate forms, such as TP (Banner et al. 2009, Hoffman et al. 2009b). Therefore, the length of riparian wetlands may act as a barrier reducing overland transport. This agrees with previous studies which show that overland flow is reduced considerably by wetlands (Weller et al. 1996) and that the spatial arrangement of riparian zones, particularly gaps between patches that function to reduce flow, is important for material flux (Weller et al. 1998). In contrast to phosphorus, our analyses of nitrogen show that all three measures of wetland presence had similar model explanatory power. Groundwater pathways allow for the transport of dissolved substances over long distances. Because nitrogen is transported in dissolved forms primarily through subsurface pathways (Peterjohn & Correll 1984, Walsh and Kunapo 2009), wetland proximity may be less important for nitrogen concentrations in lakes. Previous studies report that nitrogen removal in wetlands varies with hydrologic flushing, showing low removal in periods of high flow (Johnston et al. 1990, Jordan et al. 2003). Because our samples were taken in the spring, it is possible that frequent rain events flushed nitrogen through the system (Inamdar et al. 2009). 104 In addition to nutrients, which are commonly included in analyses of wetland processing, we compare the ability of wetland measures to explain patterns in a variety of other water chemistry variables. Our results show no appreciable difference in predictive capacity for models of reactive and conservative ions when wetlands are measured by extent, connected area or connected length. These results are in agreement with those from Johnston et al. (1990), who found that specific conductance, a surrogate for dissolved ions, was related to measures of both wetland extent and proximity. Riparian zones Riparian zone dynamics are another mechanism affecting the relationship between LULC and ecosystem responses. Riparian LULC has been shown to be important for both surface and subsurface transport of materials to streams and lakes (Peterjohn and Correll 1984, Gregory et al. 1991, Groffman et al. 2002, Groffman et al. 2004). We show that some lake water chemistry variables were predicted as well from either catchment scale LULC or riparian LULC, and in some cases, riparian LULC had higher explanatory power than catchment LULC. Specifically, the riparian only model of SRP showed greater explanatory power over the GW legacy model. Because SRP shows a stronger relationship with wetland area than wetland length, we believe that this is more an indication of wetland effects on SRP than riparian effects in general. On the other hand, the riparian only model of Si shows a riparian signal not present in the GW legacy models investigating wetlands. Recent work shows that biological processing of Si by plants is an important source of Si to freshwater systems (Derry et al. 2005, Struyf et al. 2009).Therefore, it is possible that the relationship between riparian LULC and Si we observe is due to biogeochemical cycling by riparian vegetation (Struyf et al. 2009). 105 Riparian zones are recognized primarily for reducing nonpoint source pollution (Gregory et al. 1991, Groffman et al. 2002, Craig et al. 2008). Results from our comparison of individual mechanisms (GW legacy, wetland processing, and riparian LULC) with models combining these mechanisms provides further evidence that riparian zones are important to relate LULC to aquatic ecosystem characteristics. We show that when mechanisms are combined to offer a more holistic view, the explanatory power of the models increased. Our study also provides a comparison across water chemistry variables not often included in studies relating LULC to aquatic ecosystem dynamics. In fact, we only found two other studies reporting results linking conservative ions to landscape features (Ryszkowski et al. 1999, Wayland et al. 2002). In a study of biogeochemical barriers within agriculture fields, termed shelterbelts, Ryszkowski et al. (1999) show the utility of using Ca and Mg as conservative tracers for comparison with concentrations of nutrients under trees that have been planted as windbreaks. Our results relating conservative ions to catchment and riparian LULC show that diffuse transport of these chemicals through groundwater pathways overwhelms a weaker relationship found with riparian zones. Future studies of landscape connections to aquatic ecosystems can benefit from the use of naturally occurring conservative ions as tracers. To date, there have been numerous methods proposed to incorporate riparian zones, that mainly focus on defining the most appropriate spatial extent of measurement (Baker et al. 2006, van Sickle and Johnson 2008). In contrast, our emphasis was to model the impact of the riparian zone (however it is defined) jointly with groundwater and wetland influences. Our results support the early sentiment of Gregory et al. (1991) who say: 106 “The importance of riparian zones far exceeds their minor proportion of the land base because of their prominent location within the landscape and the intricate linkages between terrestrial and aquatic ecosystems.” We built our models hierarchically by adding riparian processes to the model of GW legacies in an attempt to mathematically express the important role of riparian zones. Conclusions This study investigated the roles of groundwater pathways, wetland processing and riparian zone dynamics for relationships with lake water chemistry. We offer insights about how groundwater interacts with land transformation to create legacy effects. We combine groundwater flow dynamics with other mechanisms known to influence relationships between LULC and lake water chemistry, namely wetland and riparian zone processing, to offer a more complete model of land-water interfaces. Our results show that incorporating all of these mechanisms is important for modeling lake water chemistry and that looking at them individually can lead to misinterpretations (e.g., SRP with riparian versus wetland) and lower predictive power overall. Through the use of naturally occurring conservative tracers, we provide a basis for comparison against nutrient relationships to the landscape. By categorizing the chemistry variables by their key characteristics of solubility and reactivity, we are better equipped to explore other mechanisms that are important for the physical transport and biogeochemical transformations of these chemicals. 107 Table 3.1. Size and water chemistry characteristics of study lakes, including minimum, maximum, mean, and coefficient of variation. Phosphorus species include total phosphorus (TP), and soluble reactive phosphorus (SRP). Nitrogen species include total nitrogen (TN), nitrate (NO3), and ammonia (NH4). Reactive ions include calcium (Ca), silica (Si), and sulfate (SO4). Conservative ions include chloride (Cl), potassium (K), magnesium (Mg), and sodium (Na). Lake characteristic 2 Lake area (km ) Phosphorus TP (µg/L) SRP (µg/L) Nitrogen TN (mg/L) NO3 (mg/L) NH4 (µg/L) Reactive Ions Ca (mg/L) Si (mg/L) SO4 (mg/L) Conservative Ions Cl (mg/L) K (mg/L) Mg (mg/L) Na (mg/L) Minimum Maximum Mean CV 0.06 2.60 0.65 0.93 5.4 0.33 45.0 3.45 24.1 1.15 0.37 0.56 0.53 1.85 0.94 0.30 0.00 1.03 0.20 1.20 3 143 36 0.96 16 0.00 92 5.61 52 1.53 0.38 0.88 1.9 88.4 26.1 0.83 2 0.58 3.9 0.9 242 7.53 24.7 77.3 66 2.09 15.4 24.0 0.90 0.57 0.26 0.88 108 Figure 3.1. Map showing the Huron River Watershed within the state of Michigan, including detailed hydrography features along with outlines for the cities of Ann Arbor and Ypsilanti. Ann Arbor Ypsilanti 109 Figure 3.2. Map showing land use/cover (LULC) from 1938, 1955, 1968, 1978, and 1996 within the Huron River Watershed. 1938 1955 1968 1978 1996 110 Figure 3.3. Surficial geology zones and optimized hydraulic conductivity (K, m/d) values used in the final model. Glacial Outwash = glacial outwash sand and gravel and postglacial alluvium. Ice Contact = ice contact outwash sand and gravel. Percent of each type within the Huron River groundwater sourceshed (HRGW) upstream of Ypsilanti, and the mean percent for the study lake groundwater sourcesheds are shown. The HRGW is dominated by glacial outwash. The study lakes had similar geology as the HRGW, in general. Surficial geology type Glacial Outwash Glacial Till - Fine Glacial Till - Coarse Ice Contact End Moraine Glacial Till - Medium Water % within HRGW 39.9 3.2 1.2 5.8 31.4 16.8 1.6 K (m/d) 7.3 8.2 9.4 9.4 14.6 14.6 500 111 Mean % within lake groundwater sourcesheds 46.9 0.9 1.0 24.9 15.6 7.7 2.9 Figure 3.4. Observed versus predicted water levels (m) recorded for 15,581 wells used to estimate values for hydraulic conductivity (K). The line of best fit and coefficient of determination are shown. 370 350 R² = 0.9024 330 Predicted (m) 310 290 270 250 230 210 190 190 240 290 Observed (m) 112 340 Figure 3.5. Steps to create a land use legacy map. Step 1: The simulated groundwater travel times in years were calculated for each model location down gradient to a discharge point (e.g. river) using Darcy’s Law, as detailed in Methods. Step 2: Groundwater travel time year was reclassified to bracket times with available LULC data (in this case, 1996, 1978, 1968, 1955, 1938, or Pre-1900). Step 3: LULC type for each cell was assigned using time categories from step 2. Step 4: LULC from 1996, showing areas that differ with the legacy LULC with black cross-hatching. 113 Figure 3.6. Groundwater travel times reclassified to represent the 5 time steps of land cover data. Time classes span from midpoint to midpoint around a specified year. Percent of each groundwater travel time category within the Huron River groundwater sourceshed (HRGW), and the mean percent for the study lake groundwater sourcesheds are shown. The HRGW is dominated by groundwater needing less than 10 years to travel from it’s source to a surface water delivery pathway. The study lakes had similar groundwater travel times as the HRGW, in general. Groundwater travel time 0-9 years 9-23 years 23-34 years 34-41 years 41-96 years >96 years Assigned to time step 1996 1978 1968 1955 1938 Pre-1900 % within HRGW 53.1 22.5 8.8 3.6 9.6 2.4 Mean % within lake groundwater sourcesheds 62.1 19.1 6.9 2.9 7.1 1.8 114 Figure 3.7. Land use/cover within the model area showing A) 1996 and B) legacy land use/cover. Percent of each LULC type within the Huron River groundwater sourceshed (HRGW), and the mean percent for the study lake groundwater sourcesheds are shown. Maximums are shown in bold. The HRGW is dominated by urban LULC type in 1996, but agriculture dominates in the legacy map. The study lakes had similar LULC representation as the HRGW, in general. A) 1996 land use/cover B) Legacy land use/cover 1996 Urban Agriculture Open Forest Water Wetland % within HRGW 27.0 25.4 17.2 17.9 5.1 7.4 Legacy Mean % within lake groundwater sourcesheds 22.2 13.7 18.6 20.9 14.4 10.3 115 % within HRGW 22.8 29.7 16.3 19.1 5.1 7.0 Mean % within lake groundwater sourcesheds 18.9 16.3 18.5 22.2 13.9 10.1 Figure 3.8. Areas within the Huron River groundwater sourceshed (HRGW) with differences between 1996 LULC and legacy LULC. The model boundary is indicated with the bold black line. Study lake groundwater sourcesheds are indicated with dark grey lines and light grey shading. Red areas indicate where legacy land cover differs from 1996 land cover (10.5% of the area within the HRGW, 8.4% of the area within the study lake groundwater sourcesheds). 116 2 Figure 3.9. Coefficients of determination (R ) from regression models of total phosphorus (TP), soluble reactive phosphorus (SRP), total nitrogen (TN), nitrate (NO3), ammonia (NH4), silica (Si), calcium (Ca), sulfate (SO4), chloride (Cl), potassium (K), magnesium (Mg), and sodium (Na) using correlative legacy models from Chapter 2, LULC characterized by the current time step (current LULC), and LULC characterized by the legacy map. Asterisks indicate model chosen by AIC weights. 1.0 * 0.9 * * * 0.8 * 0.7 * 0.6 R2 * 2 R 0.5 * 0.4 0.3 * * * 0.2 Correlative Legacy * Current LULC GW Legacy 0.1 0.0 TP SRP TN NO3 NO3 NH4 NH4 Si 117 Ca SO4 SO4 Cl K Mg Na 2 Figure 3.10. Coefficients of determination (R ) from regression models of total phosphorus (TP), soluble reactive phosphorus (SRP), total nitrogen (TN), nitrate (NO3), ammonia (NH4), silica (Si), calcium (Ca), sulfate (SO4), chloride (Cl), potassium (K), magnesium (Mg), and sodium (Na) using LULC characterized by the GW legacy models measuring wetlands in three ways: all wetlands in the watershed (GW legacy), area of wetland connected to rivers, and length of wetland connected to rivers. Wetland area includes the entire area of the wetland polygons which connect to rivers. Wetland length differs from wetland area by only accounting for the area of the connected wetland polygons within a 50m buffer of rivers. 1.0 0.9 0.8 0.7 0.6 R2 2 R 0.5 0.4 0.3 GW Legacy Wetland Area Wetland Length 0.2 0.1 0.0 TP SRP TN NO3 NO3 NH4 NH4 Si Ca 118 SO4 SO4 Cl K Mg Na 2 Figure 3.11. Coefficients of determination (R ) from regression models of total phosphorus (TP), soluble reactive phosphorus (SRP), total nitrogen (TN), nitrate (NO3), ammonia (NH4), silica (Si), calcium (Ca), sulfate (SO4), chloride (Cl), potassium (K), magnesium (Mg), and sodium (Na) using correlative legacy models from chapter 2, LULC characterized by the legacy map, LULC within the 50m riparian zone buffer only, and legacy LULC combined with LULC from the 50m riparian zone buffer. Asterisks indicate model chosen by AIC weights. The asterisk for nitrate indicates the correlative legacy model. 1.0 * 0.9 * * 0.8 0.7 2 R2 R * * * * * * 0.6 0.5 * 0.4 0.3 * * Correlative Legacy 0.2 GW Legacy Riparian Buffer Only 0.1 NO3 NH4 NO3 NH4 SO4 GW Legacy + Riparian Buffer 0.0 TP SRP TN Si Ca 119 SO4 Cl K Mg Na REFERENCES 120 REFERENCES Baker, M.E., D.E. Weller, and T.E. Jordan. 2006. Improved methods for quantifying potential nutrient interception by riparian buffers. Landscape Ecology 21: 1327-1345. Banner, E.B.K., A.J. Stahl, and W.K. Dodds. 2009. Stream discharge and riparian land use influence in-stream concentrations and loads of phosphorus from Central Plains watersheds. Environmental Management 44:552-565. Boutt, D.F., D.W. Hyndman, B.C. Pijanowski, and D.T. Long. 2001. Identifying potential land use-derived solute sources to stream baseflow using ground water models and GIS. Ground Water 39(1): 24-34. Brazner, J. C., N. P. Danz, A. S. Trebitz, G. J. Niemi, R. R. Regal, T. P. Hollenhorst, G. E. Host, E. D. Reavie, T. N. Brown, J. M. Hanowski, C. A. Johnston, L. B. Johnson, R.W. Howe, And J. J. H. Ciborowski. 2007. Responsiveness of Great Lakes wetland indicators to human disturbances at multiple spatial scales: a multi-assemblage assessment. Journal of Great Lakes Research 33:42-66. Burcher, C.L. 2009. Using simplified watershed hydrology to define spatially explicit ‘zones of influence’. Hydrobiologia 618: 149-160. Chauvat, M., V. Wolters, and J. Dauber. 2007. Response of collembolan communities to landuse change and grassland succession. Ecography 30: 183-192. Craig, L.S., M.A. Palmer, D.C. Richardson, S. Filoso, E.S. Bernhardt, B.P. Bledsoe, M.W. Doyle, P.M. Groffman, B.A. Hassett, S.S. Kaushal, P.M. Mayer, S.M. Smith, and P.R. Wilcock. 2008. Stream restoration strategies for reducing river nitrogen loads. Frontiers in Ecology and the Environment 6(10): 529-538. Derry, L.A., A.C. Kurtz, K. Ziegler, and O. Chadwick. 2005. Biological control of terrestrial silica cycling and export fluxes to watersheds. Nature 433:728-731. Detenbeck, N.E., C.A. Johnston, and G.J. Niemi. 1993. Wetland effects on lake water quality in the Minneapolis/St. Paul metropolitan area. Landscape Ecology 8(1): 39-61. Doherty, J. 2004. PEST Model-Independent parameter estimation User’s Manual: 5th Edition. Watermark Numerical Computing, Australia. Available from: http://www.sspa.com/pest Floyd, W.C., S.H. Schoenholtz, S.M. Griffith, P.J. Wigington Jr., and J.J. Steiner. 2009. Nitratenitrogen, land use/land cover, and soil drainage associations at multiple spatial scales. Journal of Environmental Quality 38:1473-1482. 121 Foster, D.R., G. Motzkin, and B. Slater. 1998. Land-Use History as Long-Term Broad-Scale Disturbance: Regional Forest Dynamics in Central New England. Ecosystems 1(1): 96119. Foster, D.R., F. Swanson, J. Aber, I. Burke, N. Brokaw, D. Tilman, and A. Knapp. 2003. The importance of land-use legacies to ecology and conservation. BioScience 53(1): 77-88. Fraterrigo, J.M., and J.A. Downing. 2008. The Influence of land use on lake nutrients varies with watershed transport capacity. Ecosystems 11:1021-1034. Frimpong, E.A., T.M. Sutton, K.J. Lim, P.J. Hrodey, B.A. Engel, T.P. Simon, J.G. Lee, and D.C LeMaster. 2005. Determination of optimal riparian forest buffer dimensions for stream biota–landscape association models using multimetric and multivariate responses. Canadian Journal of Fisheries and Aquatic Sciences 62:1-6. Gergel, S.E., M.G. Turner, and T.K. Kratz. 1999. Dissolved organic carbon as an indicator of the scale of watershed influence on lakes and rivers. Ecological Applications 9(4): 13771390. Gregory S.V., Swanson F.J., McKee A. & Cummins K.W. (1991) Ecosystem perspectives of riparian zones. BioScience, 41, 540–551. Groffman, P.M., N.J. Boulware, W.C. Zipperer, R.V. Pouyat, L.E. Band, and M.F. Colosimo. 2002. Soil nitrogen cycle processes in urban riparian zones. Environmental Science and Technology 36: 4547-4552. Groffman, P.M., N.L. Law, K.T. Belt, L.E. Band, G.T. Fisher. 2004. Nitrogen fluxes and retention in urban watershed ecosystems. Ecosystems 7(4): 393-403. Harding, J.S., E.F. Benfield, P.V. Bolstad, G.S. Helfman, and E.B.D. Jones III. 1998. Stream biodiversity: The ghost of land use past. Proceedings of the National Academy of Sciences of the United State of America 95(25): 14843-14847. Hay-Chmielewski, E. M., P. W. Seelbach, G. E. Whelan, and D. B. Jester, Jr. 1995. Huron River assessment. Fisheries Special Report No 16. Michigan Department of Natural Resources. Hoffman, C.C., C. Kjaergaard, J. Uuski-Kamppa, H.C.B. Hansen, and B. Kronvang. 2009b. Phosphorus retention in riparian buffers: Review of their efficiency. Journal of Environmental Quality 38: 1942-1955. Hollister, J.W., P.V. August, and J.F. Paul. 2008. Effects of spatial extent on landscape structure and sediment metal concentration relationships in small estuarine systems of the United States’ Mid-Atlantic Coast. Landscape Ecology 23: 91-106. 122 Howard-Williams, C. 1985. Cycling and retention of nitrogen and phosphorus in wetlands: a theoretical and applied perspective. Freshwater Biology 15:391-431. Hunsaker, C.T., and D.A. Levine. 1995. Hierarchical approaches to the study of water quality in rivers. BioScience 45(3): 193-203. Inamdar, S., J. Rupp, and M. Mitchell. 2009. Groundwater flushing of solutes at wetland and hillslope positions during storm events in a small glaciated catchment in western New York, USA. Hydrological Processes 23:1912-1926. Johnston, C.A., N.E. Detenbeck, and G.J. Niemi. 1990. The cumulative effect of wetlands on stream water quality and quantity. A landscape approach. Biogeochemistry 10(2): 105141. Jordan, T.E., D.F. Whigham, K.H. Hofmockel, and M.A. Pittek. 2003. Nutrient and Sediment Removal by a Restored Wetland Receiving Agricultural Runoff. Journal of Environmental Quality 32:1534-1547. Kelly, V.R., G.M. Lovett, K.C. Weathers, S.E.G. Findlay, D.L. Strayer, D.J. Burns, and G.E. Likens. 2008. Long-term sodium chloride retention in a rural watershed: legacy effects of road salt on streamwater concentration. Environmental Science and Technology 42: 410-415. King, R.S., M.E. Baker, D.F. Whigham, D.E. Weller, T.E. Jordan, P.F. Kazyak, and M.K. Hurd. 2005. Spatial considerations for linking watershed land cover to ecological indicators in streams. Ecological Applications 15(1): 137-153. McTammany, M.E., E.F., Benfield, and J.R. Webster. 2007. Recovery of stream ecosystem metabolism from historical agriculture. Journal of the North American Benthological Society 26(3): 532-545. Moreno-Mateos, D., U. Mander, F.A. Comin, C. Pedrocchi, and E. Uuemaa. 2008. Relationships between landscape pattern, wetland characteristics, and water quality in agricultural catchments. Journal of Environmental Quality 37: 2170-2180. Osborne, L.L. and M.J. Wiley. 1988. Empirical relationships between land use/cover and stream water quality in an agricultural watershed. Journal of Environmental Management 26(1): 9-27. Peterjohn W.T., and D.L. Correll. 1984. Nutrient dynamics in an agricultural watershed: observations on the role of a riparian forest. Ecology 65:1466–1475. Pint, C.D., R.J. Hunt, and M.P. Anderson. 2003. Flowpath delineation and ground water age, Allequash Basin, Wisconsin. Ground Water 41(7): 895-902. 123 Pijanowski, B., D. K. Ray, A. D. Kendall, J. M. Duckles, and D. W. Hyndman. 2007. Using backcast landuse change and groundwater travel-time models to generate land-use legacy maps for watershed management. Ecology and Society 12(2): 25. Rutledge, D.T. 2001. Changes in land cover and wildlife habitats in two watersheds in the lower peninsula of Michigan. Ph.D. dissertation, Michigan State University, East Lansing, MI. Rutledge, D.T., and C.A. Lepczyk. 2002. Landscape change: patterns, effects, and implications for adaptive management of wildlife resources. In: Integrating Landscape Ecology into Natural Resource Management. J. Liu, and W.W. Taylor Ed. New York, NY: Cambridge University Press. p 312-333. Ryszkowski, L., A. Bartoszewicz, and A. Kedziora. 1999. Management of matter fluxes by biogeochemical barriers at the agricultural landscape level. Landscape Ecology 14: 479492. Soranno, P.A., S.L. Hubler, S.R. Carpenter, and R.C. Lathrop. 1996. Phosphorus loads to surface waters: a simple model to account for spatial pattern of land use. Ecological Applications 6(3): 865-878. Sponseller, R.A., E.F. Benfield, and H.M. Valett. 2001. Relationships between land use, spatial scale and stream macroinvertebrate communities. Freshwater Biology 46: 1409-1424. Strayer, D.L., R.E. Beighley, L.C. Thompson, S. Brooks, C. Nilsson, G. Pinay, and R.J. Naiman. 2003. Effects of land cover on stream ecosystems: roles of empirical models and scaling issues. Ecosystems 6: 407-423. Sturyf, E., W. Opdekamp, H. Backx, S. Jacobs, D.J. Conley, and P. Meire. 2009. Vegetation and proximity to the river control amorphous silica storage in a riparian wetland (Biebrza National Park, Poland). Biogeosciences 6: 623-631. Van Sickle, J., J. Baker, A. Herlihy, P. Bayley, S. Gregory, P. Haggerty, L. Ashkenas, and J. Li. 2004. Projecting the biological condition of streams under alternative scenarios of human land use. Ecological Applications 14(2): 368-380. Van Sickle, J., and C.B. Johnson. 2008. Parametric distance weighting of landscape influence on streams. Landscape Ecology 23: 427-438. Walsh, C.J., and J. Kunapo. 2009. The importance of upland flow paths in determining urban effects on stream ecosystems. Journal of the North American Benthological Society 28(4): 977-990. Wayland, K.G., D.W. Hyndman, D. Boutt, B.C. Pijanowski, and D.T. Long. 2002. Modelling the impact of historical land uses on surface-water quality using groundwater flow and solute-transport models. Lakes & Reservoirs: Research and Management 7: 189-199. 124 Wayland, K.G., D.T. Long, D.W. Hyndman, B.C. Pijanowski, S.M. Woodhams, and S.K. Haack. 2003. Identifying relationships between baseflow geochemistry and land use with synoptic sampling and R-mode factor analysis. Journal of Environmental Quality 32: 180-190. Weller, C.M., M.C. Watzin, and D. Wang. 1996. Role of wetlands in reducing phosphorus loading to surface water in eight watersheds in the Lake Champlain basin. Environmental Management 20(5): 731-739. Weller, D.E., T.E. Jordan, and D.L. Correll. 1998. Heuristic models for material discharge from landscapes with riparian buffers. Ecological Applications 8(4): 1156-1169. Wetzel R.G. and G.E. Likens. 2000. Limnological analyses, 3rd ed. Springer-Verlag. Whigham, D.F., C. Chitterling, and B. Palmer. 1988. Impacts of freshwater wetlands on water quality: A landscape perspective. Environmental Management 12(5): 663-671. Xenopoulos, M.A., D.M. Lodge, J. Frentress, T.A. Kreps, S.D. Bridgham, E. Grossman, and C.J. Jackson. 2003. Regional Comparisons of watershed determinants of dissolved organic carbon in temperate lakes from the Upper Great Lakes region and selected regions globally. Limnology and Oceanography 48(6): 2321-2334. Yuan, Y., R.L. Binger, and M.A. Locke. 2009. A review of effectiveness of vegetative buffers on sediment trapping in agricultural areas. Ecohydrology 2: 321-336. Zhang, X., X. Liu, M. Zhang, R.A. Dahlgren, and M. Eitzel. 2009. A review of vegetated buffers and a meta-analysis of their mitigation efficacy in reducing nonpoint source pollution. Journal of Environmental Quality 39:76-84. 125