3... u .2... ... I“ s. - £111.33)! . 6.5.55... .2... .. 2...: 3 u. ....3 ..gsv.n PI . -55.}. 23:31.5. a ...: 25â€.“ {£313.}. 3 £235... . :gaï¬â€˜ :3311. .51.»..2... v ’n...‘ , .... t ......k . HS... . . I .....zt...’ 5.45.“? . i3...- y ... :ur. .. an... a .. ....‘4. ... ...! . “Lorin. . lie? . .. . a... . .r. . .. ...... ... 1:. 13.1.... .. .. ivynï¬ï¬‚ngtfllinwrmfla .3...†. . ... . . 1...... . 33...... .. .4: c . «.2... .. t y 5.. .1... . .. 2 r: r 1.. LI! 1... . . s : v7.4... :3. r... . 5.1.3.5... . . ......3. . 2.303 .vktis ‘ S w “......n . ....w , Z .5 . r ......m...:zra..s 29$}; ... .1, ...... . .. .2... , 3.1.2}. , _.“ l4. .....19 at? d. , . 1L ï¬duflhflw «s. :91“... .3 2.5a.3§\u . RHâ€... ......r 25.! 35:5... . .31.: 1.. .. )5 .I . Auk}. ‘33. 3.5.}... 13...}... a ......4.....fi~...3i..3 2-. . x .. I. . 5.1)»... l 51.1.- ‘31.. I ...... . 1...... .. :3: ...; ......z ... 1. 3X5? , I. goo“! This is to certify that the thesis entitled Modeling current and historic habitat for Canada lynx (Lynx Canadensis) in the Upper Peninsula of Michigan presented by Daniel W. Linden has been accepted towards fulï¬llment of the requirements for the MS. degree in Fisheries and Wildlife â€â€”— l. 1â€â€œ Maj? Professlor's Signature n. I / 5/ Oi... l I Date MSU is an Aï¬â€˜innative ActiorVEqual Opportunity Institution LIBRARY Michigan State University r | PLACE IN RETURN BOX to remove this checkout from your record. To AVOID FINES return on or before date due. MAY BE RECALLED with earlier due date if requested. DATE DUE DATE DUE DATE DUE Mbfltï¬as 2/05 p:/C|RC/DateDue.indd-p.1 MODELING CURRENT AND HISTORIC HABITAT FOR CANADA LYNX (LYNX CAN ADEN SIS) IN THE UPPER PENINSULA OF MICHIGAN By Daniel W. Linden A THESIS Submitted to Michigan State University in partial ï¬Jlï¬llment of the requirements for the degree of MASTER OF SCIENCE Department of Fisheries and Wildlife 2006 ABSTRACT MODELING CURRENT AND HISTORIC HABITAT FOR CANADA LYNX (LYNX CANADENSIS) IN THE UPPER PENINSULA OF MICHIGAN By Daniel W. Linden In the ruling to list Canada lynx (Lynx canadensis) as a federally threatened species, the US Fish and Wildlife Service identiï¬ed the Great Lakes region as an area of concern. While there is no current evidence of a resident population in the Upper Peninsula (UP) of Michigan, trapping and track records over the past century suggest the area was periodically invaded by lynx during population eruptions in Canada. My objectives were to quantify past and present forest conditions in the UP for lynx habitat potential, and estimate changes in habitat connectivity between Canada and the UP. I used a spatially explicit, landscape-level habitat model that required multiple layers of spatial data compiled in a GIS to describe lynx habitat components. Forest inventory plots were used to provide detailed stand attribute information. The historical range of variability in presettlement forest conditions was simulated with the landscape age-class demographics simulator (LADS). Outputs from the model indicated that potential habitat has become more widespread under current conditions than that which existed during the presettlement era. Foraging quality has increased in many areas across the UP, but remains low throughout. Non-habitat interspersion is limiting under current conditions and has increased dramatically from that of presettlement in the eastern UP. This increase has resulted in lower habitat quality and decreased connectivity of habitat in the eastern UP. The establishment of a resident lynx population in the UP is questionable. ACKNOWLEDGEMENTS Funding for this project was provided by the Michigan Department of Natural Resources (MDNR), USDA Forest Service (USFS), Boise Cascade, LLC, the Michigan Agricultural Experiment Station and the Michigan State University, Department of Fisheries and Wildlife. This project was supported by the Federal Aid in Restoration Act under Pittman-Robertson project W-147-R. Special thanks to Dr. Pat Lederle and Dr. Dean Beyer of the MDNR for making the majority of funding possible, and to the Hiawatha and Ottawa National Forests for their contributions. I would like to thank my major advisor, Dr. Rique Campa, for providing daily support and encouragement throughout my graduate experience. His open-door policy and optimistic attitude fostered a great working atmosphere. I am also thankful for the collaborative insight provided by my committee members, Dr. Gary Roloff, Dr. Dean Beyer, Dr. Kelly Millenbah and Dr. Steve Friedman. Special thanks to Dr. Roloff for contributing invaluable knowledge and experience with the lynx model. Several people provided data and expertise that made aspects of this research possible. Mr. Geoff Holden, USF S, was instrumental in the integration of conï¬dential plot locations with my spatial data. Dr. David Cleland, USFS, provided a number of useful datasets without hesitation, as well as helpful information regarding forest ecology in the Great Lakes. Dr. Michael Wimberly, South Dakota State University, forwarded his Landscape Age-class Dynamics Simulator (LADS) along with some personal guidance on the execution of the program. I want to thank my numerous friends and colleagues in the Department of Fisheries and Wildlife for their support. Dr. Michael Jones and Dr. Michael Wilberg iii provided some custom programming that was very helpful. Ty Wagner was an indispensable source of statistical information, as well as a good friend. My lab mates, Ali Felix and Dan Walsh, created an enjoyable work environment for which I am truly grateï¬il. They enhanced my development as a student and scientist, and I thank them for their friendship. I would also like to thank my friends in the Bence/Jones lab, and elsewhere in the department, for the many good times over the past 3 years. Lastly, I wish to thank my family, especially my parents, Fred and Laurie Linden, for offering endless support and encouragement throughout my life. Words cannot express my gratitude for having grown up around such good people. iv TABLE OF CONTENTS LIST OF TABLES ............................................................................................................ vii LIST OF FIGURES ......................................................................................................... viii FOREWORD ...................................................................................................................... 1 CHAPTER 1: TEMPORAL CHANGES IN HABITAT POTENTIAL FOR CANADA LYNX IN THE UPPER PENINSULA OF MICHIGAN ........................ 5 INTRODUCTION .................................................................................................. 5 STUDY AREA ..................................................................................................... 10 METHODS ........................................................................................................... 13 Lynx Habitat Model ....................................................................................... l3 Foraging Component .............................................................................. l4 Denning Component ............................................................................... 19 Interspersion Component ........................................................................ 21 Final Habitat Potential ............................................................................ 22 Ecological Land Classiï¬cation for the Upper Peninsula ............................... 23 Current Ecological Conditions ............................................................... 29 Current vegetation ........................................................................... 3O Presettlement vegetation ................................................................. 32 Soils ................................................................................................. 35 Ecoregions ....................................................................................... 37 ELU grids ........................................................................................ 38 Forest plot surveys .......................................................................... 38 Linking plots to ELU grid ............................................................... 45 Satellite imagery .............................................................................. 52 Presettlement Ecological Conditions ...................................................... 54 Simulating forest age classes ........................................................... 55 Estimating habitat quality within age classes .................................. 58 Analysis of model outputs ............................................................................. 62 RESULTS ............................................................................................................. 63 Current Ecological Conditions ...................................................................... 63 Presettlement Ecological Conditions ............................................................. 72 Differences between Current and Presettlement conditions .......................... 77 DISCUSSION ....................................................................................................... 83 CHAPTER 2: TEMPORAL CHANGES IN HABITAT CONNECTIVITY FOR CANADA LYNX IN THE UPPER PENINSULA OF MICHIGAN ...................... 91 INTRODUCTION ................................................................................................ 9 1 STUDY AREA ..................................................................................................... 93 METHODS ........................................................................................................... 93 RESULTS ............................................................................................................. 94 DISCUSSION ....................................................................................................... 96 MANAGEMENT IMPLICATIONS .............................................................................. 102 LITERATURE CITED ................................................................................................... 103 APPENDIX A: Habitat suitability model for the North American lynx: Lakes States version ........................................................................................................ 112 APPENDIX B: Schematic of processes used for the Lynx Model ................................. 139 APPENDIX C: Accuracy assessment of IFMAP landcover ........................................... 140 APPENDIX D: Calculating horizontal cover with the Stand Visualization System ............................................................................................................................. 143 APPENDIX E: Landscape-scale modeling of forest structure using Landsat TM imagery and FIA data ................................................................................. 147 APPENDIX F: Metadata for spatial layers ..................................................................... 153 vi LIST OF TABLES Table 1. Habitat types originally described by Coffman et al. (1980) for the Upper Peninsula of Michigan with associated soil properties and “habtype†groupings. ................................................................................ 26 Table 2. Dominant tree species assemblages within seral stages of each habtype for the Upper Peninsula of Michigan. ......................................... 27 Table 3. Areas of land-cover classes identiï¬ed by the 2001 IF MAP in the Upper Peninsula of Michigan. .................................................................. 31 Table 4. Area of cover types in the presettlement vegetation layer from Comer et al. (1995) located in the Upper Peninsula of Michigan. ........... 33 Table 5. Classes within spatial layers that delineated habtypes and ELU categories in the Upper Peninsula of Michigan. ....................................... 39 Table 6. Forested classes from condensed ecoregion-ELU grid for the Upper Peninsula. ....................................................................................... 49 Table C. 1. Error matrix of forest classiï¬cations between IFMAP and FIA surveys in the Upper Peninsula. .............................................................. 142 vii LIST OF FIGURES Figure 1. Location of the Upper Peninsula (black outline) and the extent of forest cover (green) within the Great Lakes region. Images in this thesis are presented in color. ..................................................................... 11 Figure 2. Production function describing the relationship between hare habitat suitability index and stand attributes for: (a) horizontal cover (y = 103/[1+1400e'0" 1"] ); and (b) understory dominance (for x <60, y = 0.8333x + 50). ................................................................... 16 Figure 3. Viability relationship for snowshoe hares developed from previous studies (from Roloff 2003). The viability threshold occurs at 12 young/year (HSI = 60), and the marginal threshold occurs at 5 young/year (HSI = 25) .............................................................................. 18 Figure 4. Production functions relating lynx habitat quality to each of the three model components within a 250 ha home range: (a) foraging (for x <90, y = 1.1x); (b) denning (for x >400 and <1750, y = -0.0741x + 129.64); and (c) interspersion (for x >300 and <1500, y = -0.0833x - 25). .................................................................................... 20 Figure 5. Spatial distribution of ecological land units (ELUs) in the Upper Peninsula of Michigan with subsection boundaries from the Cleland et al. (2005b) map of ecoregions. ................................................ 41 Figure 6. Production function of hare habitat quality as predicted by stem cover units. Hare density was scaled to index habitat quality, and the regression model from Litvaitis et al. (1985) was used to determine the slope of the relationship with stem cover units (for y >23.04 and <55.65, y = 3.067x —— 70.667). ...................................... 46 Figure 7. Spatial distribution of climate zones created from groupings of similar biophysical units (Cleland et al. 2005a) within and between ecoregion sections, along with associated estimates of natural ï¬re rotations (N F R) in years ............................................................................ 57 Figure 8. Fluctuation in the proportion of early successional forest (10-40 yrs old) over a 10,000-yr simulation of severe ï¬res in the Upper Peninsula of Michigan, with arrows identifying highest and lowest years. Fire was simulated with the Landscape Age-class Dynamics Simulator under the natural disturbance regime that occurred prior to European settlement .............................................................................. 61 viii Figure 9. Figure 10. Figure 11. Figure 12. Figure 13. Figure 14. Figure 15. Figure 16. Figure 17. Figure 18. Figure 19. Frequency distributions of quality scores for hare habitat quality and lynx foraging quality across the study area, as assessed by HSI models measuring horizontal cover with manipulated, HSI-l (a, d) and original, HSI-2 (b, e) seedling distributions, and the H81 model measuring stem cover units, HSI-3 (c, f). ................................................. 64 Hare habitat quality in the Upper Peninsula of Michigan between 2000—2003. The inset illustrates a histogram of the proportion of quality scores across the study area. ......................................................... 66 Lynx foraging quality in the Upper Peninsula of Michigan between 2000-2003. The inset illustrates a histogram of the proportion of quality scores across the study area. ......................................................... 67 Lynx denning quality in the Upper Peninsula of Michigan between 2000-2003. The inset illustrates a histogram of the proportion of quality scores across the study area. ......................................................... 69 Lynx interspersion quality in the Upper Peninsula of Michigan between 2000-2003. The inset illustrates a histogram of the proportion of quality scores across the study area. ................................... 70 Overall lynx habitat quality in the Upper Peninsula of Michigan between 2000-2003. The inset illustrates a histogram of the proportion of quality scores across the study area. ................................... 71 Locations of viable lynx home ranges according to habitat quality in the Upper Peninsula of Michigan between 2000-2003. Habitat units were aggregated into simulated home ranges using Plume’s (2005) G18 function. ................................................................................. 73 Hare habitat quality during the 2 simulation years of presettlement ï¬re disturbance with the highest (top) and lowest (bottom) amount of early successional vegetation. In the top map, a conflagration responsible for creating a large area of early successional vegetation can be seen in the eastern UP. ................................................. 74 Lynx foraging quality during the 2 simulation years of presettlement ï¬re disturbance with the highest (top) and lowest (bottom) amount of early successional vegetation .................................... 75 Lynx denning quality during the 2 simulation years of presettlement ï¬re disturbance with the highest (top) and lowest (bottom) amount of early successional vegetation .................................... 76 Lynx interspersion quality during the 2 simulation years of presettlement ï¬re disturbance with the highest (top) and lowest (bottom) amount of early successional vegetation .................................... 78 ix Figure 20. Figure 21. Figure 22. Figure 23. Figure 24. Figure 25. Figure 26. Figure 27. Figure D] Figure D.2. Figure E.1. Figure E.2. Overall lynx habitat quality during the 2 simulation years of presettlement ï¬re disturbance with the highest (top) and lowest (bottom) amount of early successional vegetation .................................... 79 Locations of viable lynx home ranges according to habitat quality for the 2 simulation years of presettlement ï¬re disturbance with the highest (top) and lowest (bottom) amount of early successional vegetation. Habitat units were aggregated into simulated home ranges using Plume’s (2005) G18 function. .............................................. 80 Number of viable home ranges simulated in each of the 3 lynx habitat quality grids, using the Plume (2001) GIS function to aggregate habitat units. Viable home ranges were deï¬ned as having an average quality >70. ................................................................. 82 Relative frequency of bobcat harvests in the Upper Peninsula of Michigan between 1985-2005. The dotted line represents the boundary below which annual snowfall is <270 cm/year ......................... 90 Location of sources (black dots) and destination (white star) for calculating the least-accumulative-cost corridor for lynx in the Upper Peninsula of Michigan. .................................................................. 95 Least-accumulative—cost values for the eastern and western corridors created for each cost surface in the Upper Peninsula of Michigan. .................................................................................................. 97 Relative travel costs for eastern corridors simulated in each cost- surface map of the Upper Peninsula of Michigan, for the high (A) and low (B) presettlement conditions and present conditions (C). ........... 98 Relative travel costs for western corridors simulated in each cost- surface map of the Upper Peninsula of Michigan, for the high (A) and low (B) presettlement conditions and present conditions (C). ........... 99 SVS drawings of a mature red pine stand with hardwood understory. .............................................................................................. 143 Proï¬le view of stand with a red box around each height stratum. .......... 146 Location of the 8 scenes of Landsat 7 imagery acquired for the Upper Peninsula, including paths 21-25 (east to west) and rows 27-29 (north to south). ............................................................................ 148 Relationship between predicted and observed values of arcsine transformed horizontal cover for each of the 5 spectral images in the Upper Peninsula. ............................................................................... 151 FOREWORD The US Fish and Wildlife Service (USFWS) made their ï¬nal ruling on 24 March 2000 to list Canada lynx (Lynx canadensis) as a threatened species in accordance with the Endangered Species Act of 1973, following an investigation of the status of lynx populations in the contiguous United States (USFWS 2000). They determined that some current land management practices had the potential to negatively affect lynx and/or lynx habitat in the absence of adequate protection for the species. The ruling distinguished four regions (the Northeast, Great Lakes, Northern Rockies & Cascades, and Southern Rockies) collectively as one “Distinct Population Segment†(USFWS 2000) based on the criteria that lynx populations in those regions were discrete and signiï¬cant from populations in northern Canada. The separation by international political boundaries made them discrete, mostly through disparities in management policies between the US and Canadian governments, while differences in vegetation types, climate, and ecology made them signiï¬cant. In light of the ï¬nal ruling, government agencies, such as the USDA Forest Service, are faced with developing and implementing management strategies that facilitate lynx populations on public lands. Attaining this management goal has been problematic because an understanding of current and historical lynx distributions and habitat use patterns in the southern regions of their range, such as the Great Lakes, is largely incomplete. Lynx are commonly associated with the boreal forest region of northern Canada, also known as the taiga (McCord and Cordoza 1982, McKelvey et al. 2000a). Their distribution in the contiguous United States is limited to southern boreal forests, consisting of subalpine conifers in the West, and mixed coniferous-deciduous forests in the East (Aubry et al. 2000). These southern boreal forests have been increasingly fragmented by human disturbance, and lynx habitat patches tend to be smaller than those in the north. Some habitat patches in the southern boreal forests are able to support resident lynx populations (McKelvey et al. 2000b), but it is uncertain how source-sink dynamics among subpopulations have been affected by human-induced fragmentation, or how this impacts lynx viability in these areas. In other parts of the southern boreal forest, at the limits of the species’ range, lynx presence is attributed “solely†to the dispersal of individuals from core populations in Canada (USFWS 2000). Habitat quality in these regions is, therefore, assumed to be poor; the temporal and spatial extents to which these regions have acted as sinks to the core populations in Canada are undeï¬ned. This complicates the formulation of management policies, considering that the overall importance of these regions (e. g. the Great Lakes) to the species is questionable. The scant amount of information concerning lynx in the Great Lakes region consists of occurrence data, as documented by trapping and/or track records from the past 150+ years (Aubry et al. 2000, Beyer et al. 2001). These records appear to illustrate the historical rarity of lynx in Michigan, though they do not provide any indication of individual ï¬tness or population viability. Also, abundance (i.e., number of records per geographic area) can be a misleading measure of habitat quality in the absence of ancillary ecological information (Van Home 1983), so the occurrence data alone are not useful for determining whether there were viable populations of lynx historically in the state. Wood and Dice (1924) warned that some of the historical lynx records from Michigan were unreliable given the potential misidentiï¬cation of bobcat (Lynx rufus) for lynx, but several veriï¬ed records dated between 1842-1912 do exist for the Upper and Lower Peninsulas. The lynx was considered “extremely rare, if not already extinct†in the ninth Biennial Report (1937-38) of the Michigan Department of Conservation, and only 6 records of lynx occurred between 1904-1939 for the Upper Peninsula (UP), with a lone record in 1917 for the Lower Peninsula (McKelvey et al. 2000a, Beyer et al. 2001). Burt (1946) illustrated a map showing “questionable†records of lynx across most of Michigan and as far south as Monroe County, but reiterated the belief that the species was most likely extirpated from the state, except for a possible remnant on Isle Royale. Harger (1965:152) believed that lynx were “deï¬nitely making a comeback†in the UP due to the influx of trappings in 1962, but McKelvey et al. (2000a) pointed to an irruption of lynx populations in Canada as the cause for this increase. McKelvey et al. (2000a) compared lynx occurrence data from Michigan, Wisconsin and Minnesota with trapping records from Ontario, Manitoba, and Saskatchewan spanning 50+ years, and found a strong correlation when the Canadian data were lagged 3 years. They also examined the potential correlation between lynx occurrence/harvest data and snowshoe hare (Lepus americanus) harvest data in Minnesota and found no relationship, thus supporting the hypothesis that migration was the primary force influencing lynx populations in the Great Lakes region. Extensive track surveys by the Michigan Department of Natural Resources (DNR) in the 1990s found no sign of lynx in the UP, prompting Beyer et al. (2001) to conclude there was no evidence of a resident population. Whether or not a resident population exists in Michigan, it is possible that individuals dispersing from Canada during a population increase may enter the UP in search of resources (McKelvey et al. 2000a). Thus, the identiï¬cation of these resources and their location within the landscape becomes an important task for state and federal agencies in Michigan. The objectives of this research project were to quantify current and historic (i.e., prior to European settlement) habitat potential for Canada lynx in the UP of Michigan, analyze temporal trends in the amount and distribution of habitat in the region, and to examine the dynamics of potential migrations from Canada into Michigan. The process involved formulating spatial data in a Geographic Information System (GIS) that could describe the landscape, from both time periods, in terms of the forest attributes that deï¬ne the current understanding of lynx habitat requirements. The project also involved simulating the dispersal routes of individuals from Canada into Michigan, given the obstacles that may be perceived in the current landscape. Major shifts in forest dynamics across the UP due to anthropogenic influences (e. g., human development, harvesting, ï¬re suppression) have most likely affected the amount and distribution of lynx habitat since presettlement, and the increased fragmentation due to roads and development in the Great Lakes region on a whole has affected the ability of lynx to disperse into and throughout the area. The results of these analyses may guide the formulation of management practices and policies for lynx conservation in the UP, and the methods used can be applied other regions in the contiguous United States. CHAPTER 1: TEMPORAL CHANGES IN HABITAT POTENTIAL FOR CANADA LYNX IN THE UPPER PENINSULA OF MICHIGAN INTRODUCTION The listing of Canada lynx (Lynx canadensis) as a threatened species in the contiguous United States by the US Fish and Wildlife Service (USFWS) in March of 2000 necessitated the implementation of management policies that facilitate the persistence of lynx on public lands (USFWS 2000). This is complicated by an incomplete understanding of lynx ecology in the southern boreal forests of the United States, which represent the southern periphery of the species’ geographic range (McCord and Cordoza 1982). In the Great Lakes, the only information regarding the species comes from historical trapping and occurrence data (Aubry et al. 2000), and the scarcity of veriï¬ed records suggests the species was traditionally rare throughout Michigan, with nearly all occurrences in the Upper Peninsula (UP). Currently, there is no evidence of a resident population in the UP (Beyer et al. 2001), though transient lynx from neighboring states and Canada may disperse into the area in search of resources (McKelvey et al. 2000a). Accordingly, the USFWS called for “accurate mapping of lynx habitat in the Great Lakes Region†(USFWS 2000: 16057) to provide the location and distribution of resources necessary for the species’ persistence and identify potential areas for conservation and management. One resource that is critical for determining lynx survival is its primary prey, the snowshoe hare (Lepus americanus). The close relationship between fluctuations of lynx and hare populations in the boreal forest ecosystem is well documented, and patterns of habitat use exhibited by lynx are assumed to be highly correlated with those of hare (Keith 1963, Nellis et al. 1972, Brand et al. 1976, McCord and Cordoza 1982, Koehler and Aubry 1994). Lynx are known to prey upon squirrels and other small mammals, as well as ungulates and grouse during the summer and periods of low hare density (see diet summary in Koehler and Aubry 1994), but the availability of hare seems to limit the persistence of lynx, regardless of alternative prey abundance (Ward and Krebs 1985, Mowat et a1. 2000). There is some debate among biologists as to whether or not southern populations of snowshoe hare exhibit the same predator-prey cycle found in the north (see review in Hodges 2000). In the southern boreal forests, habitat conditions are thought to be suboptimal in comparison to the northern boreal forest, and may not allow hare populations to reach the high densities found in the north (Koehler 1990). Marginal forest conditions and the presence of other predators (e.g., coyote [Canis lupus], bobcat [Lynx ruï¬ls]) not abundant in the north could be keeping hare populations at a stable low density (Wolff 1981, Koehler and Aubry 1994), resulting in a cycle with reduced amplitude or no cycle at all (Hodges 2000). Wolff (1980, 1981) hypothesized that hare populations in the south were unable to successfully disperse among habitat patches due to the higher densities of predators, resulting in the absence of a cycle. In the northern boreal forest, when the pressures of predation and competition on hares are both low, individuals are able to disperse from optimal habitat and persist in suboptimal areas as populations increase. The stability of predator populations in the southern boreal forest, along with increased herbivore competition, preclude hare populations from reaching the high densities needed to initiate the predator-prey cycle (Wolff 1980, 1981). This “refugium†model (Wolff 1980, 1981) has been tested in other studies, which have found that habitat patches in some southern regions are too small (i.e., <10 ha) to provide signiï¬cant refugia for hare, resulting in similarly high predation rates between optimal and suboptimal habitat conditions (Keith et al. 1993, Wirsing et al. 2002). Thus, the fragmentation of habitat patches decreases the ability of hares to successfully disperse, as well as increasing the vulnerability of philopatric individuals to predation. Patches of sufï¬cient size (e.g., ~50 ha [Litvaitis et al. 1985]) are necessary for hare densities in southern boreal forests to rival those in the north (Wirsing et al. 2002). One distinct difference between northern and southern boreal forests is the relatively high frequency and extent of ï¬re disturbance in the north, which can create widespread areas of early successional vegetation types important to hares, interspersed with a mosaic of mature forest patches (Fox 197 8, Keith et al. 1993, Agee 2000). Snowshoe hares select habitat based on the density of understory vegetation providing security cover and winter browse (Brocke 1975, Wolff 1980, Litvaitis et al. 1985, Parker 1986, Hodges 2000), which is generally most abundant in forests during early seral stages. The combination of ï¬re suppression and exclusion in the contiguous United States, along with naturally longer ï¬re return intervals across boreal regions of the Northeast and Intermountain West, results in less frequent disturbances of a lower intensity than those found in the north. Koehler (1990:849) explained that the early successional patches he found in Washington were “small, isolated islands†created by wind-throw and lightning, typically <1 ha, and unlikely to support sufï¬cient numbers of hares. Fire disturbance is among the mechanisms believed to drive the lynx-hare cycle in northern Canada because of its high frequency and widespread impact, resulting in the constant re—initiation of secondary succession (and consequently, abundant cover) across large tracts of forest (Fox 1978). In northern Michigan, the current lengths of ï¬re rotations across the landscape are estimated to have increased tenfold over those that existed prior to European settlement (Cleland et al. 2004); this has undoubtedly altered the dynamic interactions between forest disturbances and hare populations in the region. While it is apparent that hare populations require an adequate amount of early successional vegetation types to persist (Hodges 2000), lynx have been shown to beneï¬t from an interspersion of mature forest (O’Donoghue et al. 1998, Mowat et al. 2000). Lynx are unable to effectively hunt in vegetation that is too dense for them to traverse, making early seral stands a refuge for snowshoe hare (Wolff 1980); an interspersion of mature forest would create edges that are navigable for lynx to hunt. Mature forest stands may also provide necessary denning conditions for lynx, as the most common characteristic found to be an indicator of denning habitat is the amount of downed woody debris (Koehler 1990, Koehler and Brittell 1990, Slough 1999, see review in Mowat et al. 2000), which is frequently scattered throughout older mesic forests. An old forest dominated by conifers with a sparse canopy late in succession could potentially contain a dense understory and an adequate array of woody debris, thus providing both foraging and denning opportunities for lynx (Buskirk et al. 2000). The absence of quantitative studies on lynx-habitat relationships in the UP would require the use of theoretical habitat models to evaluate the landscape for lynx. A framework proposed by Roloff and Haufler (1997) integrates the use of habitat suitability index (HSI) modeling (USFWS 1981) with a Geographic Information System (GIS) to allow a landscape-level analysis of habitat at multiple spatial scales that could address the viability of an individual organism, sub-population, or whole population in a given area (Roloff and Haufler 1997). This framework was similar in concept to habitat evaluation procedures (USFWS 1980), and was demonstrated using Canada lynx as the focal species (Roloff and Haufler 1997). The process involved quantifying habitat suitability at the home-range scale across a grid map of the landscape using HSI models, which assume that a species will use a certain area more frequently when it contains the necessary life requisites (USFWS 1981). The resulting grid map is used to simulate the spatial distribution of hypothetical lynx home ranges, based on quality and quantity thresholds of habitat suitability described by previous studies, which then provides a measure of habitat potential across the landscape (Roloff and Haufler 1997). An important step in the framework is the validation of model assumptions with observational data. Nylen- Nemetchek (1999) used this procedure with a lynx HSI model (Roloff 2001) to assess lynx habitat potential within the Riding Mountain National Park in Manitoba and found a signiï¬cant positive correlation between model outputs and surveyed lynx presence. Our study objectives were to determine the amount and distribution of potential Canada lynx habitat in the UP of Michigan, and examine differences between the current and historic (i.e., prior to European settlement) ranges of variability in forest conditions and subsequent habitat potential. Habitat potential for each time period was quantiï¬ed using the modeling framework proposed by Roloff and Haufler (1997) with the lynx HSI model developed by Roloff (2001). Results of this analysis could help guide natural resource managers with planning objectives for lynx conservation and management in the Great Lakes region, as well as providing a context for inferences about historic lynx distributions in the UP. STUDY AREA The UP of Michigan contains approximately 42,610 km’, located in the northern portion of the Great Lakes region. It is bounded by Lake Superior to the north, Lakes Michigan and Huron to the south, and the state of Wisconsin to the west (Figure 1). Over 80% of the land area is forested (Leatherberry and Spencer 1996), and the region falls within the Laurentian Mixed Forest Province, representing a transition from broadleaf deciduous to boreal coniferous forests (Albert 1995). The climate is humid continental, with signiï¬cant influence from Lake Superior on the north side. Winters are moderately long and cold, with mean temperatures of -9°C along the shorelines and -13°C inland; summers are short and cool, with mean temperatures of 19°C along the shorelines and 17°C inland (Eichenlaub et al. 1990). The length of growing season ranges from 180 days along the shorelines to 130 days inland (Eichenlaub et al. 1990). Annual precipitation averages 76-96 cm depending on location and is evenly distributed throughout the year. Snowfall occurs most regularly between November and April, and average seasonal totals range from 150 cm across the southern UP to >500 cm along the Keweenaw Peninsula and northern portions of the UP receiving lake effect snow off Lake Superior (Eichenlaub et al. 1990). The physiography of the UP is deï¬ned by numerous glacial landforms (e.g., outwash and lacustrine plains, ground and end moraines) that were created during the last retreat of glaciers >10,000 years before present (Dickman and Leefers 2003). A boundary of Paleozoic and Precambrian bedrock divides the UP roughly in half: the eastern section is characterized by a lower elevation (155-195 m) and relatively flat topography with sandy lake plains that are poorly drained, while the western section 10 -94° -92° -90° -88° -86° -84° I 48°- -48" 46°- —46° 44°- ‘ - 44° -94° ~92° -90° -88° -86° -84° Figure 1. Location of the Upper Peninsula (black outline) and the extent of forest cover (green) within the Great Lakes region. Images in this thesis are presented in color. (at elevations ranging 184-604 m) contains varying depths of glacial drift forming ground and end moraines over igneous and metamorphic bedrock, with exposed bedrock ridges along the coast of Lake Superior (Albert 1995). Each section encompasses a variety of smaller ecoregions deï¬ned by similar landforms, soils and climate regimes; these factors interact to influence natural disturbance and successional patterns among forest communities in the UP (Albert 1995, Frelich 2002). Human influences after European settlement have altered the structure and composition of forests in the UP. Extensive logging and catastrophic ï¬res during the late 19th and early 20th centuries have resulted in second growth forest throughout most of the region, including an increased presence of aspen (Populus spp.) (Dickman and Leefers 2003). Species such as eastern hemlock (T suga canadensis) and eastern white pine (Pinus strobus) have a minor presence in the landscape compared to that of the presettlement era due to heavy selective logging (Dickman and Leefers 2003), and ï¬re suppression and exclusion over the past century have favored an increased dominance of shade tolerant hardwoods (e. g., maple [A cer spp.]) (Zhang et al. 2000). Despite these changes, the species compositions of forests in the UP remain heavily dependent on physiographic conditions. Hydric lowlands are populated by conifers such as (in order of predominance) northern white-cedar (T huja occidentalis), balsam ï¬r (Abies balsamea), black spruce (Picea mariana), tamarack (Larix laricina) and hemlock, with hardwoods including black ash (F raxinus nigra), red maple (A cer rubrum) and balsam poplar (Populus balsamifera). Mesic uplands are populated by hardwoods such as sugar maple (Acer saccharum), aspen (Populus tremuloides, Populus grandidentata), red maple, yellow birch (Betula 12 alleghaniensis), paper birch (Betula papyrifera), beech (F agus grandifolia), basswood ( T ilia americana), and white ash (F raxinus americana) with conifers including ï¬r, cedar, hemlock, and white spruce (Picea glauca). Xeric uplands are populated by conifers such as jack pine (Pinus banksiana), red pine (Pinus resinosa), ï¬r, white pine and spruce, with hardwoods such as red maple, aspen, and northern red oak (Quercus rubra). Land use in the UP is primarily recreation and timber harvest. Out of the 33,900 km2 of forested land, 40% belongs to state and federal governments, including two national forests (the Ottawa and Hiawatha), and 18% is owned by timber industry (Leatherberry and Spencer 1996). There are approximately 328,000 residents; over 75% of the land area has a population density <4 persons/kmz, and 95% has a population density <15 persons/km2 (US. Census Bureau 2000). Wildlife species that are relevant to lynx ecology in the UP include coyote, bobcat, gray wolf (Canis lupus), and ï¬sher (Martes pennanti), as well as snowshoe hare, white-tailed deer (Odocoileus virginianus) and moose (A Ices alces). METHODS LYNX HABITAT MODEL Roloff’s (2001) lynx habitat model uses a limiting factor approach to characterize three components (life history traits — foraging, denning, non-habitat interspersion) demonstrated as necessary for lynx survival. The hypothesized requirements of foraging, denning, and interspersion were related to those suggested by Koehler and Aubry (1994). The input for the model consists of an ecological land classiï¬cation, in the form of a grid- based GIS, that can adequately stratify the variability in forest vegetation attributes relating to each habitat component (Roloff 2001). These habitat components are assessed l3 within the allometric home range for lynx (250 ha), a theoretical scale at which habitat selection is hypothesized to be strongest (Roloff and Haufler 1997). Several calibrations were necessary to account for differences in biogeoclimatic conditions between the Great Lakes and the Intermountain West, for which the model was originally developed (Roloff 2001). A similar approach was taken by Nylen- Nemetchek (1999) in Manitoba. Roloff (2003) (Appendix A) incorporated these changes into a new model speciï¬cally designed for the Great Lakes region, based on information obtained from the literature and recommendations from local forest and wildlife biologists. The changes involved reductions in the magnitude of variables related to snow accumulation, topographic influence and forest stature (Roloff 2003). The indices of habitat quality and quantity that were developed for each component in the original validated model were unchanged in the new model. All references to the lynx model, hereafter, will be to this new model developed in 2003. A schematic of the modeling process is illustrated in Appendix B. Foraging Component The availability of prey in the winter has been shown to affect the size of a lynx home range, and the degree to which home ranges overlap (McCord and Cardoza 1982, Ward and Krebs 1985); therefore, the foraging component of the model focuses on ecological conditions during the winter season. The foraging component consists of a snowshoe hare sub-model that assesses hare habitat quality within vegetation types based on the availability of palatable browse and winter security cover in the understory (Roloff 2003). For this study, I assumed that horizontal understory cover was the most limiting factor, based on evidence from the literature regarding the demonstrated importance of 14 horizontal cover for hare (Brocke 1975, Wolff 1980, Litvaitis et al. 1985, Parker 1986, Hodges 2000) and the high variability in deï¬nitions of palatable browse (see review in Hodges 2000). Therefore, hare habitat quality was indexed by 2 measures: (1) horizontal cover, deï¬ned by the percentage of vegetation cover measured perpendicular from the ground within two vertical height strata, 0-1 m and 1-2 m (accounting for variable snow depths throughout winter) at a standing distance of 10 m; and (2) understory dominance, deï¬ned by the percentage of coniferous species among all trees with a height to crown <3 m. Habitat quality for each variable (portrayed by HSI score on a scale of 0-100) is determined by a production function (Figure 2); quality increases logistically with the percentage of horizontal cover (Figure 2a), and linearly with the percentage of conifers (Figure 2b). These production functions were based on relationships demonstrated in the literature concerning perceived thresholds of horizontal cover, below which habitat is not provided for hare (Brocke 1975, Wolfe et al. 1982, Parker 1986, F erron and Ouellet 1992), and the additional importance of a conifer component in the understory for providing thermal cover during winter (Buehler and Keith 1982, Orr and Dodds 1982, Monthey 1986). A total HSI score for horizontal cover is calculated by combining the H81 scores for each height stratum with an arithmetic mean, which allows the horizontal cover in a particular height stratum to contribute to habitat quality regardless of the cover availability in the other height stratum. The overall hare habitat quality is calculated as the geometric mean of H81 scores for horizontal cover and understory dominance. This model structure insures that a vegetation type will not qualify as suitable hare habitat if horizontal cover is inadequate. 15 (a) Horizontal Cover O S — O _ oo 0 _. \O :73 z 9, 4 o __ N o ——< T l l T l l 20 40 6O 80 ' 100 Cover(%) (b) Understory Dominance O S ——-4 o L 00 G __ \O a :I: O _. V O _ N o _ l l l l T l 0 20 40 60 80 100 Conifer (9%) Figure 2. Production function describing the relationship between hare habitat suitability index and stand attributes for: (a) horizontal cover (y = 103/[1+1400e'0'l 1x] ); and (b) understory dominance (for x <60, y = 0.8333x + 50). The snowshoe hare sub-model produces a grid map of hare habitat quality with values ranging ï¬om 0-100, describing non-suitable to Optimal conditions, respectively. Habitat units (= quality X quantity) are preferentially aggregated by a raster-based GIS function (Plume 2005) into hare home ranges of differing size. The function places seeds randomly on the landscape, and habitat units with the highest value are selected within the search window until a habitat unit goal is achieved (Plume 2005). The habitat unit goal is based on the allometric home range for hare (4.5 ha) (Boutin et al. 1986), which describes the minimum area required by a mammal according to its mass. A home range is created for each seed as long as the habitat unit goal can be achieved before the home range reaches a maximum area (20 ha for hares). As the quality of hare habitat decreases, the number of pixels required to meet the habitat unit goal, and thus, the size of the simulated hare home range, increases. For example, one hypothetical home range would have a total area equal to the minimum of 4.5 ha if the average quality of the contributing habitat units was 100, while another would be 6.0 ha if the average quality was 75. Home ranges are considered viable (i.e., contributing to population ï¬tness) at an average quality 260 and marginal (i.e., not contributing to population ï¬tness) at an average quality 225 and <60. Home ranges with an average quality <25 are deemed non-viable and precluded from further analysis, under the assumption that these areas do not provide habitat for hares of a high enough quality to survive and reproduce. The thresholds are based on the inverse relationship between reproductive output and home range size of hares observed in previous studies (Behrend 1962, Dolbeer and Clark 1975, Cary and Keith 1979, Sievert and Keith 1985) (Figure 3), and the premise of these thresholds relates to forage availability for lynx (Roloff 2003). Areas with viable hare home ranges l7 Viability Relationship for Snowshoe Hares 0 Behrend(l962) 2 _ E Sievert and Keith (1985) o O _ g I—‘ g Dolbeer and Clark ( 1975) E o :l: Cary and Keith (1979) m _ o —J l l l l l 0 5 10 15 20 Offspring/Year l l l l l l 0 20 40 HSI 60 80 100 Figure 3. Viability relationship for snowshoe hares developed from previous studies (from Roloff 2003). The viability threshold occurs at 12 young/year (HSI = 60), and the marginal threshold occurs at 5 young/year (HSI = 25) would potentially support higher reproduction of hares, and thus, contribute to lynx foraging more than areas with marginal hare home ranges. Areas with non-viable hare home ranges are presumably void of snowshoe hares that can survive and reproduce. A production function describes the relationship between lynx foraging quality and the number of hypothetical hare home ranges within an allometric lynx home range (Figure 4a). The function is based on the estimate that lynx need to consume 1 hare every 2 days during the winter (Brand et al. 1976, as cited in Roloff 2003), and the winter is assumed to last approximately 180 days in Michigan (Eichenlaub et al. 1990). Therefore, lynx foraging quality is highest at 290 hare home ranges and lowest at 0 hare home ranges, with a linear trend between. The grid map is generated by counting the hare home ranges within a 250 ha moving window around each pixel. Hare home ranges that are viable contribute twice as much as those that are marginal, to account for the increased reproduction and survival of hares in higher quality habitats. The resulting grid map describes lynx foraging availability at the scale of a lynx home range. Denning Component The denning component deï¬nes quality based on thresholds of forest structure and the degree of juxtaposition with summer foraging habitat for lynx (Roloff 2003). In the Great Lakes region, vegetation types are assumed to provide potential denning sites when they contain 211.49 rn2 ha'1 of trees with an average DBH of 20 cm, 250% canopy closure, and mesic soil conditions (Roloff 2003). These values are intended to represent mature, productive forest conditions, which are hypothesized to provide the necessary requirements for lynx denning habitat (Koehler 1990, Koehler and Aubry 1994, Mowat et al. 2000). The patches of these vegetation types are required to be at least 2.0 ha in size, 19 (a) Lynx Foraging Quality (b) Lynx Deming Quality § - § i o _ O _ 00 cc 9 _ o _. 0 \6 E73 5 :1: I o _ o _ 1' <1- 0 _ o _ N (‘1 o - o — l l l l l l l l l l 0 20 40 60 80 100 0 500 1000 1500 # Hare Home Ranges Average Distance to Denning (m) (c) Non-lynx Habitat Interspersion O E -1 g _ o —1 \G B :l'.‘ o —1 ‘Q 8 _. o —1 I l l l 0 500 1000 1500 Average Distance to Non-lynxhabitat (m) Figure 4. Production functions relating lynx habitat quality to each of the three model components within a 250 ha home range: (a) foraging (for x <90, y = 1.1x); (b) denning (for x >400 and <1750, y = -0.0741x + 129.64); and (c) interspersion (for x >300 and <1500, y = -0.0833x - 25). 20 adjoined by other lynx habitat for 250% of their perimeter, and within 0.8 km of suitable summer foraging habitat (deï¬ned as vegetation types providing 220% understory cover) (Roloff 2003). A mosaic of lynx summer foraging and denning habitat would provide female lynx the ability to move kittens between denning sites to avoid predation (Koehler and Aubry 1994). Denning quality is, therefore, determined by a function of the average distance to lynx denning habitat from a 100 x 100 m grid of points within a lynx allometric home range (Figure 4b). Average distances to denning habitat <400 m provide the most suitable conditions and those >1750 m provide non-suitable conditions. The resulting grid map describes lynx denning potential at the scale of a lynx home range. Interspersion Component The interspersion component addresses travel requirements of lynx, and is a measure of the average distance to “non-lynx†habitat within a lynx allometric home range (Figure 4c). Throughout North America, daily movements of lynx have averaged 5-10 km depending on hare densities and season (Nellis and Keith 1968, Brand et al. 1976, Parker et al. 1983, Koehler 1990), and in Washington, they have demonstrated a reluctance to cross openings (e.g., ï¬elds, clear-cuts, developments) with little cover that are >91 m across (Koehler 1990, Koehler and Brittell 1990). These potential “barriers†could disrupt movements between foraging and denning sites (Koehler and Aubry 1994), so this habitat model component relates the amount of “non-lynx†habitat interspersion throughout the landscape to habitat quality. An average distance of >1500 m to the nearest non-lynx habitat within a lynx allometric home range provides the most suitable conditions, while an average distance <300 m provides non-suitable conditions. Pixels are assigned non-lynx habitat if they are identiï¬ed as having no suitability for denning or 21 foraging and do not provide adequate cover for travel. Inadequate cover includes human developments, water bodies, permanent openings and vegetation types having vegetation <2 m tall that are >91 m from suitable habitat (Roloff 2003). Forest patches with <440 trees ha'1 and <50% horizontal cover that are >91 m in length also qualify as non-lynx habitat. The hypothesized effect of habitat interspersion on lynx should be interpreted carefully, since this relationship is unknown (Roloff 2003).The resulting grid map of habitat interspersion describes the general contiguity of lynx habitat at the scale of a lynx home range. Final Habitat Potential The three components previously described are combined using a geometric mean to provide an index to the overall habitat suitability for lynx. The use of a geometric mean results in zero habitat suitability for areas that are unsuitable in any one or more of the habitat components, since each component is determined to be important for lynx survival (Koehler and Aubry 1994, Roloff 2003). The grid map of lynx habitat suitability is used to create hypothetical home ranges for lynx with the Plume (2001) G18 function, in the same way home ranges are created for hare. A habitat unit goal based on the allometric home range for lynx (250 ha) supports the simulation of viable, marginal, and non-viable lynx home ranges across the landscape (Roloff and Haufler 1997). Viable lynx home ranges meet a minimum quality threshold of 70, which is based on a relationship between lynx home range sizes and habitat quality, as indexed by ï¬tness indicators (e. g., survival, litter size, pregnancy rate) (see Figure 3: Roloff and Haufler 1997) demonstrated in previous ï¬eld studies. The threshold represents the point at which lynx home range size no longer increases rapidly with decreased habitat quality, and 22 ï¬tness indicators shift from good to bad (Roloff and Haufler 1997). Marginal lynx home ranges have an average quality 225 and <70, and represent areas where lynx may persist and reproduce, but would be vulnerable to temporal fluctuations in resource availability (Roloff and Haufler 1997). The number and spatial conï¬guration of these hypothetical home ranges would provide insight into the habitat potential for lynx in the region. Prior to evaluating habitat potential with the lynx model, the composition and structure of vegetation types needed to be quantiï¬ed and mapped across the landscape for the assessment of habitat suitability. Roloff (2003) suggested the use of an ecological land classiï¬cation that described biotic and abiotic factors of the landscape to properly stratify the variation in vegetation attributes deï¬ning habitat quality in the H81 models. ECOLOGICAL LAND CLASSIFICATION FOR THE UPPER PENINSULA The ecological land classiï¬cation used in this study was based upon the “Habitat Classiï¬cation System†constructed by Coffman et al. (1984), describing “habitat types†(Daubenmire 1966) of the UP that were deï¬ned by associations of tree species occurring within similar ecological conditions that had the potential to support similar successional trajectories and climax vegetation communities. The ecological conditions were characterized by a speciï¬c range of variation in abiotic factors (e.g., climate, landform, soils) that affect the composition and structure of vegetation within the various seral stages that are initiated by disturbance or developed through succession (Kotar and Burger 2000). The habitat type concept has been used in the Great Lakes region (Coffman et al. 1984, Kotar et al. 1988, Kotar and Burger 2000, Felix et al. 2004) to guide forest and wildlife management strategies by improving land-cover and forest-type classiï¬cations with metrics of site potential and an understanding of community 23 dynamics. Using the Coffman et al. (1984) habitat types as a basis for the ecological land classiï¬cation facilitated the stratiï¬cation of vegetation attributes necessary for the lynx model into ecologically meaningful and homogenous classes that were developed speciï¬cally for the forest communities of the UP. Coffman et al. (1984) described 21 habitat types for the UP, covering a gradient of soil texture and drainage classes that can influence the composition and dominance of tree species throughout seral stages and into late succession. Habitat types were named by combining the genera of the dominant climax species with that of an understory indicator species. Cofï¬nan et a1. (1984) provided a detailed description for each habitat type, including the landforrn and soil characteristics on which it occurs, the effects of silvicultural practices on species composition at the start of secondary succession, site indices for the primary species within each seral stage, and one or more speciï¬c trajectories outlining the canopy replacement of seral species by shade-tolerant species. Coffman et al. (1984) aggregated habitat types into more general categories denoted as “series†that were deï¬ned and labeled by the genera of the dominant climax species; these categories were further aggregated (though not mutually exclusively) into series groups based on soil moisture, as determined by texture and drainage properties. In summary, each of the 21 original habitat types belonged to l of 8 series, and each series belonged to 1 or more of 3 series groups, depending on the range of soil conditions that could be dominated by a particular series. Group I incorporated 3 series associated with sandy soils retaining little moisture and experiencing periodic drought: the Pinus series, Acer-Quercus series, and T suga series; Group II incorporated 4 series associated with soils with adequate moisture and ï¬ne textures ranging from loamy sand to clay loam: the 24 Acer—Quercus series, T suga series, Acer— T suga series, and Acer series; Group III incorporated 5 series associated with wet mineral or organic soils with impeded or variable drainage: the Tsuga series, Acer-Tsuga series, T suga- Thuja series, Fraxinus series, and Picea series (Coffman et al. 1984). Each of the 3 categorizations (habitat types, series, series groups) was meant to simplify a continuum of ecological conditions. I condensed the 21 original habitat types into 8 classes (hereafter, “habtypesâ€), based on combinations of similar series and series groups (Table 1). Each of the 8 habtypes contained several early successional stages with a variety of potential overstory compositions (Table 2). The PiVa habtype consisted of habitat types identiï¬ed by Coffman et al. (1984) that were dominated by jack pine and red pine on sandy outwash soils with excessive drainage. The QuAc habtype consisted of red oak/red maple dominated types that had potential early stages of aspen/birch, pine, or spruce/ï¬r occurring on sandy soils with good drainage. The TsMa habtype consisted of hemlock/maple dominated types with early stages of aspen/birch, pine, or spruce/ï¬r, occurring on well-drained sands and sandy loams. The Acer habtype consisted of sugar maple dominated types with early stages of aspen and northern hardwoods (i.e., beech, basswood) occurring on rich loamy soils with good drainage. The TsTh-dry habtype consisted of hemlock/cedar dominated types with early stages of aspen/birch, spruce/ï¬r and maples occurring on clays and loams with moderate to somewhat poor drainage. The Frax habtype consisted of black ash dominated types with early stages of seral hardwoods occurring on clays and loams with poor drainage. The TsTh-wet habtype consisted of cedar/hemlock dominated types with early stages of seral hardwoods and spruce/ï¬r occurring on a variety of poorly drained soils. Finally, the Picea habtype consisted of 25 Table 1. Habitat types originally described by Cofï¬nan et al. (1980) for the Upper Peninsula of Michigan with associated soil properties and “habtype†groupings. . . . b Habitat Typea Sorl (drainage / texture) Habtype Pinus- Vaccinium-Deschampsia E sand PiVa Pinus- Vaccinium-Carex E sand Quercus—A cer-Epigaea W sand QuAc A cer-Quercus- V accim'um W sand Tsuga-Maianthemum- Vaccinium W loamy-sand TsMa Tsuga-Maianthemum W sandy-loam Acer- T suga-Dryopteris W loam Acer Acer- V iola-Osmorhiza W loam Acer-Osmorhiza-Caulophyllum W silt loam Tsuga-A cer-Mitchella MW clay TsTh-dry Tsuga— Thuja-Lonicera MW clay T suga- Thuja-Petasites SP clay Tsuga-Maianthemum- C optis SP loam F raxinus—Impatiens P loam/clay F rax F raxinus-Mentha— C arex P loam F raxinus—E upatorium P clay Tsuga- Thuja-Mitella VP sand/ loam TsTh-wet T suga- Thuja-Sphagnum VP organic Picea-Osmunda VP deep organic Picea Picea- Chamadaphne-Sphagnum VP deep organic a The Acer—Quercus- Viburnum habitat type was not included, due to its minor extent in the Upper Peninsula. Soil drainage codes: E = excessively drained; W = well-drained; MW = moderately well-drained; P = poorly drained; VP = very poorly drained. 26 Table 2. Dominant tree species assemblages within seral stages of each habtype for the Upper Peninsula of Michigan. Habtype Early—serala Mid-seral Late PiVa jack pine jack pine jack pine, red pine QuAc (l) aspen, white birch ( 1) white pine red oak, red maple (2) jack, red pine (2) white spruce, ï¬r TsMa (l) aspen, white birch (1) white pine hemlock, maple, (2) jack, red pine (2) white spruce, ï¬r beech, yellow birch Acer aspen sugar maple, basswood, sugar maple, hemlock white ash, beech TsTh-dry aspen, white birch, (1) white spruce, ï¬r hemlock, cedar, ï¬r, balsam poplar (2) maple maple Frax aspen, white birch, balsam black ash, red maple black ash, red maple poplar TsTh-wet white birch, balsam poplar, cedar, ï¬r, black spruce, cedar, hemlock, ï¬r red maple tamarack Picea lowland brush, white birch, black spruce, tamarack, cedar, black spruce, tamarack red maple ï¬r, jack pine a Seral stages are not equally distributed through time across habtypes. Numbered species groups within seral stages represent different potential trajectories of succession. 27 black spruce/tamarack dominated types with early stages of cedar and ï¬r occurring on very poorly drained organic soils and peatlands. The combination of a habtype and seral stage (as deï¬ned by overstory species) resulted in an ecological land unit (ELU), which provided the basis of the ecological land classiï¬cation necessary for the lynx model. The classiï¬cation system developed by Coffman et al. (1984) was intended as a ï¬eld guide for in situ classiï¬cation of forest stands that could improve the predictions of forest management prescriptions; it did not involve mapping habitat types across the landscape. Mapping the spatial conï¬guration and extent of the ELUs across the UP required the integration of multiple spatial data layers within a GIS. Each data layer represented components of the ecological conditions that could identify the occurrence of a habtype, and within each habtype, a seral stage. My objectives included the assessment of current and presettlement habitat conditions for lynx, to examine how forest composition and structure effected lynx habitat potential within each time period. Theoretically, habtypes should be similar between time periods where certain ecological conditions (e.g., soils) had not been measurably altered; however, at any given point in time the distribution of seral stages within habtypes will be dictated by forest disturbances (Kotar and Burger 2000). The current distribution of seral stages results from a disturbance history with considerable human influence (e.g., ï¬re suppression and exclusion, development, harvesting), and is very different from the “shifting mosaic†of seral stages that historically occurred in the Great Lakes region due to natural disturbances (Zhang et al. 1999, F relich 2002). Therefore, it was necessary to incorporate different data layers and develop separate methods to map the location and quantify the vegetation of ELUs for current and presettlement landscapes. 28 Current Ecological Conditions A combination of digital spatial data and region-wide forest inventory surveys was used to delineate ELUs and quantify the current species composition and structure of forest types in the UP. Thematic spatial data included: (1) current vegetation, to separate forested from non-forested areas and describe the species composition of the overstory; (2) presettlement vegetation, which indicated potential natural communities in the absence of human disturbance; (3) soils, to provide a context for identifying habitat types; and (4) ecoregions, which deï¬ned broad ecological units having similar climate and physiography. The intersection of these thematic GIS layers created the framework for identifying ELUs on the landscape. Survey data provided plot-level measurements of stand composition and structure within each ELU to quantify speciï¬c vegetation attributes required by the lynx model. A ï¬nal data layer consisted of unclassiï¬ed satellite imagery that was used with the survey data to model forest structure as a continuous variable across the landscape and account for variation within ELUs. Spatial data were compiled and manipulated using GIS software which included ARC/INFO 8.1 and ArcView 3.2 (Environmental Systems Research Institute, Redlands, Calif.) and ERDAS Imagine 8.7 (ERDAS, Inc). All thematic spatial data were projected in the Michigan GeoRef Coordinate System. The ecological land classiï¬cation was created in a raster- based format to facilitate the calculation of habitat units within a GIS; habitat quality for hare was assessed with a grid resolution of 30 m, while that for lynx was assessed with a grid resolution of 90 m. The multiple data sources that were collected and the methods that were employed to utilize them as input for the lynx model included the following: 29 Current vegetation — I obtained a land-cover dataset publicly available from the Michigan Department of Natural Resources (MDNR) which had been created for their Integrated Forest Monitoring, Assessment and Prescription (IF MAP) project in cooperation with the GAP analysis program (MDNR 2001, Donovan 2005). The IFMAP land cover was developed using supervised and unsupervised classiï¬cation of 3 seasons (spring, summer, fall) of Landsat satellite imagery collected at a 30 m resolution between 1999-2001. Forested pixels were deï¬ned to have 225% canopy cover; forest composition was deï¬ned according to the dominant species (260% cover) in the canopy, with “mixed†categories for stands without dominant species. Land cover was classiï¬ed to Anderson Level II with an estimated accuracy of 87%, and for some classes Level III, with an estimated accuracy between 37-87% (Donovan 2005). Subdedi (2005) found large inaccuracies in the IF MAP prediction of some Anderson Level 111 categories, especially with aspen and oak, across public lands in the UP. My own accuracy assessment of the map revealed similar results (Appendix C); therefore, most Anderson Level 111 categories were condensed to Level II for my analysis (Table 3). Because mixed forests in the UP generally succeed to hardwoods in the uplands and conifers in the lowlands, each of the mixed forest cover types was grouped into the appropriate condensed category. The decision to group the mixed forest cover types in this manner was supported by ï¬ndings from the accuracy assessment. Although the loss of detail made it difï¬cult to distinguish canopy species, the IFMAP layer represented the most recent land-cover data for the UP, and the Anderson Level 11 categories that provided broad classes of overstory composition were accurate. The accuracy and resolution of this dataset made it the most reliable spatial layer for 30 Table 3. Areas of land-cover classes identiï¬ed by the 2001 IFMAP in the Upper Peninsula of Michigan. IFMAP land cover Low/High Intensity Urban Airports Road/Parking Lot Non-Vegetated Farmland Row Crops Forage Crops/Non-Tilled OrchardsNineyards/N ursery Parks/Golf Courses Sand, Soil, Exposed Rock Mud Flats Other Bare/Sparsely Vegetated Herbaceous Openland Upland Shrub/Low Density Trees Northern Hardwood Association Oak Association Aspen Association Mixed Upland Deciduous Upland Mixed Forest Pines Other Upland Conifers Mixed Upland Conifers Water Lowland Deciduous Forest Lowland Coniferous Forest Lowland Mixed Forest Lowland Shrub Floating Aquatic Emergent Wetland Mixed Non-Forest Wetland Condensed Developed Agriculture Non-vegetated Openland Upland shrub Upland deciduous Upland coniferous Water Lowland deciduous Lowland coniferous Lowland shrub Wetland Area (hQ 76,603 126,628 29,271 160,842 84,067 1,929,651 509,131 239,319 1,562 708,289 250,428 206,608 identifying landscape attributes; for this reason it was used to describe current vegetation and correct any inconsistencies in the descriptions of ecological conditions resulting from the intersection of other spatial layers. Presettlement vegetation — I obtained a spatial layer of presettlement vegetation that was constructed by Comer et al. (1995) using GLO surveyor notes recorded between 1816 and 1856, before the extensive logging that took place during the late 18005. This data layer provided information on forest types that occurred prior to European settlement, and thus, described potential species assemblages determined by natural disturbance and physiography, in the absence of signiï¬cant human influence (Table 4). The map was created by interpretation of recorded landscape descriptions for each 2.56 km2 section and occurrences of tree species tallied at and between 1.6 km section comers (and 0.8 km quarter comers) (Comer et al. 1995). Comer et al. (1995) combined this information with ancillary GIS data (e.g., landforms, wetlands, topography) to delineate forest cover types. The exact spatial accuracy of the map is unknown and interpolated boundaries between section comers and quarter comers were expected to contain some error where ancillary data were unable to distinguish actual borders between cover types. Also, the potential bias of surveyors for long-lived species (e.g., beech) (Manies et a1. 2001), combined with the scale of survey data which can be expected to miss small (< 20 ha) wetlands and wind-throw disturbances, has most likely resulted in the under- representation of seral communities (Comer et al. 1995). Regardless of these potential inaccuracies, the cover classes delineated by the presettlement layer identiï¬ed the dominant species within many of the forest types (Table 4), and provided more detailed 32 Table 4. Area of cover types in the presettlement vegetation layer from Comer et al. (1995) located in the Upper Peninsula of Michigan. Cover Type , Area (ha) Upland Non-forested Exposed bedrock 3,650 Sand dune 1,266 Grassland 37 Oak/Pine Barrens 387 Pine Barrens 25,556 Upland Forested Aspen/Birch 72,085 Beech/Sugar Maple/Hemlock 533,844 Sugar Maple/Hemlock 939,430 Sugar Maple/Basswood 86,216 Sugar Maple/Yellow Birch 383,908 Hemlock/Yellow Birch 118,957 Hemlock/White Pine 215,463 White Pine/Mixed Hardwood 16,430 Mixed Pine/Oak 6 White Pine/Red Pine 129,552 Jack Pine/Red Pine 81,117 Spruce/Fir/Cedar 360,584 Lmiland Non-forested Wet Prairie 3,909 Shrub swamp/Emergent marsh 9,290 Muskeg/Bog 127,946 Lake/ River 985,918 Lowland Forested Black Ash 243 Mixed Hardwood 55,598 Cedar 104,552 Mixed Conifer 94,182 33 information concerning potential species assemblages than that of the IF MAP layer. The GLO data have proven useful throughout the Great Lakes region as a means to reconstruct natural disturbance regimes and to compare current forest conditions with those of the presettlement era (F relich 1995, Zhang et al. 1999, Schulte and Mladenoff 2001, Cleland et al. 2004, Friedman and Reich 2005). The presettlement layer was combined with the IF MAP layer after being converted to a raster format, and 2 iterations of reclassiï¬cation were used to provide the most informative description of the vegetation. The ï¬rst iteration involved reclassifying pixels that described inaccurate ecological conditions (e. g., upland versus lowland), using a nearest-neighbor approach. I assumed the IF MAP layer provided greater accuracy because it was constructed from remotely sensed imagery; therefore, it was given precedence in determining the ï¬nal attributes for a pixel. For example, if a pixel described current vegetation as upland deciduous and presettlement vegetation as cedar swamp, its presettlement description was reassigned the nearest presettlement cover type that was an acceptable match for upland deciduous. Acceptable matches were based on the distinction between upland and lowland for each layer, with the exception of 2 presettlement cover types: spruce-ï¬r-cedar and aspen-birch. These 2 cover types contained tree species that could tolerate a range of ecological conditions; therefore, pixels belonging to those presettlement cover types that were located on a lowland cover type for the IF MAP layer were not adjusted, in spite of the fact that their presettlement classiï¬cation was originally upland (Table 4) (Comer et al. 1995). Non-forested presettlement pixels were assigned the nearest forested pixel value if the IF MAP cover type was forested. 34 Visual inspection of grid overlays indicated that pixels located between upland V and lowland cover type boundaries were responsible for the majority of inconsistencies, and 35% of misclassiï¬ed pixels, representing 6% of all forested areas, involved a mismatch between upland hardwoods and lowland conifers. This was expected due to the interpolation technique used to create the presettlement layer (Comer et a1. 1995) and the heterogeneous landscape of the UP, especially in the central and western portions where a mosaic of drumlins and lowland depressions occur. The second iteration of reclassiï¬cation involved reassigning pixels that were described as having a presettlement cover type of aspen-birch. This presettlement cover type would not provide useful information for identifying an ELU because aspen and/or birch forest types could occur as a seral stage in nearly every habtype (Table 2). Therefore, a nearest-neighbor approach was again utilized to assign a new presettlement cover type based on acceptable upland/lowland matches within close proximity. The 2 iterations of reclassifying pixels resulted in the most informative combination of data layers, given that I assumed: (1) any contradictions in vegetation were reflecting spatial inaccuracies in the map layers and not actual conditions, and (2) the aspen-birch cover types represented early seral stages and would have been replaced through succession by surrounding forest types. Soils — The spatial data that were used to describe soils were obtained from the State Soil Geographic (STATSGO) database for Michigan, developed by the Natural Resources Conservation Service (USDA, NRCS 1994). The boundaries of map units in the STATSGO layer are broad groupings of soil associations, which contain the most common soil series that occur together on the landscape. Soil properties were the key factors in the determination of habitat types by Coffman et al. (1984), and aside from 35 regional differences in climate, would account for most variation in vegetation composition and structure within otherwise similar forest types, given the generally small differences in local topography throughout the Great Lakes region (Kotar and Burger 2000). The information provided by the STATSGO database that was used to identify habtypes for each soil series (component) within a map unit included percent coverage within the map unit, surface texture, drainage, depth to water table, depth to bedrock, common woodland species (with site indices), and understory indicator species (both woody and non-woody). Additional characteristics of each soil series were obtained from the NRCS Ofï¬cial Soil Series Descriptions (OSD), which provided a description of common landforms, typical pedons, and forest vegetation. Unfortunately, the STATSGO layer was mapped at a coarse scale of 1:250,000, and the NRCS has suggested that its use be limited to interpretation at broad geographic regions (USDA, NRCS 1994). The more detailed data provided by the Soil Survey Geographic (SSURGO) database are currently unavailable for greater than half of the UP (USDA, NRCS 2006), thus, the STATSGO database was the only available source of soil information with associated spatial data that could be integrated into a GIS for my extensive study area. The limitations in the STATSGO soils data were mediated with several approaches. Each STATSGO map unit contained multiple soil components within its boundary, and each soil component represented a soil series. There were 84 unique map units for the UP with 181 unique soil series; the extents of 15 soil series covered ~50% of the UP, while that of 35 soil series covered ~75%. Map unit polygons ranged in size from 200 ha to 370,000 ha, and consequently, large map units encompassed a variety of 36 soil series (including both upland and lowland) that were not spatially delineated. The STATSGO layer was converted to raster format and combined with a reclassiï¬ed IFMAP layer that delineated upland from lowland cover types. The resulting grid delineated STATSGO map units into upland and lowland soils, and a relational database was used to assign properties from hydric soil components to the lowland soil pixels, and that of mesic/xeric soil components to the upland soil pixels. The soil information from the map unit components was interpreted alongside the vegetation described by the IFMAP and presettlement layers to identify the most probable habtype. Similar to the adjustment of presettlement cover types, this process resulted in a logical assignment of soil information to the grid pixels, and avoided the inconsistent descriptions of ecological conditions that would have occurred if map units were assigned information from only the most dominant soil component (e.g., upland deciduous on very poorly drained organic muck). Ecoregions — A spatial layer of ecoregions was used to separate otherwise similar habtypes that existed within different climatic and geophysical regions of the UP. Several ecoregion delineations exist for the Great Lakes region, including those of Omemik (1987), Albert (1995), Bailey (1995), and Cleland et al. (2005b), and some of these delineations are continually being reï¬ned with increases in the amount and resolution of available spatial data (Cleland et al. 2005b). 1 used the ecoregions layer developed by Cleland et al. (2005b), which depicted the 2 units of the subregion planning scale (sections and subsections) as outlined by the Forest Service national hierarchical framework of ecological units (Cleland et al. 1997). The UP contained 21 subsections that were contained within 7 sections, ranging in size from 20-670 km2 and 580- 37 1,538 kmz, respectively. Ecoregions were grouped into 2 “subregions†based on their broad-scale location within the UP; sections in the eastern subregion included 212R and 212T, while those in the western subregion included 212.1, 212L, 212S, 212X, and 212Y. This ecoregions map was chosen over the others because it was being actively used by the Forest Service (Cleland et al. 2005b) and, consequently, coordinated with other data collected for my study (see following). Also, the locations of ecoregion boundaries were mostly consistent with previously developed maps. ELU grids — The intersection of soils and vegetation data resulted in 4,781 unique combinations of forest cover categories. The high number of unique combinations was a result of the soil layer, which was not consolidated prior to the grid intersection in order to preserve all available information from the STATSGO database. Only one seral stage and ELU could be described for most habtypes due to a lack of decisive information from the vegetation layers, and one ELU (TsTh-Picea) was created to describe areas that could not be split between the two lowland coniferous habtypes. The base ELU grid contained 14 ELUs, with 2 non-forest classes and 12 forest classes (Table 5). The ecoregions layer was used to stratify the forested ELUs into smaller units and create an ecoregion-ELU grid with 220 ecoregion-ELUs. Each grid value contained the ELU and subsection identiï¬cation (e.g., PiVa-212Rb). The forested classes in the ecoregion-ELU grid ranged in total area from <1 ha to 158,715 ha, with a median area of 4,004 ha (Figure 5). Forest plot surveys — Vegetation attributes necessary for the lynx model were obtained from plot surveys collected by the Forest Inventory and Analysis (F IA) program of the Forest Service (USFS 2006). The FIA program provides detailed measurements of the nation’s forests on 5-10 year cycles (5 years in Michigan), using a double sampling 38 23:55 80 23:55 30 32.20% 22% £155 233. 95 SN .wom $23 .595 .595 .0om .Nom den .92 .373: .1: £53 .053 .053 .053 .503 .3935 .57: .m>-:m .95 J03 .ofl .m2 .m3 .03-N3 .03 .02 .02 .m2 .3: .3: -Im .03-mm-mm .Em-2m-m< maceotcoo 053: 20-5%... 3.00-0.5... :5 SN .0842†den .33 .53 .02 .33 .03 .w53 $5793 m>-2m .o: .003 .0. .3070? :2 .m3 .0333 62-0.2 .mfl .Im-2m .BmEzm .mm-2m-m< £63280 0.8—m: 89¢. 83.. £323. 83 23.6.5 080 mac—52060 0:39. 22$. Nev/3m... Add? . A500 mom .58 08308033 .m>-2m .EmEzm .0368 £9.03 :3: .NM: .22 .553 .55_-m0_ .0702 .3935 max: .m>-mm d? .NS -03 £3-52 .02 62 ._m_ .wï¬ .52 .33 .5: -mm .03-mm-mm .Im-2m-m< 380088 053: «Emu. 3323. 23.30 080 23.50 33 283280 053: 350 No.35 A09 man .53 .mom-mom .wa .0262 .33 .53 .mw3 .m3 .2: .053 .m53 .35_ .053 £03 .503 .03 .N03 .003 .m2 .02 Add? £2.80 A: :52 .53 .03 .33 .03 .92 .02 .3: .23 .wï¬ .53 .52 .OM .Om-mZE .53-“: £2.30 a 288088 ESE: 3.50 EEO 8N0 mom .08 .m3 .02 .NS .33 .33 .2: .9: .503 .N0_ .33 .m2 .53 .03 .v2 .32 .wQ .52 .5: 5:3 LTOM dag: £8.80 a: 32868 ESE: «>5 «>5 032m <2 <2 2232 .925. 2.2% <2 02% 93:03 .833 .05— -:oqo .0o§owo>-:oc $2 $2 .obd—sotwa .032050 32 328-28: seaflowo> “5823085 sozSowo> 6050 250a: Dam :58 .590 an: 2232 0855 M 283022 mo Esmï¬aom Saab 05 E motowï¬mo Dam 05 39302 083530 35 me??— Etmam E53, momma—U .0 205. 39 630800 003.35 cmeocz 20 5 00:8 320 .005 288608 Em and 05:32 0 EmamESEmzu £333 :85 30:?» u m; .3565. 2355 2:0 823 n 95 £08203.er 0.525: 308.82? 520.5: u 03 AEEeaooE. x85 2qu Swsm n Em .3323? 3:5: 05¢ 08 " n2 A953 mauxmzmv x00 02 n OM .Aaxmbxamum £503 :85 .593 H mm .Auaemï¬tg 9.353 059 xofl. u E. .Emnmheng 0.05.5 x0280: 5830 n Em .Auzuutmï¬a 6.5.50 083an M 3m A3525 08.25 83% 0.003 H mm .Acmï¬aég 8.3 5 E 8803 H mm 60me 9.32.3855 :3 0.003 n <m .QmEeEEmh “Snack 6082000202 “Eamon: comma n m< .Aeaohhzexw wzwukv 0080 5.6005 n m< m 28232 «2.53 .213; .93.: .9: .5 .§ .5 ....2 .2: A~0m0_.m¢_ at 832 .3: .3: £32 $731322 .5 .52 ....2 .5: .2: 35:2: .332 .5 .22 .52 .3: .212: 9:95 00950.3: .90an 3:8 .03-.._m-mm .953 <m QEQBm 000308: dEaBm 328 63-3-8 deg... <m 03 £893 000308: .9526. 50:8 63-2-8 deg... <m 03 dEaBm 0830.8: .9536 3.008 03-3-8 .955 <m .mm-m< 285.3050 00232 380.3080 05—32 38888 05232 3050600 05052 32m 805 63-5.5 535$ 52..— «8:. “85-55; 835$ 52..— 958 390 an: Ea“: oomhï¬m «cozï¬owo> 80820395 =00<Swu> €2.50 030mm Dam .6203 m 2.05 40 H N .mcommouooo mo 58 3253 5 go @530 05 89a 35953 53833 .33 5mE8—2 mo flaw-55¢ 5&3 2: E GOA-n: SE: c5_ 505208 mo 5:35va 559m .m 053 < \ macaommoam 8 mzoamxoom -.. .395 I .... . m. . “GA 32.25 I .5 - ... 7.51:: n-Im 5c... D _. ,, 8.: SEE: D ./. .r.....< D , T 8:22.: a“ .... \ a. _ . U<5 l . «>5 I 7.29.7.5: I , 1 34m ,. 4l for stratiï¬cation approach (Smith 2002). Remote sensing imagery is used to identify forested land across a state, which is then systematically ï¬eld sampled within a hexagonal grid of 5 interlocking panels, with one plot survey for every 2,428 ha (6,000 ac). The state of Michigan provided funding to increase the intensity of ï¬eld sampling and essentially triple the number of plot surveys conducted, reducing the area represented by each plot survey to 809 ha (2,000 ac). The plot surveys consist of a 4-subplot cluster, each subplot having a ï¬xed radius of 7.32 m, within which trees having a diameter at breast height (dbh) 212.7 cm (5.0 in) are individually measured. A11 trees with a dbh <12.7 cm are counted within a smaller ï¬xed radius (2.07 m) microplot, and only those with a dbh 22.54 cm (1 in) are individually measured. Data that are recorded for individually measured trees include species, dbh, height, crown ratio, and crown class (e.g., dominant, co-dominant), while seedlings (dbh <2.54 cm) are simply tallied by species. Landscape attributes that deï¬ne one or more “conditions†(e.g., land use, forest type, physiographic class, stand density, stand size) are also recorded for each plot. The plot data that were available for Michigan at the beginning of my study included information for 80% of the F IA plots in the 6th cycle (4/5 sub-cycles), collected between 2000-2003. Data from the ï¬nal sub-cycle (2004) became available in December 2005, resulting in a total of 4,741 plots with recorded tree data for the UP. Plot data were downloaded from the FIA website (http://ï¬a.fs.fed.us/) in two different formats: (1) spreadsheets that could be imported to a relational database; and (2) text ï¬les that were formatted for use with the Forest Vegetation Simulator (FVS). The FVS is a framework used by the Forest Service to standardize forest growth and yield modeling, with variant growth models calibrated for speciï¬c regions (Dixon 2002). I used F VS to calculate the 42 per hectare stand attributes required for the H81 models (e. g., basal area, average dbh, stem density, canopy closure) from the F IA data. The F IA data were ï¬ltered to remove plots that could not be used with my mapping methodology. Plots with <25% canopy closure were removed to ensure that vegetation attributes from the FIA data corresponded to the IF MAP deï¬nition of forested cover. Plots with >1 landscape condition were removed to ensure that survey data were describing homogenous samples that could be linked to single ELUs. This ï¬ltering process reduced the number of available plots by ~30%, leaving 3,256 plots for my analysis. The frequency distributions of stand attributes necessary for the model were compared between the original and ï¬ltered datasets to verify that no bias was introduced through the ï¬ltering process. The stand attributes required for the snowshoe hare sub-model included horizontal cover and understory dominance. Understory dominance was calculated by the percentage of conifer stems among stem densities of trees with a height to crown <3 m. Horizontal cover was not measured in the FIA surveys, so a new method was developed to extract an estimate from the data that were gathered (Appendix D). A post-processing function in F VS was used to create tree-lists for the Stand Visualization System (SVS) (McGaughey 1997), a program that creates 3-dimensional diagrams of 0.40 ha forest plots based on survey data, with each tree rendered according to its attributes (e. g., dbh, height, species). Plant form deï¬nitions are used to deï¬ne appearances for different species and can be controlled by the user; I altered plant form deï¬nitions in SVS by removing leaves from deciduous species to simulate winter conditions. For each plot, the proï¬le-view diagram in SVS was saved and converted to a 1-bit (black and white) 43 bitmap. The proportion of black to white pixels within each height strata determined the horizontal cover for a given plot (Appendix D). The method of calculating horizontal cover with SVS was complicated by the absence of size measurements for seedlings, which would have a signiï¬cant influence on estimates of horizontal cover. The Forest Service assigns a common dbh value of 0.25cm (0.1 inch) to all seedlings in the FVS-ready ï¬les, for the purposes of growth modeling. I felt that a common dbh of 0.25 cm did not accurately reflect the true variation in seedling sizes, given that seedlings were deï¬ned as any tree with a dbh <2.54 cm. Large seedlings (e. g., dbh ~ 2.0 cm) had the potential to provide more horizontal cover than small ones (e. g., dbh ~ 0.25 cm), depending on species, and the absence of this distinction could result in conservative estimates of horizontal cover. Therefore, I assigned each plot a frequency distribution of seedling dbh values that ï¬t a negative exponential trend and ranged 0.25-2.29 cm (0.1-0.9 inches) in increments of 0.25 cm (0.1 inches). This procedure resulted in a high ratio of small to large seedlings, which I hypothesized was a more adequate representation of the vegetation. Original plot data were retained to examine the effect of manipulating the seedling size distribution on measurements of horizontal cover. Seedling heights were determined by F VS using dbh-height models that were speciï¬c to each species. The uncertainty in my estimates of horizontal cover warranted an additional measure of hare habitat quality from the plot data. Litvaitis et al. (1985) described a strong relationship between estimated stand densities of snowshoe hares (indexed by live- trapping and pellet counts) and “stem cover units†(SCUs), calculated as the count of trees and shrubs <7.5cm dbh and >05 m tall, with coniferous stems weighted by a factor 44 of 3. The weighting factor for conifers was based on the estimated difference in visual obstruction provided by coniferous stems and deciduous stems (Litvaitis et al. 1985). I used Litvaitis et al.’s (1985) linear regression model of hare density as predicted by SCUs to create a production function for hare habitat quality (Figure 6). Hare habitat quality was highest at 255,652 SCUs/ha and lowest at <23,043 SCUs/ha, corresponding to the maximum and minimum densities of hares in the regression model (Litvaitis et al. 1985). This alternative index deï¬ned hare habitat quality in a similar manner to the indices in the original hare sub-model (Roloff 2003); horizontal cover and SCUs were essentially measurements of the same stand attribute, and the presence of conifers increased quality in both indices. In summary, hare habitat quality was measured in 3 separate ways: an HSI similar to that of Roloff (2003) which combined understory dominance with SVS estimates of horizontal cover using manipulated seedling distributions (HSI-1) and original seedling distributions (HSI-2); and an HSI using the regression equation from Litvaitis et al. (1985) measuring stem cover units (HSI—3). Linking plots to ELU grid — The actual coordinates of F IA plot locations were mandated conï¬dential by Congress in 2000, in accordance with amendments to the Food Security Act of 1985, to ensure that private landowners are not identiï¬ed with plots and that long-term integrity of plot data is maintained. The use of FIA plot locations by outside researchers is facilitated through the Spatial Data Services (SDS) for each F IA regional unit. I cooperated with Geoff Holden at the North Central Research Station SDS (USFS, St. Paul, Minnesota) to have the FIA plot locations overlaid with the spatial data layers. Conï¬dentiality concerns over the amount of detail and the multiple iterations of the ELU grid restricted the amount of information that could be disclosed. Also, the ï¬nal 45 Hare Habitat Quality o o— _ V? 0.. x o— o a .1: a R :1: E’. :2: II. 0.. v _ V! o o_ N o— — c: I l l 1 l I l 0 10 20 30 40 50 60 Stem Cover Units (thousand/ha) Figure 6. Production function of hare habitat quality as predicted by stem cover units. Hare density was scaled to index habitat quality, and the regression model from Litvaitis et al. (1985) was used to determine the slope of the relationship with stem cover units (fory >23.04 and <55.65,y = 3.067x — 70.667). 46 20% of the FIA survey data was released after SDS had completed the analysis. For these reasons, the overlay of FIA plot locations by SDS was used only to examine the accuracy of the land cover map (Appendix C), and in modeling forest structure using satellite imagery (see following section). The link between F IA plots and the ELU grid was established by matching the ecological conditions described by each data source. Determinant variables for the plot data included physiographic class and species composition of the canopy and sub—canopy, which were used to identify the most probable of the 12 ELUs (Table 5) for each plot. Species composition was deï¬ned by FVS calculations of species crown-cover within 3 canopy strata that I hypothesized would provide the best indication of ELU: (1) open grown, dominant, and co-dominant trees; (2) intermediate and overtopped trees; and (3) seedlings. Plot estimates of cover were analyzed with TWINSPAN (Hill 1979, Hill and Smilauer 2005), a program that conducts 2-way indicator species analysis using ordination and provides a hierarchical classiï¬cation of sample plots according to species abundance and site ï¬delity. Species compositions were summarized by TWINSPAN group, and each plot was assigned an ELU based on an interpretation of its TWINSPAN group and physiographic class. Ecoregion subsection (identiï¬ed in the F [A database) was used to further stratify the plots to match the ecoregion—ELU grid, and subsection ELUs that did not contain 210 F IA plots were aggregated by section or subregion, depending on the number of available FIA plots. For example, the PiVa ELU had >10 FIA plots within section 212R, but <10 F IA plots within subsection 212Rb. Therefore, all grid values belonging to PiVa-212Rb were aggregated at the section level to become PiVa-212R. The PiVa ELU also had <10 47 F IA plots within section 212S, resulting in those grid values being aggregated at the subregion level to become PiVa-212west. In this way, ecoregion—ELUs were required to contain 210 samples for adequate measurements of stand attributes used in the lynx model. The ï¬nal grid of condensed ecoregion-ELUS contained 2 non-forested classes (non-forest and shrub) and 134 forested classes (Table 6); 88% of forested land was described by ELUs aggregated at the subsection level, while 75% and 50% was described by subsection ELUs with >20 and >60 FIA plots, respectively. The transfer of stand attributes from the plot data to the grid classes involved the calculation and assignment of mean values, and the random simulation of a value’s spatial distribution within ecoregion-ELUS, to account for the local variation that would otherwise have been masked by the assignment of a mean. The attributes that were spatially simulated included each of the HSI measurements for hare habitat quality, since the scale of habitat selection by hares (4.5 ha) would be most affected by local heterogeneity. The local variation of other stand attributes necessary for the lynx model, including the denning and interspersion variables (e. g., basal area, average dbh, canopy cover, stem density), were assumed to have little effect on lynx habitat selection due to the larger scale of analysis (250 ha). Consequently, these attributes were averaged for each ecoregion-ELU. I simulated the spatial variation of hare habitat quality for each of the 3 HSI estimates by randomly assigning pixels within each ecoregion-ELU into 1 of 8 quantile classes. Two of the quantile classes represented the upper and lower 5% of the distribution of HSI values from the F IA plots for an ecoregion-ELU, and 6 represented 15% intervals between. The median values were calculated for each quantile class and 48 Table 6. Forested classes from condensed ecoregion-ELU grid for the Upper Peninsula. ELU Subregion Section Subsectiona Grid code F 1A plots Area (ha) PiVa east _ 111 14 507 212R 3011 13 2,728 a 31 1 1 10 7,156 west 211 10 12,151 QuAcl east 121 86 1,807 212R 3021 80 6,612 a 3121 53 33,708 b 3221 14 8,829 west 221 20 1,680 212S 4021 19 9,222 c 4221 12 7,310 QuAc2 east 122 18 2,710 212R 3022 16 7,198 west 222 26 848 2128 4022 25 1,833 c 4222 18 1,494 TsMal east 212R 3031 82 15,811 a 3131 32 53,148 b 3231 19 37,948 e 3531 24 37,100 212T 5031 30 12,437 b 5131 23 42,686 west 231 48 22,248 212S 4031 32 15,607 c 4231 13 27,803 q 4431 12 1 1,560 TsMa2 east 212R 3032 330 6,606 a 3132 194 158,716 b 3232 39 34,970 c 3332 22 17,888 e 3532 69 63,243 212T 5032 176 4,372 b 5132 149 84,573 e 5332 23 17,962 west 232 264 7,477 212.1 1032 22 6,848 o 1332 15 9,403 212S b 4132 14 12,216 c 4232 84 75,030 n 4332 41 28,434 q 4432 65 48,720 212K 6032 28 3,132 _ q 6332 18 14,832 Acer east 212R 3042 83 6,260 3142 35 15,490 _ 3542 37 22,499 49 Table 6 (cont’d). ELU Subregion Section Subsectiona Grid code F [A plots Area (ha) Acer east 212T ‘ 5042 82 1,315 b 5142 79 156,677 west 212.1 1042 221 5,819 b 1142 149 152,201 c 1242 64 61,610 212S b 4142 40 48,224 c 4242 112 123,150 n 4342 100 103,361 q 4442 68 81,785 212x 6042 86 205 c 6242 70 70,691 q 6342 16 20,509 212Y a 7142 16 12,945 TsTh-dryl east 212R 305 1 45 19,187 e 3551 24 37,828 212T 5051 15 2,523 b 5151 12 5,469 west 251 72 24,004 2121 1051 13 42,636 212L b 2151 11 6,626 212S 4051 34 10,939 c 4251 15 23,305 n 4351 10 16,895 TsTh-dry2 east 212R 3052 88 14,267 d 3452 25 25,927 e 3552 43 35,994 212T 5057- 17 2,698 b 5152 13 18,170 west 252 340 4,520 212] b 1152 41 31,959 c 1252 30 22,673 0 1352 54 63,733 2121. b 2152 26 31,869 2128 4052 98 2,409 b 4152 20 15,603 c 4252 14 10,914 n 4352 60 50,962 212Y a 7152 85 90,029 Frax east 212R a 3160 13 2,891 b 3260 27 10,905 c 3360 10 7,503 (1 3460 14 957 e 3560 35 6,576 212T 5060 61 16,284 b 5160 52 70,565 L west 260 73 9,409 50 Table 6 (cont’d). ELU Subregion Section Subsectiona Grid code F IA plots Area (ha) Frax west 212] 1060 29 2,757 b 1160 19 4,783 212S 4060 36 1,277 c 4260 1 1 25,447 11 4360 l 1 4,544 q 4460 12 8,393 TsTh-wet east 212R 3 3170 28 13,521 b 3270 44 11,145 c 3370 13 5,611 (1 3470 10 1,061 e 3570 81 23,348 212T b 5170 89 34,939 e 5370 18 11,551 west 270 80 2,689 212.1 1070 13 2,213 2123 4070 55 2,784 c 1 4270 24 2,905 n l 4370 13 2,809 q 4470 15 280 TsTh-Picea east 212R a 3180 46 51,165 b 3280 78 76,444 c 3380 15 19,622 d 3480 15 11,074 e 3580 90 66,507 212T 5080 138 613 b 5180 l 18 120,468 e 5380 20 16,567 west 2121 1080 13 27,163 b 1180 10 24,204 212L 2080 10 7,449 212S 4080 85 9,341 c 4280 47 56,549 11 4380 14 22,783 q 4480 21 20,712 212x 6080 12 19,263 Picea east 21 2R 3090 63 1,405 a 3190 18 819 b 3290 34 492 212T b 5190 29 306 west 290 37 8,934 212S g 4090 24 1,656 c 1 4290 23 8.445 alSubsection codes are combined with section codes to be consistent with those found in the ecoregions map constructed by Cleland et al. (2005b). 51 assigned to the corresponding pixels to create the map of hare habitat quality. In this way, grid pixels would have a similar proportion of HSI values on the landscape to that as measured by the plot data, within each ecoregion-ELU. I examined the effects of the random pixel assignment on lynx model output with a sensitivity analysis. Hare habitat quality (indexed by SCUs) was randomly assigned to pixels within 20 separate grids of the UP. Mean quality within a 4.5-ha moving window was calculated for each, and a grid of the standard deviation for each pixel was created to identify a sample area (1,200 kmz) with high variability. This sample area was input to the lynx model for each of the 20 grids, and 3 grids were chosen to be assessed by the lynx model across the entire study area. These 3 grids represented the range in outputs calculated for the smaller sample area, and covered the variability in hare habitat quality and subsequent lynx foraging quality that resulted from the random pixel assignment. Non-forested ELUs were not represented by FIA plot data and were assumed to provide no habitat to hare and lynx, with the exception of upland and lowland shrub cover types, which were assigned a mid-range HSI score of 50 for hare. These areas could provide refugia for hare, especially in lowland thickets of alder and willow, but they would not qualify as habitat for lynx, based on the lack of overhead cover. Therefore, the effect of the arbitrary assignment of quality for shrub cover types was assumed to have little effect on the estimation of lynx habitat potential. Satellite imagery — I obtained spectral imagery taken by the Enhanced Thematic Mapper Plus (ETM+) of the Landsat 7 satellite, with the intentions of predicting a continuous surface of stand attributes to improve upon the mapping of local variation provided by the random simulations (Appendix E). Eight images were required to cover 52 the entire UP, with dates of acquisition ranging between late summer and early fall of 2000-2002. Images were rectiï¬ed and georeferenced to zone 16 of the UTM system, and mosaics were created where possible. The FIA data that was available at the time of analysis included surveys conducted between 2000-2003 for the 6th inventory of Michigan, resulting in 3,809 forested plots. A portion of the modeling was conducted by Geoff Holden at the North Central SDS due to the conï¬dentiality of F IA data. The F IA plot locations were overlaid with the spectral imagery and a k-nearest neighbor (KNN) analysis was used to impute plot attributes continuously across the landscape. The KNN method assigns attribute values to non-sampled pixels from those that are sampled (i.e., located underneath a ground plot) based on the distance between the pixels in a multi-dimensional feature space deï¬ned by the combination of spectral values from each wavelength band. The analysis was performed using functions within ERDAS Imagine; a 3x3 mean ï¬lter was applied to the spectral imagery, of which bands 1-5 and 7 were used with k = 5. The prediction accuracy of the KNN model was estimated by holding out 40% of the FIA plots from the training set, and comparing the observed attribute values with those predicted by the model. The estimated accuracy of the KNN model was poor throughout the UP (Appendix E); as a result, I did not use the stand attribute maps for modeling lynx habitat. The proportion of explained variation (R2) in the predictions ranged from 0.01 to 0.22, with most of the UP having an R2 below 0.10 (Appendix E). The poor accuracy could have been caused by several factors, including limitations in the availability of cloud-free Landsat 7 imagery for the entire region, the use of one season of imagery as opposed to multiple seasons, the sampling intensity of FIA plots, and the lack of a relationship 53 between spectral values and stand attributes. Several alternative approaches to the KNN method, including the use of generalized additive models (Frescino et al. 2001), may have yielded better predictions of forest structure. Constraints on the access to F IA plot locations and the amount of time available prevented any further analyses. Presettlement Ecological Conditions The range of conditions that existed in the UP during the centuries prior to European settlement was estimated using spatial data that described ecological characteristics of the landscape and a simulation model that reproduced the effects of major disturbance on stand conditions. The presettlement forests exhibited a “shifting mosaic†of seral stages that was driven by disturbance and succession (Frelich 2002); therefore, it was my objective to capture a range of forest conditions, and resultant lynx habitat potential, provided by the landscape under a natural disturbance regime. Habtypes were mapped across the UP by combining data from the STATSGO soil layer and the original presettlement vegetation map from Comer et al. (1995), using similar methods to those applied for current conditions (Table 5), aside from the exclusion of the IFMAP layer. Seral stages were simulated across the UP with the landscape age-class demographics simulator (LADS), a stochastic simulation program that mimics the temporal and spatial patterns of ï¬re disturbance and subsequent forest age classes across a landscape through time (Wimberly et a1. 2000). Historical ï¬re regimes in the UP were estimated with a map developed by the Forest Service, similar to one that was created for northern lower Michigan (Cleland et al. 2004). The simulation output provided maps of age classes within habtypes, and formed the basis of ELUs for the presettlement era. Habitat quality of the ELUs was estimated and the lynx model was 54 used to predict the historical variability of lynx habitat potential prior to European settlement. The methods used to generate input for the LADS program and estimate habitat quality for ELUs included the following: Simulating forest age classes — A detailed description of the LADS program can be found in Wimberly et al. (2000), and I obtained the program directly from Michael Wimberly (Associate Professor, Department of Geography, South Dakota State University). The program simulates the spatial variability in ï¬re regimes across large- scale regions of macroclimate and small-scale changes in topography and vegetation (Wimberly et al. 2000). The required input involves 4 raster layers of landscape characteristics and a set of parameters for determining the properties of the simulation and the behavior of ï¬re within and across scales. The raster layers include: (1) climate zones, which delineate regions with speciï¬c ï¬re rotations; (2) topographic zones, which portray local differences in topography that affect susceptibility to ï¬re; (3) summary zones, which identify optional subsets of the study area for which summaries are desired; and (4) a buffer zone, which identiï¬es areas from which ï¬res can propagate but landscape summaries are not calculated. The input parameters include the total length of the simulation, length of bum-in period, length of summary steps, ï¬re shape calibration, wind direction and intensity, ï¬re susceptibility within topographic zones, age-class deï¬nitions, length of ï¬re rotation, and the mean and standard deviation of ï¬re size. Wimberly et a1. (2000) designed the program for analysis of large landscapes (>1,000,000 ha) with a coarse pixel resolution (>4 ha); the raster layers were aggregated to a resolution of 180 m as a compromise. Also, islands were masked out of the maps 55 due to the cellular automata approach of ï¬re spread in the program, which did not allow for disjunct pieces of land. The most critical input map was that of the climate zones, which determined the ï¬re regime parameters for a given location. The input map of climate zones for the UP was based on a spatial layer obtained from the Forest Service (Dave Cleland, USFS, Rhinelander, Wisconsin), which delineated the landscape into “biophysical units†(BPUs) based on ecological characteristics that affected ï¬re susceptibility and spread (i.e., landform, lake density, soil texture and drainage) (Cleland et al. 2005a). Cleland et al. (2005a) assigned each BPU to 1 of 6 stand-replacing ï¬re rotation categories, and used GLO survey notes to quantify the natural frequency of stand-replacing ï¬res within each BPU. This work was an expansion on the analysis of northern lower Michigan (Cleland et al. 2004), and BPUs were mapped across the entire Laurentian Mixed Forest Province, with natural ï¬re rotations summarized by BPU and ecoregion section (Cleland et al. 2005a). The accuracy of ï¬re rotation estimates decreases with the size of the landscape unit (Cleland et al. 2005a); therefore, I condensed BPU-section groups that had similar ï¬re ï¬equencies, and recalculated ï¬re rotations based on the ï¬re frequency and area of each new BPU (Figure 7). The mean and standard deviation of ï¬re sizes in the UP were estimated to be 342 ha and 1,347 ha, respectively (Cleland et al. 2005a). The LADS program used rotation length and ï¬re size parameters to compute the mean number of ï¬res per decade within a climate zone, which was simulated for each decade as a Poisson random variable. Individual ï¬res were randomly sized according to a lognorrnal probability distribution, calculated from the input ï¬re size parameters. 56 9| N N \ \ 98> E QEZV 25:88 Bu ESE: mo 886ch @8383.“ 53> woo—m .2558 no Ammoom an no cum—20v 3.8: 306335 3:86 .«o @392on Goa 3:88 mason 088:0 .«o "~23“:me 35QO .N. 0.59m E280 nooéoo can £55 «and I 33. I 22 I 2.2 a :3 I to D can U S. n ma E 8. I 9. I 8 I an... 5.2 Twc \ Exch— 8- 9. mm “ x A? a . 9 \m c-"b ,5 57 I created the topographic zones with a 30—m US. Geological Survey digital elevation model (DEM) by calculating grids of flow accumulation and aspect within ARC/INFO 8.1. These functions identiï¬ed local ridges (having a flow accumulation = 0) and areas that were flat or had southern aspects, which were assumed to have a higher susceptibility to ï¬re (Zhang et al. 1999). The combination of the 2 grids resulted in the following 4 topographic zones (in order of increasing susceptibility): local valleys with slopes not facing south, local ridges with slopes not facing south, flat or south-facing local valleys, and flat or south-facing local ridges. The 30 m grid was aggregated to 180 m resolution using a majority ï¬lter to match the other input maps. The summary zone was set to all forested areas within the contiguous boundary of the UP (excluding islands). The buffer zone included a 20 km strip along the border of Wisconsin to accommodate the propagation and spread of ï¬res from outside the UP. Age classes were deï¬ned by the amount of time since the last stand-replacing ï¬re; I deï¬ned 11 intervals with the following endpoints: 10, 20, 30, 40, 60, 80, 100, 150, 200, 250, and >250 years. Fire susceptibility was assumed to be uniform across age classes, and ï¬re severity was set at a constant high level to simulate only stand-replacing ï¬res. Wind direction was set at 135° (northwesterly) with a variability of 90° and a moderately low intensity (3 on a scale of 1-10). Simulation length was 10,000 yrs, with a 1000-yr bum-in period, a 10-yr run interval and a 50-yr summary interval. The output from the LADS program included maps illustrating the distribution of age classes across the UP for each of the 200 summary intervals, as well as a table of age-class proportions for each. Estimating habitat quality within age classes — The vegetation attributes necessary to quantify habitat quality in the lynx model could not be directly measured for 58 the presettlement conditions. Therefore, estimates of quality for age classes within habtypes were developed using information from the literature and simpliï¬ed assumptions about the structure and composition of vegetation. Hare habitat quality was influenced by age class and habtype by assigning high scores to early age classes (<40 yrs) and habtypes with a conifer component, both characteristics that have been consistently described as important for snowshoe hare (see review in Hodges 2000). For the 0-10 year age class, all habtypes other than PiVa were assigned a bare HSI score of 50, to account for variability in regeneration of aspen-birch within the ï¬rst 10 years of a severe ï¬re. The PiVa habtype was assigned a score of 100 for the 0-10 year age class, assuming that jack pine regeneration would be dense after several years (Frelich 2002). The 10-20 and 20-30 year age classes were assigned a score of 100, under the assumption that these age classes would provide the greatest amount of understory cover for hare. The 30-40 year age class was assumed to be the last suitable stage before stem exclusion would result in a sparse understory; it was assigned a score of75. For the QuAc, TsMa, and TsTh-dry habtypes, the age classes between 40-150 yrs were assigned a score of 25; these habtypes contained a potential spruce-ï¬r and/or cedar component that might have provided cover for hare. The PiVa was no longer suitable after 40 yrs due to the assumption that poor soil quality would result in a relatively sparse mature forest. The Acer and Frax habtypes were also unsuitable after 40 yrs, since dominance by hardwoods was expected to provide inadequate understory cover during the winter, even in areas where regeneration of shade tolerant species was high (i.e., sugar maple in the Acer habtype). The QuAc, TsMa, and TsTh-dray habtypes became 59 unsuitable after 150 yrs, when shade tolerant hardwoods would likely attain dominance in the absence of ï¬re (F relich 2002). The TsTh-wet and Picea habtypes were assigned a score of 50 for all age classes >40 yrs, since they were dominated by spruce-ï¬r and cedar; it was assumed that multi-aged mature stands within these habtypes would provide pockets of suitable cover for hare. Non-forested areas provided no suitable hare habitat, aside from shrub types, which were assigned an HSI score of 50. The assignment of HSI scores for the denning and interspersion components of the lynx model were even more simpliï¬ed than that for hare quality. Habtypes that provided potential denning sites based on mesic soil conditions included QuAc, TsMa, Acer, and TsTh-dry; within those habtypes, age classes >60 yrs qualiï¬ed as denning sites. Travel cover was assumed to be provided by age classes >10 yrs within all habtypes. All non-forested areas were considered unsuitable for lynx. Based on my methods of assigning HSI scores, habitat suitability for lynx during the presettlement era would be most influenced by fluctuations in the amount of early age classes providing highly suitable hare habitat. Therefore, the total proportion of age classes between 10-40 yrs on the landscape was calculated for each of the 200 simulated maps (Figure 8), and maps representing the highest and lowest proportions were selected for analysis with the lynx model. This allowed us to examine the range in habitat suitability that might have occurred under the presettlement disturbance regime. The resolution of the grid maps was increased to 30 m after assessing hare habitat quality so that habitat units could be sufï¬ciently accumulated into hare home ranges using Plume’s (2005) GIS ï¬mction, though I recognize that this procedure did not increase the spatial accuracy of habitat components. 60 g 0.09- Ln "cg 0.08- <—-high (0.081) .2 a 0.074 8 g 0.06- >3 ":3 0.054 LL] ‘8 0.04 8 g 0'03†<—low (0.028) e 0.02- Q.‘ 1 1 1 l 1 l 1 I I 0 1000 3000 5000 7000 9000 Sirmlation Year Figure 8. Fluctuation in the proportion of early successional forest (10-40 yrs old) over a 10,000-yr simulation of severe ï¬res in the Upper Peninsula of Michigan, with arrows identifying highest and lowest years. Fire was simulated with the Landscape Age-class Dynamics Simulator under the natural disturbance regime that occurred prior to European settlement. 61 ANALYSIS OF MODEL OUTPUTS The lynx model generated 5 maps of habitat quality for each version of the ecological land classiï¬cation, representing indices for snowshoe hare habitat, the 3 life history components for lynx (foraging, denning, interspersion), and overall lynx habitat. Each map described habitat quality on a scale of 0-100, ranging from unsuitable to optimal conditions. There were 3 versions of the ecological land classiï¬cation: 1 for current conditions and 2 for presettlement conditions. The assessment of current conditions involved several variations for describing snowshoe hare habitat on the landscape, though each utilized the same ecoregion-ELU grid. Three versions of hare habitat quality were created from the separate HSI models (HSI-1, HSI-2, HSI-3), and another 3 were created for testing the sensitivity of random pixel assignment. The outputs of lynx foraging quality were compared within each set and one representative map was chosen to calculate overall lynx habitat and illustrate current conditions. Presettlement conditions were represented by 2 simulated landscapes of forest structure resulting from a natural disturbance regime. Each version had a separate ELU grid describing the arrangement of age-classes resulting from ï¬re within each habtype, and subsequently, separate inputs and outputs for each of the habitat quality indices. Pairwise comparisons were made within and between the outputs from current and presettlement conditions using several statistics. The pixel by pixel differences between paired maps was quantiï¬ed by the Volume of Intersection Index (VOI) (Seidel 1992, Millspaugh et al. 2000). The VOI is a measure of the overlap between utilization distributions, with values ranging from 0 (no overlap) to 1 (direct overlap). This index 62 provided a measure of the similarity in the location of suitable habitat between a pair of maps. Each habitat quality map was converted into a utilization distribution, with HSI values weighted so that the total volume of the utilization distribution was equal to 1. This procedure removes the magnitude of difference in habitat quality between paired maps, so the mean and standard deviation of the differences were also calculated. In this way, if a pair of maps had a similar spatial arrangement of suitable habitat, but one suggested a much higher level of suitability, the V01 score might be close to 1, while the mean difference would be relatively high. Likewise, if a pair of maps described suitable habitat at vastly different locations, but with a similar frequency of occurrence, the V01 score would be close to 0 and the mean difference would be low. Comparisons were also made among the outputs of simulated lynx home ranges created by Plume’s (2005) GIS function, in regards to the number and location of viable home ranges across the study area resulting from each habitat quality map. Lynx home ranges were simulated 5 times for each map; the mean number of home ranges was calculated for comparison, and an average map of home ranges was chosen for illustration. Prior to the comparative analyses, each output map of habitat quality was aggregated to a 90-m grid resolution and clipped to the same extent. All processing was completed within ARC/INFO 8.1, aside from the simulation of lynx home ranges, which used a custom program written by Plume (2001). RESULTS CURRENT ECOLOGICAL CONDITIONS The assessment of bare habitat quality, and subsequent lynx foraging quality, was affected by the choice of HSI model for current conditions (Figure 9). Horizontal cover 63 Hare Habitat Qualitygl l Lynx Foraging Quality (a) (d) 0 20 40 60 80 100 0 20 40 60 80 100 0)) Map Proportion 0 20 40 60 80 100 (C) 00 ulâ€._r_E_ . 0 20 40 60 80 100 0 20 40 60 80 100 Habitat Quality Figure 9. Frequency distributions of quality scores for hare habitat quality and lynx foraging quality across the study area, as assessed by HSI models measuring horizontal cover with manipulated, HSI-1 (a, d) and original, HSI-2 (b, e) seedling distributions, and the HSI model measuring stem cover units, HSI-3 (c, i). 64 within a stand was likely overestimated when the seedling size distributions were manipulated, resulting in an overly optimistic portrayal of habitat quality by the ï¬rst model (HSI-l). The description of habitat quality was more similar between the HSI model measuring horizontal cover with original size distributions (HSI-2) and the one measuring stem cover units (HSI-3) (Figure 9). Volume of Intersection (VOI) indices for lynx foraging quality ranged from 0.79 to 0.83, indicating that the HSI models were locating suitable habitat in similar areas. Mean differences in lynx foraging quality were high when comparing output from HSI-1 to both HSI-2 (mean = 43.84, SD = 18.07) and HSI-3 (mean = 39.77, SD = 18.3), and low when comparing HSI-2 to HSI-3 (mean = - 4.07, SD = 12.42). I chose to retain the output from HSI-3 for further analysis, given the high uncertainty in the measurements of horizontal cover, but also, the relative agreement between output from HSI-3 and HSI-2. The random pixel assignment had little effect on the assessment of lynx foraging quality throughout the study area. The VOI index was 0.94 for each comparison, and mean differences in foraging quality among the 3 grids ranged from -0.60 to 0.40, with standard deviations <4.35, suggesting a high degree of similarity. Thus, local differences in hare habitat quality among the simulated grids did not translate into signiï¬cant differences in lynx foraging quality. The grid with a mean quality most typical of the 3 was chosen to represent lynx foraging quality for the current conditions. The amount and distribution of suitable habitat differed greatly between each output of the lynx habitat components. Hare habitat quality and thus, lynx foraging quality, was mostly poor throughout the study area (Figures 10, 11), with the highest suitability (HSI = 40-60) in lynx foraging occurring in the eastern portion of the UP. The 65 _ E \ ‘1 .88 >95 05 $98 $88 5:35 we 8:588 2: «0 EEon: a 8353: 82: 2F .mooméoom 5033: 53:22 mo 3335: 5.53 05 E 3:56 ESE 8mm .0: Eswï¬ 87;- 21;- 27:- 3-....- 8-3- ELVD Stan 2-5! 85“ 25' a: ed wd o.— 66 \ H .~ .88 Exam 2: £88 388 aims—u mo GOP—omen 2: mo EEonmE a 8:232: coma: 251 .mooméoom 563$: newEBE mo 233:8: Saab 05 E 3:33 wï¬wï¬o: x83 .: 85Eâ€: 8.5! 2..-: I E-.. - 5: ... t ed wd o.— 67 poor suitability in hare habitat was evident from the F IA data, as 74% of surveyed plots had an HSI score of 0 for hare, while 93% had an HSI score <60. The only habtype with >10% of its FIA plots scoring >60 for hare suitability throughout the UP was the TsTh-wet habtype (23%). Within the eastern ecoregion subsections 212Ra and 212Rb, 42% of FIA plots representing the TsTh-wet habtype scored >60 for hare suitability; the location of these subsections corresponds with the areas of the eastern UP that contained the highest suitability in lynx foraging across the study area. Lynx denning quality was highly suitable throughout most of the UP (Figure 12), as only a few areas (<10% total area) were at an unsuitable distance (>400 m) to the vegetation types that provided denning areas. The high suitability for the denning component was due to the fact that over 53% of FIA plots met the requirements for den sites according to the lynx model, and combined with the high heterogeneity in habitat types across the UP, den sites were always within a suitable distance from a given point on the landscape. The interspersion of non-lynx habitat was highly variable in the study area (Figure 13), with 20% of the total area having an HSI score <10 and 36% having an HSI score = 100. The extensive areas of non-forested wetlands, agricultural land, and human developments resulted in low levels of interspersion suitability across the eastern UP, while the contiguous forests of the western UP resulted in a high suitability with pockets of low suitability scattered evenly throughout. Overall lynx habitat quality was mostly mid-range across the UP, with 49% of the total area having an HSI between 40 and 70 (Figure 14). Quality was limited by the high amount of non-suitable (H S1 = 0) habitat according to the interspersion component, and 68 k H F .83 >95 05 $82" $88 5:26 mo 5anon 2: mo EEmonE m 8352: 62: 2:- .mOCN-OOON auction 5&522 mo EsmEcom $95 2: E 3:125 wEESu ES: N: 059: 2: - 3- 21;- :x- E I E.- 3' cc- .4. I 3:3.qu EtzmuHU 3-2! :7: H 25' Em: 1 ed r wd re.— .4 69 ~ b‘ ‘ don-n >95 05 $88 $83 3:36 mo coma-Sacha 2: mo ESwSE: a 8:52;: BE: 2: .mooméoom 5256: :meEE mo Esmanm 6.53 2: E 3:36 €63.5an x83 .2 car: 2: - 3 I .. 21; I \ ex- : I III'IIII E.- _c I .3 =6- 7- I 3. - .4 D to s.- :. WU 9o 3.- _m a su- : H no 2 -c I Q. 5: - 70 \ .\ :. a 3353: 8%: EH .mooméoom 50389 53:22 mo 255:8 6an 2: E 5:85 55s: a. =Eo>O .3 “uh—mi .88 >95 05 $88 880m @626 we aoEomoa 05 mo EEwSmE a . v E; c2 8. cm mm o 71 the generally low suitability of foraging habitat. A mean of 603.8 viable lynx home ranges (SD = 26.3) were simulated for the current conditions, with locations throughout most portions of the UP (excluding the islands) (Figure 15). Viable home ranges were concentrated in 6 of the 21 ecoregion subsections (212Jb, 212Ra/b/e, 212Sn/q) covering nearly half of the UP; only 5% of viable home ranges were located in other subsections. PRESETTLEMENT ECOLOGICAL CONDITIONS The small size and relatively infrequent occurrence of ï¬re throughout most of the UP resulted in a mature forest with little early successional vegetation on average. Over 97% of the study area had an average ï¬re frequency >100 yrs, and 75% had an average frequency >500 yrs. The greatest amount of early successional vegetation (IO-40 yrs old) in a simulation year was still only a small percentage (8.1%) of the total study area (Figure 8), but nearly 3 times the magnitude of the lowest amount. The main difference between the 2 simulation years (high, low) that were chosen to capture the range of presettlement forest conditions was a large conflagration that occurred in the eastern portion of the UP, which created >150,000 ha of early successional vegetation. As a consequence, hare habitat quality and lynx foraging quality fluctuated signiï¬cantly between the high and low simulation years (Figures 16, 17). The large differences in lynx foraging quality were not exhibited by the other lynx habitat components. Denning quality was only slightly decreased by the large expanse of early successional vegetation created by the conflagration in the high simulation year (Figure 18). The conflagration did not result in one contiguous patch of early successional vegetation; it contained many islands of mature forest that were saved from the burn, and thus, did not present a complete loss of areas for denning. Interspersion 72 N _ N N & .onn bzmsc a wEBï¬ me @2526 05>» mownï¬ qun 2%?» .325 553 wqunwocwwm com corona m5 $88 08:3 0% was: Govï¬oxo mums—m5 ng22 mo Esmï¬aom Saab 05 E 22:28 Beta .28 toga—=86 muwï¬u 080: 55 05mg .2 Emmi .~ 73 8* I w, «is: â€I I, HSI . {/7 Eo-Io “‘ ,rf ) n -20 { 5., 53:21-30 L N [:]3I-4o A Izc—Z Qu-So o 25 so 100 l50km - 5l-60 -6l-70 — -71-so -x|-90 I -9l-l()() .U _ _ ~ We?) LOW Figure 16. Hare habitat quality during the 2 simulation years of presettlement ï¬re disturbance with the highest (top) and lowest (bottom) amount of early successional vegetation. In the top map, a conflagration responsible for creating a large area of early successional vegetation can be seen in the eastern UP. 74 Illn-_-___ "“i-‘ «574 q HSI k -0-|0 ‘ 0Q fl? 3 11-20 { r†{51 2-1-30 _ N -3l-40 -:-_=1 -4l-5() A o 25 50 mo ISOkm -51-60 -6l-70 — 4 -7l-80 -8l-9() -9l-l()0 “ (â€D Figure 17. Lynx foraging quality during the 2 simulation years of presettlement ï¬re disturbance with the highest (top) and lowest (bottom) amount of early successional vegetation. A 75 N [;:_ [31-40 A .12—21 C] 4| -50 0 25 so 100 l50km Figure 18. Lynx denning quality during the 2 simulation years of presettlement ï¬re disturbance with the highest (top) and lowest (bottom) amount of early successional vegetation. 76 quality was nearly the same between simulation years (Figure 19), as only non-forested vegetation types, water bodies, and the relatively few recent burns (<10 yrs) contributed obstacles to lynx in the presettlement landscapes. Over 90% of the study area had an HSI score >80 for denning and interspersion quality in each of the 2 simulation years. Similar to the current conditions, overall lynx habitat quality during the presettlement era was limited by foraging quality. The majority of the study area had an HSI score between 40 and 70 for the high (61%) and low (64%) simulation years (Figure 20); HSI scores >70 occurred across 7% of the study area for the high year, and 3% for the low year. In each simulation year, the highest amount of suitable habitat occurred in the eastern portion of the UP, where average ï¬re frequencies where 3 times the magnitude of those in the west. A mean of 677.8 viable lynx home ranges (SD = 8.6) were simulated for the high year, while 264.8 (SD = 7.1) were simulated for the low year (Figure 21). Home ranges were mostly located in the eastern portion of the UP, and the major difference between simulation years occurred at the conflagration, where a high concentration of home ranges appeared during the high year. DIFFERENCES BETWEEN CURRENT AND PRESETTLEMENT CONDITIONS The differences between current and presettlement conditions were unique to each lynx habitat component. These differences reflected the spatial properties of vegetation types across the landscape, as well as the structural attributes within them. The VOI for lynx foraging quality between the current and presettlement maps ranged from 0.56—0.59, with a mean pixel difference of 2.87 (SD = 23.01) between the current and high presettlement year, and 6.98 (SD = 18.16) between the current and low presettlement year. Although the average foraging quality was slightly higher for the 77 . (I HSI Ff“ _ o - I0 N III—2 A 02550 I00 150km LOW 33‘ -_‘_ \- _ V â€\ 1?“) (V ft“? l. 0 -5r f- a} n f < I" h. 3 2'1! 4 =ll-20 @2140 @3140 EM -50 -51-60 -6l-70 -7l-80 -81-90 -91-1(m Figure 19. Lynx interspersion quality during the 2 simulation years of presettlement ï¬re disturbance with the highest (top) and lowest (bottom) amount of early successional vegetation. 78 In It 13..-- ‘ HSI (f7 =0- Io N" 0 if} â€-20 { O -. air-2140 _ N |:j3I-40 A .1—21-41-50 0 25 50 mo ISOkm -51-60 -6l-70 ~— -7|-80 LOW ‘ -m-9o a -9I-I00 I_.I_l‘J-__ gr . _ . V I“. P? if? UHâ€) iï¬/ t Figure 20. Overall lynx habitat quality during the 2 simulation years of presettlement ï¬re disturbance with the highest (top) and lowest (bottom) amount of early successional vegetation. 79 at u "lair Figure 21. Locations of viable lynx home ranges according to habitat quality for the 2 simulation years of presettlement ï¬re disturbance with the highest (top) and lowest (bottom) amount of early successional vegetation. Habitat units were aggregated into simulated home ranges using Plume’s (2005) GIS ï¬mction. 80 current conditions than those during presettlement, the V01 scores indicate a large difference in the location of suitable habitat for each time period. The VOI between each simulation year for presettlement was 0.72, with a mean difference of 4.1 (SD = 19.88), illustrating the changes caused by the conflagration. The VOI for lynx denning quality between the current and presettlement maps ranged from 0.94-0.96, with mean differences of 4.63 (SD = 20.99) and 2.74 (SD = 19.09) between the current map and the high and low maps, respectively. The VOI for denning quality between the presettlement simulation years was 0.97, with a mean difference of -1.89 (SD = 15.79). Lynx interspersion quality was the same between presettlement simulation years, but very different when compared to current conditions. The VOI for interspersion was 0.70, and the mean difference between current and presettlement maps was -35.55 (SD = 39.13), indicating a large overall decrease in interspersion quality from the presettlement era to the present. The VOI for overall lynx habitat quality ranged from 0.68-0.69 for current and presettlement conditions, with mean differences of 0. 14 (SD = 33.25) and 1.51 (SD = 31.66) between the current map and the high and low maps, respectively. The number of simulated lynx home ranges for the current conditions were within the range exhibited by presettlement conditions (Figure 22), though the locations of the home ranges were very different (Figures 15, 21). Under present conditions, home ranges were grouped in several locations throughout the UP, including the far western region near the border of Wisconsin. This region in the west is void of lynx home ranges under presettlement conditions, due the lack of foraging habitat. 81 Home Range Count l—Current—l I P “' ‘ A Figure 22. Number of viable home ranges simulated in each of the 3 lynx habitat quality grids, using the Plume (2001) GIS function to aggregate habitat units. Viable home ranges were deï¬ned as having an average quality >70. 82 DISCUSSION Roloff’ s (2003) model quantiï¬ed habitat suitability for Canada lynx from the scale of an allometric home range (250 ha), at which habitat selection for a species is hypothesized to be strongest (Roloff and Haufler 1997). The model consisted of 3 habitat components representing life requisites that influenced the ability of a given landscape to support an individual. I integrated multiple forms of data within a GIS to quantify the spatial arrangement and structural attributes of vegetation types for assessing lynx habitat suitability under current and presettlement forest conditions in the UP. For the analysis, the large size of the UP made it a logistical obstacle to collect data that could adequately describe the vegetation attributes necessary for the lynx model. I used multiple spatial layers portraying biotic and abiotic characteristics of the landscape to delineate the study area into ecological land units (ELUs), representing vegetation types that had similar composition and structure. The ELUs were based on a previous ï¬eld study which illustrated relationships between tree species assemblages through succession and soil/landform properties in the UP (Coffman et al. 1983). The use of “habitat type†mapping increased my ability to stratify forest conditions beyond what could be provided by simple land-cover maps. The structural attributes for ELUs were obtained in 2 different ways: for current conditions, a region-wide inventory of forest stands conducted by the F IA program provided stand-level measurements of vegetation; for presettlement conditions, a simulation program developed for modeling the effects of ï¬re disturbance on the spatial distribution of forest age-classes was used to estimate the range in forest conditions that might have occurred prior to European settlement. The synthesis of these data provided the best available estimates of both present and past 83 stand conditions in the UP, and allowed us to determine the nature of lynx habitat during both time periods. The current distribution of lynx foraging habitat according to the model is very different from that which existed prior to European settlement. The output from both time periods indicated mostly poor quality foraging throughout the UP, due to a lack of early successional (10-40 yrs) vegetation types providing the understory cover required for a high abundance of snowshoe hares. According to the 2000-2003 FIA data, the ages of forest stands throughout the UP were normally distributed with a mean of ~60 yrs; approximately 15% of forest stands were at an age between 10-40 yrs. Using a 10,000-yr simulation of natural ï¬re under a presettlement disturbance regime in the UP, the greatest amount of early successional vegetation in a given year was estimated at 8%, and forest stands had a mean ï¬re return interval of 1,250 yrs. While the current landscape has been drastically altered by human actions, it would appear that logging has actually improved conditions for snowshoe hare and lynx in the UP, given that presettlement disturbances may have been too infrequent to provide adequate amounts of high quality habitat for hare. Both timber harvesting and severe ï¬res have the ability to create adequate early successional vegetation types for hare, and each has been shown to beneï¬t lynx in the long term (Fox 1978, Litvaitis et al. 1985, Koehler 1990). It would appear that current timber harvesting strategies do not occur at a large enough scale and high enough frequency to support high densities of hare and a viable lynx population in the UP, and neither did the presettlement ï¬re regime. One example of a noticeable difference between current and presettlement foraging quality occurred in subsection 2121b of the western UP, where quality was much higher under current conditions, especially for the 84 Ottawa National Forest. In this ecoregion subsection, sugar maple forest types had twice the proportion of F IA plots with a hare HSI >60 than found elsewhere in the UP. Reasons for the increased suitability are unclear, though nearly half of the highly suitable F IA plots had a record of silvicultural treatment in this ecoregion. The forests of the presettlement era within this ecoregion were dominated by late successional sugar maple, which was assumed to be of poor suitability for hare. Denning quality was shown to be adequate throughout the UP, both during the presettlement era and in the present. The denning requirements for lynx were once thought to be strictly associated with late successional forest (Koehler and Aubry 1994), but Slough (1999) found a wider range of conditions, including those found in regenerating forests, provided suitable denning sites for female lynx. The deï¬nition of suitable habitat in the denning component of the model reflected the need for older forests, but not necessarily late successional, and therefore nearly half the forest stands according to the FIA data met the model requirements. Forest conditions according to the presettlement simulation were expected to support an abundance of denning opportunities, given the low frequency of disturbance and subsequently, the relatively high amount of mature forest across the landscape. Interspersion quality was the habitat component that changed most signiï¬cantly between time periods. Given that the presettlement landscape would only contain non- habitat where natural features such as water bodies and non-forested wetlands occurred, as well as recently burned forest stands, it was reasonable to assume that conditions for traveling would not be limiting during that time period. The amount of non-habitat has increased considerably since presettlement due to the conversion of forest into 85 agricultural land throughout large portions of the eastern and southern UP, as well as the development in and around cities such as Sault St. Marie and Marquette. Another region of non-habitat for lynx in the eastern UP is the Seney National Wildlife Refuge, an area once dominated by forested swamps that was mostly converted into a managed wetland after extensive logging in the late 19th century. The combination of these factors has resulted in an overall decrease in the contiguity of forested vegetation types that are important for lynx traveling needs, most notably in the eastern UP. A number of assumptions were required for modeling current and historic forest conditions in the UP that may have effected the results. For current conditions, the F IA data represented a point sampling of forest structure, and therefore, the mapped ELUs did not have attribute information that could be considered accurate at the pixel level, especially for the ï¬ne scale structural measurements. This was evident from my method of simulating a random pixel distribution, which inferred that variation in structure was occurring at a distance of a 30 m (pixel size) with no regard to spatial autocorrelation. Given the systematic sampling which results in F IA plots being separated by potentially long distances, and my inability to obtain accurate plot locations due to conï¬dentiality, the scale at which forest structure varied could not be validated. I caution that the prediction of stand variables not be used to assess forest conditions at scales less than that of townships (~92 kmz). The ï¬ltering of F IA plots based on condition and canopy cover requirements reduced the dataset by 30%; an examination of the distributions of stand variables did not ï¬nd signiï¬cant differences. Even so, by classifying ELUs into discrete mapping units I eliminated the potential importance of ecotones, which can effect the assessment of forest conditions in a heterogeneous landscape. This might have had the 86 most impact in sparsely forested spruce bogs, where pixels of forest (>25% canopy cover) are scattered throughout lowland shrub vegetation, and any FIA plots located within would most likely have had <25% canopy cover and been ï¬ltered out. A more intensive analytic approach to classifying the landscape would be necessary to account for ecotones. The varying accuracy of the digital maps, especially that of soil types, may also have produced errors in assessing habitat types across the landscape. The use of STATSGO soils data should be limited to large scale analyses, and given the near completion of the SSURGO dataset for Michigan (USDA, NRCS 2006), its utility may become obsolete in the future. The accuracy of the land cover dataset made it nearly impossible to distinguish canopy species with any conï¬dence (Appendix C). New approaches to classifying land cover using combinations of high and low resolution imagery collected across multiple seasons would improve the ability to predict more detailed vegetation attributes. For the presettlement conditions, many assumptions were made when modeling natural ï¬re disturbance in the UP. Fire rotations were estimated using 10+ years of data collected during the GLO surveys of the 18003; this small sample may not have been adequate to accurately document ï¬re disturbances in the region. Also, I assumed that ï¬re disturbance regimes did not change over time, an assumption that could be violated with broad-scale changes in climate (F relich 2002, Cleland et al. 2005a). The interaction of many factors results in the propagation and persistence of tree species assemblages in a given area, not all of which can be captured by modeling the physiography of a region (Frelich 2002). I did not model disturbances other than those created by severe ï¬res, but 87 while small surface ï¬res, wind, and insects have all contributed to changes in forest composition and structure in the past, none have been of the same magnitude as that of severe ï¬res (Cleland et al. 2005a). Finally, I assumed that stand age could be correlated with habitat quality for snowshoe hare, which is not always true (Hodges 2000), as understory structure can be influenced by factors other than stand age. Despite these caveats, I feel that my methods provided an adequate way to estimate habitat quality for lynx over a large study area and across 2 time periods. The current absence of lynx in the UP is contrary to the output, which suggested there were several areas throughout the UP and close to the eastern and western borders that were capable of supporting multiple viable home ranges for lynx (Figure 15). The model assessed habitat quality in terms of the vegetation attributes on the landscape that would affect the availability of resources necessary for survival; it did not incorporate other factors such as competition, which will invariably affect habitat quality for a given species. The presence of coyote (Canis latrans) and bobcat (F elis rufus) in the UP can reduce the availability of snowshoe hare for lynx, even in areas where hare habitat quality is high (Koehler and Aubry 1994). These predators are not found in high abundance across the northern boreal forests of Canada, and their presence, combined with the generally low densities of hare across southern boreal forests (Koehler 1990, Koehler and Aubry 1994, Hodges 2000), may ultimately create unsuitable conditions for lynx in the United States. It is difï¬cult to quantify the effects of these competitors on lynx foraging quality in the UP. Hoving et al. (2005) found that lynx occurrence in the Northeast was best predicted by annual snowfall, and suggested a critical threshold of 270 cm/year, below which the competitive advantage of lynx over other predators in the snow is lost. 88 Annual snowfall throughout more than half of the northern UP is >270 cm/year, especially within the snow belts (Eichenlaub et al. 1990), and 75% of the viable lynx home ranges that were simulated fell north of the snowfall threshold. Bobcat harvest records from the past 20 yrs indicate a close agreement with the snowfall threshold, as most harvests have been in areas of the southern UP with <270 cm/year (Figure 23). The trend of a warming climate may affect lynx habitat potential in the UP and beyond if snowfall decreases and the northward expansion of species such as bobcats and coyotes is facilitated (Hoving et al. 2005). Another factor influencing the absence of lynx in the UP may relate to their ability to disperse from Canada. Lynx are known to be capable of long-distance dispersals (>1,000 km) (Nellis et al. 1972, Mech 1977, Ward and Krebs 1985), and historical records indicate that lynx were once able to enter the UP during population eruptions in Canada (McKelvey et al. 2000a). In the early winter of 2003, a lynx was incidentally trapped in Hiawatha National Forest near the city of St. Ignace, which marked the ï¬rst sighting/trapping of an individual in more than 20 years for the UP (Stephen Sjogren, USF S biologist, personal communication); the origins of the trapped lynx were undetermined. In the eastern UP, habitat quality was poor throughout the northern end and annual snowfall was below the critical threshold in the southern end (Figure 23), suggesting this area could be a potential barrier to lynx dispersal from Ontario. If the lynx trapped in 2003 was a disperser from Ontario, it is obvious that migration across this area of the UP is not impossible. An examination of the connectivity of habitat between Canada and the areas of suitable habitat within the UP may provide insight into the current status of lynx in the region. 89 K x x a c.8380 chmv m. zflaocm 3:53 £053 323 @2258 2: 95858 on: peace 2:. .38-me 52209 cmeBE mo 235:8 Saab 05 E 8332 Boson mo 55:62.“ gum—om .mm oSwE Wu M .M\ U m1. .M/ .532 8. on n o .. . ¢ v V l .. ob. . .ntt < Z . . II . “mu ... . __. \\. - 4 o . 3 a . o ..\ f x m? 1g. h? ..\ xa ... t . t .\ . Bo .x. . 1,. . \ ..\ — I a f.» I k: \aa. 3.“ \\ ‘Qflm : . . T I.‘ 84¢“; _ . \KWV \ b. a. ... u \ .. .1 .. . I' a. .. I Kin \x t (3.11“ .w \\\\i £ME I .\ . 1).. ho=o=c°hm 1" .1 1“ ‘ l/ f ‘V\ \‘Kw k‘ \\ wigs—om x... v.t\\..\lt . 51 ./ ‘Mm\ K 7111* ((er / \ QCK‘ “If ... 90 CHAPTER 2: TEMPORAL CHANGES IN HABITAT CONNECTIVITY FOR CANADA LYNX IN THE UPPER PENINSULA OF MICHIGAN INTRODUCTION In listing the Canada lynx (Lynx canadensis) as a threatened species in the contiguous United States, the US Fish and Wildlife Service (USFWS) called for “accurate mapping of lynx habitat in the Great Lakes Region†(USFWS 2000:16057) to provide the location and distribution of resources necessary for the species’ persistence and identify potential areas for conservation and management. The habitat modeling framework proposed by Roloff and Haufler (1997) was used with a habitat suitability index (HSI) model (Roloff 2003) to quantify the Upper Peninsula (UP) of Michigan in terms of its ability to provide the necessary life requisites for lynx (Chapter 1). The analysis identiï¬ed areas of the UP that were potentially important for lynx persistence in the region, though no evidence of a resident lynx population is currently present (Beyer et al. 2000), and provided an historical context by assessing forest conditions and lynx habitat suitability in a presettlement-era landscape. The results of the analysis indicated that lynx habitat potential was high in several areas of the UP, both in the past and present, and it was hypothesized that factors other than the condition of forest communities were preventing a resident population from being established. One of the competing hypotheses was that habitat connectivity between the core lynx population in Canada and suitable habitat in the UP may be poor enough to preclude the successful migrations necessary for lynx establishment in the region. 91 Lynx are known to travel great distances in search of habitat that can support adequate foraging opportunities, especially during times when the local abundance of their primary prey, snowshoe hare (Lepus americanus), is extremely low (Nellis et al. 1972, Mech 1977, Ward and Krebs 1985). Mech (1977) reported a dispersal distance of 483 km in Minnesota, and long-distance movements of lynx in Canada have been estimated to exceed 1,000 km (Slough and Mowat 1996, Poole 1997). It is clear that distance alone would not prevent lynx from migrating out of the core population in Canada to areas of the Great Lakes such as the UP. Historical accounts of lynx presence in the UP include an influx of occurrences during the 19608, when an irruption of populations in Canada caused a large migration of individuals into the Great Lakes region (McKelvey et al. 2000a). The obstacle presented by the St. Mary’s River in the eastern UP cannot be considered insurmountable to lynx, given that migrations have been documented in the past, as well as the fact that winter ice cover on the river has averaged >90% over the past 30 years (Assel 2003). The ability, or apparent inability, of lynx to migrate into the UP may be more related to habitat quality than physical obstacles. The objective of this analysis was to examine the connectivity of lynx habitat across the UP during the past and present with a Geographic Information System (GIS). The outputs from a prior assessment of lynx habitat suitability (Chapter 1) were used to provide input for quantifying the “cost†associated with traveling from outside the UP to areas within that could potentially support lynx. Cost would be deï¬ned by the amount of energy expended on a dispersal, and/or the risk of mortality associated with conditions on the landscape (e. g., starvation, incidental trapping). Comparing the nature of least-cost travel corridors between current and presettlement landscapes might contribute insight 92 regarding the current absence of lynx in the UP, and identify areas that would be important for lynx conservation and management. STUDY AREA (see Chapter 1) METHODS I assumed that dispersal by lynx would be associated with areas of good habitat, and that the least-cost corridor between 2 points on the landscape would indicate the travel route providing the highest probability of survival by an individual. Determining the least-cost corridor required a map of habitat quality that could be converted into a cost surface, over which an accumulated cost could be calculated between grid pixels to identify corridors with the least resistance. Multiple cost surfaces were compared to quantify differences in the number and location of least-cost corridors, as well as the magnitude of cost calculated for each corridor. Maps of lynx habitat quality were created using the modeling framework proposed by Roloff and Haufler (1997) with an HSI model (Roloff 2003) constructed for lynx. A detailed description of the data and processes used to model lynx habitat can be found in Chapter 1. Lynx habitat quality was quantiï¬ed for 2 time periods in the UP, with 1 map for current conditions and 2 maps for presettlement conditions. The 2 maps for presettlement conditions represented simulation years with the highest and lowest amounts of early successional forest, capturing the range in lynx habitat quality exhibited under a natural disturbance regime prior to European settlement. The maps of habitat quality ranged in HSI score from 0-100, representing unsuitable to optimal habitat conditions, respectively. These maps were converted into cost surfaces by assigning a 93 relative cost value based on the HSI score. I assigned a cost of 20 to HSI scores of 0 to ensure that non-habitat was heavily weighted against in the cost surface. The remaining HSI scores were grouped into 4 classes at intervals of 25 (1-25, 26-50, 51-75, 76-100), and each class was assigned a cost based on its rank, with the most suitable class (e.g., 76-100) receiving a cost of 1. In this way, poorer quality habitat had a higher cost of traversing, with non-habitat having a very high cost. Connectivity was assessed between 3 points in the UP (Figure 25), with source locations near the border of Canada on the eastern side and the border of Wisconsin near Lake Superior on the western side, and a centralized destination located in the middle of the UP. This destination represented an area where viable home ranges had been simulated under current conditions. I used 2 functions within the GRID extension of ARC/INFO 8.1 to create the least-cost corridors for each map. The costdistance function calculated the least- accumulative-cost distance for each pixel on the cost surface grid to pixels from each of the 3 points of interest. The costdistance grids were paired using the corridor function to calculate the sum of accumulative costs between each of the 2 source locations and the destination; pixels with the lowest accumulative costs identiï¬ed the least-cost corridors between points. The least-cost corridors were compared by year and point pair to quantify differences in their locations and magnitudes. I calculated the average habitat quality and total amount of non-suitable habitat crossed by each path. RESULTS The least-cost corridors illustrated some large differences in lynx habitat connectivity between current and presettlement forest conditions, and the direction of difference depended on the point pair. The change from presettlement conditions to 94 Figure 24. Location of sources (black dots) and destination (white star) for calculating the least-accumulative-cost corridor for lynx in the Upper Peninsula of Michigan. 95 current conditions was relatively minor in the west and more pronounced in the east (Figure 25). These differences reflect the changes in land use across each side of the UP. The lowest accumulated cost for the eastern corridor was 36% and 12% higher for the current map than that for the high and low presettlement maps, respectively (Figure 25). Only the current corridor required the crossing of non-suitable habitat (1.08 km in length), which occurred directly next to the source point in the east. The mean HSI score for the current eastern corridor (63.5) was between those of the presettlement (high = 76.0, low = 62.5). In mapping the corridors under a similar scale of accumulative cost, the two presettlement maps provided the same route, though the high presettlement had less relative cost (Figure 26). The current corridor provided a different route from that of the presettlement maps and the relative cost was higher (Figure 26). The lowest accumulated cost for the western corridor was 4.6% and 2.8% lower for the current map than for the high and low presettlement maps, respectively (Figure 25). None of the corridors required crossing non—suitable habitat. Similar to the eastern corridors, the mean HSI score for the current western corridor (60.0) was bounded by that of the high (61.3) and low (58.2) presettlement corridors. The current corridor provided a different route than that of each presettlement map (Figure 27), but the difference was minor. The small differences in accumulative cost between all 3 maps are reflected in the similarity between the size of each corridor (Figure 27). DISCUSSION The least-cost corridor currently available to lynx in the western UP was similar to those provided under presettlement forest conditions, while in the eastern UP, the current corridor changed signiï¬cantly in a negative way. Contiguous forest patches in 96 ‘ IEastem ‘ 1:] Western Cost (thousands) 1———Current—-I I— Presettlement—4 Figure 25. Least-accumulative—cost values for the eastern and western corridors created for each cost surface in the Upper Peninsula of Michigan. 97 Relative Cost high Figure 26. Relative travel costs for eastern corridors simulated in each cost-surface map of the Upper Peninsula of Michigan, for the high (A) and low (B) presettlement conditions and present conditions (C). 98 , Relative Cost V' high Figure 27. Relative travel costs for western corridors simulated in each cost-surface map of the Upper Peninsula of Michigan, for the high (A) and low (B) presettlement conditions and present conditions (C). 99 the west provided lynx with an adequate array of travel habitat under current conditions, as they did under a natural disturbance regime prior to European settlement. Habitat quality was lower in the west under presettlement conditions, but this did not signiï¬cantly alter the corridor. Habitat quality in the east was higher under presettlement conditions; large patches of unsuitable habitat caused the corridor under current conditions to be located along a different route and accrue a higher cost. The existence of unsuitable habitat is a reflection of the widespread agricultural land and human developments in the eastern UP that were not present during the presettlement era. The only entry point for lynx from Ontario is located directly within an area of unsuitable habitat for lynx near the city of Sault St. Marie. Whether a lynx would survive a migration through this area is unknown, but it represents the only major obstacle to dispersal across the eastern UP. The western UP offers a number of travel corridors, and its proximity to lynx populations in Minnesota suggests that successful dispersals would be more probable from that region. Relating habitat quality to cost is a difï¬cult task, since the relative cost of traveling through an area for a species cannot be readily quantiï¬ed. Little information exists in the literature for determining how lynx might view a landscape in terms of travel cost. Intuitively, areas with low habitat quality would be traversed at greater speeds than those of higher quality; how this would affect the probability of a successful dispersal is unknown. If cost is deï¬ned by the time spent or distance traveled while dispersing, good quality habitat could theoretically have a higher “cost†than that of poor quality habitat, if the good quality habitat is small in area and does not provide enough resources for home range establishment. Individuals that spend time hunting in small patches of good quality 100 habitat could have an increased risk of mortality compared to those that travel quickly to large contiguous patches. The time and distance over which a lynx could travel without foraging might indicate the importance of habitat quality during a dispersal event. More information on the nature of lynx dispersals in the Great Lakes region is necessary to better understand the probability of a population being established in the UP. 101 MANAGEMENT IMPLICATIONS Our results indicated that lynx habitat potential in the UP has changed from that which existed prior to European settlement. Human alterations on the landscape have both increased and decreased habitat quality. The development of human-made structures and agricultural practices in the eastern UP has signiï¬cantly decreased habitat potential in that region, while forestry practices throughout the UP have resulted in a wider distribution of potential habitat. Conversion of forests to other land-uses is usually not reversible, and therefore, the low habitat quality in many regions of the UP, especially in the east, cannot be changed. While it appears that timber harvesting has increased the quality of lynx habitat, this increase is realistically small. Snowshoe hares require widespread areas of early successional vegetation, the size of which can only be created by relatively large disturbances such as severe ï¬re, or clearcutting. The negative social connotation that follows clearcut harvesting can make the practice undesirable and/or impossible across large areas, especially where public opinion seeks to prevent it. Without a consistent cycle of forest disturbance that could allow hare populations to thrive, the chances of a lynx population being established in the UP are low. A further investigation of the areas I qualiï¬ed as suitable habitat is necessary to validate the model outputs, but also to determine whether other factors (i.e., interspeciï¬c competition, snowfall) that were not modeled might have an effect on habitat quality. Given that high snowfall has been hypothesized to separate lynx from other competitors, the current trend of a warming climate will presumably increase the northern expansion of these species may further decrease habitat potential for lynx in southern boreal forests, including those of the UP. 102 LITERATURE CITED Agee, J. K. 2000. Disturbance ecology of North American boreal forests and associated northern mixed/subalpine forests. Pages 39-82 in L. F. Ruggiero, K. B. Aubry, S. W. Buskirk, G. M. Koehler, C. J. Krebs, K. S. McKelvey, and J. R. Squires, editors. Ecology and Conservation of Lynx in the United States. University Press of Colorado, Boulder, USA. Albert, D. A. 1995. Regional landscape ecosystems of Michigan, Minnesota, and Wisconsin: a working map and classiï¬cation. US. Forest Service General Technical Report NC-l78, North Central Forest Experiment Station, St. Paul, Minnesota, USA. Assel, R. A. 2003. Great Lakes Ice Cover, First Ice, Last Ice, and Ice Duration: Winters 1973-2002. NOAA TM GLERL-125. Great Lakes Environmental Research Laboratory, Ann Arbor, MI. Aubry, K. B., G. M. Koehler and J. R. Squires. 2000. Ecology of Canada lynx in southern boreal forests. Pages 373-396 in L. F. Ruggiero, K. B. Aubry, S. W. Buskirk, G. M. Koehler, C. J. Krebs, K. S. McKelvey, and J. R. Squires, editors. Ecology and Conservation of Lynx in the United States. University Press of Colorado, Boulder, USA. Behrend, D. F. 1962. An analysis of snowshoe hare habitat on marginal range. Proceedings of the Northeast Section of The Wildlife Society. Monticello, New York, USA. Beyer, D. E., Jr., B. Roell, J. Hammill, and R. Earle. 2001. Records of Canada lynx in Michigan’s Upper Peninsula. Canadian Field Naturalist 1152234-240. Boutin, S., C. J. Krebs, A. R. E. Sinclair, and J. N. M. Smith. 1986. Proximate causes of losses in a snowshoe hare population. Canadian Journal of Zoology 64:606-610. Brand, C. J ., L. B. Keith, and C. A. Fischer. 1976. Lynx responses to changing snowshoe hare densities in Alberta. Journal of Wildlife Management 40:416-428. Brooke, R. H. 1975. Preliminary guidelines for managing snowshoe hare habitat in the Adirondacks. Transactions of the Northeast Section, The Wildlife Society, Fish and Wildlife Conference. 32:46-66. Buehler, D. A. and L. B. Keith. 1982. Snowshoe hare distribution and habitat use in Wisconsin. Canadian Field Naturalist 96: 19-29. Burt, W. H. 1946. The mammals of Michigan. The University of Michigan Press, Ann Arbor, USA. 288pp. 103 Buskirk, S. W., L. F. Ruggiero, K. B. Aubry, D. E. Pearson, J. R. Squires and K. S. McKelvey. 2000. Comparative ecology of lynx in North America. Pages 397- 418 in L. F. Ruggiero, K. B. Aubry, S. W. Buskirk, G. M. Koehler, C. J. Krebs, K. S. McKelvey, and J. R. Squires, editors. Ecology and Conservation of Lynx in the United States. University Press of Colorado, Boulder, USA. Cary, J. R. and L. B. Keith. 1979. Reproductive change in the 10-year cycle of snowshoe hares. Canadian Journal of Zoology 57:375-390. Cleland, D. T., P. E. Avers, W. H. McNab, M. E. Jensen, R. G. Bailey, King, T., and Russell, W. E. 1997. National Hierarchical Framework of Ecological Units. Pages 181-200 in Ecosystem Management Applications for Sustainable Forest and Wildlife Resources. Yale University, New Haven, Connecticut, USA. Cleland, D. T., T. R. Crow, S. C. Saunders, D. I. Dickmann, A. L. Maclean, J. K. Jordan, R. L. Watson, A. M. Sloan, and K. D. Brosofske. 2004. Characterizing historical and modern ï¬re regimes in Michigan (USA): A landscape ecosystem approach. Landscape Ecology 19:311-325. Cleland, D.T., T. R. Crow, S.C. Saunders, A. L. McLean, and D. I. Dickman. 2005a. Characterizing Historic and Contemporary Fires Regimes in the Lake States. Final Report, Joint Fire Science Program, US. Forest Service. Cleland, D. T., J. A. Freeouf, Keys, Jr., J. E., G. J. Nowacki, C. A. Carpenter, and W. H. McNab. 2005b. Subregions of the conterminous United States. US. Forest Service, Washington, DC, USA. Digital map. Coffman, M. S., E. Alyanak, J. Kotar, and J. E. Ferris. 1984. Field guide, habitat classiï¬cation system for the Upper Peninsula of Michigan and Northeastern Wisconsin. Third Edition. CROF S, School of Forestry and Wood Products, Michigan Technological University, Houghton, USA. Collins, W. B., and E. F. Becker. 2001. Estimation of horizontal cover. Journal of Range Management 54:67-70. Comer, P. J., D. A. Albert, H. A. Wells, B. L. Hart, J. B. Raab, D. L. Price, D. M. Kashian, R. A. Corner, and D. W. Schuen. 1995. Michigan presettlement vegetation, as interpreted from the General Land Ofï¬ce Surveys 1816-1856. Michigan Natural Features Inventory, Lansing, USA. Digital map and database. Crookston, N. L.; and A. R. Stage. 1999. Percent canopy cover and stand structure statistics from the Forest Vegetation Simulator. US. Forest Service General Technical Report RMRS-24, Rocky Mountain Research Station, Ogden, Utah, USA. Daubenmire, R. 1966. Vegetation: Identiï¬cation of typal communities. Science 151:291-298. 104 Dolbeer, R. A. and W. R. Clark. 1975. Population ecology of snowshoe hares in the central Rocky Mountains. Journal of Wildlife Management 39:535-549. Donovan, M. 2005. Michigan GAP analysis project. Pages 47-49 in Maxwell et al., editors. Gap Analysis Bulletin No. 13 USGS/BRD/Gap Analysis Program. Moscow, Idaho, USA. Dickman, D. I. and L. A. Leefers. 2003. The forests of Michigan. University of Michigan Press, Ann Arbor, USA. Dixon, G. E. comp. 2002. Essential F VS: A user’s guide to the Forest Vegetation Simulator. Internal Report, US. Forest Service, Forest Management Service Center, Fort Collins, Colorado, USA. Eichenlaub, V., J. Harman, F. Numberger, and H. Stolle. 1990. The Climatic Atlas of Michigan. University of Notre Dame Press, Indiana, USA. F azakas, Z., and M. Nilsson. 1996. Volume and forest cover estimation over southern Sweden using AVHRR data calibrated with TM data. International Journal of Remote Sensing 17:1701—1709. Felix, A. B., H. Campa, 111, K. F. Millenbah, S. R. Winterstein,, and W. E. Moritz. 2004. Development of landscape-scale habitat-potential models for forest wildlife planning and management. Wildlife Society Bulletin 32:795-806. Ferron, J ., and J. P. Ouellet. 1992. Daily partitioning of summer habitat and use of space by the snowshoe hare in southern boreal forest. Canadian Journal of Zoology 70:2178-2183. Fox, J. F. 1978. Forest ï¬res and the snowshoe hare-Canada lynx cycle. Oecologia 31:349-374. Franco-Lopez, H., A. R. Ek, and M. E. Bauer. 2001. Estimation and mapping of forest stand density, volume, and cover type using the k-nearest neighbors method. Remote Sensing of Environment. 77:251-274. F relich, L. E. 1995. Old forest in the Lake States today and before European settlement. Natural Areas Journal 15: 157-167. F relich, L. E. 2002. Forest dynamics and disturbance regimes: studies from temperate evergreen-deciduous forests. Cambridge University Press, United Kingdom. F relich, L. E., and C. G. Lorimer. 1991. Natural disturbance regimes in hemlock- hardwood forests of the upper Great Lakes region. Ecological Monographs 61 :145-164. 105 Frescino, T. S., T. C. Edwards, Jr., and G. G. Moisen. 2001. Modeling spatially explicit forest structural variables using generalized additive models. Journal of Vegetation Science 12: 15-26. Friedman, S. K. and P. B. Reich. 2005. Regional legacies of logging: departure from presettlement forest conditions in northern Minnesota. Ecological Applications 15:726-744. Haapanen, R., A. R. Ek, M. E. Bauer, and A. O. Finley. 2004. Delineation of forest/nonforest land use classes using nearest neighbor methods. Remote Sensing of Environment 89:265-271. Harger, E. M. 1965. The status of the Canada lynx in Michigan. The Jack-Pine Warbler 43:150-153. Hill, M. O. 1979. TWINSPAN - a FORTRAN program for arranging multivariate data in an ordered two-way table by classiï¬cation of the individuals and attributes. Cornell University, Ithaca, New York, USA. Hill, M. O. and Smilauer, P. 2005. TWINSPAN for Windows version 2.3. Centre for Ecology and Hydrology, Huntingdon, England & University of South Bohemia, Czech Republic. Hodges, K. E. 2000. Ecology of snowshoe hares in southern boreal and montane forests. Pages 163-206 in L. F. Ruggiero, K. B. Aubry, S. W. Buskirk, G. M. Koehler, C. J. Krebs, K. S. McKelvey, and J. R. Squires, editors. Ecology and Conservation of Lynx in the United States. University Press of Colorado, Boulder, USA. Hoving, C. L., D. J. Harrison, W. B. Krohn, R. A. Joseph, and M. O. O’Brien. 2005. Broad-scale predictors of Canada lynx occurrence in eastern North America. Journal of Wildlife Management 69:739-751. Keith, L. B. 1963. Wildlife’s ten-year cycle. University of Wisconsin Press, Madison, USA. Keith, L. B., S. E. M. Bloomer, and T. Willebrand. 1993. Dynamics of a snowshoe hare population in fragmented habitat. Canadian Journal of Zoology 71 :1385-1392. Koehler, G. M. 1990. Population and habitat characteristics of lynx and snowshoe hares in north central Washington. Canadian Journal of Zoology 68:845-851. Koehler, G. M., and J. D. Brittell. 1990. Managing spruce-ï¬r habitat for lynx and snowshoe hares. Journal of Forestry 88: 10-14. 106 Koehler, G. M. and K. B. Aubry. 1994. Lynx. Pages 74-98 in L. F. Ruggiero, K. B. Aubry, S. W. Buskirk, L. J. Lyon, and W. J. Zielinski, tech. eds. The scientiï¬c basis for conserving forest carnivores: American marten, ï¬sher, lynx, and wolverine. US. Forest Service General Technical Report RM-254, Rocky Mountain Research Station, Fort Collins, Colorado, USA. Kotar, J ., J. A. Kovach, and C. T. Locey. 1988. Field guide to forest habitat types of northern Wisconsin. Department of Forestry, University of Wisconsin, Madison, USA. Kotar J ., and T. L. Burger. 2000. Field guide to forest habitat type classiï¬cation for north central Minnesota. Terra Silva Consultants, Madison, Wisconsin, USA. Leatherberry, E. C., and J. S. Spencer, Jr. 1996. Michigan forest statistics, 1993. US. Forest Service, Resource Bulletin NC-170, North Central Forest Experiment Station, St. Paul, Minnesota, USA. Lefsky, M. A., W. B. Cohen, A. T. Hudak, S. A. Acker, and J. L. Ohmann. 1999. Integration of lidar, Landsat ETM+ and forest inventory data for regional forest mapping. International Archives of Photogrammetry and Remote Sensing 32:119-26. Litvaitis, J. A., J. A. Sherbume, and J. A. Bissonette. 1985. Influence of understory characteristics on snowshoe hare habitat use and density. Journal of Wildlife Management 49:866-873. Manies K. L., D. J. Mladenoff, and E. V. Nordheim. 2001. Assessing large-scale surveyor variability in the historic forest data of the original US. Public Land Survey. Canadian Journal of Forest Research 31:1719—1730. McCord, C. M., and J. E. Cardoza. 1982. Bobcat and lynx. Pages 728-766 in J. A. Chapman and G. A. Feldhamer, eds. Wild Mammals of North America. Johns Hopkins University Press, Baltimore, Maryland, USA. McGaughey, R. J. 1997. Visualizing forest stand dynamics using the stand visualization system. Proceedings of the ACSM/ASPRS Annual Convention and Exposition, American Society for Photogrammetry and Remote Sensing 4:248-257. McKelvey, K. S., K. B. Aubry, and Y. K. Ortega. 2000a. History and distribution of lynx in the contiguous United States. Pages 207-264 in L. F. Ruggiero, K. B. Aubry, S. W. Buskirk, G. M. Koehler, C. J. Krebs, K. S. McKelvey, and J. R. Squires, editors. Ecology and Conservation of Lynx in the United States. University Press of Colorado, Boulder, USA. 107 McKelvey, K. S., S. W. Buskirk, and C. J. Krebs. 2000b. Theoretical insights into the population viability of lynx. Pages 21-38 in L. F. Ruggiero, K. B. Aubry, S. W. Buskirk, G. M. Koehler, C. J. Krebs, K. S. McKelvey, and J. R. Squires, editors. Ecology and Conservation of Lynx in the United States. University Press of Colorado, Boulder, USA. McRoberts, R. E., M. D. Nelson, and D. G. Wendt. 2002. Stratiï¬ed estimation of forest area using satellite imagery, inventory data, and the k-Nearest Neighbors technique. Remote Sensing of Environment 82:457-468. Mech, L. D. 1977. Record movement of a Canadian lynx. Journal of Mammalogy 58:676-677. Michigan Department of Conservation. 1938. Ninth biennial report. Michigan Department of Natural Resources - Forest, Mineral and Fire Management Division. 2001. IF MAP/GAP Upper Peninsula Land Cover. Michigan Department of Natural Resources, Lansing, USA. Digital map. Millspaugh, J. J ., G. C. Brundige, R. A. Gitzen, and K. J. Raedeke. 2000. Elk and hunter space-use sharing in South Dakota. Journal of Wildlife Management 64:994- 1003. Monthey, R. W. 1986. Responses of snowshoe hares, Lepus americanus, to timber harvesting in northern Maine. Canadian Field Naturalist 100:568-570. Mowat, G., K. G. Poole, and M. O’Donoghue. 2000. Ecology of lynx in Northern Canada and Alaska. Pages 265-306 in L. F. Ruggiero, K. B. Aubry, S. W. Buskirk, G. M. Koehler, C. J. Krebs, K. S. McKelvey, and J. R. Squires, editors. Ecology and Conservation of Lynx in the United States. University Press of Colorado, Boulder, USA. Nellis, C. H. and L. B. Keith. 1968. Hunting activities and success of lynxes in Alberta. Journal of Wildlife Management 32:718-722. Nellis, C. H., S. P. Wetrnore, and L. B. Keith. 1972. Lynx-prey interactions in central Alberta. Journal of Wildlife Management 36:320-329. Nylen-Nemetchek, M. 1999. Lynx (F elis lynx) in Riding Mountain National Park: An assessment of habitat availability and population viability. MS. Thesis, University of Manitoba, Winnipeg, Canada. O’Donoghue, M., S. Boutin, C. J. Krebs, D. L. Murray, and E. J. Hofer. 1998. Behavioral responses of coyotes and lynx to the snowshoe hare cycle. Oikos 82: 169-183. Omemik, J. M. 1987. Ecoregions of the conterminous United States. Annals of the Association of American Geographers 77:118-125. 108 Orr, C. D., and D. G. Dodds. 1982. Snowshoe hare habitat preferences in Nova Scotia spruce-ï¬r forests. Wildlife Society Bulletin 10:147-150. Parker, G. R. 1986. The importance of cover on use of conifer plantations by snowshoe hares in northern New Brunswick. The Forestry Chronicle 62: 159-163. Parker, G. R., J. W. Maxwell, L. D. Morton, and G. E. J. Smith. 1983. The ecology of the lynx (Lynx canadensis) on Cape Breton Island. Canadian Journal of Zoology 61 :770-786. Plume, D. A. 2005. HomeGrower — Version 3. Forest Capital Partners, LLC, Portland, Oregon, USA. Poole, K. G. 1997. Dispersal patterns of lynx in the Northwest Territories. Journal of Wildlife Management 61:497-505. Rogowitz, G. L. 1988. Forage quality and use of reforested habitats by snowshoe hares. Canadian Journal of Zoology 66:2080-2083. Roloff, G. J. 2001. Habitat suitability model for the North American lynx. Unpublished Report. Timberland Resources, Boise Cascade Corporation, Boise, Idaho, USA Roloff, G. J. 2003. Habitat suitability model for the North American lynx: Lake States version. Unpublished Report. Timberland Resources, Boise Cascade Corporation. Boise, Idaho, USA. Roloff, G. J ., and J. B. Haufler. 1997. Establishing population viability planning objectives based on habitat potential. Wildlife Society Bulletin 25:895-904. Schulte, L. A., and D. J. Mladenoff. 2001. The original US. Public Land Survey Records: their use and limitations in reconstructing pre-European settlement vegetation. Journal of Forestry 99:5-10. Seidel, K. S. 1992. Statistical properties and applications of a new measure of joint space use for wildlife. MS. Thesis, University of Washington, Seattle, USA. Sievert, P. R., and L. B. Keith. 1985. Survival of snowshoe hares at a geographic range boundary. Journal of Wildlife Management 49:854-866. Slough, B. G. 1999. Characteristics of Canada lynx, Lynx canadensis, maternal dens and denning habitat. Canadian F ield-Naturalist 113:605-608. Slough, B. G. and G. Mowat. 1996. Population dynamics of lynx in a refuge and interactions between harvested and un-harvested populations. Journal of Wildlife Management 60:946-961. Smith, W. B. 2002. Forest inventory and analysis: a national inventory and monitoring program. Environmental Pollution 11628233-8242. 109 Subdedi, N. 2005. Comparative analysis of forest classiï¬cation in forest management information databases in Michigan. MS. Thesis, Michigan State University, East Lansing, USA. Tokola, T., J. Pitkanen, S. Partinen, E. Muinonen. 1996. Point accuracy of a non- parametric method in estimation of forest characteristics with different satellite materials. International Journal of Remote Sensing 17: 2333-2351. Tomppo, E. 1991. Satellite image-based national forest inventory of Finland. International Archives of Photogrammetry and Remote Sensing 28:419-424. US. Census Bureau. 2000. Census 2000 - Michigan total population by block group. <http://factï¬nder.census.gov>. Accessed 3 June 2006. US. Department of Agriculture, Natural Resources Conservation Service. 1994. State Soil Geographic (STATSGO) database for Michigan. USDA, NRCS, Fort Worth, Texas, USA. Digital map. US. Department of Agriculture, Natural Resources Conservation Service. 2006. Soil Survey Geographic (SSURGO) database for Michigan. <http://soildatamart.nrcs.usda.gov> Accessed 8 January 2006. US. Fish and Wildlife Service. 1980. Habitat Evaluation Procedures (HEP). US. Fish and Wildlife Service ESM-102. US. Fish and Wildlife Service. 1981. Standards for the development of habitat suitability index models. US. Fish and Wildlife Service ESM-103. US. Fish and Wildlife Service. 2000. Endangered and threatened wildlife and plants: Determination of threatened status for the contiguous US. distinct population segment of the Canada lynx and related, ï¬nal rule. Federal Register 65:16052- 16086. U.S. Forest Service. 2006. Forest Inventory and Analysis Program. <http://www.ï¬a.fs.fed.us>. Accessed 8 January 2006. Van Horne, B. 1983. Density as a misleading indicator of habitat quality. Journal of Wildlife Management 47:893—90 1. Ward, R., and C. J. Krebs. 1985. Behavioral responses of lynx to declining snowshoe hare abundance. Canadian Journal of Zoology 63:2817-2824. Wimberly, M. C., T. A. Spies, C. J. Long, and C. Whitlock. 2000. Simulating historical variability in the amount of old forests in the Oregon Coast Range. Conservation Biology 14:167-180. 110 Wirsing, A. J ., T. D. Steury, and D. L. Murray. 2002. Demographic analysis of a southern snowshoe hare population in a fragmented habitat. Canadian Journal of Zoology 80: 169-177 Wolfe, M. L., N. V. DeByle, C. S. Winchell, and T. R. McCabe. 1982. Snowshoe hare cover relationships in northern Utah. Journal of Wildlife Management 46:662- 670. Wolff, J. O. 1980. The role of habitat patchiness in the population dynamics of snowshoe hares. Ecology Monograph 50:111-130. Wolff, J. O. 1981. Refugia, dispersal, predation, and geographic variation in snowshoe hare cycles. Pages 441-449 in K. Myers and C. D. Maclnnes, editors. Proceedings of the World Lagomorph Conference. University of Guelph, Ontario, Canada. Wood, N. A. and L. R. Dice. 1924. Records of the distribution of Michigan mammals. Papers of the Michigan Academy of Science, Arts, and Letters 3:425-445. Zhang, Q.F, K. S. Pregitzer, and D. D. Reed. 1999. Catastrophic disturbance in the presettlement forests of the Upper Peninsula of Michigan. Canadian Journal of Forest Research 29: 106-1 14. Zhang, Q.F, K. S. Pregitzer, and D. D. Reed. 2000. Historical changes in the forests of the Luce District of the Upper Peninsula of Michigan. American Midland Naturalist 143294-110. 111 APPENDIX A: HABITAT SUITABILITY MODEL FOR THE NORTH AMERICAN LYNX: LAKES STATES VERSION Gary J. Roloff Timberland Resources, Boise Cascade Corporation 1564 Tomlinson Road Mason, MI 48854 GaryRoloff@bc.com This habitat model was originally written for the intermountain west region of the United States. In an eflort to calibrate the model for use in the Lake States, modiï¬cations to the mathematical equations were implemented in a manner consistent with how the model was applied in a veriï¬cation study in Manitoba, Canada (Nylen-Nemetchek 1999). These modiï¬cations were based on the fact that most portions of the Lake States tend to have less accumulated snow, less topographic influences, and an overall smaller forest stature than the intermountain west forest communities. Although outputs from the calibrated lynx model in Manitoba corresponded to the occurrence of lynx as expected (Nylen- Nemetchek I 999), this Lake States version should be viewed as a hypothesis on lynx habitat relationships in the Region. Throughout this document, the citations in support of model variables remained unchanged from the original intermountain west documentation. Introduction In 1995, concerns over lynx (Lynx canadensis) population viability in the conterminous United States prompted several resource groups to petition the species for listing under the Endangered Species Act. Since that time, lynx conservation has received considerable attention (e.g., Paquet and Hackman 1995, Washington Department of Natural Resources 1996, Ruggiero et al. 2000). Lynx were listed as federally threatened in 2000 throughout the conterminous United States range (US. Fish and Wildlife Service 2000). Lynx are also listed on Appendix II of the Convention on International Trade of Endangered Species (CITES) and are currently identiï¬ed as a sensitive species by several US. Forest Service regions (Macfarlane 1994). Concerns surrounding the effects of land management activities on lynx populations in the conterminous United States necessitated development of a model that quantitatively assessed these effects on habitat suitability. Although few data exist on lynx in the United States from which to build habitat models, considerable research has been conducted in Canada and Alaska (see review in Ruggiero et al. 2000). Lynx data from Canadian and Alaskan studies must be applied cautiously to the lower United States because of the unique features associated with the ' 112 southern range of lynx distribution. These unique features of southern lynx range, with an emphasis on western studies, include: o The inherently peninsular and disjunct distribution of suitable habitats (Koehler and Aubry 1994). o The lack of dramatic fluctuations in both lynx and snowshoe hare (Lepus americanus) populations (Dolbeer and Clark 1975, Sievert and Keith 1985, Koehler 1990, Koehler and Aubry 1994). 0 More consistently low hare densities, comparable to hare population lows in Canada and Alaska (Koehler and Aubry 1994). 0 Consistently low lynx densities making the effects of fur-harvests on populations in some areas additive rather than compensatory (Koehler and Aubry 1994). 0 Higher human densities coupled with low lynx densities potentially causing both direct (e.g., fur harvest) and indirect (e. g., land development causing displacement) anthropogenic effects on population persistence to be of greater magnitude. 0 The potential importance of immigration from the north for short-term population persistence (Koehler and Aubry 1994). o The range overlaps of lynx, bobcat (F elis ruï¬is), and coyote (Canis latrans) (Koehler and Aubry 1994), and the propagation of bobcat and coyote range extensions that typically accompany anthropogenic development. The habitat suitability index (HSI) modeling concept (U .8. Fish and Wildlife Service 1981) was used to develop a lynx habitat model for the southern Rocky Mountains that was ï¬eld veriï¬ed (Nylen-Nemetchek 1999) and as such, formed the basis for a Lake States model. In this model, habitat potential for a lynx home range was divided into foraging, denning, and travel requisites. The lynx model uses a limiting factor approach (a concept founded in ecological theory) in that the most limiting resource(s) is assumed to have the greatest impact(s) on the population. A premise of HSI models is that limiting factors can be portrayed using mathematical relationships between vegetation structures, spatial metrics, and indices of habitat quality (US. Fish and Wildlife Service 1981). Theoretically, these limiting factors can be expressed as an index to animal ï¬tness. Overview Critical considerations prior to running the lynx model are the resolution, accuracy, and precision of the land classiï¬cation system and associated vegetation attribute information for the 113 planning landscape. The land classiï¬cation must be capable of characterizing vegetation structures and spatial arrangements at a resolution compatible with lynx and snowshoe hare habitat use. The ideal stratiï¬cation is a stand-based (minimum mapping unit around 2 ha) ecological classiï¬cation system that integrates existing vegetation conditions and site potentials (e.g., geology, soils). An ecological classiï¬cation system is recommended to reduce the variability in quantifying understory vegetation attributes since these attributes are primary components of the lynx model. Deviations from these baseline recommendations for land stratiï¬cation will reduce the robustness and utility of the model output. The lynx model is divided into three components: 1) foraging, 2) denning, and 3) interspersion. The lynx model was speciï¬cally developed for the mountainous habitats of Washington, Idaho, and Montana, corresponding to the southern extension of lynx range in the Rocky Mountains, however, the model framework may be applied to other regions, in this case the Lake States region. In applying the model to the Lake States, each input variable was calibrated to the biogeoclimatic conditions characteristic of the region. For example, model variables that index winter browse availability for snowshoe hares will differ across regions depending on average snow depths. These types of differences must be accounted for when applying this framework to other regions. Lynx habitat in the southern Rocky Mountains consists primarily of 2 structurally different forest types occurring at opposite ends of the forest seral gradient: 1) early successional forest structures that contain high numbers of prey (especially snowshoe bare), and 2) late- successional forest structures for denning (Koehler and Aubry 1994). Second-growth forests with dense understories also may support abundant hare populations (John Weaver, Northern Rockies Conservation Cooperative, Missoula, MT, pers. comm). Intermediate seral stages with sparse understories serve as travel cover, functioning to provide connectivity between foraging and denning patches (Koehler and Aubry 1994). Recent data from Montana and Maine suggest that late-successional forests may not be as important as originally thought for denning, but may serve as snowshoe hare refugia thereby contributing to foraging habitat. Literature reviews and consultation with experts on lynx and snowshoe hare ecology were used to develop this lynx model. The model is not stand-based, but rather, it is designed to evaluate habitat quality in an area that corresponds to the allometric home range of lynx (250 ha; Roloff and Haufler 1997). Within a 250 ha area, habitat quality is expressed on a scale of 0-100, denoting "poor to good" habitat, respectively. Subsequently, habitat units from each allometric home range are aggregated into viable, marginal, or non-viable areas, the size of which depends on habitat quality (Roloff and Haufler 1997). 114 Lynx Foraging Component Forage availability during the winter months appears to be the most important criterion in the determination of lynx home range size and degree of home range overlap (McCord and Cardoza 1982, Ward and Krebs 1985). Lynx populations covary with snowshoe hare numbers (Brand et al. 1976, Brand and Keith 1979, Parker 1981, Bailey et al. 1986), and lynx tend to choose habitats where hares are most numerous (Murray et al. 1994). Although prey switching has been documented in the southern Rocky Mountains, the underlying determinant of lynx ï¬tness appears to be related to winter snowshoe hare abundance. Thus, the foraging component of the lynx model is based on winter snowshoe hare habitat quality. Snowshoe hare habitat is assessed using an HSI model, and the results of the hare model are incorporated into a lynx foraging assessment. HABITAT SUITABILITY MODEL FOR SNOWSHOE HARE Overview Important components of hare habitat have been reported for different vegetation types and include dense woody vegetation (Adams 1959, Monthey 1986, Koehler 1990, Keith et al. 1993), stem diameter of browse (Keith 1974), continuity of coniferous cover (Brocke 1975), habitat interspersion (Keith et al. 1993), the distance to lowland forest cover (Conroy et al. 1979), and patch size (Thomas et al. 1997). The snowshoe hare model is divided into two primary components: 1) foraging, and 2) security cover. These components are mathematically combined into an overall index of winter hare habitat quality at the map-polygon and home range levels. Winter foraging and security cover conditions are assumed to be limiting to hares (Hart et a1. 1965, Dolbeer 1972, Keith and Windberg 1978, Pease et al. 1979, Keith et al. 1984, Boutin et al. 1986, Keith et al. 1993). In this model, summer habitats are not considered a limiting factor. To index the quality of snowshoe hare habitat, it is assumed that measures of understory cover and species composition in different height strata can be used (summarized by Ferron and Ouellet 1992). In support of this assumption, Thomas et al. (1997) demonstrated signiï¬cant relationships between hare population indices and the horizontal and vertical cover of understory vegetation. Since few snowshoe hare studies have been conducted in the Paciï¬c Northwest and Lake States, the vegetation-hare relationships depicted in this model were inferred from Thomas et al. (1997). Studies conducted across North America were used to supplement Thomas et al.’s (1997) work. Snowshoe Hare Winter Food Component Winter availability of palatable browse is believed to be a limiting factor of snowshoe hare populations (e.g., Windberg and Keith 1976, Pease et al. 1979, Vaughan and Keith 1981, 115 Sinclair et al. 1988, Sullivan and Sullivan 1988, O’Donoghue and Krebs 1992). In this model, the amount of winter browse for snowshoe hares is assessed using two different measures: 1) the amount of live horizontal (or lateral) cover, and 2) the amount of live vertical cover. Both measures are used to represent the "thickness" of forage for hares. Horizontal and vertical cover correlate with understory stem density (Gysel and Lyon 1980, Litvaitis et al. 1985), although this relationship may be weak (Thomas et al. 1997). In the southern Rocky Mountains, forage for hares is quantiï¬ed in three height strata; 0-1, 1-2, and 2-3 m to account for variations in availability as a result of changing snow depths and the ability of hares to "clip" down vegetation from unreachable heights (Keith et al. 1984, Sullivan and Moses 1986). In the Lake States region, however, snow depths rarely exceed 2 m in height, particularly for areas inland from the Great Lakes. Thus, this Lake States lynx model is based on 2 height strata, 0-1 m and 1-2 m. Horizontal cover is measured along the geometric plane that corresponds to the ground (i.e., the thickness if one stands and tries to look through a vegetation type) whereas vertical cover is measured along a geometric plane perpendicular to the ground (i.e., the thickness if one looks up). The woody browse component in this model includes all live plants since hares have been documented to feed on a variety of species (Table 1). However, it is recognized that some browse is unpalatable or of higher quality to hares and if possible, this model component should include only those species used by hare for foraging. Thomas et al. (1997) found that highest browse use occurred in vegetation types with 30 to 40% horizontal cover of live vegetation. Use of vegetation types for foraging declined as woody cover approached <20% (F erron and Ouellet 1992, Thomas et al. 1997). These ï¬ndings roughly correspond to other studies that found highest hare use during winter in vegetation types with 250% horizontal cover (Carreker 1985, Parker 1986). Thus, optimal foraging habitat for snowshoe hares is provided by vegetation types with 235% horizontal cover of live vegetation (Fig. l). Hare winter foraging habitat quality declines as horizontal cover decreases, and habitat is unsuitable when 510% horizontal cover of live vegetation is provided (Fig. 1). Horizontal cover for foraging habitat is measured for the 0-1 and 1-2 m height strata. Thomas et a1. (1997) also associated vertical cover with the intensity of snowshoe hare browsing. Highest browse levels corresponded to about 80% vertical cover. Browse use approached zero as vertical cover declined to about 20%. In this model, vertical cover of live vegetation is optimum at 280% and provides no foraging habitat at 520% (Fig. 2). Similar to horizontal cover, vertical cover is measured for two height strata. For both horizontal and vertical cover relative to snowshoe hare browsing potential, overall habitat quality is assessed independently for each strata (i.e., an increase in browse in one 116 stratum cannot offset a decrease in another stratum). The rationale behind this logic is that snow levels dictate the heights at which hares can access browse, thus, the different strata cannot compensate for each other (i.e., if the 0-1 m strata is unavailable, the quality does not matter because hares cannot access it). Two HSI scores are calculated from ï¬gure 1: 1) horizontal cover 0-1 m tall (HSI/1am, winnfood, hcov, 0-1), and 2) horizontal cover 1-2 m tall (HSI/tare,Wilt!,f00d,hC0V,1—2)° Similarly, two HSI values are derived from ï¬gure 2: 1) vertical cover 0-1 m tall (HSIhare,wint, food, vcov, 0- 1) and 2) vertical cover 1-2 m tall (HSIhare,wint ,fooa', vcov, 1- 2). The equation for calculating hare foraging HSIs for the Lake States from horizontal and vertical cover is presented in equations 1 and 2, respectively. Equation 1: (HSIhare, wint, food, hcov, 0 - I + HSIhare, wint, food. hcov] - 2) 2 = HSIhare, wint, food, hcov Equation 2: (HSIhare, wint, food, vcov, 0 - 1 + HSIhare, wint, food, vcovI - 2) = HSI , 2 hare, wmt, food, vcov The foraging habitat quality for snowshoe hare is based on the arithmetic mean of HSIhare,wint,food,hcov and HSIhare,wintfood vcov (Equation 3). An arithmetic mean was selected because some foraging habitat can be provided (i.e., the foraging HSI > 0) if only horizontal or vertical foraging cover is present. For example, densely-stocked woody species with single-stem growth forms that do not have spreading crowns [e.g., aspen (Populus tremuloides)] will tend to exhibit suitable horizontal cover during winter months whereas the vertical cover provided by this vegetation community may be marginal. Using the arithmetic relationship in Equation 3, horizontal or vertical foraging cover can equal 0 and the foraging HSI can still be >0. Both horizontal and vertical foraging cover is weighted equally in the winter food component (Equation 3). 117 Equation 3: (HSIhare, wint, food, hcov + HSIhare, wint, food, vcov) = HSI , 2 hare, wmt, food Snowshoe Hare Winter Security Cover Component The presence of adequate winter security cover has been recognized as the primary determinant of hare habitat quality (Buehler and Keith 1982, Wolfe et al. 1982, Sievert and Keith 1985, Rogowitz 1988). In this model, cover is deï¬ned as any structure (live or dead) that provides security for snowshoe hares. Winter security cover for hares is assessed using three measures of structure and composition: 1) understory vegetation composition (Severaid 1942, deVois 1962, Bookhout 1965a, b, Buehler and Keith 1982, Orr and Dodds 1982), 2) horizontal cover in three height strata (Brocke 1975, Wolfe et al. 1982), and 3) vertical cover in three height strata (Wolff 1980). Understory Vegetation Composition Snowshoe hares appear to select habitats based on vegetation structure as opposed to species composition (Litvaitis et al. 1985, Ferron and Ouellet 1992) and will use most forest cover types if adequate understory vegetation exists. However, some researchers have demonstrated that vegetated stands <3 m tall dominated by conifer species provide better habitat as opposed to stands dominated by deciduous species (deVos 1962, Buehler and Keith 1982, Orr and Dodds 1982, Monthey 1986). A subjective evaluation of the dominant understory vegetation type 53 m tall is used to index winter cover composition (HSI/rare,wint,cov,d0m)- The following criteria were developed to calculate HSIhare,wint,cov,dom3 . a Understory Cover Dominance Class Deciduous Mixed Coniferous None HSI=50 HSI=75 HSI=100 HSI=0 a . . . . Based on a subjective evaluation of understory cover S3 m tall. "Decrduous" = >60% understory in deciduous species; "Mixed" = 40% S Deciduous/Coniferous Cover 560%; "Coniferous" = >60% understory in coniferous species. If no understory cover exists, the HSI defaults to 0. 118 Horizontal Security Cover Brooke (1975) suggested that horizontal cover is the single most important stimulus in selecting cover to avoid predation. Parker (1986) found that snowshoe hare population indices were related to horizontal cover in the 1-2 and 2-3 m height strata. Winter horizontal cover (HSIhare,wint,cov,hcov) includes both live and dead vegetation and inanimate objects (e.g., rocks, root wads). Optimum horizontal cover is provided at 290%, and horizontal cover 340% is deemed unsuitable winter habitat (Carreker 1985, Ferron and Ouellet 1992) (Fig. 3). Separate horizontal cover HSIs are calculated for height strata 0-1 and 1-2 m and these are subsequently combined using an arithmetic mean to produce HSIhare,wint,cov, hcov (Equation 4). Equation 4: (HSIhare, wint, cov,hcov,0 - I + HSIhare, wint, cov, hcov] - 2) = HSI , 2 hare, wmt, cov, hcov Vertical Security Cover Vertical vegetation cover is also considered an important component of hare habitat quality (Wolff 1980). Vertical cover is deï¬ned as the percent cover of live or dead material. Again, multiple strata are used to account for variations in cover availability as a result of changing snow depths. Optimal vertical cover occurs at 290%, and vertical cover 340% provides unsuitable habitat (HSIhare,wint,cov,vcov) (Fig. 4). Separate vertical cover indices are calculated for height strata 0‘1 (HSI hare, wint, cov,vcov, 0-1 ) and 1’2 (HSIhare,wint,cov,vcov,1-2) m and are subsequently combined using an arithmetic mean to produce HSIhare,wint,c0v,vcov (Equation 5). An arithmetic mean was selected because if vertical cover is provided in one stratum, the vegetation type provides functional security cover for at least a portion of the winter (i.e., until snow covers it). Both vertical cover strata are weighted equally in the winter vertical cover component (Equation 5). 119 Equation 5: [HSIhara wint, cov, vcov,0 - 1 + HSIhare, wint, cov,vcov1 - 2) _ HSI 2 _ hare, wint, cov, vcov Winter security cover for snowshoe hares (HSIhare,wint,cov) is computed by ï¬rst establishing whether suitable security cover exists (as the arithmetic mean of HSI hare, wint, cov, hcow and HSIhare,wint,cov,vcov) (Equation 6). Subsequently, the arithmetic mean from the cover calculation is geometrically combined with the understory composition component (HSIhare,wint,cov,dom) (Equation 6). This mathematical relationship will cause HSI/1am, wint,cov to score as unsuitable if appropriate cover conditions are not provided. Equation 6: h . t h + HSI h . 0'5 are, wzn, cov are, wm, vcov x HSI . : HSI ' 2 hare, wmt, dom hare, wmt, cov Cachlating the Snowshoe Hare Winter HSI Winter habitat conditions are assumed to limit snowshoe hare populations, and thus, winter HSI values drive the ï¬nal HSI calculation. Winter habitat components (forage and cover) are integrated into one habitat value using a geometric mean. If the winter HSI for one habitat component equals 0, the ï¬nal HSI equals 0 (i.e., suitable forage and cover must be present to provide hare habitat). Equation 7 is used to calculate the snowshoe hare winter HSI (H SI hare, wint)- Equation 7: (HSI )05 = HSI . x HSI . . hare, wmt, food hare, wmt, cov hare, wmt 120 Calculating the Lynx Forage Component The index HSIhare,wint provides a polygon-level assessment of snowshoe hare habitat quality. The next step in the modeling process for lynx is to relate the polygon-level depiction of hare habitat quality to the allometric home range of lynx (250 ha). It is assumed that for each allometric home area to support lynx, some minimum level of foraging habitat (i.e., snowshoe hare habitat) is required. These habitats must themselves be of sufï¬cient quality to support consistent and abundant numbers of snowshoe hares. Applying the methodology of Roloff and Haufler (1997), home range functionality thresholds were established for snowshoe hares based on an evaluation of hare studies. Similar to relationships demonstrated for other wildlife species, the home range of Lagomorphs appears negatively associated with habitat quality (Boutin 1984, Hulbert et al. 1996). It is assumed that hares will exhibit smallest home ranges when habitat conditions are optimum and that hares have largest home ranges or become nomadic in unsuitable habitats (see Roloff and Haufler 1997). Although few home range studies have quantiï¬ed habitat quality and estimates of annual home range sizes for hares are uncommon, existing literature and allometric theory were used to estimate home range functionality thresholds for snowshoe hares (Roloff and Haufler 1997). The smallest documented home range for snowshoe hares is 1.4 ha for females (mid- summer to fall)(Ferron and Ouellet 1992). Ferron and Ouellet’s (1992) estimate is smaller than the allometric home range for snowshoe hares [4.5 ha, assuming an average mass of 1,400 g (Boutin et al. 1986) and using the equation for primary consumers from Harestad and Bunnell 1979], but note that Ferron and Ouellet (1992) did not estimate an annual range. Studies conducted over longer time periods have demonstrated larger home ranges. For example, Dolbeer and Clark (1975) estimated a home range of 8.1 ha using mark-recapture techniques from mid-April to early September in Colorado. Similarly, Sievert and Keith (1985) documented annual home ranges >10 ha in Wisconsin. Neither of these studies occurred in what would be considered optimum habitat conditions (Dolbeer and Clark 1975, Sievert and Keith 1985), thus, for an annual estimate of snowshoe hare home range in optimal habitat, the allometric scale (4.5 ha) seems to be a reasonable minimum area threshold (Fig. 5). The maximum documented home range (excluding nomadic individuals) is 16 ha (Behrend 1962). Behrend’s (1962) study occurred at the southern edge of snowshoe hare range in presumably fragmented and thus sub-optimal habitat (Sievert and Keith 1985). Sievert and Keith (1985), working in fragmented habitats in Wisconsin, also documented home ranges >10 ha in size. Based on assumptions between home range size and expected productivity (see Fig. 5), this 121 model assumes that hares exhibiting home ranges of 16 ha or larger are not contributing to the viability of the population (see Roloff and Haufler 1997). Habitat quality thresholds (Roloff and Haufler 1997) were inferred by comparing home range size to hare productivity under the assumption that larger home ranges correspond to lower quality habitats and thus lower productivity (Fig. 5). The maximum annual productivity of snowshoe hares (18 young/female) has been recorded from the center of their geographic range (i.e., central Alberta) in what many assume to be an area of high quality habitats (Cary and Keith 1979). This model assumes that maximum reproductive output corresponds to optimum habitat conditions, that habitat quality scores scale linearly with reproductive output, and that the maximum documented home range corresponds to habitat quality in which a female only replaces herself annually (i.e., 2 offspring per year assuming a 50-50 sex ratio)(Fig. 5). Using documented productivity rates and estimates of home range area, a viability relationship was established for snowshoe hare (Fig. 5). A 4.5 ha (corresponding to the allometric home range of hares) area under optimal habitat conditions (i.e., HSI score = 100) provides 4.5 habitat units for hares (consistent with the US. Fish and Wildlife Service’s process for calculating habitat units). The accumulation of habitat units to the 4.5 target occurs in a grid- environment and is based on the snowshoe hare HSI score (see Roloff and Haufler 1997). The higher the quality of the habitat, the fewer number of pixels required to achieve the habitat unit target. As habitat quality declines, more pixels are required and home range size increases (Roloff and Haufler 1997). Using the viability relationship developed for snowshoe hares (Fig. 5) and the output from the snowshoe hare HSI model, the forage potential of each lynx home range is scored according to the number of viable, marginal, and non-viable hare ranges it encompasses (Fig. 6). Hare home ranges above the viability threshold (60 HSI, Fig. 5) count double towards the home range tally whereas marginal home ranges (between 25 and 60 HSI, Fig. 5) count one each. Non-viable home ranges do not contribute towards the lynx forage score. It is estimated that lynx require about 600 g of food/day (or a hare every 2 days) to subsist during winter (Brand et al. 1976). Assuming that the winter season starts in November and extends through April (about 180 days), this would imply that 90 hares would support a lynx through winter. Thus, 90 hare home ranges in a lynx home range was considered optimum (Fig. 6). The resulting HSI score from tallying hare home ranges and applying the sum to ï¬gure 6 is the foraging score for the 250 ha lynx home range. 122 Lynx Denning Component Delineation of potential lynx denning habitat is a 3-phase process: I) identify vegetation types that provide vegetative structure and size deemed suitable for denning, II) identify vegetation types that are properly arranged within a home range area, and III) identify vegetation types that provide suitable denning micro-sites. Components of suitable lynx denning habitat include: 1) vegetation cover type, 2) mesic site conditions, 3) canopy closure, 4) the area of the vegetation type, 4) juxtaposition and interspersion, and 5) the amount and arrangement of downed woody debris. These stand- and site-based components are integrated into a single estimate of denning habitat quality (HSIlynx,den) for the home range area. Management for denning habitat should also emphasize minimizing human disturbance. PHASE I: (Vegetation type, site potential, canopy closure, and the area of the vegetation type) Potential denning sites are initially delineated by vegetation cover type and site conditions. Vegetation types classiï¬ed as forested with an average diameter of 20 cm providing 2 11.49 m2/ha (50 ft2/ac) of basal area on mesic sites are assumed to satisfy denning requirements in the Lake States region. These areas must also have >50% canopy closure (where "canopy" is deï¬ned as trees >5 min height) and be a minimum of 2 ha in size. Also, rock crevices, caves, and overhanging banks may be used for denning sites (Hoover and Willis 1987). PHASE II: Juxtaposition and Interspersion Denning sites (the 2 ha patch, rocks, crevices, or banks identiï¬ed in Phase I) must be in close proximity to forage cover (Koehler and Brittell 1990). At least 50% of the perimeter of 2 ha patches identiï¬ed as potential denning sites in “Phase I†must be adjacent to some form of "suitable lynx habitat" (i.e., habitat identiï¬ed as denning, foraging, or travel). Also, 30% of the land within 0.8 km of the potential denning sites must be in suitable summer foraging habitat. Suitable summer foraging habitat is based on the habitat potential score calculated for the 0-1 m vertical cover measurement in the snowshoe hare model. Snowshoe hares forage on a variety of herbaceous vegetation during the spring and summer months (Wolff 1978), and thus, hare forage is not limiting. It is assumed, however, that snowshoe hares are more vulnerable to predation in open areas (Adams 1959, Dolbeer and Clark 1975, Sievert and Keith 1985), and thus, vegetation cover for security is the limiting factor during spring and summer. Snowshoe hare summer security cover can be estimated based on the amount of cover (both live and dead) 0-1 m tall. The HSI for summer cover (HSIhare,sum,cov) is derived 123 from a measure of vertical cover provided by all objects within the 0-1 m height strata. Optimum summer cover for hares exists in stands providing 260% vertical cover, and summer cover habitat quality is 0 when 520% vertical cover exists (Fig. 7). Map polygons with a summer forage HSI value >10 satisfy the lynx forage requirements. A map polygon may provide both suitable forage and denning, in which case the denning area is counted towards the 30% foraging. If these criteria are satisï¬ed, the map polygon is a potential denning site and an assessment of downed woody debris is performed. Map polygons not identiï¬ed as potential denning habitat through the ï¬rst 2-phases are assigned a HSIIynx, den, stand value of 0. These sites are not evaluated for Phase 111 attributes. Also, denning structures can be constructed in association with sites that satisfy the Phase I and 11 criteria. PHASE III: Dead and Downed Woody Debris Dead and downed woody debris include logs, stumps, and upturned root masses (Koehler 1990, Koehler and Brittell 1990). Potential lynx dens generally consist of inter-tangled woody material with interstitial spaces large enough to provide lynx cover. Lynx dens have been described as having a high density (40 pieces per 50 m) of downed woody debris that were vertically structured 0.3-1.2 m above the ground (Brittell et al. 1989, Koehler 1990). These types of structures are often dependent on micro-site characteristics (e. g., areas susceptible to wind throw; drainages) and are often uncommon across entire forest stands, thus, it was deemed impractical to systematically sample this attribute within a stand and base the estimate on a mean value. Additionally, recent data from Montana and Maine have demonstrated that lynx exhibit greater variability in den site selection than previously thought. These studies have documented lynx denning in young forest plantations in slash piles and jackstrawed tree stems. The common theme among all lynx denning data is that some minimal amount of cover is required and that these areas have access to abundant hare populations. The lynx model assumes that habitat patches identiï¬ed in Phase I and Phase II, with some occurrence of downed wood or slash, provide suitable denning sites for lynx. Within lynx home ranges, multiple denning sites are important. Females may move kittens to better foraging areas or to avoid disturbance (Koehler and Brittell 1990). In low quality habitat, the inability of females to move kittens may increase kitten mortality (Koehler and Aubry 1994). Assuming that the majority of denning stands contain suitable micro-sites (veriï¬ed by walk-through inventories), the quantity and spatial extent of denning stands is used to index denning habitat quality. An optimal home range is assumed to contain a mosaic of vegetation 124 types that include foraging and denning habitat (Koehler and Aubry 1994). In the lynx habitat model, the denning score in a home range is based on the average distance in a home range to denning sites. The variable HSI lynx, den is calculated on the premise that multiple, interspersed denning sites in a home range is of better habitat quality than a home range containing few, blocked sites. To assess each home range area, a 100x100 m grid of points is overlaid and the average nearest distance to a suitable denning site from all points is calculated. Optimum denning habitat is provided when the average distance to denning sites is £400 m and denning habitat is deemed unsuitable if average distance is 21,750 m (Fig. 8). Under these parameters, optimum conditions roughly correspond to a denning site every 16 ha. Lynx Interspersion Component The interspersion component is designed to address the "travel corridor" needs of lynx (Washington Department of Wildlife 1993). Lynx travel through and on a variety of vegetation cover types and landscape features including thinned and un-thinned forested stands, regeneration, open meadows (S 100 m in width), ridges just above timberline, roads, and forest trails (Taylor and Shaw 1927, Parker 1981, Brittell et al. 1989, Koehler 1990). This model assumes lynx will traverse most cover types except open or sparsely-stocked stands >100 m in width. The interspersion component of this model uses a 2-step process: 1) identify areas of "non-lynx" cover, and 2) index the amount and spatial distribution of "non-lynx" habitat in the home range. "Non-lynx" habitat is deï¬ned as map polygons (or portions of map polygons) with a winter foraging or denning HSI of 0 that are: a) permanent "openings" >91 m in width (e.g., meadows), b) map polygons with perennial vegetation <2 m tall and >91 m in width, and c) map polygons with <440 trees/ha (178 trees/ac) having a 2-3 m understory providing <50% visual obstruction. It is important to note that some map polygons may be split during this process, (i.e., a portion of the polygon is >91 m in width and the other portion is <91 m in width). These portions need to be segregated during the analysis to reduce assessment error. Suitable travel cover is subsequently delineated as map polygons not identiï¬ed as forage, denning, or "non-lynx" habitat. The interspersion index (HSI/ynx' inter) is based on the average nearest distance within a home range to "non-lynx" polygons. A systematic grid (100 x 100 m) is used to estimate the average distance to “non-lynx†polygons (Fig. 9). The HSIbmx’ inter is based on the premise that a lower average distance to "non-lynx" polygons equates to a more interspersed conï¬guration of 125 habitats, and thus, to a greater probability of lynx encountering travel barriers (Fig. 10). The 100 x 100 m grid corresponds to the maximum hypothetical distance lynx will traverse without sufï¬cient cover. Model simulations conducted on z5,000 ha in potential lynx habitat in northeast Washington demonstrated that the size of the sample grid had negligible impact on the average nearest distance to “non-lynx†habitats (Table 2). Of more importance is the relationship depicted in ï¬gure 10. Lynx will traverse long distances to fulï¬ll their life requisites. For example, Brand et al. (1976) and Nellis and Keith (1968) found that lynx traveled 8.8 km hunting during hare population lows and 4.7 km when hares were plentiful. Parker et al. (1983) calculated daily cruising distances of 6.5-8.8 km in winter and 73-101 km during summer in Nova Scotia. Koehler (1990) documented females foraging 6-7 km from their den sites. The habitat model for lynx penalizes landscapes that restrict these movements. Figure 10 attempts to quantify the effects of barriers to movement on habitat quality (i.e., how often can lynx encounter movement barriers without detracting from habitat quality?) A low average nearest distance to “non-lynx†habitat in a home range (i.e., the chances of encountering a “non-lynx†polygon are high) equates to a poor habitat quality rating (Fig. 10). As with all of the relationships depicted in this model, the distances in ï¬gure 10 are believed to be conservative approximations and should be reï¬ned with empirical data. ComputzLion of Overall Lvnx HSI The 3 primary components of the lynx habitat model; foraging (HSI lynx ,foodia denning (HSI lynx, den), and interspersion (HSIIynx, inter) are combined into one index value (HSI lynx) depicting overall habitat suitability for lynx in the 250 ha area. All components of the lynx model are weighted equally and deemed critical for a functional home range therefore a geometric mean was used to represent the ï¬nal HSI (Equation 8). The geometric mean causes the ï¬nal HSI to equal 0 if any of the components equal 0. These 250 ha areas can be subsequently aggregated into home ranges of differing functionality and used for resource planning and modeling (Roloff and Haufler 1997). Equation 8: )0.33 _ HSI (HSI _ lynx lynx, food x HSIlynx, den x HSIlynx, inter 126 Literature Review Adams, L. 1959. An analysis of a population of snowshoe hares in northwestern Montana. Ecological Monographs 29:141-170. Bailey, T.N., E.E. Bangs, M.F. Portner, J .C. Malloy, and R.J. McAvinchey. 1986. An apparent overexploited lynx population on the Kenai Peninsula, Alaska. Journal of Wildlife Management 50:279-290. Behrend, DP. 1962. An analysis of snowshoe hare habitat on marginal range. Proceedings of the Northeast Section of the Wildlife Society, Monticello, NY. 7pp. Bookhout, T. A. 1965a. The snowshoe hare in upper Michigan-its biology and feeding co- actions with white-tailed deer. Michigan Dept. of Conservation, Resource Division Report 38. l9lpp. Bookhout, T. A. 1965b. Feeding co-actions between snowshoe hares and white-tailed deer in northern Michigan. Transactions of the North American Wildlife Natural Resources Conference 30:321-335. Boutin, S. 1984. Effect of late winter food addition on numbers and movements of snowshoe hares. Oecologia 62:393-400. Boutin, 8., OJ. Krebs, A.R.E. Sinclair, and J .N.M. Smith. 1986. Proximate causes of losses in a snowshoe hare population. Canadian Journal of Zoology 64:606-610. Brand, C.J., L.B. Keith, and CA. Fischer. 1976. Lynx responses to changing snowshoe hare densities in central Alberta. Journal of Wildlife Management 40:416-428. Brand, C.J., and LB. Keith. 1979. Lynx demography during a snowshoe hare decline in Alberta. Journal of Wildlife Management 43:827-849. Brittell, J .D., R.J. Poelker, S.J. Sweeney, and GM. Koehler. 1989. Native cats of Washington. Unpublished Report, Washington Department of Wildlife, Olympia. Brooke, RH. 1975. Preliminary guidelines for managing snowshoe hare habitat in the Adirondacks. Trans. Northeast Seo., The Wildl. Soc. Fish Wildl. Conf. 32:46-66. Buehler, D.A., and LB. Keith. 1982. Snowshoe hare distribution and habitat use in Wisconsin. Canadian Field Naturalist 96:19-29. Cary, J .R., and LB. Keith. 1979. Reproductive change in the lO-year cycle of snowshoe hares. Canadian Journal of Zoology 57:375-390. Carreker, R.G. 1985. Habitat suitability index models: snowshoe hare. US. Fish and Wildlife Service, Biological Report 82. 21pp. Conroy, M.J., L.W. Gysel, and GR. Dudderar. 1979. Habitat components of clear-out areas for snowshoe hares in Michigan. Journal of Wildlife Management 43:680-690. 127 deVos, A. 1962. Changes in the distribution of the snowshoe hare in southern Ontario. Canadian Field Naturalist 76: 183-189. Dolbeer, RA. 1972. Population dynamics of the snowshoe hare in Colorado. Ph.D. Dissertation, Colorado State Univ., Fort Collins. 210pp. Dolbeer, R.A., and W.R. Clark. 1975. Population ecology of snowshoe hares in the central Rocky Mountains. Journal of Wildlife Management 39:535-549. F erron, J ., and J.-P. Ouellet. 1992. Daily partitioning of summer habitat and use of space by the snowshoe hare in southern boreal forest. Canadian Journal of Zoology 70:2178-2183. Gysel, L.W., and L.J. Lyon. 1980. Habitat analysis and evaluation. Pages 305-327 in SD. Sohemnitz, ed. Wildlife management techniques manual. The Wildlife Society, Washington, DC. Harestad, A.S., and FL. Bunnell. 1979. Home range and body weight - A reevaluation. Ecology 60:389-402. Hart, J .S., H. Pohl, and J .S. Tener. 1965. Seasonal acclimatization in varying hare (Lepus americanus). Canadian Journal of Zoology 43:731-744. Hoover, R.L., and D.L. Wills. (eds.). 1987. Managing forested lands for wildlife. Colorado Division of Wildlife, Denver. 459pp. Hulbert, I.A.R., G.R. Iason, D.A. Elston, and RA. Raoey. 1996. Home-range sizes in a stratiï¬ed upland landscape of two lagomorphs with different feeding strategies. Journal of Applied Ecology 33:1479-1488. Jones, K.J., R.S. Hoffman, D.L. Rice, et al. 1992. Revised checklist of North American mammals north of Mexico, 1991. Occ. Papers of Texas Tech Univ. 146: 1-23. Keith, LB. 1974. Some feature of population dynamics in mammals. Trans. Intern. Congress Game Biol. 11:17-67. Keith, L.B., and LA. Windberg. 1978. A demographic analysis of the snowshoe hare cycle. Wildlife Monographs 58. 79pp. Keith, L.B., J .R. Cary, O.J. Rongstad, and MC. Brittingham. 1984. Demography and ecology of a declining snowshoe hare population. Wildlife Monographs 90. 43pp. Keith, L.B., S.E.M. Bloomer, and T. Willebrand. 1993. Dynamics of a snowshoe hare population in fragmented habitat. Canadian Journal of Zoology 71 :1385-1392. Koehler, G.M. 1990. Population and habitat characteristics of lynx and snowshoe hares. Canadian Journal of Zoology 68:845-851. Koehler, G.M., and J .D. Brittell. 1990. Managing spruce-ï¬r habitat for lynx and snowshoe hares. Journal of Forestry 88:10-14. 128 Koehler, G.M., and KB. Aubry. 1994. Lynx. Pages 74-98 in LR Ruggiero, K.B. Aubry, S.W. Buskirk, L.J. Lyon, and W.J. Zielinski, tech. eds. The scientiï¬c basis for conserving forest carnivores: American marten, ï¬sher, lynx, and wolverine. USDA Forest Service GTR-RM-254. Litvaitis, J .A., J .A. Sherbume, and J.A. Bissonette. 1985. Influence of understory characteristics on snowshoe hare habitat use and density. Journal of Wildlife Management 49:866-873. Macfarlane, D. 1994. Appendix C: National Forest system status information. Pages 176-184 in L.F. Ruggiero, K.B. Aubry, S.W. Buskirk, L.J. Lyon, and W.J. Zielinski, tech. eds. The scientiï¬c basis for conserving forest carnivores: American marten, ï¬sher, lynx, and wolverine. USDA Forest Service GTR-RM-254. Martin, A.C., H.S. Zim, and AL. Nelson. 1951. American wildlife and plants. A guide to wildlife food habits. Dover Publ., Inc., New York. 500pp. McCord, C.M., and J .E. Cardoza. 1982. Bobcat and lynx. Pages 728-766 in J .A. Chapman and GA. Feldhamer, eds. Wild mammals of North America. John Hopkins University Press, Baltimore, MD. Monthey, R.W. 1986. Responses of snowshoe hares, Lepus americanus, to timber harvesting in northern Maine. Canadian Field Naturalist 100:568-570. Murray, D.L., S. Boutin, and M. O’Donoghue. 1994. Winter habitat selection by lynx and coyotes in relation to snowshoe hare abundance. Canadian Journal of Zoology 72: 1444- 1451. Nellis, CH, and LB. Keith. 1968. Hunting activities and success of lynxes in Alberta. Journal of Wildlife Management 32:718-722. Nylon-Nemetchek, M. 1999. Lynx (F elis lynx) of Riding Mountain National Park: An assessment of habitat availability and population viability. MS. Thesis, University of Manitoba, Winnipeg, Manitoba, Canada. O’Donoghue, M., and OJ. Krebs. 1992. Effects of supplemental food on snowshoe hare reproduction and juvenile growth at a cyclic population peak. Journal of Animal Ecology 61 :63 1-641. Orr, CD, and D.G. Dodds. 1982. Snowshoe hare habitat preferences in Nova Scotia spruce-ï¬r forests. Wildlife Society Bulletin 10:147-150. Paquet, P., and A. Hackman. 1995. Large oamivore conservation in the Rocky Mountains. A long-term strategy for maintaining free-ranging and self-sustaining populations of carnivores. World Wildlife Fund Canada, Toronto. 52pp. Parker, GR. 1981. Winter habitat use and hunting activities of lynx (Lynx canadensis ) on Cape Breton Island. Canadian Journal of Zoology 61:770-786. Parker, G.R., J.W. Maxwell, L.D. Morton, and GE]. Smith. 1983. The ecology of the lynx (Lynx canadensis) on Cape Breton Island. Canadian Journal of Zoology 61 :770-786. 129 Parker, GR. 1986. The importance of cover on use of conifer plantations by snowshoe hares in northern New Brunswick. The Forestry Chronicle 62:159-163. Pease, J.L., R.H. Vowles, and LB. Keith. 1979. Interaction of snowshoe hares and woody vegetation. Journal of Wildlife Management 43:43-60. Rogowitz, G.L. 1988. Forage quality and use of reforested habitats by snowshoe hares. Canadian Journal of Zoology 66:2080-2083. Roloff, G.J., and IE. Haufler. 1997. Establishing population viability planning objectives based on habitat potentials. Wildlife Society Bulletin 25:895-904. Ruggiero, L. F., K. B. Aubry, S. W. Buskirk, G. M. Koehler, C. J. Krebs, K. S. McKelvey, and J. R. Squires. 2000. Ecology and conservation of lynx in the United States. University Press of Colorado, Boulder, Colorado. Scott, DP, and RH. Yahner. 1989. Winter habitat and browse use by snowshoe hares, Lepus americanus, in a marginal habitat in Pennsylvania. Canadian Field Naturalist 103:560- 563. Severaid, J.H. 1942. The snowshoe hare, its life history and artiï¬cial propagation. Maine Dept. Inland Fish Game. 95pp. Sievert, PR, and LB. Keith. 1985. Survival of snowshoe hares at a geographic range boundary. Journal of Wildlife Management 492854-866. Sinclair, ARE, 0]. Krebs, J.N.M. Smith, and S. Boutin. 1988. Population biology of snowshoe hares. III. Nutrition, plant secondary compounds and food limitation. Journal of Animal Ecology 57:787-806. Sullivan, T.P., and RA. Moses. 1986. Demographic and feeding responses of a snowshoe hare population to habitat alteration. Journal of Applied Ecology 23:53-63. Sullivan, T.P., and D.S. Sullivan. 1988. Influence of stand thinning on snowshoe hare population dynamics and feeding damage in lodgepole pine forest. Journal of Applied Ecology 25:791-805. Taylor, WP, and W.T. Shaw. 1927. Mammals and birds of Mount Rainier National Park. US. Government Printing Ofï¬ce, Washington, DC. Telfer, ES. 1972. Browse selection by deer and hares. Journal of Wildlife Management 36:1344-1349. Thomas, J.A., J.G. Hallett, and MA. O’Connell. 1997. Habitat use by snowshoe hares in managed landscapes or northeastern Washington. Washington Dept. Fish and Wildlife, Unpublished Report 58pp. US. Fish and Wildlife Service. 1981. Standards for the development of habitat suitability index models. US. Fish and Wildlife Service ESM 103. 130 Vaughan, M.R., and LB. Keith. 1981. Demographic response of experimental snowshoe hare populations to over-winter food shortage. Journal of Wildlife Management 45:354-380. Ward, R., and OJ. Krebs. 1985. Behavioral responses of lynx to declining snowshoe hare abundance. Canadian Journal of Zoology 63:2817-2824. Washington Department of Natural Resources. 1996. Lynx habitat management plan. Washington Dept. Natural Resources, Olympia. 180pp. Washington Department of Wildlife. 1993. Status of the North American lynx (Lynx canadensis) in Washington. Unpublished Report, Washington Dept. Wildlife, Olympia. Windberg, L.A., and LB. Keith. 1976. Snowshoe hare population response to artiï¬cial high densities. Journal of Mammalogy 57:523-553. Wolfe, M.L., N.V. DeByle, C.S. Winchell, and TR. McCabe. 1982. Snowshoe hare cover relationships in northern Utah. Journal of Wildlife Management 46:662-670. Wolff, J .O. 1978. Food habits of snowshoe hares in interior Alaska. Journal of Wildlife Management 46:662-670. Wolff, J .O.. 1980. The role of habitat patchiness in the population dynamics of snowshoe hares. Ecological Monographs 50:1 1 1-130. 131 Table l. Documented browsed and un-browsed species and associated geographic region for the snowshoe hare (Martin et al.1951, Adams 1959, Telfer 1972, Keith et al. 1984, Rogowitz 1988, Sinclair et al. 1988, Sullivan and Sullivan 1988, Scott and Yahner 1989, Thomas et al. 1997). Snowshoe Hare Food Type Browsed Un-browsed or used < expected Tree Species: Tree Species: Scotch pine (Pinussylvestris) (NY) sugar maple (Acer saccharum) (NY, PA) striped maple (Acer pensylvanicum) (PA) beech (F agus grandifolia) (NY, PA) yellow birch (Betula alleghaniensis) (PA, Nova Scotia) black cherry (Prunus serotina) (NY, PA) grey willow (Salix glauca) (Yukon) Norway spruce (Picea abies) (NY) bog birch (Betula glandulosa) (Yukon) pin cherry (Prunus pensylvanica) (PA) lodgepole pine (Pinus contorta) (B.C., WA) red maple (Acer rubrum) (PA) red pine (Pinus resinosa) (NY) western white pine (Pinus monticola) white pine (Pinus strobus) (NY) (WA) white spruce (Picea alba) (NY) grand ï¬r (A bies grandis) (WA) paper birch (Betula papyrifera) (NY, Nova Scotia) western larch (Larix occidentalis) (WA) aspen (Populus spp.) (NY, MN) aspen (Populus spp.) (WA) willow (Salix spp.) (MN, UT) alder (Alnus spp.) (WA) Douglas-ï¬r (Pseudotsuga menziesii) (UT, WA, MT) red alder (Alnus rubra) (WA) Shrub Species: ponderosa pine (Pinus ponderosa) (MT) ocean spray western serviceberry (Amelanchier alnifolia) (MT, WA) ninebark (Physocarpus spp.) (WA) white cedar (Thuja occidentalis) (Nova Scotia) squaw berry balsam ï¬r (A bies balsamifera) (Nova Scotia) rose (Rosa spp.) (WA) bearberry (A rctostaphylus uva-ursi) (WA) Shrub Species: Oregon grape (Berbis repens) (WA) southern arrowwood (Viburnum dentatum) (NY) boxwood (Pachystima myrsinites) (WA) bristly dewberry (Rubus hispidus) (NY, PA) blackberry (Rubus allegheniensis) (NY, PA) beaked hazelnut (Corylus comuta) (MN) snowberry (Symphoricarpos spp.) (UT) Oregon grape (Berberis repens) (WA, MT) ceanothus (Ceanothus spp.) (WA) bearberry (Arctostaphylos uva-ursi) (WA, MT) birchleaf spiraea (Spiraea betulifolia) (MT) rose (Rosa spp.) (WA) huckleberry (Vaccinium spp.) Table 2. Effect of grid cell size on the computation of mean nearest distance to non-lynx habitat. No. Points Mean Nearest Grid Cell Size (113 Generated Distance (m) 50x50 19,939 1,388 100x100 4,977 1,384 500x500 197 1,356 1000x1000 47 1,303 132 HSI Score Browse Horizontal Cover 0 11111 . . r 0 10 20 30 40 50 60 70 80 90 100 Percent Cover HSI Score Browse Vertical Cover 0 10 20 30 40 50 60 70 80 90 100 Percent Cover 133 Figure 1. Relationship between horizontal cover of browse and HSI score for 0-1 and 1-2 m height strata. Line equation between 10 and 35% is y = 4x—40. Figure 2. Relationship between vertical cover of browse and HSI score for 0-1 and 1-2 m height strata. Line equation between 20 and 80% is y = l.666x-33.33. Horizontal C over HSI Score 10 20 30 40 50 60 70 80 90 Percent Cover 100 HSI Score Vertical Cover 10 20 3O 40 50 60 70 80 90 100 Percent Cover 134 Figure 3. Relationship between horizontal cover and HSI score for 0-1 and 1-2 m height strata. Line equation between 40 and 90% is y = 2x-80. Figure 4. Relationship between vertical cover and HSI score for 0-1 and 1-2 m height strata. Line equation between 40 and 90% is y = 2x-80. Viability Relationship for Snowshoe Hares 18 16 '1 ehrend (1962) a 14 a go 12 - g 0 Sievert and Keith (1985) b M 10 d D E: E8 _ Dolbeer and Clark (1975) 6 ‘ Cary and Keith (1979) c 4 .. 2 .. 0 I l I 0 5 10 15 20 I I T l 1 0 25 50 75 100 Offspring/Y ear Habitat Quality Score a Annual productivity assumed to be 2 in that females are only replacing themselves. b Home range presented as >10 ha in literature. Here assumed as 11 ha. c Inferred home range size based on allometric theory. Figure 5. Viability relationship for snowshoe hare developed for the lynx habitat model. Viability threshold was established at 12 young/year (60 score). The marginal threshold was established at 5 young/year (25 score). 135 Lynx Foraging HSI Score 100 90 801 70* 60“ 50 40‘ 30 20m 10* HSI Score 0 7 ï¬T 0 10 20 30 40 50 60 70 80 90 100 # of Snowshoe Hare Home Ranges in 250 ha Vertical Security Cover 100 90 80‘ 70‘ 60“ 50‘ HSI Score 40 30‘ 20’ 10‘ 0 0 10 20 30 40 50 60 70 80 90 100 Vertical Cover 0-1 m (%) 136 Figure 6. Relationship between the number of snowshoe hare home ranges in a lynx allometric home range and the lynx forage HSI score. Line equation between 0 and 90 home ranges is y = 1.1x. Figure 7. Relationship between summer security cover and HSI score for the 0-1 m height strata. Line equation between 20 and 60% is y = 0.25x-50. Denning Interspersion 100 Figure 8. Relationship between 90- the average distance to suitable denning sites in a 250 ha lynx 80 home range and HSI score. Line 70’ equation between 400 and g 60 1,750 m is y = -.O74lx + 129.64. 0 2 50 m :5 40 30 20 10 0 0 400 800 1200 1600 Average Distance to Suitable Den Site (m) 250 ha Lynx Home Range 100m [:I = Non-lynx habitat Figure 9. Calculating the average distance to non-lynx habitat using a 100x100 m sample grid. Distance from each grid intersection to the nearest non-lynx habitat is measured. 137 HSI Score .1) O Non-lynx Habitat Interspersion UrCh CO r—Nw Ooco 400 800 1200 1600 0 Average Distance to Non-lynx Habitat (In) 138 Figure 10. Relationship between the average distance to “non-lynx" habitat in a 250 ha lynx home range and HSI score. Line equation between 300 and 1500 m is y = 0.08333x-25. APPENDIX B: SCHEMATIC OF PROCESSES USED FOR THE LYNX MODEL Combine spatial layers describing _ attributes of vegetation and physiography Ecologlcal Land Classiï¬cation /urrent land cover / / : historic land cover / soils I . Quantify habitat quality across ecological land / ecoregl‘ms classes for each lynx model component FORAGING I DENNING I I INTERSPERSION I Hare Habitat Hare Home Ranges V W Foraging Quality Denning Quality Interspersion Quality /|// //./ I 1 Combine grids to compute ï¬nal grid of lynx habitat quality Lynx Habitat Quality APPENDIX C: ACCURACY ASSESSMENT OF IFMAP LANDCOVER A land-cover dataset for Michigan was created through the 2001 IFMAP project of the Michigan DNR, and used in the Michigan GAP Analysis Project (Donovan 2005). The map is widely used as the source of land-cover information for the state, with an estimated accuracy of 87% at Anderson Level II and 37-87% at level III, and can be obtained from the Geographic Data Library of the Michigan DNR (MiGDL) on their website (<http://www.mcgi.state.mi.us/mgdl/>). The version that is released by the MiGDL contains Anderson Level III cover classes for deciduous (e.g., northern hardwoods, oak, aspen, other upland deciduous, mixed upland deciduous, lowland deciduous) and coniferous (e.g., pines, other upland conifers, mixed upland conifers, lowland coniferous, lowland mixed) forests. A problem with using these cover classes for habitat assessment is the variability in their accuracy, which could lead to erroneous conclusions about the distribution of forest resources on the landscape. Subdedi (2005) examined the accuracy of the IF MAP forest classiï¬cations on public lands using Forest Inventory and Analysis (F IA) plots collected by the Forest Service as ground references. He found overall accuracies of 63.6% on state forests and 64.8% on national forests, when using the Anderson Level 111 categories provided by the IFMAP layer. User’s accuracies ranged from 17.6-60.0% for aspen and 23.5-45.7% for oak on public lands. These results suggest that the map predictions of certain forest types may not be very useful at the Level III category. I conducted a similar analysis to that of Subdedi (2005), using F IA plots on all ownerships as ground references for the IF MAP forest classes in the entire UP. At the 140 time of the analysis, the plots that were available included 80% of the 6th cycle for Michigan. Plots were ï¬ltered to include only those that were forested (>25% tree cover) and containing a single condition, to ensure that canopy information was consistent throughout the 4 subplots. This resulted in 2,328 plots spread across the UP. The tree information for each plot was used to calculate the percent cover of canopy trees for each species, and plots were labeled according to the classiï¬cation rules used by the IF MAP project (MDNR 2001). The FIA plots were overlaid with a grid of the land-cover dataset by the North Central Spatial Data Services (SDS) of the Forest Service, and the cover type for each plot was determined by the majority of pixels within a 3x3 window around each plot location. An error matrix was constructed to compare the IF MAP classiï¬cations for each plot, as determined by the map and tree list data (Table CI). The overall accuracy of the IF MAP forest classiï¬cation was estimated at 45.9%, with relatively poor user’s accuracies for most forest classes other than northern hardwoods (79.8%) (Table C.1). After grouping IF MAP classes according broader categories of forest type (Table 3), overall accuracy increased to 74.3%. Our results ï¬irther illustrate the potential problems with using the IF MAP land- cover dataset for mapping forest types in the UP. The heterogeneity of forest types across the landscape of the UP might explain the difï¬culty in accurately predicting their occurrence using satellite imagery. Aspen forest types that are being succeeded by conifers (e.g., balsam ï¬r) may also present problems during interpretation. My analysis utilized a ï¬ltered set of F IA plots, which may not have provided the best representation of reference points. Using all FIA plots may have yielded different results, though plots with multiple conditions would likely have decreased the accuracy of IFMAP predictions. 141 .88:on mmooseoa n $5 59508 Meow: I. <3 a mo 0.? 0.: n: 3 .1 _ v. 6 ed 3N od Re 22 DE 2 2 3N N: a: a. 3 mm NoN 5 v as 33 0.0 a: ON NN M: N N N S 2 _ a 2 mm season a: _.: a _ N N N _ use Bee 2332 N New 3. cm oNN NM «N w NN a. m N 2 assess 2232 N 3N we N m N e N 2 N 2 32.20% Baas N a: mNN a v 3 z N m N an N NN use Exec use: NN 9o 2 m N N N N aces. 22% Bees _N 2.: K N v a. 2 w 2 m 2 N 2 V eees 22% 55o 8 NNN 32 N2 2 _ N ON 2 NN X z o _ 2 85a a 3. N N N 8828.0 22% 33%. M: 3: EN 4 N 2 N N e 4 E a. 3 can 2 2c ... _ N v.8 2 NE N3 _ 2 MN _ a. 2 5 on _ :3 . €853: 8058: E .202 $5 =29 N N N NN _N ON a M: 8 2 E He,8 $25 88 , - aoa <E ‘ ‘ .Esmchm 5&3 05 E 982:5 <5 98 n22": c0253 gonna—mama? 388 we x53: Sum 4.0 053. 142 APPENDIX D: CALCULATING HORIZONTAL COVER WITH THE STAND VISUALIZATION SYSTEM Horizontal cover has been demonstrated as an important habitat component for snowshoe hare (Brocke 1975, Wolff 1980, Litvaitis et al. 1985, Parker 1986, Hodges 2000), and many other wildlife, as it can provide a range of life history requirements, including thermal and escape cover, and available browse. Common measurements of horizontal cover involve counting the proportion of covered grid squares on a proï¬le board held perpendicular to the ground and viewed from a speciï¬c distance. This type of sampling is subject to large measurement error on the part of the observer (Collins and Becker 2001), as the decision to count a grid square as “covered†is subjective at best. Forest inventories do not typically measure horizontal cover as a stand attribute, which limits their use in assessing habitat for species that require it. The Stand Visualization System (SVS) (McGaughey 1997) is a software program that creates 3—dimensional renderings of 0.40 ha stands (Figure D. 1), using forest inventory data for input. SVS has the ability to interpret a number of common sampling Figure D.1 SVS drawings of a mature red pine stand with hardwood understory. 143 designs that are used for forest inventory, and when incorporated with the Forest Vegetation Simulator (F VS), stands can be grown and subjected to harvest scenarios over time. SVS contains a number of tree form ï¬les that outline parameters used to create a realistic diagram for each tree species. Parameters can be manipulated by the user and include the shape of the crown, the number, positioning, and color of branches and leaves, and the manner in which these parameters change with height and diameter. For example, the leaves on deciduous stems can be removed to simulate the appearance of leaf-off season. The proï¬le view illustrates all trees within a user-speciï¬ed strip width, which can simulate the horizontal cover provided at the speciï¬ed distance. I developed a method to estimate horizontal cover (within 0-1 m and 1-2 m from the ground) using the proï¬le view of an image in order to assess the quality of forest stands for winter hare habitat. SVS images were created for a series of FIA plots using the post-processing ï¬inction of FVS, and tree forms for each deciduous species were altered to remove leaf cover. The view azimuth was changed to 0° for the perspective view, which placed the range poles at either end of the proï¬le image. Images were saved as 1280 x 960 Windows bitmaps in SVS. A single image was examined to determine the pixel coordinates of a box encompassing each height stratum (Figure D2). The 3 m increments on the range poles measured 28 pixels, which meant the height of each 1 m stratum would approximate to 9 pixels, with a length between range poles of 642 pixels. These dimensions should be consistent for all SVS images that are saved as 1280 x 960 bitmaps. Each image was converted to a monochrome portable bitmap and read into R as a text ï¬le. The ï¬le contained a matrix of binary values which represented the black (1) 144 and white (0) pixels from the bitmap. Dividing the sum of values within the 9 rows and 642 columns for each stratum by the total pixels (57 78), resulted in the estimate of cover. The following R code was used for the calculations: svs <- read.table("standl.pbm", skip=32189, nrows=943) stratuml <- sum(sum(svs[64l,24:36]), sum(svs[642:658,l:36]), sum(svs[659,l:l7]), sum(svs[677,8:36]), sum(svs[678:694,l:36]), sum(svs[695,l]), sum(svs[712,28:36]), sum(svs[713z729, l :36]), sum(svs[730,l :21]), sum(svs[748,12:36]), sum(svs[749:765,1:36]), sum(svs[766,1:5]), sum(svs[783,32:36]), sum(svs[784:800,l:36]), sum(svs[801,l :25]), sum(svs[8l9,l6:36]), sum(svs[820:836,1:36]), sum(svs[837,1:9]), sum(svs[854,36]), sum(svs[855:87l,l :36]), sum(svs[872,l:29]), sum(svs[890,20:36]), sum(svs[89l :907,1:36]), sum(svs[908,l : 13]), sum(svs[926,4:36]), sum(svs[927:942,1:36]), sum(svs[943,1:33])) / 5778 stratum2 <- sum(sum(svs[l,24:36]), sum(svs[2: 18,1:36]), sum(svs[l9,l:l7]), sum(svs[37,8:36]), sum(svs[38:54,l :36]), sum(svs[55,l]), sum(svs[72,28:36]), sum(svs[73z89, l :36]), sum(svs[90,1:21]), sum(svs[108,12:36]), sum(svs[109: 125,1 :36]), sum(svs[126,l :5]), sum(svs[ l43,32:36]), sum(svs[144: 160,1 :36]), sum(svs[ 16 1 ,1 :25]), sum(svs[l79,l6:36]), sum(svs[l80: l96,1:36]), sum(svs[l97,1:9]), sum(svs[214,36]), sum(svs[215:23l,1:36]), sum(svs[232,l :29]), sum(svs[250,20:36]), sum(svs[251:267,1:36]), sum(svs[268,1:13]), sum(svs[286,4:36]), sum(svs[287:302,1:36]), sum(svs[303,l :33])) / 5778 stratuml = 0.5896504 (59%) stratum2 = 0.5705434 (57%) 145 \ . r ( .lv ‘ . ‘ . r : . y . . .. .' ..11II,IllI... E . ‘_ 1 IE! . .ï¬ I F ."1 .lwt. .I .. Figure D.2. Proï¬le view of stand with a red box around each height stratum. 146 APPENDIX E: LANDSCAPE-SCALE MODELING OF FOREST STRUCTURE USING LANDSAT TM IMAGERY AND FIA DATA The mapping of forested land using remotely-sensed imagery has become a standard procedure in the management of natural resources at the landscape scale. These classiï¬ed maps of land cover contain varying levels of detail and accuracy, depending on the extent of the mapped landscape and the resolution of the imagery, but the amount of information they can relay is limited. Evaluations of wildlife habitat often require ï¬ne- scale measurements of vegetation structure, which are not provided by land-cover datasets. Obtaining this type of information for forest inventories can be costly for large landscapes, involving a combination of high resolution aerial photography and ï¬eld- based measurements at the stand level. Even when this information becomes available, the temporal resolution of data acquisition may be too coarse to capture alterations in vegetation structure due to growth and disturbance through time. Our objective was to map the local variation in forest structure across the UP, which required the integration of remotely-sensed imagery with ï¬eld surveys to utilize relationships between vegetation structure and spectral reflectance from the ground. The k-nearest neighbor (KNN) method assigns attribute values to non-sampled pixels from those that are sampled (i.e., have an associated ground plot) based on the distance between the pixels in a multi-dimensional feature space, as deï¬ned by the combination of spectral values in each wavelength band of a remotely-sensed image. This method has been used to predict forest attributes across large landscapes (Tomppo 1991, F azakas and Nilsson 1996, Tokola et al. 1996, Franco-Lopez et al. 2001), and in the Great Lakes 147 region with Forest Inventory and Analysis (FIA) plots serving as ground references (McRoberts et al. 2002, Haapanen et al. 2004). I acquired Enhanced Thematic Mapper Plus (ETM+) imagery from the Landsat 7 satellite to model forest structure using the KNN method, with ground references provided by the FIA plots collected between 2000- 2003 for the 6th inventory of Michigan. The UP required 8 scenes of Landsat imagery for complete land coverage, including row 28 of paths 21-22, rows 28-29 of path 23, and rows 27-28 of paths 24-25 (Figure E. 1). Images were acquired from a range of dates between 2000-2002 due to the limited availability of cloud-free imagery; images 5 and 6 were from late August (8/20), 0 25 50 100 150km \ Figure E.1. Location of the 8 scenes of Landsat 7 imagery acquired for the Upper Peninsula, including paths 21-25 (east to west) and rows 27-29 (north to south). ---1 l 1' ‘V l 148 images 3 and 4 from early September (9/08), images 1, 2 and 8 from mid to late September (9/17-9/29), and image 7 from early October (10/1 1). Each image was rectiï¬ed and geo-referenced to zone 16 of the UTM system, using the GRS 1980 spheroid and NAD83 datum. A vector layer of roads from the Michigan Center for Geographic Information (MCGI) was used for geo-referencing. Small clouds or areas of haze were masked out and each image was clipped to the boundaries of the UP. Three pairs of images were successfully mosaicked together, and the ï¬nal 5 images were named as follows: EAST = images 1 + 2, MIDEAST = images 3 + 4, MIDWEST = images 5 + 6, WESTl = image 7, WEST2 = image 8. The images were transferred to Geoff Holden (USF S) at the Spatial Data Services (SDS) of the North Central Research Station, where FIA plot locations were overlaid with the spectral imagery and values of plot attributes were imputed across the landscape using the KNN method. Spectral values from bands 1-5 and 7 were used to deï¬ne the feature space, and a 3 x 3 mean ï¬lter was applied to each image to match the 4-subplot sampling design of the F IA plots. The number of nearest neighbors (k) was set to 5 and the distance metric was Euclidean. The prediction accuracy of the model was estimated by holding out 40% of the FIA plots from the training set and comparing the observed with predicted pixel values using error matrices and linear regression. In the ï¬rst iteration, the plot attribute was categorical, describing the structural stage of each plot as deï¬ned by the Forest Vegetation Simulator (FVS). The stages were based on the size of trees and number of height strata, with dbh thresholds for sapling/pole and pole/large set at 12.7 cm (5 in) and 30.48 cm (12 in), respectively, and percent cover of valid height stratum set to 20 (Crookston and Stage 1999). Several 149 combinations of structural stages were used to vary the number of attribute categories (i.e., 4 stages versus 6 stages), but the KNN procedure was unable to produce an overall classiï¬cation accuracy >40%. The second iteration involved plot attributes that were continuous, including horizontal cover (as measured by FVS) and basal area. Horizontal cover was corrected for normality using the arcsine transformation, aï¬er converting it to a proportion, while basal area was determined to be normally distributed. Initial examination of basal area predictions on a single image were comparable to that of horizontal cover; therefore, basal area was not assessed any further. Predictions of horizontal cover were developed for each of the 5 images; all had poor relationships between predicted and observed values, with most having an adjusted R2 <0.10 (Figure E.2). Our inability to accurately predict horizontal cover using the Landsat imagery could be attributed to a number of different reasons. The data provided by the 6 ETM+ bands were most likely insufï¬cient for capturing any differences in reflectance caused by vegetation structure, especially that of the understory. The analysis was limited to one date for each region due to the availability of satellite imagery (determined by acquisition cost and cloud cover), and the range of acquisition dates (20 August — 11 October) could arguably have spanned multiple seasons. Using 3 or more dates of imagery that spanned the growing season of vegetation, as well as ancillary GIS data that stratiï¬ed the landscape based on multiple ecological characteristics (i.e., topographical position, ecoregion), may have improved the predictions (Franco-Lopez et a1. 2001). The KNN method assumes homogeneity in forest composition across the analysis area, which was not the case for most of the UP. In spite of this, the greatest prediction accuracy came 150 1.5 observed 1.0 0.5 r r 1 0.0 0.5 1.0 1.5 predicted 1.5 observed 1.0 0.5 Figure E.2. 0.5 predicted observed MIDEAST R2 = 0.08 P < 0.001 o. _ O O . o 1 1 1 1 0.0 0.5 1.0 1.5 predicted WESTI 2— 112-0.01 . P = 0.15 ° 5.1537, G. .. o ' —‘ o... o g‘ . . o C :0 V! _. . . . C q .. o I 1 I 1 0.0 0.5 1.0 1.5 predicted l j I 0.0 0.5 1.0 1.5 predicted Relationship between predicted and observed values of arcsine transformed horizontal cover for each of the 5 spectral images in the Upper Peninsula. 151 from the EAST image, which had a far greater contrast in forest types than the areas captured by the western images. A stratiï¬cation of the landscape, other than that created by the Landsat scenes, was not practical due to the low sampling intensity of the F IA plots. A potential source of error in the methodology was in matching F IA plots to pixels within the satellite image, as plots cannot be expected to lie perfectly centered within the 3 x 3 pixel windows, which may effect the assignment of spectral values. Finally, the measure of horizontal cover used as the dependent variable was an estimated value with an unknown accuracy; any hypothesized relationship between horizontal cover and the spectral reflectance of the vegetation would rely on an accurate estimate. Earlier examinations of the KNN method for classifying forest structure have had some success (Franco-Lopez 2001), and the method has proven useful in predicting compositional attributes across large landscapes (MoRoberts et al. 2002). The application of this method to predictions of understory structure using Landsat imagery is dependent on a relationship between understory vegetation and the reflectance of the canopy. The use of high resolution imagery and/or light detection and ranging (lidar) data in conjunction with Landsat imagery may provide a more accurate means of predicting forest structure across large landscapes (Lefsky et al. 1999). 152 APPENDIX F: METADATA FOR SPATIAL LAYERS Type Data Layer Description Source input IFMAP/GAP Current land cover for Michigan Department of Natural Resources Land Cover Michigan using 1997-2001 (MN DR) Geographic Data Library: satellite imggery. <http://www.mcgi.state.mi.us/mgdl/> Michigan’s Vegetation type boundaries Presettlement interpreted from GLO Vegetation Surveys 1816-1856 for Michigan STATSGO soils Map units of soil for Michigan associations for major soil groups in Michigan Ecological Delineations of regions with David T. Cleland Subregions of the similar climate and USDA Forest Service United States physiography across the Southern Research Station United States 5985 County Road K Fire regime Delineations of regions with Rhinelander, WI 54501 categories for the a similar risk of ï¬re under email: dcleland@fs.fed.us UP natural ï¬re disturbance <http://www.ncrs.fs.fed.us/gla/> regime for the Upper Peninsula of Michigan output ecoregion-ELUs Grid of ecological land units Daniel W. Linden for the UP describing similar habitat types and tree canopy species, with ecoregion identiï¬cation in the Upper Peninsula of Michigan lynx foraging quality in the UP Map of habitat suitability index (HSI) scores for lynx foraging quality in the Upper Peninsula of Michigan lynx denning quality in the UP Map of HSI scores for lynx denning quality in the Upper Peninsula of Michigan lynx interspersion Map of HSI scores for lynx quality in the UP interspersion quality in the Upper Peninsula of Michigg lynx habitat Map of HSI scores for quality in the UP overall lynx habitat quality in the Upper Peninsula of Michigan Dept. of Fisheries and Wildlife 13 Natural Resources Bldg. Michigan State University East Lansing, MI 48824 email: lindend 1 @msu.edu Henry (Rique) Campa, III Dept. of Fisheries and Wildlife 13 Natural Resources Bldg. Michigan State University East Lansing, MI 48824 email: campa@msu.edu Dean E. Beyer, Jr. Wildlife Division Michigan Dept. of Natural Resources 213 West Science Building Northern Michigan University Marquette, MI 49855 email: dbeyer@nmu.edu 153 IIIIIIIILIIIIIIIIIIIIII