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... . :gafi‘
: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.... ..
..

ivynfiflngtfllinwrmfla
.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
fiduflhflw
«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 fulfillment
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 Afi‘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

 

 

Mbfltfias

 

 

 

 

 

 

 

 

 

 

 

 

 

 

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 fiJlfillment 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 identified 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 confidential
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
gratefiil. 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 Classification 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 identified 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 classifications 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 fire
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 fires 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

fire 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 fire disturbance with the highest (top) and lowest
(bottom) amount of early successional vegetation .................................... 75

Lynx denning quality during the 2 simulation years of
presettlement fire disturbance with the highest (top) and lowest
(bottom) amount of early successional vegetation .................................... 76

Lynx interspersion quality during the 2 simulation years of
presettlement fire 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 fire 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 fire 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 defined 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

Profile 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 final 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 significant 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 significant. In light of the final 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 undefined. 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 fitness 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 misidentification of bobcat (Lynx rufus) for

lynx, but several verified 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 “definitely 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 identification 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 define
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, fire
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 verified 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 rufils]) 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 significant 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 sufficient 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 fire 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 fire suppression and exclusion in the contiguous United
States, along with naturally longer fire 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 sufficient 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 fire
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 benefit
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
significant 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 quantified
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 significant 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 defined 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 defined 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 fires 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 fire
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 fir (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 fir, cedar,
hemlock, and white spruce (Picea glauca). Xeric uplands are populated by conifers such
as jack pine (Pinus banksiana), red pine (Pinus resinosa), fir, 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 fisher
(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 classification, 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 specifically 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 definitions of palatable browse (see review in
Hodges 2000). Therefore, hare habitat quality was indexed by 2 measures: (1) horizontal
cover, defined 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,
defined 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 fiom 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 fitness) at an average
quality 260 and marginal (i.e., not contributing to population fitness) 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 defines 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 (defined 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., fields, 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 identified 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 fitness
indicators (e. g., survival, litter size, pregnancy rate) (see Figure 3: Roloff and Haufler
1997) demonstrated in previous field studies. The threshold represents the point at which

lynx home range size no longer increases rapidly with decreased habitat quality, and

22

fitness 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 configuration 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 quantified and mapped across the landscape for
the assessment of habitat suitability. Roloff (2003) suggested the use of an ecological
land classification that described biotic and abiotic factors of the landscape to properly

stratify the variation in vegetation attributes defining habitat quality in the H81 models.

ECOLOGICAL LAND CLASSIFICATION FOR THE UPPER PENINSULA

The ecological land classification used in this study was based upon the “Habitat
Classification System” constructed by Coffman et al. (1984), describing “habitat types”
(Daubenmire 1966) of the UP that were defined 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 specific 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

classifications 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
classification facilitated the stratification of vegetation attributes necessary for the lynx
model into ecologically meaningful and homogenous classes that were developed
specifically 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. Coffinan 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 specific
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 defined 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 fine 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 identified 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/fir
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/fir,
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/fir
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/fir

occurring on a variety of poorly drained soils. Finally, the Picea habtype consisted of

25

Table 1. Habitat types originally described by Coffinan 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, fir
TsMa (l) aspen, white birch (1) white pine hemlock, maple,
(2) jack, red pine (2) white spruce, fir beech, yellow birch
Acer aspen sugar maple, basswood, sugar maple, hemlock
white ash, beech
TsTh-dry aspen, white birch, (1) white spruce, fir hemlock, cedar, fir,
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, fir, black spruce, cedar, hemlock, fir
red maple tamarack
Picea lowland brush, white birch, black spruce, tamarack, cedar, black spruce, tamarack

red maple

fir, 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 fir occurring on
very poorly drained organic soils and peatlands. The combination of a habtype and seral
stage (as defined by overstory species) resulted in an ecological land unit (ELU), which
provided the basis of the ecological land classification necessary for the lynx model.
The classification system developed by Coffman et al. (1984) was intended as a
field guide for in situ classification 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 configuration 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., fire 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 defined 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 specific vegetation
attributes required by the lynx model. A final data layer consisted of unclassified 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 classification 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 classification of 3 seasons
(spring, summer, fall) of Landsat satellite imagery collected at a 30 m resolution between
1999-2001. Forested pixels were defined to have 225% canopy cover; forest
composition was defined according to the dominant species (260% cover) in the canopy,
with “mixed” categories for stands without dominant species. Land cover was classified
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 findings from the accuracy assessment.

Although the loss of detail made it difficult 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 identified 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 significant 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 identified 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 reclassification were used to provide the
most informative description of the vegetation. The first 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 final 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-fir-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
classification 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 misclassified 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
reclassification 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 Official 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 reclassified 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 refined 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
identification (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_ .wfi .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 .wfi .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 Esmfiaom Saab 05 E motowfimo 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 .Auaemfitg 9.353 059 xofl. u E. .Emnmheng

0.05.5 x0280: 5830 n Em .Auzuutmfia 6.5.50 083an M 3m A3525 08.25 83% 0.003 H mm .Acmfiaé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“: oomhfim

«cozfiowo> 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 stratification approach (Smith 2002). Remote sensing imagery is used to identify
forested land across a state, which is then systematically field 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 field
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 fixed 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 fixed 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 define 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 final 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://fia.fs.fed.us/) in two different formats: (1)
spreadsheets that could be imported to a relational database; and (2) text files 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 specific 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 filtered 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 definition 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 filtering
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 filtered datasets to verify that no bias was introduced
through the filtering 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 definitions are used to define appearances for different
species and can be controlled by the user; I altered plant form definitions in SVS by
removing leaves from deciduous species to simulate winter conditions. For each plot, the

profile-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 significant 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 files, 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 defined 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 fit 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 specific 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 defined 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 confidential by Congress in 2000, in accordance with amendments to the Food
Security Act of 1985, to ensure that private landowners are not identified 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. Confidentiality 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 final

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 defined 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 classification of sample plots according to species
abundance and site fidelity. 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 (identified 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 final 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 rectified 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 confidentiality 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 defined by the combination of spectral
values from each wavelength band. The analysis was performed using functions within
ERDAS Imagine; a 3x3 mean filter 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 fire disturbance and subsequent forest
age classes across a landscape through time (Wimberly et a1. 2000). Historical fire
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 fire 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 fire within and across scales. The raster layers include: (1) climate
zones, which delineate regions with specific fire rotations; (2) topographic zones, which
portray local differences in topography that affect susceptibility to fire; (3) summary
zones, which identify optional subsets of the study area for which summaries are desired;
and (4) a buffer zone, which identifies areas from which fires 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, fire shape calibration,
wind direction and intensity, fire susceptibility within topographic zones, age-class
definitions, length of fire rotation, and the mean and standard deviation of fire 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 fire 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
fire 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 fire 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 fire rotation categories, and used
GLO survey notes to quantify the natural frequency of stand-replacing fires 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 fire rotations summarized by BPU and ecoregion section (Cleland et al.
2005a). The accuracy of fire rotation estimates decreases with the size of the landscape
unit (Cleland et al. 2005a); therefore, I condensed BPU-section groups that had similar
fire fiequencies, and recalculated fire rotations based on the fire frequency and area of
each new BPU (Figure 7). The mean and standard deviation of fire 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 fire size parameters to compute the mean number of
fires per decade within a climate zone, which was simulated for each decade as a Poisson
random variable. Individual fires were randomly sized according to a lognorrnal

probability distribution, calculated from the input fire 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 identified 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 fire (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 filter 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 fires from outside the UP. Age
classes were defined by the amount of time since the last stand-replacing fire; I defined
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 fire
severity was set at a constant high level to simulate only stand-replacing fires. 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 simplified
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 first 10 years of a severe fire. 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-fir 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 fire (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-fir 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 simplified 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 qualified 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 sufficiently accumulated into hare home ranges using Plume’s
(2005) GIS fimction, 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 fires 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 classification, 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 classification: 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 fire 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 quantified 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 first
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 significant
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: Eswfi

 

 

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 wfiwfio: 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 fire 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 fire 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 significantly
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 wEBfi me @2526 05>» mownfi qun 2%?» .325 553 wqunwocwwm com corona m5 $88 08:3 0% was:
Govfioxo mums—m5 ng22 mo Esmfiaom Saab 05 E 22:28 Beta .28 toga—=86 muwfiu 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 fire
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 fire
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 fire
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 fire 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 fire
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”) ifi/ t

Figure 20. Overall lynx habitat quality during the 2 simulation years of presettlement fire
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 fire 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 fimction.

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 defined as having an average quality >70.

82

DISCUSSION

Roloff’ s (2003) model quantified 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
field 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
fire 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 fire 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 fire 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 fires have the ability to create adequate early
successional vegetation types for hare, and each has been shown to benefit 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 fire 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 definition 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 significantly
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 fine 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 confidentiality,
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 filtering 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 find significant 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 filtered 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 confidence (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 fire 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 fire disturbances in the region. Also, I assumed that fire
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 fires, but

87

while small surface fires, 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 fires (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 difficult 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 first 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 identified 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 defined 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 quantified 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 identified 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 significantly 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 significantly 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 difficult task, since the relative cost of
traveling through an area for a species cannot be readily quantified. 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 defined 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 significantly 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 fire, 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 qualified as suitable habitat is necessary to
validate the model outputs, but also to determine whether other factors (i.e., interspecific
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 classification. 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 fire 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
classification 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 Office 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: Identification 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 fires 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 classification 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-fir 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 scientific
basis for conserving forest carnivores: American marten, fisher, 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 classification 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. Stratified 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-fir 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 classification 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://factfinder.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, final rule. Federal Register 65:16052-
16086.

U.S. Forest Service. 2006. Forest Inventory and Analysis Program.
<http://www.fia.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, modifications to the mathematical equations
were implemented in a manner consistent with how the model was applied in a verification study in
Manitoba, Canada (Nylen-Nemetchek 1999). These modifications 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
identified 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 rufiis), 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 field
verified (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 fitness.

Overview
Critical considerations prior to running the lynx model are the resolution, accuracy, and

precision of the land classification system and associated vegetation attribute information for the

113

planning landscape. The land classification must be capable of characterizing vegetation
structures and spatial arrangements at a resolution compatible with lynx and snowshoe hare
habitat use. The ideal stratification is a stand-based (minimum mapping unit around 2 ha)
ecological classification system that integrates existing vegetation conditions and site potentials
(e.g., geology, soils). An ecological classification 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
stratification 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 specifically 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
fitness 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 significant relationships
between hare population indices and the horizontal and vertical cover of understory vegetation.
Since few snowshoe hare studies have been conducted in the Pacific 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 quantified 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 findings
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 figure 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 figure 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 defined 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 defined 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 first 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 final 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 final 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 sufficient 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 quantified 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 figure 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 classified 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 defined 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 identified 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 identified as potential denning sites in “Phase I” must be adjacent to some form of
"suitable lynx habitat" (i.e., habitat identified 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 satisfied, the map polygon is a potential

denning site and an assessment of downed woody debris is performed. Map polygons not
identified as potential denning habitat through the first 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 identified 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 (verified 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 defined 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 identified 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 configuration 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
sufficient 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 figure 10. Lynx will traverse long distances to fulfill 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 figure 10 are believed to be conservative approximations and should

be refined 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 final HSI (Equation 8). The geometric mean causes the final 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 stratified
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-fir 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 scientific basis for conserving

forest carnivores: American marten, fisher, 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
scientific basis for conserving forest carnivores: American marten, fisher, 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-fir
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 artificial 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 Office, 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 artificial 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 fir (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-fir (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 fir (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 fiT

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
Classification

/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 final
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 classifications 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 filtered 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 classification 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
classifications for each plot, as determined by the map and tree list data (Table CI).

The overall accuracy of the IF MAP forest classification 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 fiirther 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 difficulty in accurately predicting their
occurrence using satellite imagery. Aspen forest types that are being succeeded by
conifers (e.g., balsam fir) may also present problems during interpretation. My analysis
utilized a filtered 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 profile
board held perpendicular to the ground and viewed from a specific 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 files 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 profile view illustrates all trees within a user-specified strip width,
which can simulate the horizontal cover provided at the specified distance.

I developed a method to estimate horizontal cover (within 0-1 m and 1-2 m from
the ground) using the profile 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 fiinction 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 profile 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 file. The file 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!
. .fi I
F
."1
.lwt.
.I ..

 

Figure D.2. Profile 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
classified 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 fine-
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 field-
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 field 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 defined 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
rectified 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 final 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 define the
feature space, and a 3 x 3 mean filter 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 first iteration, the plot attribute was categorical, describing the structural
stage of each plot as defined 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
classification 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, afier 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 insufficient 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 stratified 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 stratification 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 fire under email: dcleland@fs.fed.us
UP natural fire 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
identification 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