. u...

my»?!

bus
41.:

.7.

1:.

imv_w..1hw43r.m....um .r ..

 

. r

% 53Lw1

11...! .... .6! ..

our
.. \\\‘!

3...“ .

5f.

fi

(uni-9.?!”
- .2“ I

uirtat

rfifi‘mq.
D 1

\tln1xh1'1lll foul.
1‘: t...‘l\ll’\t‘n

2...]... till}?
I323). U
:rvvlvzz. lr’rl I
it)- 11). El
1'... «

2|“. ’3...
[xii-(39%.:
I!!.b§linh.l

t
{lily
bez‘ys.¢0.lsi'€(\t71~$1li
Ivn..2\(i..‘!|1l£ 5.: l.’ .
ill!“ '0!)

I} in.) 3 ‘30:!
.Ivt...naY.A-.00

 

 

“‘73" LIBRARY
4 ;» Michigan State
I’ 0 ‘3 University

 

 

 

This is to certify that the
thesis entitled

MULTIPLE GAP-, STAND- AND LANDSCAPE-SCALE
FACTORS AFFECT REGENERATION IN MANAGED
NORTHERN HARDWOOD FORESTS

presented by

Megan Shanahan Matonis

has been accepted towards fulfillment
of the requirements for the

MS. degree in Forestry

 

 

 

MSU is an Affirmative Action/Equal Opportunity Employer

PLACE IN RETURN Box to remove this checkout from your record.
'1 TO AVOID FINES return on or before date due.
' MAY BE RECALLED with earlier due date if requested.

 

DATE DUE DATE DUE DATE DUE

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

5/08 K:IProjIAcc&Pres/ClRC/Date0ue.indd

MULTIPLE GAP-, STAND- AND LANDSCAPE-SCALE FACTORS AFFECT
REGENERATION IN MANAGED NORTHERN HARDWOOD FORESTS

By

Megan Shanahan Matonis

A THESIS
Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of

MASTER OF SCIENCE
Forestry

2009

ABSTRACT

MULTIPLE GAP-, STAND- AND LANDSCAPE-SCALE FACTORS AFFECT
REGENERATION IN MANAGED NORTHERN HARDWOOD FORESTS

By
Megan Shanahan Matonis

The perpetuation of extensive and economically, socially and ecologically
important northern hardwood forests relies largely on natural regeneration below canopy
gaps created by selection harvesting. I explored factors affecting the abundance,
composition and survival of seedlings and saplings in harvest gaps across northern
hardwood stands in the Western Upper Peninsula of Michigan with both a natural
experiment and a mesic conifer planting experiment. Sapling abundance was highly
variable; nearly fifty percent of harvest gaps were devoid of sugar maple saplings (the
dominant overstory species), but abundance was great in other stands (up to 32,100
saplings / ha). Sugar maple sapling abundance was lower on southern stands with rich
soil nutrient conditions where estimated winter deer densities were high. Deer browse
pressure was high in the study area. Almost all planted hemlock seedlings were browsed,
demonstrating the futility of restoring browse preferred coniferous species where winter
deer densities are high. Seedling and sapling abundance were positively influenced by
canopy openness and negatively by cover of competing vegetation, but effects were
weaker than those of Habitat Type and deer density. There was only limited evidence of
seed source limitations in these stands. Results illustrate that natural and artificial
regeneration in managed northern hardwood forests are affected by multiple interacting
and spatially varying factors that need to be considered when planning timber harvests,

implementing restoration efforts and modeling regeneration dynamics.

Copyright by
MEGAN SHANAHAN MATONIS
2009

 

Dedicated to my loving and caring Grandmother, the shining queen of our family. Your

presence will be missed but never forgotten.

iv

ACKNOWLEDGEMENTS

Earning my Master’s degree has been an opportunity to utilize science to explore
an interesting ecological question with management implications, but more importantly,
it has been an opportunity to interact with many great people without whom my success
would not have been possible. I would like to acknowledge the people who inspired me
to challenge myself, who assisted me with my research, who reminded me to have fun,
who provided me with intellectual insight and who offered me comfort and support.

Thank you to my many field assistants who kept great attitudes on rainy days, on
hot days when bugs swarmed our heads, and even on a day when the tire of the truck fell
off while we were driving. I was accompanied in the field by Andrea Bianco, Amanda
Falk, Phillip Kurzeja, Alyssa Nugent, Ashlie Peterson, Erik Palm, Nick Reno, Grant
Slusher and the hardworldng crew of Tom Nolta, and I was assisted in the lab by Julia
Jones, Chad Babcock and James Bussa. Thank you so much for your dedication, sunny
personalities and encouragement.

I must extend a special thank you to Jim Cousineau from whom I rented a
comfortable lakeside cabin during summer field work. Jim, a 76 year old man with the
energy of a youth in his twenties, was a warm and welcoming host who provided me with
insights from his life in the woods.

William Brondkye, Linda Lindberg, Jim Ferris and Eric Thompson from the
Michigan Department of Natural Resources, Greg Lake from American Forest
Management and Charlie Becker from Plum Creek were instrumental in finding forest

' stands for my research. I enjoyed the opportunities to develop relationships with forest

land managers so I could learn how to tailor my research to benefit their efforts at
promoting sustainable forestry.

I would like to thank Dr. Kyung-Hwan Han, Dr. Donatien Kamden and Dr.
Laurent Matuana from Michigan State University for permitting me to use their
laboratories and equipment. My committee members, Dr. Andrew Finley, Dr. Rique
Campa, Dr. Richard Kobe, and my major adviser Dr. Mike Walters were engaged in my
research and provided me with useful insights into the process of scientific inquiry. Dr.
Walters was a great advisor who provided me with assistance when needed but also the
freedom to manage my own project. I appreciate the opportunities he gave me to present
my research to forest managers and researchers. Thank you to Dr. James Millington who
provided me with intellectual support, kept me company on multiple eight hour drives to
the UP and worked with me in the field to execute the mesic conifer planting project.

Thank you to Steven Pierce and Yun Xue from the Center for Statistical Training
and Consulting for providing me with assistance as I taught myself how to analyze
generalized linear mixed models with Bayesian inference.

I am very grateful to Dr. Jack Liu and the Center for System Integration and
Sustainability (CSIS) for the wonderful working environment and the financial support
they provided me with. Due to the financial support through the USDA CSREES grant
for the Western Upper Peninsula Economic Ecological Modeling Project, I never felt
constrained to limit the scope of my research. Through CSIS I was also fortunate enough
to interact with Dr. Bill Taylor who offered me the chance to serve as the assistant for the

Michigan Department of Natural Resources Statewide Council, who helped me think

vi

about my goals and how I can achieve them in my future, and who introduced me to
Mark Rey who helped me acquire an internship with the Forest Service.

Thank you to Michigan State University for the wonderful education and financial
assistance through an Academic Achievement Graduate Assistantship.

I want to thank the Binkley-Higgins family, my second family, for their love and
support throughout the years. Thank you to my scientific mentor and lifetime adviser, Dr.
Dan Binkley, who initially inspired me to pursue forestry and who taught me to think
critically and creatively. I also want to thank my Ferency housemates, and my other dear
friends who I met through the larger Michigan State University Student Housing
Cooperative network, for all the great memories. I want to thank Mike Heyse, my
boyfriend, for helping me to control my stress level. i

I am eternally grateful to my wonderful, loving and supportive family. Thank you
to my aunts, uncles, cousins and grandparents for making it a point to let me know how
loved I am. Thank you to my Canadian relatives who invited me to their homes during
the holidays and sporadically throughout the year for much needed vacations. Thank you
to my sister, my best fi‘iend, who can always be counted on for a funny story or a kind
word. Thank you also to my fantastic parents. My mother taught me to let passion and
compassion run my life, and my father taught me self-motivation and work ethic. I want
to especially thank my dearest Grandma who passed away last year. She was the most
wonderful grandma and helped to keep me going by calling or sending letters when I

needed comfort the most. I am so thankful to have known her.

vii

TABLE OF CONTENTS

LIST OF TABLES ............................................................................................................ ix
LIST OF FIGURES ............................................................... . .......................................... xv
INTRODUCTION ............................................................................................................ 1
CHAPTER 1

GAP-, STAND- AND LANDSCAPE-SCALE FACTORS AFFECT NORTHERN
HARDWOOD REGENERATION FOLLOWING SELECTION HARVESTING

Abstract .................................................................................................................. 4

1.1 Introduction ...................................................................................................... 5

1.2 Methods ............................................................................................................ 9

1.3 Results ............................................................................................................. 21

1.4 Discussion ....................................................................................................... 28

1.5 Conclusion ...................................................................................................... 39
CHAPTER 2

SNOW DEPTH, DECIDUOUS-LOWLAND CONIFER FOREST EDGES AND
DEER DENSITY AFFECT SEEDLIN G BROWSE DAMAGE

Abstract ......................... , ........................................................................................ 77
2.1 Introduction ..................................................................................................... 78
2.2 Methods ........................................................................................................... 81
2.3 Results ............................................................................................................. 88
2.4 Discussion ....................................................................................................... 91
2.5 Conclusion ...................................................................................................... 96
APPENDICES
Appendix A: Description of Study Stands .......................................................... 111
Appendix B: Additional Details on Methodology .............................................. 117
Appendix C: Sample R2WinBUGS code and WinBUGS models ..................... 132
Appendix D: Assessing Convergence ................................................................ 139
Appendix E: Model Parameters ......................................................................... 149
LITERATURE CITED ................................................................................................... 155

viii

LIST OF TABLES

Table 1.1. Northern hardwood Habitat Type descriptions (Burger and Kotar
2003) ................................................................................................................................. 43

Table 1.2. Model specification for spatial trend, full and null generalized linear
multilevel models (GLMM), linear multilevel models (LMM) and linear models

(LM). uij is the predicted value for gap i at stand j. H’=Shannon diversity

index .................................................................................................................................. 45

Table 1.3. Abundance of saplings 1-7 m tall per 154 m2 gap plot (per hectare of

gap area): Average (Ave), median (Med), standard deviation (Stdev), minimum

(Min) and maximum (Max) and percent of gap plots where present (% gaps)

(n=347 gaps) ..................................................................................................................... 46

Table 1.4. Sapling type by species and height. Gap colonizers germinated after
gap formation and advanced regeneration saplings were present before harvest ............. 47

Table 1.5. Summary of gap and stand-level independent variables: Average
(Ave), median (Med), standard deviation (Stdev), minimum (Min), maximum
(Max) and number of gap or stand-level observations (N obs) ........................................ 47

Table 1.6. Kendall tau rank correlations between gap and stand-level variables.
Comparisons between gap-level and stand-level variables use stand-averages for
gap-level variables to comply with independence assumptions ....................................... 48

Table 1.7. Differences in stand-level and stand average gap-level variables by

Habitat Type based on Kruskal-Wallis rank sum test. Inter-quartile ranges are

listed for those variables that significantly differ between Habitat Types.

Significant comparisons (p < 0.05) are based on Wilcoxon signed-rank test with
Bonferroni adjusted p-values ............................................................................................ 49

Table 1.8. Spatial patterns in sapling abundance 1-2 m tall, stocking levels and
Shannon Diversity Index (H’): mean parameter estimates and 95% credible

interval from posterior distributions. Positive values indicate increases in that

variable with Northing or Easting (UTM position rescaled from 0 to 1).

Parameter estimates are on the log or square root scale due to the use of a log-

link or transformations (where indicated). Null models include the overall

intercept and random stand-level intercepts for multilevel models and the overall
intercept for growth rate linear models. A decrease in DIC of 5 or more indicates
improved performance for the spatial model over the null model .................................... 50

ix

Table 1.9. Spatial patterns in gap- and stand-level variables: bootstrap parameter
estimate and 95% confidence intervals for stand-level variables or stand averages

of gap-level variables. Positive values indicate increases in that variable with

Northing or Easting (UTM position rescaled from 0 to 1). Gap- and stand-level

variables were not standardized for this analysis .............................................................. 54

Table 1.10. Deviance Information Criterion for full models (all gap- and stand—

level covariates) and null models (overall intercept and random stand-level

intercepts for multilevel models and the overall intercept for growth rate linear

models) and number of gaps and stands with observations. A decrease in DIC of

5 or more indicates improved performance for fiill model over the null model ............... 62

Table 2.1. Multilevel logistic regression models used to explore the effects of

gap- and stand-level variables on percentage (pg) of white pine or hemlock
seedlings browsed in gap i at stand j in 2008. The number of seedlings browsed

is yij, the number of seedlings relocated is nij, and the random stand-level
intercept is aj ..................................................................................................................... 99

Table 2.2. Summary of gap- and stand-level variables: Average (Ave), median
(Med), standard deviation (Stdev), minimum (Min), maximum (Max) and
number of gap or stand-level observations (N obs) ......................................................... 100

Table 2.3. Kendall’s tau rank correlation between gap- and stand-level variables.
Comparisons between gap-level and stand-level variables use stand-averages for
non-centered gap-level variables to comply with independence assumptions.

UTM Easting and Northing values were taken from the center of the deer pellet

transects ............................................................................................................................ 100

Table 2.4. Average height (stdev), number seedlings found (N .found), percent

found (of 2040 seedlings / species planted), browsed and severely browsed

(browse category 3 or 4 for hemlock only) for planted seedlings across all gaps

and stands in spring 2008 ................................................................................................. 102

Table 2.5. Number of seedlings found (N .found) and percent found, browsed,

severely browsed (hemlock only), desiccated and percentage of desiccated

seedlings showing evidence of browse for hemlock at the same 11 stands and for

pine at the same 19 stands visited in both spring 2008 and 2009. Hemlock

seedlings were only resurveyed in 2009 at stands where more than 50% of the

seedlings had not been severely browsed in 2008 ........................................................... 102

Table 2.6. Multilevel logistic regression models predicting the gap-level

percentage of a hemlock seedlings severely browsed or white pine seedlings

browsed in 2008: mean parameter estimates and 95% credible interval from

posterior distributions. Deer and DSW were grand mean centered and pellets

were group (i.e. stand) mean centered. Estimates are on the logistic scale, so a

change of one unit in the independent variables corresponds to a change of at

most [3/4 on the probability scale. A decrease in DIC of 5 or more indicates

improved performance for full model over the null model .............................................. 106

Table 2.7. Multilevel logistic regression models predicting the gap-level

percentage of a hemlock seedlings severely browsed or white pine seedlings

browsed in 2008 for a subset of thirteen stands where distance to lowland conifer

was determined: mean parameter estimates and 95% credible interval from

posterior distribution. DSW was grand mean centered and pellets and DLC were

group (i.e. stand) mean centered. Estimates are on the logistic scale, so a change

of one unit in the independent variables corresponds to a change of at most [3/4

on the probability scale. A decrease in DIC of 5 or more indicates improved

performance for full model over the null model .............................................................. 108

Table 2.8. Mean predicted probability (Pr.) (95% prediction quantiles) of severe

browse for hemlock seedlings and probability of browse for white pine seedlings

under different gap pellet count (stand centered) and deep snow week conditions.
Predictions are for an unmeasured stand based on one-to-one simulations using

the posterior distribution for Model 3 parameters ........................................................... 109

Table A.1. Description of 59 study stands including average snow depth

November 2007 — April 2008 (cm), county, ownership, Habitat Type, relative

densities (RD) of the three most abundant overstory tree species and average

stand diameter at breast height (DBH) of mature trees with dbh >20 cm, and year

of harvest entry (YOE) ..................................................................................................... 1 11

Table A2. Geological characteristics for each stand, including landform type,

soil type, upper layer soil texture and sugar maple site index at age 50 ($150) in
meters estimated for each soil type and soil series .......................................................... 114

Table 3.1. Species-specific parameter values for “b” in equation estimating dbh

from stump diameter (equation 2) and standard error of estimated dbh

(Demaerschalk and Omule 1978). “Species” refers to the species from this study

and “Equivalent species” refers to the species from Demaerschalk and Omule

1978 .................................................................................................................................. 118

Table 8.2. Parameter estimates ([3), standard errors (SE), AIC and adjusted R2

values from linear regression models to predict log transformed basal area (BAu)
of trees with dbh < 20 cm ................................................................................................ 121

xi

Table 3.3. Parameter estimates ([3), standard errors (SE), AIC and adjusted R2

values from linear regression models to predict quadratic mean diameter (QMDt)
of trees with dbh > 2.54 cm ......................................................................................

Table B.4. Sensitivity analysis: percent change in diameter grth (DG) with
given percent change in BAu or QMDt and sensitivity index (absolute percent
change in diameter grth divided by absolute percent change in BAu or QMDt)

Table 3.5. Number of section (N. sections) centers separated by different
distances. The average of the semivariance values from the points in the different
distance bins was used to develop the semivariogram (Figure B.3) ........................

Table D.l. Convergence in models of spatial trends in seedling and sapling
abundance and log transformed stocking-level: sample size (N obs.), burn-in
length (N .burnin), iterations after burn-in (N .iter), thinning factor (T.factor) and
maximum Raftery-Lewis diagnostic N (number of iterations required to have a
95% probability of estimating the 2.5% quantile with a tolerance of +/-0.5%) and

....... 122

........ 124

........ 130

I (Dependence Factor) from three parallel MCMC chains .............................................. 140

Table D2. Convergence in LMM of spatial trends in square root transformed
Shannon Diversity Index (H’): Sample size (N obs.), burn-in length (N .bumin),
iterations after burn-in (N .iter), thinning factor (T.factor) and maximum Raftery-
Lewis diagnostic N (number of iterations required to have a 95% probability of
estimating the 2.5% quantile with a tolerance of +/-0.5%) and I (Dependence
Factor) from three parallel MCMC chains ...............................................................

Table D.3. Convergence in LM of spatial trends in log transformed sapling
grth rate after harvest: Sample size (N obs.), burn-in length (N .bumin),
iterations after burn-in (N .iter), thinning factor (T.factor) and maximum Raftery-
Lewis diagnostic N (number of iterations required to have a 95% probability of
estimating the 2.5% quantile with a tolerance of +/-0.5%) and 1 (Dependence
Factor) from three parallel MCMC chains ...............................................................

Table D4. Convergence in GLMM of seedling abundance with gap- and stand-
level covariates: Sample size (N obs.), burn-in length (N .burnin), iterations after
burn-in (N .iter), thinning factor (T.factor) and maximum Raftery-Lewis
diagnostic N (number of iterations required to have a 95% probability of
estimating the 2.5% quantile with a tolerance of +/-0.5%) and I (Dependence
Factor) from three parallel MCMC chains ...............................................................

Table D5. Convergence in GLMM of sapling (1-2 m) abundance with gap-

and stand-level covariates: Sample size (N obs.), burn-in length (N .burnin),
iterations after burn-in (N .iter), thinning factor (T.factor) and maimum Raftery-
Lewis diagnostic N (number of iterations required to have a 95% probability of
estimating the 2.5% quantile with a tolerance of +/-0.5%) and 1 (Dependence
Factor) from three parallel MCMC chains ...............................................................

xii

........ 141

........ I41

........ 142

........ 143

Table D6 Convergence in LMM of log transformed sapling 1-7 m tall stocking

with gap- and stand-level covariates: Sample size (N obs.), burn-in length

(N .burnin), iterations after burn-in (N .iter), thinning factor (T.factor) and

maximum Raftery-Lewis diagnostic N (number of iterations required to have a

95% probability of estimating the 2.5% quantile with a tolerance of +/-0.5%) and

I (Dependence Factor) from three parallel MCMC chains .............................................. 144

Table D7 Convergence in LMM of square root transformed of Shannon

Diversity Index (H’) with gap- and stand-level covariates: Sample size (N obs.),

burn-in length (N .burnin), iterations after burn-in (N .iter), thinning factor

(T.factor) and Raftery-Lewis diagnostic N (number of iterations required to have

a 95% probability of estimating the 2.5% quantile with a tolerance of +/-0.5%)

and 1 (Dependence Factor) from three parallel MCMC chains ....................................... 145

Table D8. Convergence in LM of log transformed sapling growth rates with

gap- and stand-level covariates: Sample size (N obs.), burn-in length (N .burnin),
iterations after burn-in (N .iter), thinning factor (T.factor) and Rafiery-Lewis

diagnostic N (number of iterations required to have a 95% probability of

estimating the 2.5% quantile with a tolerance of +/-0.5%) and 1 (Dependence

Factor) from three parallel MCMC chains ....................................................................... 146

Table D9. Convergence in multilevel logistic models of gap-level percentage

of seedlings browsed with gap- and stand-level covariates: Sample size (N obs.),

burn-in length (N .burnin), iterations after burn-in (N .iter), thinning factor

(T.factor) and Raftery-Lewis diagnostic N (number of iterations required to have

a 95% probability of estimating the 2.5% quantile with a tolerance of +/-0.5%)

and 1 (Dependence Factor) from three parallel MCMC chains ....................................... 147

Table D.10. Convergence in multilevel logistic models of gap-level percentage

of seedlings browsed with gap- and stand-level covariates for stands where DLC

could be determined: Sample size (N obs.), burn-in length (N .bumin), iterations

after burn-in (N .iter), thinning factor (T.factor) and Raftery-Lewis diagnostic N

(number of iterations required to have a 95% probability of estimating the 2.5%

quantile with a tolerance of +/-0.5%) and 1 (Dependence Factor) from three

parallel MCMC chains ..................................................................................................... 148

Table E.1. Effects of gap and stand-level variables on seedling abundance <1 m

tall: mean parameter estimates and 95% credible interval from posterior

distributions. Positive values indicate increases in seedling abundance with a one
standard deviation increase in that covariate on the log scale. AOCa is the

reference Habitat Type. A decrease in DIC of 5 or more indicates improved

performance for full model over the null model .............................................................. 149

xiii

Table E.2. Effects of gap and stand-level variables on sapling abundance 1-2 m

tall: mean parameter estimates and 95% credible interval from posterior

distributions. Positive values indicate increases in sapling abundance with a one

standard deviation increase in that covariate on the log scale. AOCa is the

reference Habitat Type for ironwood and other species. ATD-Hp and AOCa were
combined as the reference Habitat Types to improve convergence of the sugar

maple model. A decrease in DIC of 5 or more indicates improved performance

for full model over the null model ................................................................................... 150

Table E.3. Effects of gap and stand-level variables on log transformed sapling

1-7m tall: mean parameter estimates and 95% credible interval from posterior
distributions. Only gaps with stocking-levels > 0 by species were used in these

analyses. Positive values indicate increases in stocking-levels on the log scale

with a one standard deviation increase in that covariate. AOCa is the reference

Habitat Type. A decrease in DIC of 5 or more indicates improved performance

for full model over the null model ................................................................................... 151

Table E.4. Effects of gap and stand-level variables on square root transformed

Shannon Diversity Index (H’) for seedlings, saplings 1-2 m and saplings 1-7 m

tall: mean parameter estimates and 95% credible interval from posterior

distributions. Only gaps with H’ > 0 by size class were used in these analyses.

Positive values indicate increases in H’ on the square root scale with a one

standard deviation increase in that covariate. AOCa is the reference Habitat Type.

A decrease in DIC of 5 or more indicates improved performance for full model

over the null model .......................................................................................................... 152

Table E.5. Effects of gap and stand—level variables on log transformed sapling

after harvest grth rates (m/yr): mean parameter estimates and 95% credible

interval from posterior distributions. Positive values indicate increases in growth

rate on the log scale with a one standard deviation increase in that covariate.

AOCa is the reference Habitat Type. A decrease in DIC of 5 or more indicates

improved performance for full model over the null model .............................................. 153

xiv

LIST OF FIGURES

Figure 1.1. Hypothesized factors affecting northern hardwood regeneration

throughout seedling and sapling ontogeny and the relationships between them.

This is not an all inclusive diagram, but highlights many of the factors the

literature and field observations emphasize as important ................................................. 41

Figure 1.2. Map of study area. MI=Michigan, WI= Wisconsin, ON: Ontario,
Canada ............................................................................................................................... 42

Figure 1.3. Relationship between maximum observed seedling and sapling

density in a 154 m2 gap plot among 347 study gaps versus height. This

relationship was used to develop stocking estimates. Loge(Density)= -2.316 X
Loge(Height) + 5.738 ........................................................................................................ 44

Figure 1.4. Site average sugar maple sapling (1-2 m tall) abundance per gap
plot across the study area. For reference, 10 saplings / gap plot is 650 saplings /
hectare and 227 saplings / gap plot is 14,746 saplings / hectare of gap area .................... 51

Figure 1.5. Site average ironwood sapling (1-2 m tall) abundance per gap plot
across the study area. For reference, 10 saplings / gap plot is 650 saplings /
hectare and 82 saplings / gap plot is 5,327 saplings / hectare of gap area ........................ 52

Figure 1.6. Site average other species sapling (1-2 m tall) abundance per gap
plot across the study area. For reference, 10 saplings / gap plot is 650 saplings /
hectare and 208 saplings / gap plot is 13,512 saplings / hectare of gap area .................... 53

Figure 1.7. Spatial distributions of Habitat Types by UTM Easting and

Northing (scaled to 0-1). Letters indicate significant differences between Habitat

Types based on Wilcoxon signed-rank test with Bonferroni adjusted p-value.

Violin plots show distribution of values and median (white dot), interquartile

range (thick vertical line) and range of values 1.5 * IQR (thin vertical line) ................... 55

Figure 1.8. Effects of gap and stand-level variables on seedling (<1 m tall)
abundance: mean parameter estimates, 95% credible interval (thick line), and
90% credible interval (thin line) from posterior distributions. Positive values
indicate increases in seedling abundance with a one standard deviation increase
in that covariate on the log scale. AOCa is the reference Habitat Type. SM=sugar

maple, IW=ironwood, Other=other species, SP2003= seedling production

potential (Zdiameter/distancez) for 2008, Random Int. Stdev= standard deviation

for the random stand-level intercept. Exact parameter values can be found in
Appendix B, Table E.l ..................................................................................................... 57

XV

Figure 1.9. Effects of gap and stand-level variables on small sapling (1-2 m tall)
abundance: mean parameter estimates, 95% credible interval (thick line), and
90% credible interval (thin line) from posterior distributions. Positive values
indicate increases in sapling abundance with a one standard deviation increase in
that covariate on the log scale. AOCa and ATD-Hp are reference Habitat Types
for sugar maple and AOCa is the reference Habitat Type for ironwood and other

species. SM=sugar maple, IW=ironwood, Other=other species, SPHmest=

seedling production potential (Ediameter/distancez) at time of harvest, Random
Int. Stdev= standard deviation for the random stand-level intercept. Exact
parameter values can be found in Appendix B, Table E.2 ................................................ 59

Figure 1.10. Effects of gap and stand-level variables on log transformed sapling
(1-7 m tall) stocking-level: mean parameter estimates, 95% credible interval
(thick line), and 90% credible interval (thin line) from posterior distributions.
Positive values indicate increases in sapling stocking on the log scale with a one
standard deviation increase in that covariate. AOCa is the reference Habitat Type.

SM=sugar maple, IW=ironwood, Other=other species, SPHawest= seedling

production potential (Zdiameter/distancez) at time of harvest, Random Int.

Stdev= standard deviation for the random stand-level intercept. Exact parameter
values can be found in Appendix E, Table E.3 ................................................................. 61

Figure 1.11. Observed sapling (1-2 m tall) abundance and mean posterior

predicted values and 95% credible intervals per gap plot from generalized linear

mixed models with gap- and stand-level covariates (full model) and the null

model with an overall intercept and stand-level random effect intercepts ....................... 64

Figure 1.12. Effects of gap and stand-level variables on square root transformed
Shannon Diversity Index (H’) for seedling (<1 m tall) (H’seedling), small

saplings (1-2 m tall) (H’Sapl) and all saplings (1-7 m tall) (H’SapTot): mean

parameter estimates, 95% credible interval (thick line), and 90% credible interval

(thin line) from posterior distributions. Positive values indicate increases in H’ on

the square root scale with a one standard deviation increase in that covariate.

AOCa is the reference Habitat Type. H’Overstory = H’ of overstory trees for

H’Seedling and H’ of overstory trees and stumps for H’Sapl and H’SapTot.

Random Int. Stdev= standard deviation for the random stand-level intercept.

Exact parameter values can be found in Appendix B, Table E.4 ...................................... 66

Figure 1.13. Effects of gap and stand-level variables on log transformed sapling
post-harvest grth rates (m/yr mean parameter estimates, 95% credible interval

(thick line), and 90% credible interval (thin line) from posterior distributions.

Positive values indicate increases in growth rate on the log scale with a one

standard deviation increase in that covariate. AOCa is the reference Habitat Type.
SM=sugar maple, IW=ironwood, Other=other species, Advanced Regen=

bivariate variable indicating if a variable is advanced regeneration. Exact

parameter values can be found in Appendix B, Table E.5 ................................................ 68

xvi

Figure 1.14. Sugar maple and ironwood sapling growth rate after harvest for
gap colonizers and advanced regeneration (ad regeneration) ........................................... 69

Figure 1.15. Predicting sugar maple abundance / gap plot (or / hectare of gap
area) for an observed stand: median and 95% percentile predictions based on one
to one simulations drawing from the posterior coefficient distributions and
random effects intercept for an observed ATM stand where a gap had 387
saplings in a gap plot (25,140 / ha of gap area). Simulations were run with
different covariate conditions (holding all else equal): all at observed covariate
values for this gap and stand, changing Habitat Type (HT) to AOCa/ATD-Hp,

increasing deer density from 4 to 16 deer/km2 (observed to observed plus one
standard deviation), increasing canopy openness (CO) from 5 to 12% and
increasing cover of competing vegetation (CV) from 18 to 41% ..................................... 70

Figure 1.16. Predicting sugar maple abundance / gap plot (or / hectare of gap
area) for an unmeasured stand: Median and 95% percentile predictions based on
one to one simulations drawing from the posterior coefficient distributions and

randomly drawing the stand-level intercept from ~N(u, 60,) where u is defined

by the stand-level covariates and 0a is the random stand-level intercept stand
deviation. Simulations were run with different covariate conditions (holding all
else equal) for a gap at an AOCa and ATD type stand: all at average observed

covariate values, increasing deer density from 14 to 26 deer/km2 (average

observed to average observed plus one standard deviation), increasing canopy

openness (CO) from 13 to 20% and increasing cover of competing vegetation

(CV) from 42 to 65%. A break was added to the y-axis so changes in the median
prediction could be seen .................................................................................................... 71

Figure 1.17. Percentage of gaps harvested 8-15 years ago where a species was

removed (stump species) and there was at least one sapling of the same in the

gap plot. Species codes: Con=conifer (eastern hemlock, northern white cedar,

balsam fir and white spruce), YB= yellow birch, BW= American basswood,

QA= quaking aspen, SM= sugar maple, RM= red maple, BC= black cherry,
IW=ironwood. Numbers below species codes indicate the number of gaps where

that species was removed by harvesting ........................................................................... 72

Figure 1.18. Percentage of gaps sugar maple (SM), ironwood (IW) or other

species are the dominant stump species (highest relative basal area of stumps

removed from gap) and percentage of gaps these species are the dominant

sapling species (highest relative stocking level of saplings 1-7 m tall). Percentage

of gaps with no saplings is also presented. N= 211 gaps harvested 8-15 years ago ......... 73

xvii

Figure 1.19. Average estimated deer density (deer / km?) from 1996-2000 created

from universal kriging (Appendix B, Section B.5) of 113 observations collected by

the Michigan Department of Natural Resources and deer density estimates from

winter 2007-2008. Estimates for winter 2007-2008 are positively correlated with

the MDNR 1996-2000 estimates (Kendall’s t=0.27, p-value=0.002) .............................. 74

Figure 1.20. Differences in stand average seedling abundance / 2 m2, small
sapling (1-2 m) abundance / gap plot and sapling stocking levels / gap plot by

Habitat Type. Gap plots were 154 m . Letters indicate significant differences

between Habitat Types for sugar maple based on Wilcoxon signed-rank test with
Bonferroni adjusted p-value. Kruskal-Wallis rank sum test indicated that sapling
abundance and sapling stocking were not significantly different between Habitat

Types for ironwood or other species, but it was significant for seedling of both

species. However, no pairs of Habitat Types were significantly different based

on the Wilcoxon test. Violin plots show distribution of values and median (white

dot), interquartile range (thick vertical line) and range of values 1.5 * IQR (thin

vertical line) ...................................................................................................................... 76

Figure 2.1. Location of stands and estimated winter 2007-2008 deer density

overlaid the number of weeks from November 2007 to April 2008 predicted to

have snow depth > 40 cm. MI=Michigan, WI=Wisconsin, USA. ON=Ontario,

Canada ............................................................................................................................... 98

Figure 2.2. Correlation between stand-mean centered gap distance to lowland

conifer (DLC) and stand-mean centered gap pellet group counts (PGC). Lower

image shows the correlation between the non-centered estimates. Non-parametric
Kendall’s tau rank correlations (t (p-value)) are presented ............................................. 101

Figure 2.3. Percentage of seedlings browsed / stand by species based on

browsed vs not-browsed categorization. Distributions were all significantly

different from each other based on Wilcoxon signed-rank test with bonferroni

adjusted p-value. Violin plots show distribution of values and median (white

dot), interquartile range (thick vertical line) and range of values 1.5 * IQR (thin

vertical line) ..................................................................................................................... 103

Figure 2.4. Observed percentages of hemlock seedlings severely browsed /

stand in 2008 (and 2009 in the distance to lowland conifer (DLC) plot only)

versus estimated winter deer density, stand average pellet group count (PGC),

deep snow weeks (DSW), and stand average gap-DLC. Non-parametric Kendall’s

tau rank correlations (1' (p-value)) are presented. 1: was not calculated for 2009

observed percentage browsed versus DLC due to low sample size (n=4).

**p-value < 0.05, * p-value <0.10 ................................................................................... 104

xviii

Figure 2.5. Observed percentages of pine seedlings browsed / stand in 2008

(and 2009 in the distance to lowland conifer (DLC) plot only) versus estimated

winter deer density, stand average pellet group count (PGC), deep snow weeks

(DSW), and stand average gap-DLC. Non-parametric Kendall’s tau rank

correlations (r (p-value)) are presented. **p-value < 0.05, * p-value < 0.10 .................. 105

Figure 2.6. Observed percentage of hemlock seedlings severely browsed and

white pine seedlings browsed / gap in 2008 versus mean posterior predicted

percentages per gap and 95% credible intervals from logistic mixed models.

Model 3 includes deep snow weeks and stand-centered gap-level pellet group

counts as covariates and model 1 (null model) includes only an overall intercept

and stand-level random effect intercept. Observed percentage on the x-axis has

been jittered for visualization (random jitter is equivalent for Model 3 and Model

1 plots within species) ...................................................................................................... 107

Figure B. 1. Layout of pellet count transects .................................................................. 117

Figure B.2. Effect of the prior on the standard deviation (6) of the random stand-
level intercept in the sugar maple sapling abundance generalized linear multilevel
model: mean, 90% credible interval (thick black line), and 95% credible interval
(thin black line) of the posterior distribution. IG (inverse gamma) prior was on

the precision (l/o ) and half normal and uniform density priors were on o. The
half normal distribution was ~N(0, 10,000), bounded above 0 ....................................... 126

Figure 3.3. Effect of the prior on the standard deviation (6) of the random stand-
level intercept in the gap-level percentage browsed hemlock seedlings logistic
multilevel model: mean, 90% credible interval (thick black line), and 95%
credible interval (thin black line) of the posterior distribution. IG (inverse

gamma) prior was on the precision (1/02) and half normal and uniform density
priors were on o. The half normal distribution was ~N(0, 10,000), bounded
above 0 ............................................................................................................................. 126

Figure B.4. Average estimated deer density (deer/kmz) from DNR pellet count

surveys 1996-2000 within Iron, Dickinson, Marquette and Menominee counties

versus UTM Northing for the Public Land Survey section. Log transformed deer

density = 23.8518 - 4.0e-06 >< UTM Northing based on GLS-estimation. AIC is

107.7 compared to 85.7 for the null model ...................................................................... 128

Figure B.5. Semivariogram of the residuals from the GLS-estimated
relationship between average estimated deer density (deer/kmz) and UTM

Northing. A lag distance of 2,000 m and an exponential model form were used ............ 129

Figure B.6. Observed deer density estimates versus values predicted from cross
validation .......................................................................................................................... 13 l

xix

INTRODUCTION

Northern hardwood forest (maple-beech-birch type) spans over 53.7 million acres
throughout the Great Lakes, North Central and Northeastern United States (Smith et a1
2004) and is the predominant forest type in Michigan (Potter-Witter and Lacksen 1993).
These forests have high multi-use value because they produce millions of cubic meters of
timber each year for pulp, lumber and veneer, they support summer and winter recreation
activities and provide habitat for wildlife (Filip 1973; OMNR 1998).

Northern hardwood forests have undergone, and are undergoing changes in
structure and composition due in part to the interrelated disturbances of harvesting and
deer herbivory (Kraft et al 2004). During the early commercial logging era in the Lake
States from the mid 1800’s to early 1900’s, forests were clear cut and white pine and
hemlock were virtually eliminated from many forests (Stearns 1997). As early as the
1920’s, industrial forest companies began experimenting with selection harvesting as a
sustained yield approach (Stearns 1997). With the publication of long-term studies
demonstrating the benefits of this approach (Eyre and Zillgitt 1953) and the development
of marking guides (Arbogast 195 7), selection harvesting has become common practice in
northern hardwood stands. It is advocated as a silvicultural practice that meets ecological
and economic objectives by mimicking gap phase dynamics and capitalizing on the
ability of these forests to naturally regenerate (N yland 1998). However, when seedlings
and saplings fail to grow into larger size classes, it threatens the ability of these forests to
sustain long-term timber production (N yland 1998) and could destabilize the economy in
areas such as the Western Upper Peninsula of Michigan where the timber industry is one

of the main economic drivers (Racevskis and Lupi 2006).

Forest harvesting, along with increased agricultural activity, the extirpation of
predators, and legislation to regulate hunting has created landscape conditions favorable
for white-tailed deer (Rooney 2001; Cote et al 2004). Although a natural herbivore in
northern hardwood forests, white-tailed deer populations have reached unprecedented
densities and are exerting changes on understory composition in hardwood forests
(Rooney 2001). Browsing by white-tailed deer (Odocoileus virginianus) can decrease
seedling and sapling densities (Stoeckeler et a1 1957; Horsley et al 2003; Rooney and
Waller 2003; Eschtruth and Battles 2008) by reducing height increments and survivorship
(Liang and Seagle 2002; Horsley et al 2003; Eschtruth and Battles 2008). Although
intense browse by white-tailed deer (Odocoz'leus viriginanus) is often cited as the major
cause of regeneration failure after harvest treatments in hardwood forests (Stoeckeler et a1
1957; Shafer et al 1961; Marquis and Brenneman 1981; Cook 2008), seedling and sapling
layers are affected by a suite of abiotic and biotic factors throughout their ontogeny, so
studying deer density in isolation of other factors produces limited results (Didier and
Porter 2003).

Seed source, light availability, competing vegetation, deer density and soil
nutrient and moisture conditions were identified from the literature and field observations
as key factors limiting regeneration in northern hardwood forests. Chapter 1 describes a
study that assesses the relative effect of these factors on seedling and sapling abundance
and composition below harvest gaps in 59 northern hardwood stands across the Western
Upper Peninsula of Michigan. This natural experiment took advantage of the north to
south gradient in winter deer density across this landscape and the variation in gap

sizes/light availability within stands created by harvest practices. Although causal

relationships cannot be established from natural experiments, this approach can reveal
patterns and associations between variables across their natural range of variation and can
be implemented in a more spatially intensive and extensive manner.

Restoring the white pine and eastern hemlock components that were once part of
northern hardwood forests are conservation goals across the Great Lakes region in order
to increase biodiversity and alter white-tailed deer winter and summer range availability
(Herman et a1 2004). Chapter 2 describes a study that assesses the efficacy of efforts to
plant mesic conifer seedlings in an area of the Western Upper Peninsula with a gradient
in snow depth, variation in winter deer densities and variation in landscape contexts
related to winter range suitability for deer (i.e. distance to lowland conifer stands used by
deer for winter shelter). Results from this study demonstrate a relationship between
seedling browse damage and deer density estimated by the fecal pellet technique and can
assist with directing restoration efforts.

The two studies in this thesis were performed in managed forests and focus on
seedling and sapling populations (natural or artificially planted) in harvest gaps and
explore impacts of deer herbivory, as well other factors on seedling and sapling
abundance and composition (Chapter 1) and survival (Chapter 2). Results demonstrate
that relatively simple mechanisms, such as differences in shade tolerance, cannot be used
to explain and predict forest dynamics because various interacting and spatially varying
factors act as regeneration filters. Insight into the relative contributions of regeneration
limiting factors, provided by my multivariate assessment at multiple stands across a
broad geographic area, can assist with predicting future forest structure and composition

and with managing northern hardwood forests for the many services they provided.

Chapter 1: Gap-, stand- and landscape-scale factors affect northern hardwood
regeneration following selection harvesting

Abstract

The abundance and diversity of saplings regenerating below canopy gaps are
affected by a suite of abiotic and biotic variables including seed source, light-availability,
soil moisture and nutrient regime, competition from understory vegetation and herbivory.
Multivariate approaches that explore gap-, stand- and landscape-scale factors across
broad geographic areas are necessary to assess the relative effects of these factors and to
understand spatial patterns in regeneration success. To this end, I measured seedling and
sapling layers and quantified the previously listed factors within 347 canopy gaps created
by selection harvesting across 59 northern hardwood stands in the Western Upper
Peninsula where winter deer density increases with latitude. A multilevel modeling
approach was taken to account for the hierarchical structure of the data (i.e. gaps nested
within stands), and parameter estimates from these models were used to simulate sapling
abundances under various gap- and stand-level conditions. Sapling abundance of sugar
maple (Acer saccharum Marsh.) was higher on stands with lower estimated deer
densities, higher snow depths, and less nutrient rich Habitat Types underlain with coarser
textured substrates. Habitat Type (a proxy for stand soil moisture and nutrient regime)
had the greatest effect on sugar maple sapling abundance, followed by deer density, light
availability and competing vegetation. Northern hardwood stands have the potential to
undergo compositional changes if deer densities are not reduced because
unpalatable/browse-tolerant species may replace highly browsed species. Results
demonstrate that various factors can interrupt the regeneration of even shade tolerant

species following canopy gap formation.

1.1 Introduction

In forests where stand clearing disturbances occur infrequently, sapling
regeneration relies on pulses of resource provided by canopy gaps (Runkle 1985;
Yamamoto 2000). Canopy gaps form naturally due to mature tree mortality from old age,
windfall, lightning, disease and insects (Krasny and Whitrnore 1992), but tree removal
from selection harvesting has become the primary mode by which these gaps are created
in managed northern hardwood forests. The aim of selection harvesting is to sustain the
production of high quality timber by mimicking natural gap dynamics to foster
regeneration and by maintaining an uneven-aged structure (Arbogast 1957; Guldin 1996;
OMNR 1998). This harvest practice is recommended for northern hardwood forests
(Arbogast 1957) because of evidence that it promotes the regeneration of economically
valuable shade tolerant species such as sugar maple (Acer saccharum Marsh.) (Tubbs
1968) and can facilitate regeneration of shade-intolerant and intermediate species if group
selection is used to create larger gaps (Leak 1999).

Abundant natural regeneration of desirable species is essential to the success of
the selection system (Nyland 1998). Unimodal size-density stand structures considered
unsustainable under selection harvesting (Smith et a1 1997) may develop when there are
insufficient seedling and sapling populations to recruit into larger size classes and replace
the pole and saw-log sized trees removed by harvesting (N yland 1998). Failure of
hardwood species such as sugar maple, yellow birch (Betula alleghaniensis Britton),
black cherry (Prunus serotina Ehrh.) and northern red oak (Quercus rubra L.) to
germinate or grow into taller sapling layers has been observed in many managed

hardwood forests for over half a century in the Midwestern and Northeastern United

States (Stoeckeler et a1 1957; Marquis and Brenneman 1981; Jenkins 1997; Aldrich et al
2005; Belden and Pallardy 2009), including in the Western Upper Peninsula of Michigan
(Graham 1954; Miller 2004; Donovan 2005). Sugar maple is the dominant species in
northern hardwood forests, and loss of this species from the understory could
dramatically alter future stand composition and structure (Jenkins 1997). Reduced
regeneration of commercially important species following selection harvesting causes
concern for forest managers who seek to ensure long-term timber production and
maintain certification under the Sustainable Forestry Initiative (Donovan 2005).
Regeneration failures could jeopardize the sustainability of forest productivity and thus
destabilize the economy in areas like the Western Upper Peninsula of Michigan where
the timber industry is one of the main economic drivers (Racevskis and Lupi 2006). In
addition, loss of seedling and sapling layers in hardwood forests have ecological
consequences because they serve as forage for species such as white-tailed deer
(Odocoileus viriginanus Zimmermann) (Blouch 1984) and provide vertical stand
structure necessary for some bird species (Webb et a1 1977; DeGraaf et a1 1998).

The abundance and species composition of seedling and sapling layers in harvest
gaps are affected by many abiotic and biotic factors that reduce or enhance seed
availability, germination, grth and survival. Seed production, light availability, soil
nutrients and moisture regime, competition from non-tree understory vegetation, and deer
herbivory have been identified as key factors influencing regeneration in northern
hardwood forests based on previous studies and field observations (Figure 1.1). This list
is not all-inclusive, but highlights those factors hypothesized to have the greatest effects

on sugar maple and other common northern hardwood species present in the understory.

Densities of seedlings are initially constrained by the production of seeds.
Generally positive relationships between tree size, seed production, and seedling densities
(Marquis 1975; Tubbs 1977; Hughes and F ahey 1988; Ribbens et a1 1994; Garret and
Graber 1995; Fei and Steiner 2008) imply that management decisions regarding the
retention of large mature seed trees can impact future regeneration. Once seedlings
germinate, availability of soil moisture, nutrients and light affect their growth (Canham
and Marks 1985; Canham 1988; Walters and Reich 1997; Finzi and Canham 2000;
Bigelow and Canham 2002; Webster and Lorimer 2002; Schreeg et al 2005; Kobe 2006),
survival (Caspersen and Kobe 2001; Schreeg et a1 2005) and species richness (Tubbs
1977; Runkle 1982; Runkle 1984; Burger and Kotar 2003; Schumann et a] 2003).
Increasing the size of harvest gaps increases understory light availability (Gray et al
2002; Canham et al 1990), seedling and sapling growth rates (Canham and Marks 1985;
Kobe 2006) and cover of non-tree understory vegetation (Collins et al 1985; Schumann et
al 2003; Shields and Webster 2007). Shrubs, ferns and graminoids can reduce seedling
and sapling densities, survival and height growth (Horsley and Marquis 1983; De Steven
1991; Romagosa and Robison 2003; Fei and Steiner 2008) by altering understory
microhabitat conditions (e. g. temperature and light availability) and competing for
nutrients (George and Bazzaz 1999) and water (Randall 2007).

Intense browse by white-tailed deer is often cited as the major cause of
regeneration failure after harvesting treatments in hardwood forests (Graham 1954;
Stoeckeler et al 1957; Shafer et al 1961; Marquis and Brenneman 1981; Cook 2008).
Selective browsing and differences in species tolerance to browsing (Augustine and

McNaughton 1998) affect the abundance and species richness of understory tree layers

(Stoeckeler et al 1957; Liang and Seagle 2002; Horsley et al 2003; Rooney and Waller
2003; Eschtruth and Battles 2008), occasionally causing the near elimination of highly
preferred browse species, such as is the case with eastern hemlock (T suga canadensis
(L.) Carriére) and northern white-cedar (T huja occidentalis L.) in the Great Lakes region
(Rooney et al 2000; 2002). By removing palatable species, herbivory can permit the
dominance of a few unpalatable, browse resistant / tolerant species (Horsley et al 2003;
Cote et a1 2004) such as ironwood (Ostrya virginiana (Mill.) K. Koch) or beech (Fagus
grandifolia Ehrh.) (Sage et al 2003; Miller 2004).

This study evaluates the relative impacts of deer density and other gap- and stand-
level factors on seedling and sapling layers. It takes advantage of a decreasing gradient in
winter deer density with increasing latitude in the Western Upper Peninsula of Michigan,
caused in part by the southern seasonal migration of deer to avoid deep lake-effect snow
(V erme 1973; Doepker et a1 1994; VanDeelen et al 1998). The impacts of white-tailed
deer on tree regeneration have been the focus of many studies (see reviews- Weisberg
and Bugmann 2003; Rooney and Waller 2003; Cote et a1 2004), however, few have
assessed the relative impact of deer density in conjunction with other gap- and stand-
level factors across a large geographic area (Rooney et al 2000; 2002). None have done
so for broadly distributed and economically valuable northern hardwood species. Given
the suite of factors that affect seedlings and saplings throughout their ontogeny,
multivariate assessments at multiple stands across broad geographic areas are critical for
determining the relative contributions of different factors and to reduce the risk of

assigning significance to factors that are correlated with unmeasured variables.

This study explores the following hypotheses: H1) as winter deer density
increases across a landscape-scale gradient, sapling abundance of heavily browsed
species such as sugar maple (Swift 1949; Stoeckeler et al 1957) decreases, H2) seedling
and sapling abundances, diversity and height growth rates in gaps created by selection
harvesting are affected by gap-level seed source availability (+), light availability (+) and
cover of competing vegetation (-) and by stand-level soil moisture and nutrient conditions
(+) and winter deer density (-) and H3) current understory regeneration may not
perpetuate current overstory species composition, due in part to deer browse of palatable
species.

1.2 Methods
1.2.] Study area

Northern hardwood stands included in this study (n=59 stands) were harvested
two to fifteen years ago and located across 475,000 contiguous hectares of Dickinson,
Iron, Marquette and Menominee counties in the Western Upper Peninsula of Michigan
(Figure 1.2; Appendix A, Table A. 1). This study area encompasses that for the ongoing
Western Upper Peninsula Economic-Ecological Modeling Project that integrates ecology
and economics to better understand interactions among deer, tree harvesting and
regeneration on managed forest landscapes (Laurent et al 2005; LeBouton et a1 2005;
Racevskis and Lupi 2006; Millington et al 2009 in review). This area was selected for the
domination of forest cover, the relative absence of agriculture and urban or suburban
developments, and the natural variation in deer densities and snow depth (Doepker et a1
1994; Shi et a1 2006). All but five stands were dominated by sugar maple in the overstory

(43-100% basal area of trees with dbh >20 cm), and had various components of

American basswood (T ilia americana L.), red maple (Acer rubrum L.), eastern hemlock
(Tsuga canadensis (L.) Carriére), white ash (Fraxinus americana L.), black cherry
(Prunus serotina Ehrh.), balsam fir (Abies balsamea (L.) Mill.), yellow birch (Betula
alleghaniensis Britton), and other minor coniferous and deciduous species.

Stands ranged in elevation from 283 to 551 m with underlying soils from thirty
distinct soil types (Appendix A, Table A2). Soil drainage varies from very poor to
excessive. Annual snow fall varies from 161 cm in the south to 434 cm in the north
(National Climatic Data Center 2009) due to lake-effect snow influenced by Lake
Superior (Jerome 2006). ’

1.2.2 Study stand harvest dates

Harvest dates for all stands were estimated with radial increment analysis of stem
cross sections taken from saplings growing in harvest gaps. The harvest date was
estimated by the year of release from growth suppression (Webster and Lorimer 2005).
Harvest dates estimated from release were strongly correlated with the harvest records
available for a subset of seventeen stands (Pearson’s r= 0.87), with a difference of d: 2
years (average = -0.4 years).

1.2.3 Field methods
1.2.3.a Stand-level variables: Soil moisture and nutrient regime (Habitat Type), deer
density, browse pressure and snow depth

Habitat Type, a proxy for soil moisture and nutrient regime, was determined from
diagnostic assemblages of understory vegetation at each stand (Burger and Kotar 2003).
Northern hardwood stands in the study area fall into five different Habitat Types (listed

from lowest to highest soil nutrient availability): TMC, ATM, ATD, ATD-Hp and AOCa

10

(Table 1.1). Stand classification into one of these types is not based on sapling abundance
or diversity.

Winter deer density was estimated with fecal pellet surveys at each stand in spring
2008 (Hill 2001). Fecal pellet surveys are a fairly simple method of assessing relative
local deer density (Hill 2001) with reasonable accuracy (Neff 1968; Forsyth et a1 2007,
but see Fuller 1991). However, it has been criticized because defecation rates are
correlated with forage intake, forage moisture content and percentage of young in the
herd and estimates can be biased by sampling design and observer error (Neff 1968).
Despite these limitations, no other method can be as readily utilized to determine stand-
scale variation in deer density across a large geographic area. Pellet groups were counted
along ten 50 m by 4 m transects oriented in an hour glass shape (Appendix B, Section
B.1). All transects were double counted with different observers to improve accuracy
(Jenkins and Manly 2008). The average pellet group counts from all ten transects were
used to calculate deer density following the methods of Hill 2001 (Appendix B, Section
B.1). Visible pellets were those above leaf litter deposited in November 2007, so the
estimate of deer density corresponds to deer presence in the late fall, winter and early
spring months. 1 term the estimate “winter deer density” although deer depositing fecal
pellets in northern hardwood stands likely did so while migrating between summer and
winter ranges (i.e. in November and April).

Stand-level sugar maple browse index (SMBI), a surrogate for browse pressure,
was estimated by counting the number of browsed and unbrowsed terminal twigs 30—200
cm above the ground on 5- 24 sugar maple saplings 1 to 2.5 m tall (Frelich and Lorimer

1985; Rooney et a1 2000). Each sapling was assigned to one of five browse categories

11

based on the percentage of total twigs browsed (0= no browsing, 1 = 1-25%, 2 = 26-50%,
3= 51-75%, and 4 = 76-100%). Average SMBI was calculated for each of the 25 stands
where there were more than 5 sugar maple saplings for which browse index could be
determined.

Average daily snow depth (cm) for November 2007 to April 2008 was determined
from the National Snow and Ice Data Center’s Snow Data Assimilation System
(SNODAS) snow mass and energy balance model (Barrett 2003; NOHRSC 2004).

1.2. 3. b Gap-level variables: gap size, light availability, seed production potential,
competing vegetation, seedling and sapling abundance and growth rates

Gap-level variables were quantified within six gaps at each stand, with the
exception of five stands where fewer gaps were measured because fewer existed or very
high sapling densities made data collection time-intensive. Gaps were defined as canopy
openings created by the removal of one or more trees, with the shortest inter-bole
diameter of at least ten meters. Gap sizes were estimated using the method of Runkle
1981, i.e. calculated as an ellipse using inter-bole distances between gap edge trees. Gaps
were sampled from three size strata (shortest diameter <12.5 m, shortest diameter 12.5-15
m shortest diameter >15 m) when present.

Light availability was estimated from hemispherical photographs (Canham et a1
1990; Domke et a1 2007) taken at the gap center as well as three meters from the center in
each cardinal direction to capture variability in canopy openness. Canopy openness (CO-
a proxy for light availability) was estimated from photos converted to binary images
using an automatic threshold value (SideLook v. 1.1.01 , Nobis 2005) with Gap Light

Analyzer software (Frazer et a1 1999). The average CO from the usable photos for each

12

gap was used in analyses. CO was not estimated for ten gaps where three or more of the
five photographs were unusable due to sun glare or image obstruction.

Seed production potential (SP) was estimated for each gap by summing the
quotient of diameter at breast height (dbh-1.4 m) and the squared distance for each
mature tree by species. This method assumes that a tree’s contribution to SP increases
with diameter and decreases with distance (Ribbens et a1 1994). Trees and stumps within
20 meters of the gap plot and with dbh > 5 cm for ironwood (a smaller species at
maturity) and > 20 cm for other species were included in this estimate. If the species of a
stump could not be determined in the field, it was determined in the laboratory fi'om
wood samples (Marx 2005). Stump basal diameter was converted to dbh using
relationships developed by Demaerschalk and Omule 1978 (Appendix B, Section B.2).
Stumps were not measured at one stand due to the level of decay which made species

identification and diameter measurements difficult.

Seed production potential estimates were developed for 2008 (szoog) based on

the current dbh of living trees and for the time of harvest (SPHaWest) based on stump dbh

and “grown back” dbh of living trees. Mature tree dbh was “grown back” to dbh at the
time of harvest using radial growth equations from the USDA Forest Vegetation

Simulator Lake States Variant (Bush and Brand 1993) (Appendix B, Section B.3).
SPHamst increases rapidly with the inclusion of large stumps near the gap center.

Percent ground cover of fern, grass and Carex spp, Rubus spp, and other shrubs

was visually estimated in five 1 m2 quadrats located within each gap. Numbers of tree

seedlings (<1 m tall) were tallied by species and height (to the nearest 0.25 m) in two of

13

 

the 1 m2 competing vegetation quadrats located near the gap center. Tree saplings 1-7 m

tall were tallied by species and height (to the nearest 0.25 m) in one 154 m2 gap-centered

circular plot (7 m radius). This plot encompassed the entire gap area for smaller gaps, and
a large portion of the area for larger gaps. In gaps with extreme sapling densities (> 200
saplings / gap plot), only one forth to three fourths of the gap plot was sampled to
improve efficiency, and estimates were scaled up to the entire gap plot.

If saplings were present, stem cross section were taken at 5 cm and breast height
from one sugar maple and one ironwood sapling within three different height strata (1-2
m, 2-4 m and >4 m). These species were selected because they were more consistently
representated across stands than other species. Sapling age at 5 cm and 1.4 m height was
determined (number of rings counted on the cross section plus 0.5 years to account for
the current year’s growth) in order to estimate height grth rates. Saplings taller than
1.4 m at the time of harvest were excluded from analysis because post-harvest grth
rates could not be estimated for them. The average annual height grth rate was
determined from 5 cm to the saplings height in 2008 for gap colonizers (saplings with age
at 5 cm < time since harvest) and from 1.4 m to the sapling’s height for advanced
regeneration. A different method was used for advanced regeneration because some
height growth below 1.4 m occurred during the period of suppression before harvest.
Height growth estimates were adjusted using the method of Carrnean 1972 because cross
sections did not correspond with terminal bud scars.
1.2.4 Data analysis: independent variables

I sought to explore spatial patterns in seedling and sapling abundances, stocking-

levels, growth rates, and gap- and stand-level predictor variables, and to develop models

14

to explain the impacts of these variables on regeneration. I focus analysis on sugar maple
and ironwood because they were the most abundant and consistently represented species
in the understory. All other species are pooled together as “other” because infrequent
observations within species prohibited separate analyses by species. The “other” category
is heavily dominated by red maple, white ash and black cherry with these species
representing 55%, 26% and 6% of seedlings in the “other category”, 30%, 49% and 12%
of saplings 1-2 m tall and 30%, 44% and 15% of saplings 1-7 m tall. Four size categories
of regeneration are analyzed: 1) seedlings (< l m tall), 2) small saplings that are within
the range of deer browse (1-2 m tall), 3) all saplings (1-7 m tall) and 4) a measure of
sapling stocking that takes into account densities and heights of saplings.

Comparing sapling abundance alone does not estimate the potential of the sapling
pool in a gap to recruit into the overstory because self thinning causes fewer larger
saplings to persist over time (Tubbs 1968; Runkle 1984; Kittredge and Ashton 1995; Ray
et a1 1999). To estimate the stocking level of each gap, I weighted each sapling by the
inverse of the maximum observed gap-level density for saplings of that height across all
gaps (Figure 1.3) and summed the weights of all saplings in a gap by species. If a gap
was at the maximum density for a sapling size class, it had a stocking level of 1 (i.e.
100%), however, many gaps had multiple cohorts of saplings so the stocking level was >
1. One gap with an outlying stocking level (>14) was removed from analyses for sugar
maple and another from the analyses of other species (stocking level > 7). Analyses for
sapling stocking by species were performed for gaps where the stocking level for that

species was > 0.

15

Shannon Diversity Index (H’) was calculated for seedlings, small saplings (1-2 m)
and all saplings (1 -7 m). Shannon Diversity Index analyses were only performed for gaps
where H’ was > 0 for that size class.

To improve normality, sapling stocking and sapling post-harvest growth rates
were log transformed and H’ was square root transformed.

1.2.5 Data analysis: statistical analysis

The data were hierarchically structured, with gaps nested within stands, so a
hierarchical multilevel modeling approach was taken for analysis. Ignoring the
relatedness between gaps within stands would increase the risk of type I errors because
the effective sample size is really less than the total number of gaps (Goldstein 1995;
Congdon 2003). Multilevel models incorporate information about the clustering to
produce estimates of standard errors that account for non-independence (Goldstein 1995).
I used Bayesian inference with Markov chain Monte Carlo (MCMC) sampling to
estimate parameter values and obtain credible intervals that average over the uncertainty
in fixed- and random-effect parameters (Zhao et a1 2006; Bolker et al 2009).

Multilevel models were not used for sapling growth rate analyses because
preliminary analyses using the 1me4 package (Bates and Maechler 2009) in R 2.9.1 (R
Development Core Team 2009) demonstrated that inclusion of a stand-level random
intercept did not improve model fit (AIC = 1.68 vs 52.44 for sugar maple and 12.74 vs
62.7] for ironwood linear versus multilevel linear models with all gap- and stand-level
covariates). This is likely because variation in growth rates for saplings within the same
stand was comparable to variation between all saplings (average within-stand stdev=

0.062 cm vs overall stdev=0.077 cm for sugar maple and 0.086 cm vs 0.099 for

16

ironwood) and there were many stands (4 of 27 for sugar maple and 10 of 30 for
ironwood) where growth rate after harvest could only be determined for one sapling.
Non-linear effects in the relationship between sapling height and height growth were not
included because it appeared linear over the range of observed heights.

1.2.5.a Assessing spatial patterns in regeneration (H I ) and modeling the effect of gap-
and stand-level factors (H2)

Spatial trend models (UTM Northing and Easting as covariates), full models (gap-
and stand-level covariates) and null models (overall intercept and random stand-level
intercepts only) were developed using generalized linear multilevel models for seedling
and sapling abundance, linear multilevel models for log transformed sapling stocking
level and square root transformed H’ and linear models for log transformed sapling
growth rates. Seedling and sapling abundances were modeled as negative binomial
processes arising from the distribution ~NegBin(p, r) where p = r / (r + u) and r (the
dispersion parameter) = 112/ (0’2 — u). Table 1.2 shows the formulation for the spatial
trend, full and null models and lists the covariates used for each analysis.

UTM Easting and Northing were rescaled from 0 to 1 for the spatial trend models,
and gap- and stand-level parameters for the full models were standardized to improve
convergence. Percent ground cover of fern, grass and Carex spp, Rubus spp, and other
shrubs were added together to create a composite metric, percent cover of competing

vegetation, to reduce the number of parameters and because these variables were auto-

correlated. SP2003 was used in seedling abundance models and H’ of overstory trees in

seedling H’ models because it was unlikely that trees represented by stumps were the

seed source for small seedlings except at the most recently harvested stands.

17

1.2. 5.b Model fitting, priors and convergence

Models, including the linear growth rate models, were run with WinBUGS v.1 .4.3
(Spiegelhalter et al 2003) (see code Appendix C, Section C.l- C.3) through R v.2.9.1
with the package R2WinBUGS (Sturtz et a1 2005) (see code Appendix C, Section C.5).
Three parallel chains with dispersed randomly selected starting values were run for
70,000 iterations with a burn-in of 10,000 and a thinning rate of 5 for seedling and
sapling abundance and H’ models, and for 20,000 iterations with a burn-in of 5,000 and a
thinning rate of l for stocking level and growth rate models. Some models had to be run
for longer iterations and with a higher thinning rate to improve convergence (Appendix
D, Tables D.1-D.8). AOCa and ATD-Hp Habitat Types were combined into one group
for sugar maple sapling abundance models because low variation within and between
these groups affected identifiability (the ability to estimate all parameters in the model).

Deviance Information Criterion (DIC), a method of assessing model performance
in terms of fit and complexity based on the posterior mean deviance (-2 x
log(likelihood)) and the effective number of parameters (Spiegelhalter et a1 2002), was
compared for the spatial trend models, full models and null models to determine if
predictor variables improved model performance. A decrease in DIC by 5 or more
indicates support for better model performance.

In multilevel models, a group’s (i.e. stand’s) random intercept is adjusted by the
group’s deviation from the predicted value, so predictions can closely fit observed data.
This is especially the case when few observations are made within groups and variation is
low within groups. The null model predictions could perfectly fit all gap-level

observations at a stand where variation between gaps was zero. Improvements in model

18

fit over the null model demonstrate the ability of gap— and stand-level covariates to
explain variation between and within stands given the deviation of each stand from the
mean. If DIC is lower for the null model than a model with covariates, it indicates that
relative to other unmeasured differences that might exist between stands, differences
cannot be attributable to covariates.

A uniform prior (range 0-50) was used for the standard deviation of the random
stand-level intercept for all multilevel models. The random effects standard deviation was
insensitive to prior specification (Appendix B; Section B.4). An inverse gamma prior was

used for the sapling-level precision in the sapling growth rate linear models. A

noninformative prior distribution (Normal~(u=0, 02:10,000)) was used for gap-level and

stand-level coefficients (Clark 2007; Bolker et al 2009).

Convergence was diagnosed using the R package coda (Plummer et al 2009) with
the Gelman-Rubin diagnostic and the Raferty-Lewis diagnostic. All models showed
strong evidence of convergence (Appendex D, Tables D.1-D.8).

1.2. 5. c Comparing the effect of significant gap- and stand-level variables

Simulations were used to predict the effect of a one standard deviation increase in
gap- and stand-level variables on sugar maple small sapling abundance to determine the
relative effect of significant variables, and to aid in the biological interpretation of model
parameter values. Predictions were developed using one-to-one simulation draws from
the posterior parameter distributions. Simulations for small sugar maple sapling
abundance were run for 1) the gap with the greatest observed sugar maple abundance
holding all covariates at their observed values and 2) a gap at an unmeasured stand where

the covariate values were at the mean level observed across the study area. Method 1

19

incorporates less uncertainty into the predictions because the posterior distribution for the
random stand-level effect for the observed stand can be used (Gelman and Hill 2007).
1. 2. 5. d Assessing spatial patterns in gap- and stand-level variables

Linear regression with nonparametric pairs bootstrapping was used to explore
spatial trends in stand-level variables and stand averages of gap-level variables.
Nonparametric Kruskal-Wallis rank sum test was used to determine if Easting and
Northing were significantly different between Habitat Types, and compared between
Habitat Types with nonparametric Wilcoxon signed-rank test with Bonferroni adjusted p-
values.
1.2.5.e Potential species composition changes (H3 )

Future forest composition can change if saplings regenerating in harvest gaps are
of different species than those that were removed from the gap (Runkle 1981). The
potential for species change was assessed with two methods: 1) determining the
percentage of gaps where a species was removed by harvesting (i.e. present as a stump in
the gap) and also present in the sapling layer and 2) comparing the percentage of gaps
where a species was the dominant species removed (highest relative stump basal area) to
the percentage of gaps where that species dominated the sapling layer (highest relative
sapling stocking level). These analyses only used gaps from stands that had been
harvested at least eight years before data collection because younger gap may require
more time to accumulate all possible sapling species. The first analysis only included
species represented by stumps in at least 5 gaps and the second analysis combines species

other than sugar maple and ironwood into the “other” category.

20

1.3 Results
1.3.] Characterizing regeneration abundance and composition

A total of 18 species were identified as saplings (1-7m tall) within gaps across the
study area, but regeneration was heavily dominated by sugar maple and ironwood (Table
1.3). Only six species occurred in more than 5 percent of the gap plots (ironwood, sugar
maple, black cherry, white ash, red maple and balsam fir), and only ironwood and sugar
maple occurred in more than 50% of all gap plots. Twelve percent of gap plots were
devoid of saplings of any species and 35% had less than five stems of all species
combined (< 325 saplings / ha of gap area). Sugar maple was absent from 48% of gap
plots and ironwood from 43%.

Species richness and diversity were low, with a median number of two tree
species in the seedling layer, two in the sapling layer and three in the overstory layer. The
maximum observed number of species was five in the seedling layer, seven in the sapling
layer and seven in the overstory layer. Shannon Diversity Index (H’) was 0.31 for
seedlings with a maximum of 1.46, 0.35 for saplings with a maximum of 1.48 and 0.58
for overstory trees with a maximum of 1.60. The maximum possible H’ is 1.61 with five
species and 1.95 with seven species.

Percentage of sugar maple and ironwood saplings originating as gap colonizers
decreased with increasing height (Table 1.4). The tallest gap colonizer was 2.75 m for
sugar maple (7 years old at a stand harvested 7 years ago) and 2.5 m for ironwood (9
years old at a stand harvested 9 years ago). The average age (number of growth rings at 5
cm height) of advanced regeneration at the time of gap formation for sugar maple was 16

years (stdev 14.4 years) and 9 years (stdev 8.0 years) for ironwood.

21

Average post-harvest growth rates for saplings 1-2 m tall were not significantly
different between ironwood (0.19 m / yr, stdev 0.06) and sugar maple (0.18 m / yr, 0.07),
but growth rate for ironwood (0.29 m / yr, 0.10) was significantly faster than that for
sugar maple saplings 2-4 m tall (0.24 m / yr, 0.07). The growth rate for ironwood 4-7 m
tall (0.42 m / yr , 0.07) was higher than the growth rate of the two tall sugar maple
saplings for which post-harvest growth could be calculated (0.24 m / yr).

1.3.2 Characterization of gap and stand-level variables
Harvest gaps were created by the removal of three trees (dbh > 20 cm) on

average, but as few as one and as many as thirteen, and ranged in extended gap size from

82 to 913 m2 (Table 1.5). Positive correlations were observed between canopy openness

and gap size and cover of competing vegetation (Table 1.6). Canopy openness, cover of

competing vegetation and estimated deer density were negatively correlated with time

since harvest. Estimated winter deer density was on average 13.9 deer/ kmz, ranging from

0.6 to 61.6 deer/ kmz, and was negatively correlated with snow depth. Stand average

sugar maple browse index was positively correlated with time since harvest but not with
deer density estimates.
Habitat Types were well represented (n> 10 stands) except for TMC (n= 2). Snow

depth was deeper on ATD and ATM types compared to ATD-Hp (Table 1.7). Stand

average percent cover of ferns, grass and Carex spp and SPHamst for other species were

the only other variables differing between Habitat Types. If TMC was excluded from the

Kruskal-Wallis rank sum test due to low sample size, percent cover of Rubus spp and

sugar maple SPHarvest were significantly different between the remaining four Habitat

22

Types (x2=8.27, dfi3, p-value=0.04 and x2=9.69, dfi—3, p—value=0.02 respectively).
However, no pairwise comparisons between Habitat Types were significant for cover of

Rubus spp. SPHmest for sugar maple was significantly higher on ATD type stands than

on ATD-Hp, even though SP2003 was not. This may be due to the larger number of sugar

maple stumps harvested from the area around gaps at ATD stands compared to ATD-Hp
(9.6 stumps vs 6.5 stumps, Wilcoxon rank sum test W=131.5, p=0.045).
1.3.3 Spatial patterns in regeneration and gap- and stand-level variables (H1)

Spatial trends in abundance were evident for all species groups in either seedling
or sapling size classes as well as for seedlings and saplings H’ (Table 1.8). Sugar maple
showed the most consistent and strongest spatial trends, increasing with UTM Northing
for seedling and sapling abundance, sapling stocking levels and grth rates. Sugar
maple small sapling abundance was dramatically different between 163 gap plots across
28 stands located south of latitude 46° 6’43” N and 184 gap plots across 31 stands north
of this latitude (Figure 1.4). Small sugar maple saplings were present in only 4% of gap
plots in southern stands, with a total abundance of eight saplings in the 163 gaps, but
were present in 71% of gap plots in the northern stands with abundances as great as 387
saplings in a single gap plot (25,140 saplings / ha of gap area). Consistently low
abundance of ironwood (Figure 1.5) and other species (Figure 1.6) were not observed
across stands in the southern region. Ironwood and other species small saplings were
present in 54% and 39% of southern gap plots with abundances as high as 127 and 242 /
gap plot (8,250 and 15,721 / hectare of gap area), respectively.

Although spatial patterns were evident in many models, very few showed signs of

improved performance (lower DIC) over the null model (random stand-level intercept

23

only). Due to few gap observations per stand, the stand-level intercept is able to explain
much of the variation in the observations.

Snow depth increased with Northing and marginally decreased with Easting,
whereas SMBI and estimated winter deer density decreased with Northing and Easting
(Table 1.9). Stand average percent cover of competing vegetation decreased with Easting.

Spatial trends were also observed in competing vegetation and seed production potential

for sugar maple and ironwood. Sugar maple SPHamst increased with Northing, even

though Sonog did not, because more sugar maple stumps were removed from northern

stands in the most recent harvest. Easting was significantly different between Habitat

Types (Kruskal-Willis x2: 18.07, df=4, p-value=0.001) as was Northing (Kruskal-Willis

x2: 33.63, df=4, p-value<0.001), with ATD stands generally located farther east than

AOCa, and ATD and ATM located farther north than AOCa and ATD-Hp (Figure 1.7).
1.3.4 Gap- and stand- level variables affecting seedling and sapling abundance, sapling
stocking, sapling growth rates and regeneration diversity (H2)

All measured gap- and stand-level factors significantly corresponded with at least
one species or species category in the seedling abundance (Figure 1.8), sapling
abundance (Figure 1.9) or sapling stocking levels (Figure 1.10) models.

At the gap-level, cover of competing vegetation had a consistent negative
association with seedling and sapling abundance and stocking. Canopy openness was
positively associated with small sapling abundance and stocking levels of sugar maple
and other species. Seed production potential was only positively associated with

ironwood seedling abundance.

24

Stand-level Habitat Type was consistently important in the full multilevel models,
especially for sugar maple where higher seedling and sapling abundances and stocking
levels were found in gap plots on ATD, ATM and TMC type stands compared to AOCa
and ATD-Hp. There was only one sugar maple sapling 1-2 m tall across all 85 gaps
measured on ATD-Hp stands. Differences between Habitat Types were not observed for
ironwood and other species small saplings, however specific species within the other
species group showed variation in sapling abundances between Habitat Types. No white
ash saplings (1-7 m tall) were found in gap plots on ATM type stands and only one
sapling was found across gap plots on ATD type stands. Stand-level time since harvest
was negatively associated with ironwood seedling abundance and estimated deer density
was negatively associated with small sugar maple saplings abundance.

The incorporation of gap- and stand-level variables improved model performance
(reduced DIC compared to the null model) for all sugar maple models, for ironwood
seedling and sapling stocking-level models (marginally), and for the other species sapling
abundance model (Table 1.10). Comparison of the observed versus mean posterior
predicted values for the sapling abundance models (Figure 1.11) demonstrates that the
stand-level intercept explained much of the variation, not the covariates. The stand-level
intercept provided the model with the ability to predict gaps with zero seedlings or
saplings per gap plot, but incorporation of gap- and stand-level effects improved its
ability to predict gaps with higher abundances for sugar maple and other species.

Diversity of seedlings, small saplings and all saplings was mostly insensitive to
measured gap- and stand-level variables (Figure 1.12). The Shannon Diversity Index of

overstory trees increased seeding diversity, and diversity appeared lower on ATD

25

compared to AOCa type stands across size classes. The incorporation of gap- and stand-
level covariates only improved model performance for the seedling diversity model
(Table 1.10).

Sapling post-harvest growth rates increased with height and decreased with time
since harvest for both species, but were insensitive to the other gap- and stand-level
variables except for marginal increases in ironwood growth rate with canopy openness
(Figure 1.13). Sapling growth rates post- harvest for sugar maple and ironwood were
lower for advanced regeneration saplings compared to gap colonizers at a given height
(Figure 1.14). Model performance was greatly improved with inclusion of sapling-, gap-
and stand-level variables (Table 1.10).

1.3.5 Predicting the effect of significant gap- and stand-level variables on sugar maple
abundance

Generalized linear mixed model parameter estimates were used to predict sapling
abundance for the gap with the highest observed sugar maple small sapling abundance
(3 87 saplings in a gap plot (25,140 / ha of gap area)) and monitor changes in abundance
when significant covariates were increased by one standard deviation. With covariates
held at their observed values, the median predicted abundance was 146 / gap plot (9,484 /
ha of gap area) (95% credible interval for the prediction: 11-731 /gap plot; 715 — 47,487 /
ha). The most dramatic change in predicted sapling abundance was elicited by changing
the Habitat Type classification from ATM to AOCa/ATD-Hp; predicted sapling

abundance declined to 0 / gap plot (0-2 / gap plot; 0-130 / ha) (Figure 1.15). Increasing
winter deer density from the observed value by 1 standard deviation (4 to 16 deer/kmz)

decreased the prediction to 40/ gap plot (2,598 / ha) (2-315 / gap plot; 130-22,801 / ha),

26

increasing canopy openness from 5% to 12% increased the prediction to 193 / gap plot
(12,538 / ha) (14-972 / gap plot; 909-63,142 / ha), and increasing competing vegetation
from 18% to 41% decreased the prediction to 124 / gap plot (8,055 / ha) (6-451 / gap
plot; 390 — 29,298 / ha).

Uncertainty in parameter estimates, and the importance of the random stand-level
intercept for constraining predictions, reduced the model’s ability to predict sapling
abundance for gaps at an unmeasured stand (Figure 1.16). Median predictions of sugar
maple abundance within a gap at an unmeasured stand were all below 10 / gap plot (650 /
ha). Due to similarities between the response of abundance to ATD, ATM and TMC
Habitat Types, only results from a simulated gap at an ATD type stand are presented for
comparison with predictions of abundance at an AOCa/ATD-Hp type stand. Predictions
were consistently higher for ATD compared to AOCa/ATD-Hp, and more variable for
ATD. The median predicted sapling abundance for gaps on AOCa/ATD-Hp type stands
remained at 0 regardless of changing covariate values.

1.3.6 Potential species composition changes

The percentage of gaps where a species was harvested and at least one sapling of
the same species was present varied from 7% for conifer species (n=14 gaps) to 100% for
ironwood (n=10) (Figure 1.17). About fifty percent of gaps with sugar maple (n=191),
white ash (n=6) or red maple stumps (n=21) had saplings of that species present.
Although yellow birch and basswood were removed from 20 and 24 gaps respectively,
there were only saplings of these species in about 10% of these gaps. ‘

Sugar maple was the dominant species (highest relative basal area) of stumps

removed from 74% of gaps, but only the dominant sapling species (highest relative

27

stocking level) in 37% of 211 gaps harvested 8-15 years ago (Figure 1.18). Ironwood
rarely had the highest relative stump basal area, but was the dominant sapling species in
32% of gaps. Other species were the dominant stump species and the dominant sapling
species in about 20% of gaps. Eleven percent of gaps had no regenerating saplings.
1.4 Discussion

Exploring multiple factors across multiple scales is fundamental for elucidating
the complexity of ungulate-vegetation interactions (Rooney and Waller 2003; Weisberg
and Bugmann 2003) and the causes of regeneration failure (Didier and Porter 2003;
Romagosa and Robison 2003; Sage et a1 2003). This study responded to this need by
characterizing abiotic and biotic variables operating at the landscape-, stand- and gap-
scales (Figure 1.1), and their effects on regeneration dynamics in the Western Upper
Peninsula of Michigan.

High variation was observed in the abundance of seedlings and saplings in gap
plots, with some plots full of saplings but most containing very few. Gap plots were

equivalent to the size of the growing space created by the removal of several mature
trees, and many sampled gaps were over the minimum size (100-400 m2) required for the

regeneration of shade-intolerant species (Webster and Lorimer 2005). The absence of any
sapling species from over ten percent of gap plots, the absence of shade-tolerant sugar
maple from fifty percent of gap plots, and the near complete absence of shade-intolerant
species from all gap plots implies regeneration failure in many stands.

Sugar maple sapling abundance increased with latitude where snow depth was
greater, winter deer densities were lower (H1), and less nutrient rich Habitat Types were

more abundant. All measured gap- and stand-level variables were associated with

28

changes in seedling and/or sapling abundance in selection harvest gaps (H2). 1 found
evidence that shifts in species composition are possible in some northern hardwood
stands (H3), potentially due to species sensitivities to gap- and stand-level variables.
There was less evidence of gap- and stand-level variables affecting seedling and sapling
diversity and sapling growth rates.

1.4.1. Inverse landscape-scale gradients in biotic and abiotic factors and sugar maple
sapling abundance (H1)

Snow depth increases with latitude in this region due to lake-effect snow. White-
tailed deer avoid areas with deep snow (Tierson et al 1985; Poole and Mowat 2005) and
migrate south to winter range where snow depth is lower (Verme 1973; VanDeelen et al
1998). Decreasing winter deer density with UTM Northing and Easting (Figure 1.19) has
been documented across the study area for at least half a century (Doepker et a1 1994).

Snowfall patterns also help account for the distribution of mesic forest
communities in the Great Lakes region, potentially because snow cover impacts soil
moisture, nutrient availability and fire history (Henne et a1 2007). The landscape-scale
gradient in winter deer density and Habitat Types, as influenced by snow depth, may
account for the spatial pattern of increasing sugar maple seedling and sapling abundance
with Northing and Easting. Sugar maple seed production potential at the time of harvest
also increased with Northing, but sugar maple still dominated the overstories of stands in
the southern portion of the study area. The sapling spatial trend does not reflect a spatial
trend in the range of sugar maple.

Deer density was not correlated with spatial patternsiin sugar maple sapling

abundance in northern New York State (Didier and Porter 2003), but in the Western

29

Upper Peninsula of Michigan the spatial pattern in winter deer density is inverse that for
sugar maple sapling abundance. Sugar maple saplings are heavily utilized as forage
(Swift 1949; Stoeckeler et al 195 7) and decline in abundance outside deer exclosures
(Stoeckeler et a1 1957; Miller 2004). Deer browse in the southern portion of the Western
Upper Peninsula may affect the perpetuation of sugar maple following selection
harvesting.

Confounding the north to south gradient in deer density is a gradient in Habitat
Types. Sugar maple seedling and sapling abundance and sapling stocking levels differ
between Habitat Types (Figure 1.20), being more abundant on ATD, ATM and TMC
(less-nutrient rich Habitat Types). Abundances of other species and ironwood also differ
between Habitat Types, but not to the extent of sugar maple. In the Western Upper
Peninsula of Michigan, variation in snow depth between Habitat Types, and the
associated impacts of snow on soil moisture-nutrient conditions, may account for some of
the abundance of sugar maple saplings observed on ATD, ATM and TMC stands with
coarser-textured substrates (Henne et a1 2007), but it cannot account for the paucity of
saplings on AOCa and ATD-Hp stands with finer-textured substrates. Other evidence
also points against the strength of the Habitat Type soil moisture and nutrient regime
effect on sugar maple sapling abundance: 1) differences in sapling growth rates between
Habitat Types, that could indicate differences in soil moisture and nutrient availability
(Walters and Reich 1997), were not observed, 2) the Habitat Types with lowest sugar
maple abundance are predicted to be most nutrient-rich and 3) differences in sugar maple

seedling abundance between Habitat Types were not as strong as differences in sapling

30

abundance, potentially because deer browse is interrupting the transition from seedling to
sapling.

Quantifying the availability of soil resources previously found to affect sugar
maple seedlings (i.e soil moisture, pH, nitrogen, concentrations of exchangeable
aluminum and calcium) (Walters and Reich 1997; Kobe et al 2002; Juice et al 2006)
between and within Habitat Types in this area could help disentangle the effects of
Habitat Types from that of deer. At the same time, habitat-suitability of forest stands for
deer has been found to vary among Habitat Types (Felix et a1 2004), so interactions
between deer and Habitat Type might innately interweave their impacts on regeneration.
1.4.2. Multiple gap- and stand-level factors inhibit or enhance seedling and sapling
abundance and stocking levels, but few seem to affect diversity and growth rates (H2)

Seedling abundance appears more sensitive to gap- and stand-level variables than
sapling abundance or stocking-levels. Seedlings have experienced less variation in gap
and stand-level conditions since germination compared to the older taller saplings,
potentially linking their abundance more closely to the condition of variables measured in
2008. Results from the small sapling abundance and stocking-level models were similar
because small saplings were also included in the stocking-level model, but they were
weighted less than taller saplings. Stocking-level models also did not include gaps with
zero saplings. Results of the stocking-level analysis, therefore, indicate variables that
affect the recruitment potential of a gap (i.e. Habitat Type, competing vegetation and
canopy openness) given there is some regeneration. The sapling abundance models
demonstrate impacts to smaller saplings within the range of deer browse, including

impacts contributing to instances when no small saplings are present. These instances

31

represent regeneration failure in the small size class. Many of these saplings germinated
after harvest and are directly impacted by silvicultural decisions.

Evidence was found for impacts of abiotic and biotic gap- and stand-scale factors
on seedling and sapling abundance. The relative effect of Habitat Type on predicted sugar
maple sapling abundance was the greatest, followed by deer density. Habitat Type and
deer density are spatially confounded, reducing the model’s ability to represent the
relative effect of one variable in isolation of the other. Estimated deer density did not
differ significantly between Habitat Types, but uncertainty in pellet count estimates could
contribute to this.

Gap-level canopy openness and cover of competing vegetation had weaker
impacts on sugar maple sapling abundance than stand-level factors, but effects were in
the direction predicted. Previous work corroborates my findings that increasing light
increases sapling abundance (Schumann et al 2003; Dietze and Clark in prep; Runkle
1982), and increasing cover of non-tree vegetation negatively impacts seedling and
sapling layers (Horsley and Marquis 1983; De Steven 1991; George and Bazzaz 1999;
Romagosa and Robison 2003; Fei and Steiner 2008). The negative association of cover of
competing vegetation with seedling and sapling abundance could also be contributed to
by the negative effect taller saplings have on non-tree vegetation (Collins et a1 1985).
However, controlled field experiments have found direct negative impacts of competing
vegetation cover on seedling emergence and seedling and sapling survival and growth
rates (Yawney and Carl Jr. 1970; Horsley and Marquis 1983; George and Bazzaz 1999;

Randall 2007).

32

Seed production did not appear to limit the abundance of understory regeneration,
except for ironwood seedlings. This was surprising given that seed density, and
occasionally seedling density, in forest soils show increases with the abundance, size and
proxirrrity of mature seed trees in the overstory (Marquis 1975; Ribbens et a1 1994; Garret
and Graber 1995; Shibata and Nakashizuka 1995; Clark et a1 1998). High inter-annual
variation in sugar maple seed crop production, seed viability and seed predation (Graber
and Leak 1992; Houle 1992; Garret and Graber 1995) may uncouple overstory tree size
and distances from gap-level seedling and sapling abundances. Disconnect may also be
observed between seed fall and seedling abundance due to different spatial and temporal
patterns in mortality factors (e. g. predation, pathogen attack and competition) affecting
seeds and seedlings (Houle 1992).

Seedling and sapling diversity and growth rate were not strongly associated with
gap- and stand-level factors in these northern hardwood stands. Understory tree layer
species diversity did not increase with canopy openness (based on the fill] H’ models).
Some studies report similar observation (Shields et a1 2007; Dietze and Clark in prep) but
others report diversity or richness responses to light availability (Runkle 1982; Webster
and Lorimer 2002; Schumann et a1 2003). The lack of a diversity response to light in this
study may be partly attributable to the exclusion of gaps with only 1 species (and
therefore H’=0) from analyses to comply with normality assumptions. The correlation
between canopy openness in isolation of other gap- and stand-level variables (and
ignoring autocorrelation between gaps sampled from the same stand) is positively related
to small sapling H’ (Kendall’s r=10.9, p-value=0.01) and all sapling H’ (r=14.6, p-

value<0.001) for gaps where at least one sapling was present (i.e. including those with

33

H’=0). Canopy openness remained non-significant for seedling H’ (F-0.05, p-value=
0.3). Understory seedling diversity appears more limited by overstory diversity than light
availability. However, overstory H’ did not translate into sapling H’, potentially because
some species germinate but do not survive and grow into larger size classes (Metzger and
Tubbs 1971).

Canopy openness did not affect sugar maple sapling post-harvest growth rates,
and only marginally affected ironwood growth rates. Light estimates from 2008 do not
approximate light availability at the time of gap formation due to canopy coalescence
(Domke et a1 2007), which could weaken light vs. growth rate relationships. Extended
gap size, which would not change with time since harvest, was not correlated with sugar
maple growth (r=-0.04, p-value=0.56), but it was positively correlated with ironwood
growth (r=0.18, p-value=0.01). Other studies indicate that sapling height growth
responses to light saturate beyond relatively low light levels (Canham 1988; Coates 2000;
Webster and Lorimer 2002). This causes significant differences to be observed between
height growth below a closed canopy and below canopy gaps, but not between gaps of
different sizes (Canham 1988; Webster and Lorimer 2002).

Sapling growth rates only responded strongly to time since harvest out of all gap-
and stand-level variables. Others have found sugar maple height growth to slow about 23
years after gap formation (McClure et a1 2000) but I find evidence of this effect
beginning earlier because all sampled gaps were less than 15 years old.

Gap colonizer saplings had faster growth rates than advanced regeneration of the
same height, consistent with previous findings (McClure et al 2000; Webster and Lorimer

2005). Physiological differences between gap colonizer and advanced regeneration,

34

mediated by years of suppression in advanced regeneration saplings, could lead to
differences in post-harvest growth rates. However, the differences observed in this study
may also be contributed to by differences in the growth rate estimation technique. Post-
harvest growth rates obtained using the advanced regeneration method (average growth
rate above 1.4 m tall) were slower than those using the method for gap colonizers (height
divided by age) by 0.02 to 0.08 m / year (average 0.05 m / yr) for fifteen of the nineteen
gap colonizer saplings tall enough to obtain reasonable cross sections at breast height.
Growth rates for the remaining four species were faster using the advanced regeneration
method by 0.01 to 0.11 m / year (average 0.04).
1.4.3. Potential for species composition changes due to replacement by other species
following gap creation

Shade intolerant species like yellow birch, deer browse-preferred species like
basswood, and conifers species are not regenerating in gaps where these species were
removed by harvesting. Failure of shade intolerant species to regenerate in stands
managed with selection silviculture is a common concern in northern hardwood
management (Metzger and Tubbs 1971; Shields et a1 2007). Loss of conifers (i.e. eastern
hemlock, northern white cedar, balsam fir and white spruce) from hardwood forests has
implications for forest diversity because conifers provide habitat and food for many
songbird species (Kendeigh 1945; Patmos 1995). Others have found negative effects of
deer density on sapling H’ (Horsley et al 2003; Cote et a1 2004; Miller 2004), but I found
no evidence of this in the full small sapling H’ models. However, when considered in
isolation of other variables (and including gaps with only one species), small sapling H’

is negatively correlated with estimated deer density (Kendall’s 1:=-0.17, p-value<0.001).

35

Shade tolerant sugar maple failed to regenerate in many gaps created by the
harvest removal of overstory sugar maple trees, and could be replaced by other species
especially in the southern portion of the study area. Ironwood was the most likely
successors in many gaps. Although sugar maple trees could potentially overtop ironwood
saplings after many years due to the smaller stature of mature ironwood, sugar maple
saplings (1-7 m tall) were absent from 88% of gaps where sugar maple was removed and
ironwood saplings dominated the gap plot. Ironwood is a small species at maturity, which
could account for the infrequency with which it was observed as the dominant stump
species. However, it was only removed from 10 gaps so the frequency with which it was
the most dominant sapling species is greater than would be expected if species
composition were remaining static. Ironwood saplings were unresponsive to gap- and
stand-level variables, indicating that this species has the potential to regenerate across all
Habitat Types, regardless of seed source, light availability, competing vegetation cover
and deer densities. The potential for ironwood to increase in dominance in areas with
high winter deer densities is corroborated by evidence from deer exclosure studies (Miller
2004). Replacement of commercially valuable sugar maple by the non-valuable
ironwood, and the general loss of species diversity in northern hardwood stands, are
concerning to forest managers (Miller 2004).

1.4.4. Management implications

Observations of sugar maple sapling (1-7 m tall) absence from nearly fifty percent
of gap plots challenges the notion that “securing some sort of commercially important
natural regeneration is usually a simple matter in most northern hardwood stands” (Tubbs

1977), especially for stands located in the southern portion of the Western Upper

36

Peninsula. Although early work on selection harvesting in the region supports this
statement (Eyre and Zillgitt 1953), increasing deer density across the area since the
1970’s (Doepker et a1 1994) and the potential unsuitability of selection harvesting for
securing regeneration on different Habitat Types may affect the ability to apply this
technique ubiquitously across northern hardwood stands.

The Michigan Department of Natural Resources is using the Burger and Kotar
2003 Habitat Type system to classify stands and direct management practices. I found
significant differences in regeneration abundance between Habitat Types, potentially
lending support to stratifying forest management by Habitat Type. Based on the current
northern hardwood management recommendations, ATD-Hp type stands are treated
equivalent to ATD, but my research suggests that these types should be treated separately
because sugar maple regeneration success between these Habitat Types was vastly
dissimilar. Specific recommendations may need to be developed for managing AOCa and
ATD-Hp type stands in areas with high deer density. Shelterwood cutting has been
suggested as a successful alternative to selection harvesting in areas experiencing
herbivory-mediate regeneration failure (Sage et a1 2003; Marquis and Brenneman 1981).
The development of Habitat Type-specific recommendations could facilitate ecosystem
management by coordinating appropriate silvicultural techniques with considerations of
deer habitat-potential (Felix et al 2004).

Creating larger harvest gaps may produce moderate increases in seedling and
sapling abundance, but across the range of gap sizes observed, I find limited evidence
that larger gaps increase sapling diversity. Gap size and light availability were positively

correlated with cover of competing vegetation, which in turn was negatively associated

37

with seedling and sapling abundance. I was unable to identify the ideal gap size past
which increased light may hinder regeneration by facilitating competing vegetation due
to high variation in light availability, sapling abundance and cover of competing
vegetation within gaps of the same size class.

The results suggest that removal of competing vegetation may be necessary to
ensure successful regeneration (Romagosa and Robison 2003; Horsley and Marquis
1983). Significant increases in height and radial growth were only observed in sugar
maple seedlings protected from both competing vegetation and deer herbivory (Yawney
and Carl Jr. 1970). Reducing densities of advanced regeneration ironwood may also
benefit other species by reducing competition for light and soil nutrients. Results also
highlight the importance of avoiding damage to well-formed advanced regeneration sugar
maple; few saplings > 2.5 m were gap colonizers, and advanced regeneration has a higher
chance of gap capture due to its stature (McClure et a1 2000; Cole and Lorimer 2005).

White-tailed deer browsing had stronger effects on sugar maple sapling
abundance than light or competing vegetation, so attention should be devoted to quality
deer management that maintains populations in balance with their habitat (F rawley 2005).
This management approach could strike a balance between deer and forest resources
enjoyed by people in the Western Upper Peninsula. Targeted deer harvests that reduce
local densities could allow seedlings and saplings an opportunity to outgrow the browse
line (Sage et a1 2003). Deer harvests may be more successful at improving sugar maple
regeneration in the southern portion of the Western Upper Peninsula if regulations are

altered to permit hunting during the winter months when deer yard in their winter range.

38

1.5 Conclusion

Lake Superior may play a strong role in the regeneration failure of sugar maple
saplings observed in the southern portion of the study area. Lake-effect snow leads to the
southern migration of many deer to their winter range, and may contribute to a gradient in
soil moisture and nutrient conditions. Sugar maple abundance was lower in areas with
lower snow depth, higher deer densities and less nutrient rich Habitat Types. Habitat
Types appear to have the greatest impact on sugar maple small sapling regeneration,
followed by deer density, canopy openness and cover of competing vegetation.

Previous studies have found effects of one to several of these gap- and stand-level
variables, but few have undertaken a comprehensive assessment of these variables across
a large geographic area. Even though these factors were unable to explain much of the
variability observed in seedling and sapling abundances, the results suggest that stand-
level factors have a larger impact on regeneration than gap-level factors. Temporal and
spatial variability in patterns of seed dispersal and seedling establishment (Clark et a1
1998), heterogeneity in the suitability and availability of seed bed microhabitats (Houle
1992; Marx and Walters 2008), individual-level variation in sapling growth rates (Dietze
and Clark in prep) and the chance presence of advanced regeneration saplings in the 1
location of harvest gaps (Brokaw and Busing 2000) all contribute to unpredictability in
the abundance and composition of gap-level regeneration.

Studies that focus on multiple factors will provide greater insight into the
complexities of the regeneration process, with potential implications for our
understanding of the process of succession, the mechanisms that account for species

distributions across landscapes, and the maintenance of biodiversity in forest

39

communities. An integrated, multivariate approach that matches harvest practices to
Habitat Types and local deer densities, and considers interactions between light levels
and competing vegetation, is also necessary for management practices when the goal is to
encourage the success of commercially valuable species, maintain biodiversity, and

improve wildlife habitat.

40

 

 

 

Harvest l
Regime I

Snow Depth

  

Winter Deer

   
 

  

Sapling
Population
>2m

 
   
 

2m Browse Line

 

 

 

 
 

 

Density

  

 

 

Gap Size and

 

 

Browse
Pressure

 

 

 

Competing
Vegetation

 

 

i

 

 

 

 

Density

 

 

Mature Tree

Gap Light
Availability

 

 

 

L,

 

Habitat
Type

 

 

 

 

 

Figure 1.1.

Size and Density

 

 

 

 

Local Seed
Source

 

 

 

 

Bold= operating at landscape-scale
Italics: operating at stand-scale
Underlined= operating at gap-scale

 

 

 

Growth/Survival

Sapling
Population
l-2m

Growth/Survival

Seedling
Population
<1 m

Germination

Seed
Population

Hypothesized factors affecting northern hardwood regeneration
throughout seedling and sapling ontogeny and the relationships between
them. This is not an all inclusive diagram, but highlights many of the
factors the literature and field observations emphasize as important.

41

 

 

 

 

 

 

  
 

 

 

 

 

 

Dug-:61:- Kilometer

12 18

 

 

 

 

 

 

 

Legend 0 175 350 700 1,050
0 Study stands P Kilometer
‘ ' A 1’
:_ _ _i Study area ,’ ~\\ m’ .-
[::| Counties F ’ o ‘ x i WM“,
i O \ x k / ‘
I o x I .
I ~ Wl
’ , Ml ON 4:
I I h
[I
I
I V \ x \
I, \..
,’ ° Marquette County \
I \
\
t . .
O O
I O 0
It . .. I!
I O O l
I o 9
q \ o I!
\ O Q I
\ \ o . ’ .—
\ \ O l!
\ . I
Iron County \ \ o o . 1’
\ Dickinson County 0 P 1’
\ . Q o ,’
\ \ I
\ Q [l
\ \ . . I
\ \ I
\ O i
\ I
‘ \ o O /.
\ \ \ O .' I
\ I
\ \ 0 fl

b

 

‘X

‘ , ’ Menominee County

 

 

 

Figure 1.2. Map of study area. MI=Michigan, WI= Wisconsin, ON: Ontario, Canada.

42

Table 1.1.

Northern hardwood Habitat Type descriptions (Burger and Kotar 2003).

 

 

Soil Soil
Habitat Type Primary landforms Primary soils moisture nutrient
regime regime
AOCa
(A cer saccharum / Moraines, especially Well to moderately Mesic Rich to
Osmorhiza Claytoni — ground moraines, and drained silt loam very rich
Caulophyllum thalictroides) loess deposits and loam texture
soils
ATD-Hp
(Acer saccharum- Tsuga Medium texture till Well drained and Mesic Medium
canadensis / Dryopteris plains well developed to rich
spinulosa- Hepatica variant) sandy loam soils
ATD Coarse texture Moderately well Mesic Medium
(Acer saccharum- Tsuga moraines, especially drained podzolized to rich
canadensis / Dryopteris ground moraines, and or well developed
spinulosa) loess deposits sandy and loam
textured soils
ATM . . .
End morames and Well to moderately Dry-mesrc Medium
(A cer saccharum- Tsuga . .
. . outwash covered well drained sandy to mesrc
canadensrs/ Maranthemum .
ground moraines loams and loamy
canadense)
sands
TMC Moraines, scattered Poorly drained Mesic to Medium
(Tsuga canadensis / low-lying areas, along podzolized sandy wet-mesic
Maianthemum canadense- slope bottoms and loams, loamy sands
Coptis groenlandica) drainages, and on lake and loams

and swamp borders

 

43‘

 

10000 _
(649,612)

1000 _
(64,961)

100 _
(6,496)

10 a
(650)

Seedlings or saplings I plot (Ihectare of gap area)

 

 

R2 = 0.95 ,
F-stat= 559.2 (df= 1 and 26, p< 0.001) .

 

 

Figure 1.3.

I l l l l r r I I I I
0.125 0.375 0.625 1.125 2.125 5.125

Height (m)
Relationship between maximum observed seedling and sapling density in
a 154 m2 gap plot among 347 study gaps versus height. This relationship
was used to develop stocking estimates. Loge(Density)= -2.316 X
Loge(Height) + 5.738.

44

ootfiooowoo c8533 2 wozoom a 3:865 633:; 288.6 “M3. Jaye: mom—o8 Lzomm .233:
8:7. 2:: ”mm... 5256 but MOO dob. “85¢: 0...: .8832 98 moot 232: .2 "ms:— .moo.o 0:52: .: "2.: J85: mo 2:: 05 on 33:38 55030.3
38 Homogmdm £98 5 39:83 8:335 33. ”woomdm .:o_§owo> wozooEoo ">0 .3282? Room... HOD .wotmmm Ethm $55.52 EHDHZ

 

 

 

 

 

 

 

 

 

 

 

8 "a a N6 32: 2o 6.8 seasoned
o: E as. N; 5538 Exam
dome . 9321? .8 n a: A» No 51sz ESE Amwozvoom .Etcm
$5286. meadows
u .. wozoam
1.6.322? 213;. 31;: .2 am: .
N .a C V 2 o .3392: 22.5 25.8%
lee—52 :2
He do .E .>o do as... .292 ex + on "a a No .32: 2o one £36523
:2. do .E 2.: .>o do as E as 3 $539. .522
:me do .5 2.: .>o do 3% AN» + 932...? 9x + .8 u E a No .322 22o 35.83 .523
one do .E sofizom .>o do 3:28: waammowoo
. . , 235:2 .>o do
:3. DO .5.— Ad 0 AN» + 9322.? 9x + .3 n 31va— Qm + b \ .2 no 2240 E N; 353mm
:3 do .E 88% .>o do N o aosowozi. 228%.
class. :5.
o .2 9x + on ”a a N6 .322 22 as. 532023
E 2 Ba 3 $539. .552
m .2 dome AN» + 9322.? so u :1 9 No .31in ESE Ame—uoom .mvtcm
Amos—82 wozoamvwoq
. i e \ ._ "a . E 3 acid
. As o $7. $22.8 a u 2 we A- +
m z N . C o _ o .3532: 22.6 $282
.682 26; E: m
2.3 Go 3 8
2823.5 _o>o_-u§m E08695 322.95 .286 35.65% .39: _o>o_-om0 ootoebma @0522 032...;

 

.52: .3536 ooSE:mH.I .3282 “a .N new 8m 2:9 @82on 2: fl :1 A25 2358 30:: one ASH—ad
2358 35:22: Soc: .022va 3258 65.258 28:: uoE—Eooow :2: use =£ 6:65 13on 8a oomomommooom 332 .NA 033.

45

Table 1.3. Abundance of saplings 1-7 m tall per 154 m2 gap plot (per hectare of gap
area): Average (Ave), median (Med), standard deviation (Stdev),
minimum (Min) and maximum (Max) and percent of gap plots where
present (% gaps) (n=347 gaps).

 

 

axe/£3; Med Stdev Min Max % gaps

Sugar maple 39.3 1 89.0 0 494 52
Acer saccharum Marsh. (2,553) (65) (5,782) (32,091)
Ironwood . 10.5 1 25.7 0 188 57
Ostrya virginiana (Mrll.) K. Koch (6,82) (65) (1,670) (12,213)
White ash 7.3 0 36.5 0 373 19
F raxinus americana L. (474) (2,371) (24,231)
Red maple 5.0 0 24.1 0 299 14
Acer rubrum L. (325) (1,566) (19,423)
Black cherry 2.6 0 8.5 0 94 32
Prunus serotina Ehrh. (169) (552) (6,106)
Quaking aspen . 0.7 0 7.6 0 132 4
Populus tremuloides Mrchx. (45) (494) (8,575)
Balsam fir . 0.6 0 3.6 0 42 10
Abies balsamea (L.) Mrll. (39) (234) (2,728)
White spruce 0.1 0 0.4 0 4 5
Picea glauca (Moench) Voss (6) (26) (260)
Balsam poplar 0.1 0 0.8 0 11 2
Populus balsamifera L. (6) (52) (715)
American elm 0.1 0 0.5 0 6 4
Ulmus americana L. (6) (32) (390)
Yellow birch . 0.1 0 0.4 0 5 3
Betula alleghaniensis Brrtton (6) (26) (325)
Bigtooth aspen . 0.1 0 1.1 0 20 l
Populus grandidentata Mrchx. (6) (71) (1,299)
Northern red oak 0.1 0 0.5 0 6 3
Quercus rubra L. (6) (32) (390)
American basswood <01 0 0.3 0 3 3
Tilia americana L. (<6) (19) (195)
American beech <01 0 0.2 0 2 1
F agus grandifolia Ehrh. (<6) (I 3) (130)
Black ash <01 0 0.1 0 l l
Fraxinus nigra Marsh. (<6) (6) (65)
Paper birch <01 0 0.1 O 1 1
Betula papyrifera Marshall (<6) (6) (65)
Tamarack <01 0 0.1 0 1 1
Larix laricina (Du Roi) K. Koch (<6) (6) (65)

. 66.5 16 108.8 575
A" 51’6““ (4,320) (1,039) (7,068) 0 (37,353) 88

 

46

 

 

 

 

 

Table 1.4. Sapling type by species and height. Gap colonizers germinated alter gap
formation and advanced regeneration saplings were present before harvest.
Gap Advanced Number of
colonizers (%) regeneration (%) samples

Sugar maple

l-<2mtall 53 47 118

2 - <4 m tall 9 91 93

4 - <7 m tall 0 100 81

Ironwood

l - <2 m tall 34 66 89

2 - <4 m tall 12 88 74

4 - <7 m tall 0 100 50

Table 1.5. Summary of gap and stand-level independent variables: Average (Ave),
median (Med), standard deviation (Stdev), minimum (Min), maximum
(Max) and number of gap or stand-level observations (N obs).
Ave Med Stdev Min Max N obs

Gap-level
Competing veg (%) 42 41 23 0 97 347
Gap size (m2) 191 156 103 82 913 347
Canopy openness (%) 13.4 1 1.8 7.3 1.7 54.8 337
5132008 sugar maple 6.8 6.5 2.7 0.9 15.9 347
SPHarvest sugar maple 15.9 11.4 13.6 1.7 83.6 341
spzoog ironwood 0.4 0.00 1.5 0.00 16.6 347
spHamst ironwood 1.1 0.00 5.4 0.00 67.8 341
3132008 other spp 1.8 0.6 2.8 0.00 22.7 347
spHmest other spp 4.6 1.3 8.7 0.00 66.5 341
Stand-level
Deer density (deer / kmz) 13.9 12.4 1 1.7 0.6 61.6 59
SMBI 2.2 2.3 0.9 0.4 3.6 25
Snow depth (cm) 26.6 23.3 7.4 15.6 48 59
TSHarvest (years) 9 9 3 2 15 59

 

Competing veg= % cover of competing vegetation (Rubus spp, grass/Carex spp, other shrubs and
ferns), gap size= extended gap size, SP = seedling production potential (Zdiameter/distance ) for 2008

2
and at time of harvest, deer density= estimated deer / km winter 2007-2008, SMBI= stand average
sugar maple browse index, snow depth= average snow depth from November 2007-April 2008,
TSHarvest= time since harvest.

47

 

 

Table 1.6. Kendall tau rank correlations between gap and stand-level variables.
Comparisons between gap-level and stand-level variables use stand-
averages for gap-level variables to comply with independence
assumptions.

CV Gap size CO Deer density SMBI TSHarvest
Snow depth 0.13 0.06 0.1 7* -0.21** -0.01 -0.03
TSHarvest -0.23*"' -0.05 -0.30""'I 0.16“ 0.34"
SBMI 0.06 0.17 0.01 0.22
Deer density 0.06 0.08 -0.04
CO 0.28" 0.30'Ml
Gap size 0.18"

 

*‘p-value < 0.05, ‘p—value < 0.10

CV= % cover of competing vegetation (Rubus spp, grass/Carex spp, other shrubs andzfems), gap size=
extended gap size (m ), CO = % canopy openness, deer density= estimated deer / km winter 2007-2008,
SMBI= stand average sugar maple browse index, snow depth= average snow depth from November 2007-
April 2008, TSHarvest= time since harvest.

48

2.6.235 89a 399:2 83 03:. :23» :33 afioEcwa 8: 0.33 28:889. omgbmm

 

 

3.33: :8 ~on oomv on? 038.8 32 .3 8:32. 523:: 89a BEEP. we; 02... :23» “Sacco—ma no? xm

“mug: 00min

oETomoEaT—mh smote: .3 3:5 3 com woom :8 Amoo:8m§.338m_vwv 33:33: :23:on wfizooom H mm £32: 8395 2:2: Ema owes; 38m "52m
«6.2-3 3-2. 3.65 I E... :13 on; 4 2.2 E... .656 665:2

vwmd v 86 mom 35: woommm

82. v .82.. .8382 28.51%

£3 4 a .m 6838.: 88%.

~35 v and 2:5: Lama autumn—m

Mood v 0.: 2%.: 3mg. woomdm

and: A 80... 2.2.: 92-3 02.3 2-3 32.? 2o... 6 was :5 as:
God a 2.” sea 650

~cmod v _m.a dam “Soak

not: 31.? 3:: 33.2 m. _ New Rod 4 8.2 as Em 869380

90.0 v mm.m :o_§owo> $323.50

8.3 a. 3.... 86886 Basso

R. .c v t .5 £28 6% 865:3

Sod m 42 32.6,

.23 a mom 85mg

33:22: de< m. _ T... _ m 3298 33.3 3.3 .2 4.8-6.8 52v 4 8.8 23 58.0 38m
83:29» 3379:“ omega v.33 to: 3.3.6.35

com 8: : 8:2 85 32 2322889228 28.. .z

28538 56252 o2e 2: o: Ed: 89.. 8.9.: 46 mg 632.5

 

353?: ©8353 EoboEom ~23 “moo “38-3:me :oxooz3 :0 women 2: God v 3 382388
38$:me 39¢. 333$ 533:3 3&3 bEwomea 35 83:5; 82: :8 Bow: 2: mowofl outdooéoofi .28

83 3.8 333-3533 :0 women 22,—. 3:3: .3 833.8.» 33.79% 39:33 :56 3:: 3:373:36 E moo:o:o.w5

.5; 33m...

49

o :2:0>0 .0: 0000 .0 $0: 00820:... .0 :2:0>0 .0: 0000 832:. 0.0.00.0 .83 0080.05:

 

 

 

8.3. 8.3 cm \ a: 3.0 I v. .o- 3.: 0m... I m _ .0. cm... 083:0...
3.3. :83. em \ R. 2.0 I and- o. .o- and I a. .0 30m... 2%.: 8wsm
A.» \ :. .008: :28 08: 530% w:._%mvw0u.
2.3- .08- 03 SN .0... I 8.0. .0. .0- 8... I 8...- 8... a.-. 85.80
00.0w- mndw- m: \ on. 3.0- I 3.0- 1mm..- um... I cad. 8.0 E N. mw:._%m..
3.00.- 3.8.- cm \ mm. no.0 I mmd- .50. o. .o- I 9.6. 3.3.0- mw:._000m
05:8
86$ 38$ 31m:— mod I :00 :MNN 006. I 2.0- *2... ::m :050
5.30 8.30 :0 \ X: 8.. I 3... 00.0- 0...: I mad- 3.0- 0003:8—
mmiqm and? $100. owd I .3... 1.00.... um... I 00.. 1.2.m 2%.: 885
0o... 80 Na 8. E... 85.80 05.8.0003
3.03. 3.53 «2 .mm on... I mud- 0.... N0: I mo. imam ::m :050
madam. 00.32 mm \ .mm mod I 3.0. no..- 02.. I m5». mm..- 083:8.
2.2m. $05. «3 .mm Quhlhm. Sim... 0:..NI ZN. 3.00. 2%.: 885
.2: :0w «8 0m. \ E m; 0080§0< w:._%m
otmm: no. _ m: an \ 5mm ww.m I 2.0. .00.. afim I m0..- 3.: ::m :050
3...: 00.00: :m \ 5H 00.0 I 09m- 00.0- :0. T I «0.0. 1.3.0.- 00035:.
3.32m wcdfm on \ 5mm wmd I 00.0 :0... 2.... I m..m 33mm 2%.: :85
.2: 9% NE m \ .: m. 00:00:00.: w:._000m
w .0 80E .U 802
2000:: .35 2000.: 8.8%. 008.0 \ 0%
0.0 0. 0 .z 05.3: 05532 0.8.5

 

88:88: 80:. 5:. .805 3.82.5 888$ 08 20>2 3.0.00.0 28. :. N. 00.80550 m:._%m :. 050.8: 8.8:m

2000:. 2:: 0... :0>0 .0008 8.8% 0... :0. 00:88:88: 00>0::E. 0080.0... 0:0...
:0 m .0 05 :. 080.000 < 2000:: 80:: 08: 530% :0. .:00:0.:. .8:0>0 0... 0.8 2000:. 85:23.: :0. 3:00:03. .05.

-08.: 800:8: 08 .:00:0.:. .8:0>0 0:. 0030:. 2000.: .0. Z A0080.0:. 0:0..3. 30.8.5088: :0 2:..-w0. 8 .0 0m: 0...
0. 0:0 080.: .00: 083.00 :0 mo. 0... :0 08 0088.80 8.2.88: .2 0. o 80:. 028000: :0...mo: 2RD. w:..8m :0 w:.....:0 Z
5.3 28.8.. 8... :. 00808:. 0820:. 0038.. 03.80: .m:0..:n....m.0 8.8.00: 80:. 832:. 03.008 gm: 0% 0088.80

.w.. 2.8...

50

 

 

Legend
Sugar maple 1-2 m / !
gap plot -

O O

0 1-10
0 11-50

0 51-100 ,’
. 101-227 /

—-- r

i_ _ _' Study area

l **| Counties

 

 

 

 

 

 

 

 

 

r
l
i

036

 

_ Kilometer

12 18“. 24

\ /

‘0’

 

 

l
x
\

\ l

l /
\ x

l
l

 

Figure 1.4.

Site average sugar maple sapling (1 -2 m tall) abundance per gap plot
across the study area. For reference, 10 saplings / gap plot is 650 saplings /
hectare and 227 saplings / gap plot is 14,746 saplings / hectare of gap area.

51

 

 

Legend

.7
Ironwood 1-2 m / i
gap plot a
_—____J
O O , ’
0 1-10 9‘00
0 11-25 ,I O

0 26-50
. 50-82

n- _ _' Study area

I ___._, [Counties

 

 

 

 

 

 

A \ t 0 w
\\\\i O /Q/
0 3 6 12 18 24 l\
ISIS—:2— Kilometer I \ ,
i U

 

 

 

____._J

 

 

 

Figure 1.5.

Site average ironwood sapling (1-2 m tall) abundance per gap plot across

the study area. For reference, 10 saplings / gap plot is 650 saplings /
hectare and 82 saplings / gap plot is 5,327 saplings / hectare of gap area.

52

 

Legend

Other species 1-2 m / i
gap plot

0 1-10
0 11-50 ,'

. 51-100 /
. 101-208 //

—-- l

i_ _ _' Study area '

I “-1 Counties

 

 

0 O ,.

 

 

 

o 36 12 18“] 24
lI-::_ Kilometer

 

 

:
A

 

r /°

0 I
. 0'
///

l

I
i\\
L :1

~—.._

 

 

Figure 1.6.

Site average other species sapling (1-2 m tall) abundance per gap plot
across the study area. For reference, 10 saplings / gap plot is 650 saplings /
hectare and 208 saplings / gap plot is 13,512 saplings / hectare of gap area.

53

Table 1.9. Spatial patterns in gap- and stand-level variables: bootstrap parameter
estimate and 95% confidence intervals for stand-level variables or stand
averages of gap-level variables. Positive values indicate increases in that
variable with Northing or Easting (UTM position rescaled from 0 to 1).
Gap- and stand-level variables were not standardized for this analysis.

 

 

 

Variable Northing Easting stziqnds
Estimate C1 Estimate Cl
Snow depth 22.24“ 16.29 —— 28.83 -4.75* -8.57 —- 0.08 59
Time since harvest -1.43 -4.04 — 1.45 -1.81 -4.89 - 1.39 59
SMBI -4.03** -5.26 — -3.02 -1.66** -2.20 — -1 . 10 25
Estimated deer density -17.51** -26.27 — -9.42 -10.48** -23.99 — -0.33 59
Canopy openness 4.26 -1.03 — 9.16 1.60 -3.14 — 7.26 59
Competing veg -5.29 -23.79 — 13.94 -21.28** -37.23 — -4.51 59
SP2008 sugar maple ~0.55 -2.34 — 1.47 0.68 -1.25 — 2.73 59
SPHarvest sugar maple 10.64" 3.76 — 17.56 1.43 -4.93 - 8.18 58
SP2008 ironwood -0.94** -1.86 — -0.31 -0.38** -0.95 - -0.04 59
SPHarvest ironwood -1.65** -3.40 — -0.45 -0.54 -l.86 — 0.63 58
SP2008 other spp -0.55 -2.64 — 1.32 -l .21 -3.49 — 0.60 59
SPHarvest other spp -1.51 -7.62 — 4.42 -0.69 -6.84 - 4.26 58

 

** indicates 95% C1 does not overlap 0, * indicates 90% C1 does not overlap 0
SMBI= sugar maple browse index.

54

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

1.0 -
E] Easting
I Northing
0.8 A
c»
.E
.r:
t
o
z
3 0.6 fl
0:
.5
iii
(I
m
E 0.4 " 1
co
0
in
o
a:
0.2 —
0.0 - J L—l
a A ab AB
1 T 1 1 r
AOCa ATD-Hp ATD ATM TMC
Habltat Type
Figure 1.7. Spatial distributions of Habitat Types by UTM Easting and Northing

(scaled to 0-1). Letters indicate significant differences between Habitat
Types based on Wilcoxon signed-rank test with Bonferroni adjusted p-
value. Violin plots show distribution of values and median (white dot),
interquartile range (thick vertical line) and range of values 1.5 * IQR (thin
vertical line).

55

Figure 1.8.

Effects of gap and stand-level variables on seedling (<1 m tall) abundance:
mean parameter estimates, 95% credible interval (thick line), and 90%
credible interval (thin line) from posterior distributions. Positive values
indicate increases in seedling abundance with a one standard deviation
increase in that covariate on the log scale. AOCa is the reference Habitat

Type. SM=sugar maple, IW=ironwood, Other=other species, SP2003=
seedling production potential (Zdiameter/distancez) for 2008, Random Int.

Stdev= standard deviation for the random stand-level intercept. Exact
parameter values can be found in Appendix E, Table E. l.

56

Gap-level Parameter Estimates

 

 

 

 

 

 

 

 

 

 

 

 

 

-0.6 -O.4 -0.2 0.0 0.2 0.4 0.6
l l L l I 1 I
Gap-level :
——_‘——
Canopy a 1:
Openness u _ gunman——
competing —J_:—_
Vegetation fl mu
—__’_—
SP2008 — ‘ :
1 W _
Stand-level .
Intercept - -D-
.— ” 0 SM
ID- IW
*- @, Other
‘ _ "7
fl
-‘ ,F’
“-
i_,_
-i-—_—
TMC —+ C: __
Deer '05
. d}
Densrty fl ,_ ,_
Time since _ J '0'
harVGSt q ’p
' 0
Random E , . .0-
lnt. Stdev q: _
l l 1 l l l 1
-6.0 -4.0 -2.0 0.0 0.2 0.4 0.6

Site-level Parameter Estimates
Factors Affecting Seedling Abundance

57

Figure 1.9.

Effects of gap and stand-level variables on small sapling (1 -2 m tall)
abundance: mean parameter estimates, 95% credible interval (thick line),
and 90% credible interval (thin line) from posterior distributions. Positive
values indicate increases in sapling abundance with a one standard
deviation increase in that covariate on the log scale. AOCa and ATD-Hp
are reference Habitat Types for sugar maple and AOCa is the reference
Habitat Type for ironwood and other species. SM=sugar maple,

IW=ironwood, Other=other species, SPHWM= seedling production

potential (Zdiameter/distancez) at time of harvest, Random Int. Stdev=
standard deviation for the random stand-level intercept. Exact parameter
values can be found in Appendix E, Table E2.

58

Gap-level Parameter Estimates

 

 

 

 

-O.8 -O.4 0.0 0.4 0.8 1.2
l 1 l L l 1
Gap-level
—__—
Canopy_ £1
Openness —- —-—-
Competing _ ——j——
Vegetation ....__.. ,_____D
*-
SPHarvest— ——:_—
Stand-level —-—-o—
Intercept — —D—r 0 SM
C1- IW
153' Other
ATD-Hp _,

 

 

 

 

 

 

 

—__—
ATD —1
——_—
—_—-—
TMC — 1
Deer g -°-_D_
Densrty - .—
Time since __ i?
harvest .- ,—
Random '0-
lnt. Stdev .- -
l l l l 1 1
-8.0 -4.0 0.0 4.0 8.0 12.0

Site-level Parameter Estimates
Factors Affecting Small Sapling Abundance

Figure 1.10.

Effects of gap and stand-level variables on log transformed sapling (1-7 m
tall) stocking-level: mean parameter estimates, 95% credible interval
(thick line), and 90% credible interval (thin line) from posterior
distributions. Positive values indicate increases in sapling stocking on the
log scale with a one standard deviation increase in that covariate. AOCa is
the reference Habitat Type. SM=sugar maple, IW=ironwood, Other=other
species, SPHawesfi seedling production potential (Ediameter/distancez) at
time of harvest, Random Int. Stdev= standard deviation for the random

stand-level intercept. Exact parameter values can be found in Appendix B,
Table E3.

60

Gap-level Parameter Estimates

 

 

 

 

 

 

 

 

 

 

-0.75 -0.50 -0.25 0.00 0.25 0.50 0.75
l l L I 1 I 1
Gap-level
——_—
canopy __ ——-:——
Openness m_ pun—'— ~—
Competing ...
Vegetation -i "a nun-.—
——'——
SPHarvest — “-
— _
Stand-level -—o—-
— -—:_—
lntercept 0 SM
, {:1- W
{3. Other
——j_—
ATD-Hp — “
M_*w
‘ “
ATD _. —_j_—’
1 g_«
“
ATM -— ”—
- m,_ ‘—
————
TMC _ ———:——
Deer .5 .0-
Density ... .—
. . *
Time smce _ -D-
harvest .— 3...
Random __ .0-
lnt. Stdev 9. .-
fi 1 T 1 l
-4.0 -2.0 0.0 2.0 4.0

Site-level Parameter Estimates
Factors Affecting Sapling Stocking-Levels

Table 1.10. Deviance Information Criterion for full models (all gap- and stand-level
covariates) and null models (overall intercept and random stand-level
intercepts for multilevel models and the overall intercept for grth rate
linear models) and number of gaps and stands with observations. A
decrease in DIC of 5 or more indicates improved performance for full

model over the null model.

 

 

Variable Full model Null model N. gaps / stands
Seedling abundance (GLMM)

Sugar maple 2172.61 2184.88 337 / 59
Ironwood 701.20 714.41 337 / 59
Other species 1 132.04 1132.46 337 / 59
Small sapling abundance (GLMM)

Sugar maple 1197.10 1218.17 331 /58
Ironwood 1299.96 1296.83 331 / 58
Other species 1468.80 1489.85 331 / 58
Log(Sapling stocking-levels) (LMM)

Sugar maple 540.62 544.32 164 / 45
Ironwood 622.61 625.81 184 / 49
Other species 674.85 675.80 193 / 48
Sqrt(H’) (LMM)

H’ Seedlings -112.27 -100.26 182 / 54
H’ Small Saplings -46.81 -50.06 180 / 42
H’ All Saplings -54.46 -57.47 201 /46
Log(Post-harvest growth rate) (LM)

Sugar maple 2.34 114.86 127 / 26
Ironwood 13.64 84.03 98 / 30

 

GLMM=generalized linear mixed model, LMM= linear mixed model, LM= linear model,

H’=Shannon Diversity Index

62

Figure 1.11. Observed sapling (1-2 m tall) abundance and mean posterior predicted
values and 95% credible intervals per gap plot from generalized linear
mixed models with gap- and stand-level covariates (full model) and the
null model with an overall intercept and stand-level random effect
intercepts.

63

 

 

 

 
 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

1 1 m 11 1 11111191111111 T w
1 1111 d 1
0 d 1 #11
m 1. 1 m
1 o n
1 MW m an
1 11 111111 I r o
1 1
r o
1 0 8
o
o. 11 n.
1 6
e
I. .u
m. 1 m" 4. m”
m 1 s
r l-
a m a
g h
U f.
s 1 o 1 o o
_ a _ _ _ _ _ a
o
11 1 w 1 1 11 1 W1119|111NU f m
11 11 w
111HH 1H.111 o
0 1 110 81 1 2
w 1
1 n. n 1
0 0 o
1111111111J1 3 h 10
1111 1 1 1
11111111110111
1 11111M11101111 .1 m
. - m
1 0 11011.1“. 1 2 . w
11 1 11111111111111.5111111 1 o 1 e 1 1 .1
1 11111 1 o 11 M1111111M
h 111M111n1uv11w 11 1 1 1 1 of 0
p 1101 1 o 0 1 0 1 4 p
3 11111111 ““11 1 o 1 1 1 11 $111111? 0.
m .1111. 1 1 11 11111 a. s
m 1 1u1 11111111.11 1 1 1 1 H 1 .1 . .I m r
g 11.1. 1 111 1.1.
u Hum 1 11 1.1 t
S 111 1‘91 1 0 111 1111111." 1 o O
_ _ _ _ _ _ _ _ 1-1 —
o o o o m o o o o o o o o
o o o m o o o o o o
w a 6 2 3 2 1 o 8 6

cheE. 2:520 #22 So... can \ 32.; 38.30.... ..o_._3mon cues.

100 150 200 250

50

100 1 50 200 250

50

Null Model
Observed sapling (1-2 m) abundance l gap plot

Full Model

Figure 1.12.

Effects of gap and stand-level variables on square root transformed
Shannon Diversity Index (H’) for seedling (<1 m tall) (H’seedling), small
saplings (1-2 m tall) (H’Sapl) and all saplings (1-7 m tall) (H’SapTot):
mean parameter estimates, 95% credible interval (thick line), and 90%
credible interval (thin line) from posterior distributions. Positive values
indicate increases in H’ on the square root scale with a one standard
deviation increase in that covariate. AOCa is the reference Habitat Type.
H’Overstory = H’ of overstory trees for H’Seedling and H’ of overstory
trees and stumps for H’Sapl and H’SapTot. Random Int. Stdev= standard
deviation for the random stand-level intercept. Exact parameter values can
be found in Appendix E, Table E4.

65

Gap-level Parameter Estimates
-0.06 -0.03 0.00 0.03 0.06 0.09

 

 

 

 

 

 

 

 

 

 

 

 

 

L 1 l l 1 1
Gap-level
Canopy _ —T:——
Openness -- “:1 “——
. ————
Competmg __ H
vegetation ~ .— ‘U—- n
——3——
H’Overstory — 1::
w m -
Stand-level
*-
lntercept — ———Ci-
— SF"
: .o. H’Seedling
ATD-Hp -- —_:—-— G H’sap‘l
~—-—-r plan—nu»— 16, H'SapTot
“
ATD 1 ——:_—-
"7“;
ATM _ _—-E'O.=._
*—
TMC ._ 1
Deer _ [3'0-
Density q 1,.-
Time since _ *1:
harvest .11; 1...
Random _J '0'
Int. Stdev .1 ‘D
r I l I 1
-0.3 0.0 0.3 0.6 0.9

Site-level Parameter Estimates
Factors Affecting Shannon Diversity Index

Figure 1.13.

Effects of gap and stand-level variables on log transformed sapling post-
harvest growth rates (m/yr mean parameter estimates, 95% credible
interval (thick line), and 90% credible interval (thin line) from posterior
distributions. Positive values indicate increases in growth rate on the log
scale with a one standard deviation increase in that covariate. AOCa is the
reference Habitat Type. SM=sugar maple, IW=ironwood, Other=other
species, Advanced Regen: bivariate variable indicating if a variable is
advanced regeneration. Exact parameter values can be found in Appendix
B, Table E5.

67

Intercept

Advanced
Regen

Sapling
Height (m)

Canopy
Openness

Competing
Vegetation

ATD-Hp

ATD

ATM

TMC

Deer
density

Time since
harvest

Intercept Parameter Estimates

 

 

 

 

 

 

-1.5 -1.0 -0.5 0.0
L L 1 J 0 SM
£1- IW
“F
q:-
“—
——:——
*-
#-
‘1’:
fl
4"
*-

 

“—

 

————

l

*
——:——

fl

1 l

 

1
-0.50

-0.25 0.0
Parameter Estimates

1

0.25

Factors Affecting Sapling Post-Harvest Growth

0.5

0.4

0.3

P
N

0.1

P
on

Sapling growth rate after harvest (m I yr)
0
'5

P
u

0.2

0.1

Figure 1.14.

0.... O

D DDD
[>89
D

 

e Gap colonizer
. A Ad regeneration

 

 

 

Sugar maple

 

 

—1

l T l

Ironwood

 

1.0

l l r
3.0 4.0 5.0
Height

Sugar maple and ironwood sapling growth rate after harvest for gap

colonizers and advanced regeneration (ad regeneration).

69

800

(51,909)
3.3
.9
é
a 600
*5 (30,977)
a
D.
G
D
0
O
5
g 400
a (25,994)
.D
N
2
0.
fl
E
'5
O
a 200
3 (12,992)
.6
2
O.
0

Figure 1.15.

‘

_

 

 

 

 

 

 

" 0
6
_ . T
l l l 1 I
All at HT to Deer CO CV
Observed AOCaIATD-Hp +1 Stdev +1 Stdev +1 Stdev

Covariate conditions for predictions

Predicting sugar maple abundance / gap plot (or / hectare of gap area) for
an observed stand: median and 95% percentile predictions based on one to
one simulations drawing from the posterior coefficient distributions and
random effects intercept for an observed ATM stand where a gap had 387
saplings in a gap plot (25,140 / ha of gap area). Simulations were nm with
different covariate conditions (holding all else equal): all at observed
covariate values for this gap and stand, changing Habitat Type (HT) to
AOCa/ATD-Hp, increasing deer density from 4 to 16 deer/km2 (observed
to observed plus one standard deviation), increasing canopy openness

(CO) from 5 to 12% and increasing cover of competing vegetation (CV)
from 18 to 41%.

7O

6,000

 

 

 

 

 

 

 

 

(399,797) '
E
g 5 000 + AOCa/ATD-HP
3: (324,909) ‘ "A" ATD
§
:1.
g 4,000 _
a (259,945)
8
C
3 3,000 _
'3 (194,994)
.1:
fl
2
8' 2000 _
E (129,922)
a
a)
3
a) - 1 :
'0 10 s 0 - s s 1
3 1 5 ’ ,g ’1‘ 1
'2 5 (325) = 5 5 5
'0 - : : :
g 1:1. i

0 _ l i 9 i 1 5 6
1 1 1 1
All at Deer CO CV
Average +1 Stdev +1 Stdev +1 Stdev
Covariate conditions for predictions
Figure 1.16. Predicting sugar maple abundance / gap plot (or / hectare of gap area) for

an unmeasured stand: Median and 95% percentile predictions based on
one to one simulations drawing from the posterior coefficient distributions

and randomly drawing the stand-level intercept from ~N(u, ca) where u is

defined by the stand-level covariates and 0a is the random stand-level
intercept stand deviation. Simulations were run with different covariate
conditions (holding all else equal) for a gap at an AOCa and ATD type
stand: all at average observed covariate values, increasing deer density
from 14 to 26 deer/km2 (average observed to average observed plus one
standard deviation), increasing canopy openness (CO) from 13 to 20% and
increasing cover of competing vegetation (CV) from 42 to 65%. A break
was added to the y-axis so changes in the median prediction could be seen.

71

3 8
1 1

°/o Gaps with >= 1 Same Species Sapling
8
1

20—

 

 

on .. YB 3w oA SM WA RM ec 1w

(14) (24) (20) (6) (191) (6) (21) (5) (1°)
Stump species (n gaps)

Percentage of gaps harvested 8-15 years ago where a species was removed
(stump species) and there was at least one sapling of the same in the gap
plot. Species codes: Con=conifer (eastern hemlock, northern white cedar,
balsam fir and White spruce), YB= yellow birch, BW= American
basswood, QA= quaking aspen, SM= sugar maple, RM= red maple, BC=
black cherry, IW=ironwood. Numbers below species codes indicate the
number of gaps where that species was removed by harvesting.

72

Percentage of Gaps

 

Figure 1.18.

 

 

Stumps
El Saplings

 

 

 

 

 

 

 

   

saplings
Species

Percentage of gaps sugar maple (SM), ironwood (IW) or other species are
the dominant stump species (highest relative basal area of stumps removed
from gap) and percentage of gaps these species are the dominant sapling
species (highest relative stocking level of saplings 1-7 m tall). Percentage
of gaps with no saplings is also presented. N= 211 gaps harvested 8-15
years ago.

73

 

 

 

   
  
 
 
 
 
 
   
 
 
 
 
  
 

6 12 18 24
I—

 

Legend

 

Deer / km
Winter 07-08

0 0.6-5.0
0 5.0-10.0
0 10.0-20.0
0 20.0-30.0

0 30.0-61.6

1_ _'Study area
[j Coun 'es
Deer/ km
1996-2000

54

Kilometers ' 0

   

 

 

Figure 1.19.

Average estimated deer density (deer / kmz) from 1996-2000 created from
universal kriging (Appendix B, Section B5) of 113 observations collected
by the Michigan Department of Natural Resources and deer density
estimates from winter 2007-2008. Estimates for winter 2007-2008 are
positively correlated with the MDNR 1996-2000 estimates (Kendall’s
7:0.27, p-value=0.002).

74

Figure 1.20.

Differences in stand average seedling abundance / 2 m2, small sapling (1-2
m) abundance / gap plot and sapling stocking levels / gap plot by Habitat

Type. Gap plots were 154 m2. Letters indicate significant differences
between Habitat Types for sugar maple based on Wilcoxon signed-rank
test with Bonferroni adjusted p-value. Kruskal-Wallis rank sum test
indicated that sapling abundance and sapling stocking were not
significantly different between Habitat Types for ironwood or other
species, but it was significant for seedling of both species. However, no
pairs of Habitat Types were significantly different based on the Wilcoxon
test. Violin plots show distribution of values and median (white dot),
interquartile range (thick vertical line) and range of values 1.5 * IQR (thin
vertical line).

75

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

I SM
1:1 1W
N 40 1 D Other spp
E
2 30 —
W
e:
.E
0
m
10 -
o -‘ AB -'0—-
.. 200 1
2
Q
a
8150 —
E
3.‘ 100 —
W
U)
.5
Ta. 50 -
N
.11 121 1 2L I C13
_ 213 (3:1
0 AC ‘6— B AB
§ 5 ‘
a.
1%
a: 4 7
3:
.E
x 3 —
.3
U3
E 2 -1
[s
g.
'a
N
11 21 2121 1.30:1
° “ f B B AB
1 1 *F 1 I
AOCa ATD-Hp ATD ATM TMC
(19) (15) (12) (11) (2)
More nutrient rich 5 Less nutrient rich
Habitat Type (# stands)

76

Chapter 2: Snow depth, deciduous-lowland conifer forest edges and
deer density affect seedling browse damage

Abstract

Restoration of eastern hemlock (T suga canadensis (L.) Carriere) and eastern
white pine (Pinus strobus L.) to forests in the Great Lakes region are priority
conservation goals. I planted white pine, hemlock and white spruce (Picea glauca
(Moench) Voss) seedlings in harvest gaps within northern hardwood stands in the
Western Upper Peninsula of Michigan to 1) explore the efficacy of planting efforts in
areas with high winter deer densities and 2) improve our understanding of factors
affecting deer herbivory at the stand- and gap-level. Results show that white spruce
seedlings may successfully reestablish with planting efforts because this species was
virtually un-browsed, but white pine and hemlock will not reestablish without protection
from deer. Stand-level percentages of seedlings browsed increased with estimated winter
deer density and decreased with snow depth and stand-average gap distance to lowland
conifer. Non-linearity was observed in the relationship between browse pressure and
estimated deer density, emphasizing the importance of maintaining deer densities below a
threshold level when the goal is to protect browse-preferred understory vegetation. Snow
sheltered seedlings from deer, but relationships between snow depth and estimated stand
deer density were not observed. Evidence of deer gap-use corresponded with increased
gap-level percentage of seedlings browsed, but contrary to expectations, deer-use within
a stand was higher in gaps farther from lowland conifer. The relationship between browse
pressure and estimated deer density supports the use of the pellet count technique for
estimating deer densities. However, monitoring deer browse directly on seedlings and

saplings provides better insight into impacts of herbivory.

77

2.1 Introduction

Hemlock-white pine-northem hardwood forests in the Great Lake Regions
underwent dramatic transformations in the mid 18005 to early 19003 when logging
virtually eliminated eastern hemlock (Tsuga canadensis (L.) Carriers) and eastern white
pine (Pinus slrobus L.) from the landscape (Whitney 1987; Stearns 1997; MDNR 2006).
Hemlock has not been replaced by natural regeneration due to seed source limitations
(Mladenoff and Stearns 1993), lack of suitable seed beds (Mladenoff and Stearns 1993;
Rooney et a1 2000; Marx and Walters 2008), browsing by white-tailed deer (Swift 1949;
Graham 1954; Rooney et al 2000), unsuitable climatic conditions (Mladenoff and Stearns
1993) and light limitations (Rooney et al 2000). White pine seedlings and saplings are
still found widespread throughout northern hardwood-conifer forests where seed sources
are available (Carleton et a1 1996) but declines have been observed (MDNR 2006). The
absence of a seed source greatly limits the regeneration of this species (Ahlgren 1976;
Dovéiak et a1 2003), but so does the availability of suitable seed beds (Smith 1951;
Dovéiak et a1 2003 ), light availability (Smith 1951; Krueger and Puettmann 2004),
browsing by white-tailed deer (Swift 1949; Saunders and Puettmann 1999; Krueger and
Puettmann 2004), competition from understory vegetation (Ahlgren 1976; Dovciak et a1
2003; George and Bazzaz 1999) and damage from insects and disease (Ahlgren 1976;

Krueger and Puettmann 2004).

. . . 1 . . .
Increasmg the cover of mesrc comfer stands and re-establishmg these specres as

components in hardwood forests are common conservation goals in the Great Lakes

Region (Thorns 1992; Crow et a1 1994; OMNR 2003; MDNR 2006). On State land in the

 

1 Mesic conifers include eastern hemlock, eastern white pine, balsam fir (Abies balsamea (L.) P. Mill.) and
white spruce (Picea glauca (Moench) Voss)

78

Western Upper Peninsula of Michigan, white pine and hemlock planting efforts are
underway (Herman et al 2004). It is predicted that increased cover of pure conifer and
mixed deciduous-conifer stands will improve habitat for bird species (Herman et al 2004)
like the blackbumian warbler (Dendroicafusca) which feeds and nests in hemlock trees
(Kendeigh 1945). Increasing the abundance of conifer cover could also reduce suitability
of white-tailed deer (Odocoileus viriginanus Zimmermann) summer range (Kohn and
Mooty 1971; Tierson et al 1985; Felix et a1 2004) and increase cover of winter habitat
(Graham 1954; Venue 1965; Ozoga 1968; Blouch 1984), resulting in smaller herds
dispersed over larger wintering areas (Herman et a1 2004).

Heavy deer browsing in the Great Lakes region can impact the efficacy of white
pine and hemlock restoration efforts (Rooney et al 2000). Consideration of climatic
conditions and landscape contexts that affect deer behavior and local densities could
improve restoration success. Snow depth and lowland conifer stands influence deer
migration and distribution in the winter months. Many deer in the Upper Peninsula of
Michigan migrate on average 10.9-20.3 km (depending on winter severity), and as far as
about 50 km, from summer to winter ranges where snow depth is lower and availability
of winter habitat is greater (Verme 1973; VanDeelen et a1 1998). Winter white-tailed deer
densities in northern hardwood stands have been found to decrease with increasing
1 distance to lowland conifer, likely due to the importance of conifers for shelter
(Millington et al 2009, in review). Within deciduous forests, deer browse on seedlings
and saplings may increase with proximity to lowland conifer stands.

The most direct method of assessing the impact of deer herbivory on seedlings

and saplings is to monitor browse damage, but it is also important to establish reliable

79

estimates of deer density for informing deer management and modeling the relationships
between deer density and herbivory. Fecal pellet surveys are a relatively simple and
inexpensive method of assessing relative local deer density (Hill 2001; Langdon 2001)
and positive relationships between pellet count metrics and deer density have been
empirically supported (N eff 1968; Bailey and Putman 1981; Forsyth et al 2007).
However, others have found poor performance of the pellet count method compared to
population estimates from aerial surveys (Fuller 1991) and higher variability in pellet
count population estimates compared to mark-resighting and distance sampling
techniques (Langdon 2001). Deer density estimates derived from pellet counts can be
biased because defecation rates are correlated with forage intake, forage moiSture content
and percentage of young in the herd, and they can be affected by sampling design and
observer error (N eff 1968). Regardless, they remain the cheapest and most
implementable method available for assessing relative differences between deer density
at the stand-scale across a large geographic area.

I planted mesic conifer seedlings within northern hardwood stands in the Western
Upper Peninsula of Michigan to test the following hypotheses: H1) seedling browse
damage will be higher at stands with lower snow depth, shorter distance to lowland
conifer and higher winter deer density and H2) within a stand, seedling browse damage
will increase with proximity to lowland conifer where winter deer-use is higher. The
results of this study can aid in the development of effective mesic conifer restoration
efforts and provide insights into the relationship between estimated deer densities and

browse pressure.

80

2.2 Methods
2. 2.1 Study area

This study occurred within 34 northern hardwood stands across Dickinson, Iron,
Marquette and Menominee counties in the Western Upper Peninsula of Michigan (Figure
2.1). Stands were owned by Plum Creek Timber Company Inc. (PCT), GMO Renewable
Resources Inc. (managed by American Forest Management Inc. (AFM)), and the State of
Michigan and had been selection harvested from 2003 to 2007. The study area
corresponds with that for the ongoing Western Upper Peninsula Economic-Ecological
Modeling Project (Laurent et a1 2005; LeBouton et a1 2005; Racevskis and Lupi 2006;
Millington et a1 2009 in review). This area was selected for the domination of forest
cover, the relative absence of agriculture and urban or suburban developments, and the
natural variation in deer densities and snow depth (Doepker et al 1994; Shi et a1 2006).

Northern hardwood forests are found throughout the area on upland topography,
such as drumlins, created by glacial activity and coniferous forests are prevalent in
lowland areas. Aspen and mixed upland cover types are also present. Coniferous forests
are composed of mixes of northern white cedar (Thuja occidentalis L.), black spruce
(Picea mariana (Mill.) Britten & et al.), tamarack (Larix Iaricina (Du Roi) K. Koch) and
balsam fir (A bies balsamea (L.) Mill.) and northern hardwood forests are dominated by
sugar maple (A cer saccharum Marsh), with various components of American basswood
(Tilia americana L.), yellow birch (Betula alleghaniensis Britton), red maple (Acer
rubrum L.), white ash (F raxinus americana L.) and other minor species including eastern

hemlock and white pine.

81

Stands were located within ecoregion sub-sections VIII.3.1 (loamy ground
moraines and drumlin ridges), VIII.3.2 (poorly drained outwash plains), 1X.l (steep
bedrock knobs, outwash plains and sandy ground moraines) and IX.3.2 (irregular ice-
disintegration topography) (Albert 1995). Section VIII is underlain by Cambrian-age
sandstone and Paleozoic limestone, shale and dolomite whereas section IX is underlain
by highly resistant igneous and metamorphic bedrock of the Precambrian Shield.

Annual snow fall varies from 161 cm in the southern portion of the Western
Upper Peninsula to 434 cm in the northern portion (National Climatic Data Center 2009).
Snowfall is higher in the northern parts of the Upper Peninsula due to lake-effect snow
that develops when cooler air masses from the north move over the warmer waters of
Lake Superior (Jerome 2006).

2.2.2 Field methods

Six seedlings of hemlock, white pine and white spruce were planted in ten harvest
gaps at each stand for a total of 2,040 seedlings per species in November 2007. Hemlock
seedlings were the tallest (mean 34.3 cm, stdev 7.1 cm) followed by white pine (26.7 cm,
7 cm) and white spruce (20.8 cm, 4.9 cm). Hemlock and white spruce were two year old
containerized seedlings and white pines were two year old bare root seedlings. If lowland
conifer stands could be identified in the surrounding landscape, gaps were selected at
various distances from the hardwood-conifer edge. The planted seedlings were not
individually marked so as not to attract or deter deer, but within gaps they were planted in
straight lines and maps were made of the orientation of gaps within stands to facilitate

relocation.

82

Stands were revisited to monitor browse damage and perform pellet surveys from
late April to early May 2008 after snow melted and before cover of spring ephemerals
would hinder relocation of planted seedlings. Browse damage for each hemlock seedling
was qualitatively classified into one of five categories: no browsing (0), light browsing
(1), moderate browsing (2), heavy browsing (3) and complete removal of needles (4).
White pine and white spruce seedlings were marked as browsed or unbrowsed because
natural variation in the size and branching patterns of the seedlings made it difficult to
consistently categorize browse damage. Browse damage caused by cottontail rabbits
(Sylvilagusfloridanus J. A. Allen) or snowshoe hares (Lepus americanus Erxleben) was
possible, but likely negligible. Rabbits browse minimally on hemlock (Todd 1927;
Hough 1949), white pine and spruce (Swihart and Yahner 1983), and hare browse
preference is only moderate for white pine and low for spruce and hemlock (Telfer 1972).

To estimate stand-level winter deer density, pellet groups were counted within ten
50 m by 4 m transects oriented in an hour glass shape (LeBouton et a1 2005) (Appendix
B, Section B.l). A handheld GPS unit was used to determine the UTM Northing and
Easting for the center of the pellet surveys. Transects were double counted with different
observers to improve accuracy (Jenkins and Manly 2008). The mean pellet group counts
from all ten transects were used to calculate deer density following the methods of Hill
2001 (Appendix B, Section 8.1). Visible pellets were those above leaf litter deposited in
November 2007, so the estimate of deer density corresponds to deer presence in the late
fall, winter and early spring months. 1 term the estimate “winter deer density” although
deer depositing fecal pellets in northern hardwood stands likely did so while migrating

between summer and winter ranges (i.e. in November and April).

83

Pellet groups were also counted within one 10 m by 4 m transect located parallel
to the planting lines through the center of each planted gap. Pellet groups counted along
these transects were not incorporated into the stand-level deer density estimate. Gap
pellet group counts (PGC) were used as estimates of fine-scale differences in deer-use of
the planted gaps. A deer has a higher chance of defecating in an area where it spends
more time; however, the absence of pellet groups does not necessarily correspond to the
absence of deer-use. Gap-PGC were not converted into deer densities because they were
quantified in such small areas that they could not reasonably be scaled up to a larger area.

A subset of 19 stands was revisited in late April and early May 2009 to reassess
browse damage. Eleven stands were selected because percentage of hemlock seedlings
severely browsed (category 3 or 4) was <50% in 2008 and eight stands were selected
because distance to lowland conifer (DLC) estimates could be obtained for them (DLC
could also be found at three stands selected for hemlock browse criteria). Hemlock
seedlings were only reassessed at stands with low 2008 browse pressure because severe
browse damage made relocation difficult elsewhere. Seedlings with completely
desiccated foliage were noted as was their browse status. White spruce seedling browse
status was not reassessed in 2009 because browse was infrequent the first year.

Distance to lowland conifer (DLC) could only be determined for thirteen stands

due to incomplete stand-level data across the study area. At these stands a Trimble®

Asset SurveyorTM model TSCl GPS receiver and datalogger system with built in real-

time differential correction capabilities was used to determine the location of each
planted gap. GPS locations were not determined for four gaps at one stand due to battery

failure. Distance to lowland conifer for each gap was determined using ERSI ArcMap

84

v9.2 and digitized stand-level data provided by MDNR, AFM and PCT. Although the
Integrated Forest Monitoring, Assessment and Prescription (IFMAP) project provides
forest classification for the entire Upper Peninsula (Donovan 2005), its accuracy has been
challenged even for distinguishing conifer and deciduous cover types (Linden 2006).

The number of weeks from November lSt 2007 to April 30th 2008 when the snow
depth at each stand was deeper than 40 cm was found from daily 1-km resolution maps
created by the National Snow and Ice Data Center’s Snow Data Assimilation System
(SNODAS) (Barrett 2003; NOHRSC 2004). Deer often avoid areas with snow depths
greater than 40 cm (Poole and Mowat 2005) because it hinders their movement (Kelsall
1969), and snow depths approaching this threshold can initiate their migration to winter
range (Tierson et a1 1985). This depth also corresponds closely with the 90th percentile
height of planted hemlock seedlings (41.5 cm). The estimate represents the number of
weeks when snow depth was prohibitive to deer movement and provided cover for
planted seedlings.
2.2.3 Statistical analysis

The goals of this study were to determine browse damage on planted conifer
seedlings, quantify relationships between deer densities, deep snow weeks (DSW) and
DLC and determine their effects on browse pressure. Variation in the percentage of
browsed white spruce was low within and between stands, so analysis focused on the
percentage of white pine and hemlock seedlings browsed at the stand-level and at the gap
level. Since almost all hemlock seedlings showed signs of browse, hemlock seedlings

were classified into two categories: un-browsed / minimally browsed (browse category 0,

85

1 or 2) and severely browsed (category 3 or 4). Hereafter, “browsed” is used to refer to
browsed white pine and severely browse hemlock seedlings.

To test hypothesis 1, non-parametric Kendall’s tau correlations were calculated
between stand-level percentage of seedlings browsed in 2008 and stand-mean gap-DLC,
stand-mean gap-PGC, stand deer density and stand deep snow weeks (DSW). Correlation
between percent browsed in 2009 was only found with stand-mean gap-DLC because
gap-PGC, deer density and DSW were only found for winter 2007-2008.

To address hypotheses 1 and 2, multilevel logistic regression models were
developed for hemlock and white pine to predict the gap-level percentage of seedlings
browsed in 2008 based on stand and gap-level characteristics. Models included a random
stand-level intercept to account for the hierarchical structure of the data (gaps nested
within stands) and to incorporate information about the clustering to produce estimates of
standard errors that account for non-independence (Goldstein 1995). In order to
determine the effects of differences in gap-level PGC and DLC within stands on gap-
level percentage of seedlings browsed (i.e. to isolate the effects of variation in gap-level
PGC and DLC within stands instead of between stands), these variables were centered on

the stand-mean (Enders and Tofighi 2007).

The number of seedlings browsed in gap i at stand j (yij), out of nij seedlings
relocated in 2008 was modeled as a binomial process yij~Bin(piJ-, nij). The percentage of
browsed seedlings (pij) was a function of the inverse logit of X10 + orj, where Xi’s were
gap-level covariates. The random stand-level intercept (aj) came from a normal

. . . 2 , . .
distribution ~N(yo + m, o a) where Zj s were stand-level covariates. Five models were

86

developed with different gap-level (Xi’s) and stand-level (Zj’s) covariates (Table 2.1): 1)
null model (stand-level random intercept only), 2) stand deer density and DSW, 3) stand-
mean centered gap-PGC and DSW, 4) stand-mean centered gap-PGC, gap-DLC and
DSW and 5) stand-mean centered gap-DLC only. Variables were put on more

comparable scales to improve convergence; deer density was rescaled to 10’s of deer /

kmz, gap-DLC was rescaled to 10’s of meters and stand deer density and DSW were

grand mean centered.

Models were run with WinBUGS v.1 .4.3 (Spiegelhalter et a1 2003) (see code
Appendix C, Section C.4) through R v.2.9.1 with the package R2WinBUGS (Sturtz et al
2005) (see code Appendix C, Section C.5). Three parallel chains with dispersed randomly
selected starting values were run for 40,000 iterations with a burn-in of 10,000 and a
thinning rate of 5. Models were run with a uniform prior (range 0-100) for the standard
deviation of the random stand-level intercept. The random intercept standard deviation

was insensitive to prior specification (Appendix B; Section B4). A noninformative prior
distribution (Normal~(u=0, oz=10,000)) was used for gap—level and stand-level
coefficients (Clark 2007; Bolker et a1 2009).

Gap-level percentages of seedlings browsed were overdispersed as is often the
case with count data used in logistic regression (i.e. variance is greater than explained by

the model) (Gelman and Hill 2007). However, the random stand-level intercept helped

account for the overdispersion in the data. Pearson residuals were closer to a standard

normal distribution N~(p=0, oz=1) with the inclusion of a stand-level random effect

(N~(-0.01, 1.37) vs ~N(0.00, 2.66) without stand-level random effect for pine and

N~(0.01, 1.64) vs N~(0.02, 2.99) for hemlock null models).

87

Deviance Information Criterion (DIC), a method of assessing model performance
in terms of fit and complexity based on the posterior mean deviance (-2 x log(likelihood))
and the effective number of parameters (Spiegelhalter et al 2002), was compared between
models to determine if predictor variables improved model performance. A decrease in
DIC by 5 or more indicates support for better model performance.

Convergence was diagnosed using the R package coda (Plummer et a1 2009) with
the Gelman-Rubin diagnostic and the Raferty-Lewis diagnostic. All models showed
strong evidence of convergence (Appendix D, Tables D9 and D.10).

2.3 Results
2.3.1 Summary of gap- and stand-level variables

The number of deep snow weeks (DSW) varied from 0 to 12.6 (Table 2.2) and
was negatively correlated with stand-mean gap-PGC and positively with Northing (Table
2.3). Stand-mean centered gap-PGC varied from 4.1 pellet groups less than the stand
mean to 10.9 groups more than the stand mean. Stand-mean centered gap-DLC varied
from 161 m less than the stand mean to 156 in more than the stand mean. Non-centered
gap-PGC and DLC were not significantly correlated, but stand-mean centered gap-PGC

and DLC were positively correlated (Figure 2.2). Stand-level deer density estimates
ranged from 4.1 to 64.5 deer / km2 and were not correlated with Northing or DSW.
2.3.2 Percentage browsed varied between species and increased over time

Over 97% of planted seedlings for each species were relocated in 2008 (Table

2.4). Two spruce, eight white pine and 15 hemlock seedlings were found pulled out of the

ground, so some seedlings may have been completely removed from the gaps by deer.

88

Percentage of seedlings browsed varied by species (Kruskal-Wallis x2= 86.8, p-

value < 0.001 ), with browsing lowest on white spruce, intermediate on white pine and
highest on hemlock (Figure 2.3). Over winter 2007-2008, 92% of hemlock, 30% of white
pine and 2% of white spruce seedlings were browsed (Table 2.4). Sixty percent of
hemlock seedlings were severely browsed (category 3 or 4).

In 2009, relocation of seedlings was more difficult, likely due to increased
incidence of browse and death from apparent desiccation (Table 2.5). In 2009, 9% of
hemlock and 21% of white pine relocated had completely desiccated foliage. Percentage
of seedlings browsed for both white pine and hemlock increased at all resurveyed stands
from 2008 to 2009. Percentage of hemlock seedlings severely browsed increased from
less than 50% to over 50% for all but one resurveyed stand between spring 2008 and
spring 2009, and increased to as high as 96% at one stand.

2.3.3 Deer density, deep snow weeks, deer gap-use and distance to lowland conifer
aflect stand-level percentages of browsed seedlings (H1)

Stand deer density estimates were positively correlated with the percentage of
hemlock seedlings browsed (p<0.10) (Figure 2.4), but not with percentage of pine
seedlings browsed (Figure 2.5). The impact of deer on percentage of severely browsed

hemlock was nonlinear, with percentage browsed rising rapidly with increasing deer

density. One stand was the exception with high estimated deer density (38 deer / kmz)
and low browse pressure (13%), potentially because there were 9 weeks of deep snow at
this stand. Stand average gap-PGC (related to deer gap-use) was not correlated with

hemlock browse, but showed a similar pattern to the relationship between deer density

and percent severely browsed hemlock. Stand average gap-PGC was positively correlated

89

with white pine percent browsed. Percentages of planted hemlock and white pine
browsed were negatively correlated with DSW. Hemlock browsing was greater than 45%
at all stands where DSW was less than 8, but as low as 3% where DSW was greater.
Stand-level percentages of seedlings browsed for both species were negatively correlated
with stand average gap-DLC. This relationship strengthened between 2008 and 2009 for
white pine.

2.3.4 Deep snow weeks, deer gap-use and distance to lowland conifer affect gap-level
percentages of browsed seedlings (H1 and H2)

Stand deer density did not affect gap-level percentages of seedlings browsed for
either species (Model 2, Table 2.6). In contrast, the gap-level percentage of hemlock and
white pine seedlings browsed decreased with increasing DSW and increased with
increasing stand-mean centered gap-PGC (Model 3). Compared to the null model (Model
1), inclusion of snow and gap-PGC improved model performance for hemlock and
marginally for white pine. However, differences in predicted percentages between Model
3 and the null model (Model 1) were very minor for both species (Figure 2.6).

Distance to lowland conifer and DSW for both species, and gap-PGC for pine,
were unrelated to percentage of seedlings browsed for the subset of stands where DLC
could be quantified (Model 4, Table 2.7). In the univariate model with DLC (Model 5),
increasing stand mean-centered gap-DLC corresponded to increased gap-level
percentages of seedlings browsed for hemlock, but not for pine.

Based on simulations using parameters from the model with mean-centered gap-
DLC and stand DSW (Model 3), increasing the number of DSW from 6.6 (the average

observed) to 7.6 decreased the mean predicted percentage of severely browsed hemlock

90

seedlings in a gap by 2%, and the percentage of browsed white pine by 1% (Table 2.8).
Larger decreases in the percentages were predicted (15% for hemlock and 8% for white
pine) when number of DSW increased from 6.6 to 12.6 (the maximum observed).
Increasing the gap-PGC by one over the stand-mean increased the predicted percentage
of severely browsed hemlock by 4% and the percentage of browsed white pine by 1%.
Larger increases in the percentage (24% for hemlock and 14% for pine) were predicted
when the gap-PGC increased by eleven over the stand-mean. For all simulations, the 95%
quantiles for the predictions were large due to uncertainty in the stand-level random
intercept.

Based on simulations using parameters from the model with only DLC (Model 5)
for hemlock, the mean gap-level percentage of severely browsed seedlings showed no
significant increase if gap-DLC increased 10 m from the stand-mean. If gap-DLC
increased 150 m from the stand-mean, the probability of severe browse increased from
56% (95% quantiles for the predictions: 3-99%) to 66% (5-99%).

2.4 Discussion

A high percentage of planted hemlock and white pine seedlings, both deer
browse-preferred winter foods (Blouch 1984), were browsed over a two year period in
northern hardwood forests. White spruce, a last resort forage species for deer (Blouch
1984; Dumont et a1 2005), was only minimally browsed. Percentages of browsed seedling
were negatively related to abiotic stand characteristics (i.e. length of time when snow
depth prohibits deer movement and shelters seedlings) and local landscape context (i.e.

distance to lowland conifer).

91

2. 4.1 Hypothesis 1: the percentage of browsed seedlings will be higher at stands where
snow depth and distance to lowland conifer are lower and winter deer density is greater

Stand-level percentage of hemlock seedlings severely browsed was positively, but
non-linearly, related to estimated deer density. Non-linear relationships between deer
density and browse impacts have been reported for hemlock and some hardwood species
(Rooney and Waller 2003; Eschtruth and Battles 2008). Variability in deer density
estimates from the pellet count technique (Langdon 2001) may contribute to the relative
weakness in the relationship between browse pressure and deer density. The relationship
between stand-level percentage of white pine browsed and deer appears sensitive to
pellet-count transect location and length; it was not significant with stand-level deer
density estimated from pellet counts observed along long transects radiating from a
random center but was positive with stand-average gap-pellet group counts observed
along short transects in areas with expected deer use. Both transect length and
distribution of transects are important considerations when designing fecal pellet surveys
(N eff 1968). Counting pellets in patches where deer use may be higher (i.e. in harvest
gaps with greater forage availability) would likely over-estimate stand deer densities, but
might assist with monitoring relationships between relative deer density and herbivory at
a finer spatial scale.

Stand-level percentages of white pine and hemlock seedlings browsed decreased
with increasing length of deep snow weeks and stand-average gap-distance to lowland
conifer. Browse pressure could be lower at stands where snow depth is > 40 cm for a
longer portion of the winter because: 1) snow provides cover to seedlings from deer, 2)

many deer migrate south to winter range where snow depth is lower (V erme 1973;

92

Tierson et a1 1985; Blouch 1984; VanDeelen 1998) and 3) deep snow and winter severity
cause deer to alter their foraging behavior and habitat selection to conserve energy
(Ozoga and Verme 1970; Pauley et al 1993; Dumont et al 2005; Poole and Mowat 2005).

Both deep snow weeks and deer density varied across northern hardwood stands
and had the anticipated relationships with stand-level percentage of seedlings browsed,
but were not negatively correlated with each other as was expected. I found a negative
relationship between deer density and snow depth in this region across a larger study area
that included more stands and extended farther north (Chapter 1). The latitudinal gradient
of the stands for this planting study and the number of stands measured may not be great
enough to capture the inverse relationship between deer density and snow fall, given high
variability in pellet count derived estimates of deer density (Langdon 2001).

Unlike deer density, stand-average gap-pellet group counts (i.e. deer gap-use)
were negatively correlated with deep snow weeks. It is possible that the relationship
between snow and deer was captured by counting pellet groups in gaps because deer
might avoid gaps where snow depth is greater due to lower interception of snow by tree
branches.

The observed relationship between stand-level browse and distance to lowland
conifer support previous findings that deer density is lower in stands farther from lowland
conifer stands (Millington et a1 2009 in review). The correlation between deer density
estimates and stand-average gap-distance to lowland conifer was not significant in this
study, likely because distance to lowland conifer could only be determined for thirteen

stands.

93

2.4.3 Hypothesis 2: percentage of seedlings browsed will increase with proximity to
lowland conifer where deer-use is higher

Browsing of hemlock and white pine seedlings increased with deer gap-use
(stand-mean centered gap-pellet group counts) and with gap distance to lowland conifer
(stand-mean centered gap-level distance to lowland conifer) for hemlock, but not pine.
Deer gap-use was positively correlated with gap distance to lowland conifer. The later
findings are contrary to the hypothesis that deer browse would be lower in gaps farther
from lowland conifer. Deer confinement to conifer stands increases with winter severity
(Ozoga and Gysel 1972), so deer foraging may not have been limited to gaps near conifer
stands during this fairly average winter (MDNR 2008). Also, deer disperse farther from
habitat edges when highly-browse preferred species are available (Williamson and Hirth
1985). Within stands, planted gaps may not have been located sufficiently far from each
other to infer small-scale difference in deer use of a stand. Daily winter movement for a
deer can be 341 to 1,650 meters (Heezen and Tester 1967) and the farthest distance
between gaps within a stand was 730 m (average maximum distance was 370 m). Within
one day, one deer could have visited all gaps at a stand.

The relationships between browse, gap distance to lowland conifer and deer gap-
use could also be an artifact of the uncertainty in the gap pellet group counts and the
small sample size for gap distance to lowland conifer. Gap pellet group counts were
quantified in small areas and was not a perfect estimate of deer use of a gap (i.e. there

were no gap-pellet groups in 40% of gaps where there was evidence of deer browse).

94

2.4.4 Management implications
The efficacy of efforts to restore hemlock and white pine in northern hardwood

stands with replanting may be minimal in areas with high deer densities (> 15 deer /

kmz). Hemlock seedlings will not grow into larger size classes unless protected from

white-tailed deer. White pine browse was relatively low the first year of planting, but
increased to fifty percent at re-surveyed stands after the second year. Many white pine
seedlings desiccated after two winters and one growing seasons, often independent of
browse damage, so planting containerized seedlings or taller saplings may increase white
pine survival. Planting more seedlings per acre may also satiate deer browse, however,
this comes at a higher cost and risk. Last-resort browse species, such as white spruce or
balsam fir (Blouch 1984; Dumont et a1 2005) may be the only mesic conifers that can

successfully reestablish in areas with high deer density.

Seedling browse was lower at stands with lower deer densities (< 15 deer / kmz)

that were farther from lowland conifer (>150 m)and covered by deep snow for a longer
portion of the winter (> 6 weeks). Using these criteria to select planting sites may
improve success. Snow cover will only confer protection to small seedlings, so it may
i also be necessary to plant seedlings on microsites that serve as refugia from deer, such as
slash piles or tip-up mounds (Grisez 1960; Krueger and Peterson 2006), to use plastic tree
shelters (Ward and Stephens 1995), or to plant larger seedlings with bud caps (Saunders
and Puettmann 1999). Seedlings may need protection until they reach 125 cm, past which
point terminal leader browsing by deer decreases (Saunders and Puettmann 1999).
Quality deer management that maintains populations in balance with their habitat

(Frawley 2005) could benefit mesic conifer restoration efforts and increase the abundance

95

of naturally regenerating saplings (Chapter 1). This management approach could strike a
balance between deer and forest resources enjoyed by people in the Western Upper
Peninsula, and improve habitat for other species such as birds that associate with conifer
species and require vertical stand structure provided by seedlings and saplings. The
percentage of severely browsed hemlock seedlings rose rapidly with estimated deer
densities, suggesting that if the goal is to restore browse-preferred conifer species to the

landscape, intensive hunting may be necessary to reduce deer populations below 15 deer /

kmz. Hunting may be more effective at reducing deer densities in the southern portion of

the Western Upper Peninsula if permitted during the winter when deer yard together.
2.5 Conclusion

Mesic conifers are important sources of biodiversity in northern hardwood forests,
providing food and habitat structure for birds, especially warbler species, and white-tailed
deer (Kendeigh 1945; Blouch 1984; Green 1992; Patmos 1995). Intense logging and poor
regeneration have reduced the representation of mesic conifers across the Great Lakes
Region. The acreage of hemlock-dominated mixed forests and white pine-mixed forests
are estimated to have declined over eighty five percent from circa 1800 to 2000 in the
Western Upper Peninsula of Michigan (MDNR 2006). Unfortunately, restoration of white
pine and hemlock in areas with high deer density seems unlikely.

Browse pressure on seedlings in northern hardwoods stands was affected by deer
density, local landscape context (i.e. distance to lowland conifer) and snow depth. The
percentage of hemlock seedlings browsed rose rapidly with stand deer density,
highlighting the importance of maintaining deer densities below a threshold level if the

goal is to restore browse preferred species. The relationship between browse pressure and

96

estimated deer density provides evidence that the pellet-count method can be used to
predict relative browse pressure based off relative estimates of deer density. Variability
was observed in the relationships of the pellet count estimation techniques (deer density
vs deer gap—use) with browse pressure, distance to lowland conifer and deep snow weeks,
indicating that the ability of the pellet count method to estimate deer densities, or relative
deer densities, is sensitive to sample method.

Deer herbivory can eliminate browse-preferred species from seedling and sapling
layers established by natural regeneration or planting. An understanding of factors that
affect deer distributions and browse habits are important for developing methods to
restore browse-preferred species to landscapes with high deer densities. When studying
the effects of deer on seedling and sapling layers, using fecal pellet counts in conjunction
with observations of browse on planted or naturally regenerating seedlings may improve

the inferences made.

97

  
 

 

02.55 10 15

  
  

   

20
Kilometers

v.
‘r

in: '1,-
,- 2‘ “
5’11 11 i"

O 0.1 ~‘
0 Q

  

‘Legend

Deer / km
Winter 2007-2008

e 4.1-10.0
. 10.0-200
. 20.0-300

. 30.0-40.0 3,. 1 711 i-
I. 40.0 - 64.5

  

     

Ml

1.1/113311

 

E Counties WI N
Weeks w/ snow
depth > 40 cm MI
ON

Kilometers

 

 

 

55110 220 330 440
an

 

 

Location of stands and estimated winter 2007-2008 deer density overlaid
the number of weeks from November 2007 to April 2008 predicted to
have snow depth > 40 cm. MI=Michigan, WI=Wisconsin. USA.
ON=Ontario. Canada.

Figure 2.1.

98

.30m

:85: 850:: a E 8:: 5 ov a 55:5 :50: 80:0: :80: :00: :0 8:55: 0:3 5:00 50% 8:0: n00: 8:500 :5230— 0: 00:08:: "015 5:0 ow A 5:0:

N

30:: :23 noon :54: 3 Son 8:50.62 50.: 8:003 :0 89:55 8:003 30:: 80:1113mo ”83-38 :85? Anne. \ :83 53:0: :00: :53: 88:58 u :000

 

.v 3:02 5 :585%0: 83 :m 0580: m

 

.082 a 0822.: so as am: .3: 03:88 as 01:: 38.81899 €23 :8 2 so aowmama- m m
6:: DO: :0 8:886:00 506:: :50: a 55:3 015 :0 800:: 0.: 8:295 : u :5:
Am :8 :5 80.30:: 85:88 :0 88:09.0: :0 9:003 30:0 :00: :58. :5 :55: O 0.5 :8 00: v
a 5 :3 :0 500 5:30 0: 00:5,: :0 05- a :00 0 :0: :5 0 8:0 0: 0:880 :85- :8:
a. :. a _ an :0 3 Fe 3: am e :0
Am :8 :5 :003015 mw:::00m 1:0 085020: 00:
:0 3003 305 :00: :58 :8 :53: a 55:3 051%» :00: 1:0 800:: 0:: 8:29am 3mm 88:80 5:05:53...“ m
:5 :00305 mw:::00m 0 .t. o :1 2a:
:0 085008: :0 £003 305: :00: :5: £08: :00: :58 1:0 809:: 0.: 8:0_axm 3mm :8 :009 0:0: m Adm .N + v2 .
38505805 0.519% A? +::Xv_1:w0_ 0::
5:: 8050.880 :0 83:0 0:: :0: :5 8038880 9% :0 5:5 8:208: 0:03 3:0 a» 0:0: E576 _ E: £85m}?
£085 :5: 088:0m .E 0:05 1:0 80.93500 :0: 0::08: 8 :05 8:05 :5 Z
.. QNV 833:8: 53 835:?
5:08—91: _0>0T::Sm _0>0_ 1:00 8:02 8805008 0:02

 

sq: 0: E0085: _0>0:1::::m 80:5: 0:: :8 ma: 0: 88020: 85:88 m0

:0::::: 0.: AS m: :8305 85:80.: :0 85:5: 0:1: .woom E .388 8 .: new 5 :00305 mw:::000. 0.00:8: :0 05: 0:23

1:0 33 085080: :0 83:93 38:12:80 :5 190m :0 300:8 05 0:038 0: :00: 20:08 868:8: 0:050: E8532 AN 050:.

99

Table 2.2. Summary of gap- and stand-level variables: Average (Ave), median (Med),
standard deviation (Stdev), minimum (Min), maximum (Max) and number
of gap or stand-level observations (N obs).

 

 

Ave Med Stdev Min Max N obs

Gap-level

DLC (m) 245.2 187.3 187.5 63.24 886.0 126
Stand—mean centered gap-DLC 0.0 3.7 52.7 -l61.3 156.4 126
PGC (pellets groups / gap transect) 1.5 1 2.1 O 16 338
Stand-mean centered gap-PGC 0.0 -O.2 1.7 -4.1 10.9 338
Stand-level

Deer density (deer/kmz) 22.1 17.7 14.0 4.1 64.5 34
DSW (weeks) 6.6 8.3 4.9 O 12.6 34

 

DLC= distance to lowland conifer in meters, PGC= pellet group count (the number of deer fecal pellets

2
found within a 40 m area in a planted harvest gap), DSW=deep snow weeks (number of weeks fi'om
November 2007 to April 2008 with snow depth > 40 cm).

Table 2.3. Kendall’s tau rank correlation between gap- and stand-level variables.
Comparisons between gap-level and stand-level variables use stand-
averages for non-centered gap-level variables to comply with
independence assumptions. UTM Easting and Northing values were taken
from the center of the deer pellet transects.

 

 

Northing Easting DSW (13:5; Gap-PGC
Gap-DLC 0.31 -0.51** 0.12 -0.28 0.06§
Gap-PGC -0.16 0.04 -0.32" 0.16
Deer density -0.09 O. 10 -0.07
DSW 0.56" -0.18
Easting -O.40**

 

**p-value < 0.05, *p-value < 0.10

§stand-mean centered gap-PGC and stand-mean gap-DLC are significantly positively correlated
(Kendall’s t=0.24, p-value < 0.001).

Gap-DLC= distance to lowland conifer from gap centers, DSW=deep snow weeks (number of
weeks from November 2007 to April 2008 with snow depth > 40 cm), gap-PGC= pellet group

count (the number of deer fecal pellets found within a 40 m area in a planted harvest gap).

100

Figure 2.2.

 

 

 

 

o .4 o
(D 4
“.- .
Q o o
8 . °
1, 2 ~ . ,
e o
3 o . .° '1’. I
g o o .0.0: ..000 0’ o e
2 0" . o. furiu': V3. .
8 o . .~ “:ze’o o o o .
g . . . o o ' o .
u .2 -4 O O O
C
5 o
(I)
'4 r .l I l I l l
-150 -100 -50 0 50 100 150
Stand-mean centered gap-DLC (m)
I =o.241** (< 0.001)
8 - o
6 +
O
(9
n.- o
3 4
o — 00 o o
2 - o o .000 O o
0 - “-0... o o o a. no.
200 400 600 800
Gap-DLC (m)

r = 0.057 (0.401)

Correlation between stand-mean centered gap distance to lowland conifer
(DLC) and stand-mean centered gap pellet group counts (PGC). Lower
image shows the correlation between the non-centered estimates. Non-
parametric Kendall’s tau rank correlations (t (p-value)) are presented.

101

 

 

 

Table 2.4. Average height (stdev), number seedlings found (N.found), percent found
(of 2040 seedlings / species planted), browsed and severely browsed
(browse category 3 or 4 for hemlock only) for planted seedlings across all
gaps and stands in spring 2008.

Species Average N. found % found % browsed % severely

height (cm) browsed

Hemlock 34.3 (7.1) 1982 97.2 92.4 60.3

White pine 26.7 (7.0) 2014 98.7 29.5

White spruce 20.8 (4.9) 2018 98.9 1.8

Table 2.5. Number of seedlings found (N.found) and percent found, browsed,

severely browsed (hemlock only), desiccated and percentage of desiccated
seedlings showing evidence of browse for hemlock at the same 11 stands
and for pine at the same 19 stands visited in both spring 2008 and 2009.
Hemlock seedlings were only resurveyed in 2009 at stands where more
than 50% of the seedlings had not been severely browsed in 2008.

 

 

2008 2009
Hemlock
N. found 652 609
% found 98.8 92.3
% browsed 81.7 95.6
% severely browsed 28.7 75.9
% desiccated 0 9.2
% desiccated browsed 85.7
White pine
N. found 1126 1015
% found 98.8 89.0
% browsed 19.0 51.6
% desiccated O 20.6
% desiccated browsed 30.2

 

102

 

 

 

 

 

 

 

1.0 —
0.8 ~
'0
f.‘
.‘3
In
E 0.6 —
3
2
a
(D
a
g 0.4 —
'u
o
@
tn
.\°
0.2 —
0_0 __ A
1 T F
White spruce White pine Eastern hemlock
Species

Figure 2.3.

Percentage of seedlings browsed by species / stand based on browsed vs
not-browsed categorization. Distributions were all significantly different
from each other based on Wilcoxon signed-rank test with bonferroni
adjusted p-value. Violin plots show distribution of values and median
(white dot), interquartile range (thick vertical line) and range of values 1.5
* IQR (thin vertical line).

103

 

 

 

 

 

 

 

 

 

 

 

 

. . o 9 .0 o o
O O. . . 0 _ O
13 801 . .° 0
c . ° 0 o . o
5 o. O O O
no $ 0 o o o
‘ 60‘ o 1 :
8 ° . .
5 . ° 0 .° .
2 40d 0 . ~ ° .
.0
2:
0
'5 201 . ‘ J o
> o O
a . . . °
a O O
O 0 ‘r I I I I I r I r I I T
g 10 20 30 40 50 60 2 0 1 2 3 4 5
§ Winter Deer Density (deer! km ) Stand Ave. Gap-PGC (groupslgap)
3 r = 0216* (0.073) t = 0.168 (0.167)
0
2 e
g o o o o o. ”o
c 80, . . . . o 2008
e\° o A. A , A 2009
2 . ° ° . . 3° A
O
E 60‘ o e ‘ O
3 o
0 ° 0 o O o
40- , . -
A
20- ,. -
0
O
0 . .
0 ° °
0 2 4 6 a 10 12 200 300 400 500 63) 700 800
DSW (weeks) Stand Ave. Gap-DLC (m)
r = -0.300** (0.013) «2008) = -0.564** (0.007)

Figure 2.4. Observed percentages of hemlock seedlings severely browsed / stand in
2008 (and 2009 in the distance to lowland conifer (DLC) plot only) versus
estimated winter deer density, stand average pellet group count (PGC),
deep snow weeks (DSW), and stand average gap-DLC. Non-parametric
Kendall’s tau rank correlations (1: (p-value)) are presented. I was not
calculated for 2009 observed percentage browsed versus DLC due to low
sample size (n=4).

**p-value < 0.05, * p-value < 0.10

104

 

 

 

 

 

 

 

 

e . ° .
801 ‘
O. O O
E 60- .. 1 . .
g . 9 ° 0
'8 40‘ 0 fl . 0
g °: 0 e 0 °
I. ' .0 0 . 0 .
D 20__4 . O H q . O O o O
a 0 ° 0 . e e o . . e 0 . ° ° 0
.5 ° 0 . e e 00.
a 0- O . 4 . O
0 r I I r fi i I i r l I 1
8 10 20 30 40 50 60 2 0 1 2 3 4 5
g Winter Deer Density (deerl km ) Stand AVG- Gap-PGC (groups/gap)
9- t = -0.016 (0.894) t = 0.248“ (0.042)
.3.
i . ~
.\° ° A
‘3 601 A40 .
g A A 2008
g . , - A A 2009
8 60* . . .. 0 A
22. e
40‘ e 3 0
°. 0 . e A
w ° 0 e. . o - .0 A
20 0‘ ° 0 0 ° 9 o
O O . t . 0 .
fl . e H e.
O I I I fl I I F I I I I I I I I
0 2 4 6 8 10 12 200 300 400 500 600 700 800
- DSW (weeks) Stand Ave. Gap-DLC (m)
= -0.230* (0.059) 1(2008) =-0.462** (0.030)

“[(2009) =-0.515** (0.021)

Figure 2.5. Observed percentages of pine seedlings browsed / stand in 2008 (and 2009
in the distance to lowland conifer (DLC) plot only) versus estimated
winter deer density, stand average pellet group count (PGC), deep snow
weeks (DSW), and stand average gap-DLC. Non-parametric Kendall’s tau
rank correlations (r (p-value)) are presented.

**p-value < 0.05, * p-value < 0.10

105

Table 2.6.

Multilevel logistic regression models predicting the gap-level percentage

of hemlock seedlings severely browsed or a white pine seedling browsed
in 2008: mean parameter estimates and 95% credible interval from
posterior distributions. Deer and DSW were grand mean centered and
pellets were group (i.e. stand) mean centered. Estimates are on the logistic
scale, so a change of one unit in the independent variables corresponds to
a change of at most B/4 on the probability scale. A decrease in DIC of 5 or
more indicates improved performance for full model over the null model.

 

 

 

Hemlock White pine

Mean 95% C1 DIC Mean 95% CI DIC
Model 1: Null model 1159.84 1018.37
Intercept 0.49 -0.02 — 0.99 -1.09 -1.58 — -0.62
Stand Stdev 1.43 1.09 — 1.89 1.35 1.02 — 1.79
Model 2: DSW + Deer 1159.38 1017.76
intercept 0.49“ 0.03 — 0.95 -1.10* -1.56 - -0.64
Deer 0.24 -0.08 —- 0.58 -0.13 -0.47 — 0.21
DSW -0.12* -0.22 — -0.03 -O.10"‘ -0.20 — -0.01
Stand Stdev 1.27 0.96 — 1.70 1.30 0.97 — 1.74
Model 3: DSW + PGC 1133.83 1014.59
Intercept 0.50“ 0.04 — 0.96 -1.10* -l.56 - -0.64
Pellets 0.18* 0.11 — 0.26 0.08* 0.01 — 0.14
DSW -0.14* -0.23 — -0.04 -0.10* -0.20 — -0.003
Stand Stdev 1.31 0.99 — 1.73 1.29 0.97 — 1.71
N gaps 338 338
N stands 34 34

 

Deer= estimated winter deer density in 10’s of deer / kmz, DSW=deep snow weeks (number of weeks
from November 2007 to April 2008 with snow depth > 40 cm), PGC= pellet group count (the number

2
of deer fecal pellets found within a 40 m area in a planted harvest gap).

106

Model 3 Model 1 (Null)
Eastern hemlock

 

 

 

 

 

 

 

 

 

 

   

 

 

100 7
8' 75 - I. I ~
a ..
;. l
g 50 - i J
g 25 2 ( 1 i _
°' I
1
o 1 . 2
0 25 50 75 100 0 . 25 50 75 100
Observed % seedlings severely browsed I gap
Eastern white pine
100 5 2

 

 

 

 

Predicted % I gap
0'1
0

25—

 

 

A
A
I vv A
v

 

 

       

 

 

 

 

 

0 25 50 75 100 0 25 50 75 100
Observed % seedlings browsed l gap

Figure 2.6. Observed percentage of hemlock seedlings severely browsed and white
pine seedlings browsed / gap in 2008 versus mean posterior predicted
percentages per gap and 95% credible intervals from logistic mixed
models. Model 3 includes deep snow weeks and stand-centered gap-level
pellet group counts as covariates and model 1 (null model) includes only
an overall intercept and stand-level random effect intercept. Observed
percentage on the x-axis has been jittered for visualization (random jitter
is equivalent for Model 3 and Model 1 plots within species).

107

Table 2.7. Multilevel logistic regression models predicting the probability that a
hemlock seedling was severely browsed or a white pine seedling was
browsed in 2008 for thirteen stands where distance to lowland conifer was
determined: mean parameter estimates and 95% credible interval from
posterior distribution. DSW was grand mean centered and pellets and
DLC were group (i.e. stand) mean centered. Estimates are on the logistic
scale, so a change of one unit in the independent variables corresponds to
a change of at most B/4 on the probability scale. A decrease in DIC of 5 or
more indicates improved performance for full model over the null model.

 

 

 

Hemlock White pine
Mean 95% C1 DIC Mean 95% C1 DIC

Model 1: Null model 397.18 325.29
Intercept 0.37 -0.71 — 1.42 -1.87 -2.60 - -l.19
Stand Stdev 1.84 1.15 — 2.97 1.16 0.69 — 1.95
Model 4: DSW + PGC + DLC 374.76 326.83
Intercept 0.36 -0.70 — 1.42 -1.92* -2.71 — -1.18
DLC 0.02 -0.02 - 0.05 0.02 -0.03 — 0.07
Pellets 0.42“ 0.22 — 0.63 0.06 -0.1 l - 0.22
DSW -0.18 -0.43 — 0.06 -0.05 -0.23 -— 0.12
Stand Stdev 1.81 1.11 — 2.99 1.23 0.72 — 2.12
Model 5: DLC 392.10 325.32
Intercept 0.37 -0.73 - 1.47 -1.88* -2.63 — -1.18
DLC 0.05“ 0.01 — 0.08 0.03 -0.01 — 0.07
Stand Stdev 1.87 1.16 — 3.04 1.17 0.69 — 1.96
N gaps 126 126
N stands 13 13

 

DSW=deep snow weeks (number of weeks from November 2007 to April 2008 with snow depth > 40

cm), PGC= pellet group count (the number of deer fecal pellets found within a 40 m area in a planted
harvest gap), DLC= distance to lowland conifer in 10s of m.

108

Table 2.8. Mean predicted probability (Pr.) (95% prediction quantiles) of severe
browse for hemlock seedlings and probability of browse for white pine
seedlings under different gap pellet count (stand centered) and deep snow
week conditions. Predictions are for an unmeasured stand based on one-to-
one simulations using the posterior distribution for Model 3 parameters.

 

 

Gap pellet count (stand centered) Deep snow weeks Pr. offizzlelgectrowse HQEfiEEXSC
At stand mean 6.6 weeks (ave obs) 59% (IO-96%) 30% (2-81%)
At stand mean 7.6 weeks (+1 week) 57% (9-95%) 29% (2-80%)
At stand mean 12.6 weeks (max obs) 44% (4-92%) 22% (1-72%)
+1 from stand mean 6.6 weeks (ave obs) 63% (12-97%) 31% (3-83%)
+11 from stand mean (max obs) 6.6 weeks (ave obs) 87% (44-99%) 45% (5-92%)

 

109

APPENDICES

110

APPENDIX A: Descriptions of Study Stands

 

 

Table A. 1. Description of 59 study stands including average snow depth November
2007 - April 2008 (cm), county, ownership, Habitat Type, relative
densities (RD) of the three most abundant overstory tree species and
average stand diameter at breast height (DBH) of mature trees with dbh
>20 cm, and year of harvest entry (YOE).

$1an 5;:Xa County b Owershipc “161:3? RDd’e (81:12:: (1 YOE
4 16.3 MN PCTC ATD-Hp BA 2,1,2 93:; 2% (36% 2005
20 18.7 D1 5. MI ATD-Hp 13.485132. 92‘: 2% 357.528) 2004
27 33.9 MA 5. MI ATD 13c 1802423111190 (23%) 2006
34 31.1 [R KLA ATD-Hp BA 1841:1851: f, 1% 32'9”; 2002
38 3 I .5 1R GMO ATM 17133;: 9:523 2% (269.46) 2006
803 23.3 MN 3. MI ATD-Hp 31,31): 93/; 1% (2:25) 2000
806 20.4 DI GMO AOCa BA 383: i024. (2529) 1999
807 22.6 MN GMO ATD-Hp BA 3;: 617:3 <1% (269;) 2001
808 22.2 MN GMO ATD-Hp SM 233251;? 20% (2:56) 2000
814 19.2 D1 s. M1 AOCa BA 5&733 1% (371.69) 1999
815 33.5 D1 3. MI ATM PB 18022.8ng 4% (38915) 1995
816 22.8 DI 3. M1 AOCa BA 3.573;: {ii/(i) 14% (279118) 2001
817 21.3 DI GMO AOCa BA 13:33:; 4% (379:) 2001
818 23.1 MA PCTC AOCa BA 33>: 9:33“ <1 % (3:00) 2002
819 21.2 MA PCTC AOCa SM 100% (259.59) 1999
820 48.0 MA GMO ATD 1738913282344. (257.59) 2001
821 40.5 MA GMO ATM v13 S7122 83? 4% (2:65) 2002
822 43.9 MA GMO ATM BF 132/'3 53:109.. (2572) 2002
823 40.8 MA GMO ATM BF 15:; 7r: 6% (22°79) 2003
824 31.1 MA GMO TMC YB 15;??? 4% (392;) 1999
825 22.6 MA PCTC ATD 171385940 93; 3% (279.53) 2003

 

111

Table A. 1 . continued

 

 

$3an 32;; County Owership “£511? RD (33:) YOE
826 22.6 MA PCTC ATM 53,933? (256.45) 1999
829 24.5 MA PCTC ATD BC :22 8:?6% (3:25) 1998
832 22.4 MA PCTC ATM BA it: 8:: 4% (35958) 2003
834 37.0 MA s. M1 ATD AR 15:23:13? 1% (39917) 1993
837 35.2 [R 5. MI AOCa AR 1394;??? 8% (27955) 1998
842 28.8 MA s. M1 ATD AR 280112714131 1% 83:11) 1997
844 25.7 MA PCTC ATD “322 913000 4% (2523) 2001
847 32.2 MA s. M1 ATD BC :22 732’ 8% (fig) 2002
851 24.6 DI s. MI ATD SM 5&63‘; 9% (372.01) 2000
852 22.6 D1 KLA ATD-Hp BA 194047121; 3% (389:) 2001
1042 31.9 IR 5. M1 TMC AR 3:": 31332204. (399:) 1996
1043 31.9 IR 5. M1 AOCa “331; 9;?m (369:) 1998
1044 30.4 1R 5. M1 AOCa BA :3: iii/21% (381:) 1995
1081 24.5 1R 5. MI ATM AR 1321:;73 5% (289.06) 1996
1133 22.2 D1 KLA AOCa BA figs/(7% 8% (2628) 1999
1272 17.5 D1 GMO ATD-Hp BA 3831125351 8% (3:26) 1998
1603 37.2 IR s. M1 AOCa AR 2391143318.), (31‘ 33) 2002
4024 19.4 MN PCTC ATD-Hp HM iii/07913.1: 9% 83:3) 2000
4042 20.3 MN PCTC AOCa BA 181;,83/0 4% 83:2) 1999
4052 20.5 MN PCTC ATD-Hp SM 3,2433; 9% (39903) 2002
4082 20.6 D1 3. M1 AOCa SM E3221: 1% (3636) 2001
4102 20.5 D1 GMO AOCa 1333/: 9‘3? 1% (267.55) 1996
4162 23.8 D1 s. M1 AOCa BA 3;: ?f;/°<l% (36355) 1998
4172 16.9 D1 3. MI ATD-Hp 3142;: igfim (372;) 1994
4232 29.5 MA s. MI ATM AR 18113306312704. (389; 1995

 

112

Table A. 1. continued

 

 

Stlasid 32;; County Owership H1533)? RD (33:) YOE
4324 35.1 MA 5. MI ATD R0 361;: 960410304 (22:7?) 2002
4332 36.2 MA s. MI ATM 13031;: 933% (38%?) 1998
4342 29.4 MA 3. M1 ATD “511:: 93:" 1% (362.63) 1998
4422 22.9 MN PCTC ATD-Hp “$85,131:? 6% (277.68) 2000
4423 22.8 MN PCTC ATD-Hp WA 2;: fig<m (2:69) 1997
4429 22.9 MN PCTC AOCa “545511281: 4% (289.11) 1999
4524 22.0 DI KLA ATD-Hp BA 13;,81’: 1% (32:) 1998
5032‘ 15.6 D1 s. MI ATD-Hp BA fflsfiz‘; 1% (3:56) 1997
5052 16.9 DI s. M1 AOCa BA 2;: iiifim (372.31) 1995
5104 30.3 MA s. M1 ATD AR 3): 331m (36'_'7') 1998
5292 34.6 IR s. M1 AOCa “it: 93:: 1% (2:2) 1997
5314 28.4 IR s. MI AOCa BA f£8$g 1% (2583‘; 1995
- o
5342 22.3 D1 53:: ATM BA 315322? 1% (2:65) 2001

 

aSnow Data Assimilation System (SNODAS) from the National Snow and Ice Data Center (N SIDC)
bDI=Dickinson, IR=1ron, MA=Marquette, MN=Menominee

cGMO = GMO Renewable Resources (managed by American Forest Management, Inc); KLA= Keweenaw
Land Association, Ltd.; PCTC= Plum Creek Timber Co., Inc.; S. M1 = State of Michigan

ClCalculated for all trees with dbh >20 cm within 20 m from 2-5 non—overlapping gaps at each staqbd
randomly selected from the list of all possible combinations with the most non-overlapping gaps.

eA13= American beech (Fagus grandifolia Ehrh.), AE= American elm (Ulmus americana L.), AR= red
maple (Acer rubrum L.), BA= American basswood (Tilia americana L.), BC= black cherry (Prunus
serotina Ehrh.), BF= balsam fir (A bies balsamea (L.) Mill.), BL= black ash (F rarinus nigra Marsh), BP=
balsam poplar (Populus balsamifera L.), HM= eastern hemlock (Tsuga Canadensis (L.) Carriére), IW=
ironwood (Ostrya virginiana (Mill.) K. Koch), NC= northern white-cedar (Thuja occidentalis L.), PB=
paper birch (Betula papyrifera Marshall), SM= sugar maple (A cer saccharum Marsh), QA= quaking aspen
(Populus tremuloides Michx.), RO= northern red oak (Quercus rubra L.), WA= white ash (Fraxinus
americana L.), WS= white spruce (Picea glauca (Moench) Voss), YB= yellow birch (Betula
alleghaniensis Britton)

113

Table A2.

Geological characteristics for each stand, including landform type, soil

type, upper layer soil texture and sugar maple site index at age 50 (8150) in
meters estimated for each soil type and soil series.

 

Soil

 

Stlagd Landforma Soil typeb textureb’c S150 Soil seriesb
4 Drumlinized Onaway fine sandy loam, fsal 19.8 Onaway
ground moraine 3-9% slopes
20 Ground moraine Emmet fine sandy loam, fsal 20.1 Emmet
6—18% slope
27 Disintegration Karlin sandy loam, 1-6 % sal NA Keewaydin-
moraine slopes Michigamme-Rock
Outcrop Association
34 Ground moraine Trenary fine sandy-loam, fsal 18.6 Trenary
6-18% slopes, stony
38 Disintegration Sundog very fine sandy vfsal NA Sundog
moraine loam, 6-18% slopes
803 Drumlinized Onaway fine sandy loam, fsal 19.8 Onaway
ground moraine 3-9% slopes
806 Drumlinized Emmet fine sandy loam, fsal 20.1 Emmet
ground moraine 0-6% slope
807 Bedrock-controlled Onaway fine sandy loam, fsal 19.8 Onaway
ground moraine 3-9% slopes
808 Drumlinized Onaway fine sandy loam, fsal 19.8 Onaway
ground moraine 3-9% slopes
814 Ground moraine Emmet fine sandy loam, fsal 20.1 Emmet
6-18% slope
815 Disintegration Pemene fine sandy loam, fsal 18.3 Pemene
moraine 6-18% slope
816 Ground moraine Emmet-Pemene fine fsal 20.1 / Emmet-Pemene
sandy Ioams, 0-6% slope 18.3
817 Ground moraine Emmet fine sandy loam, fsal 20.1 Emmet
6-18% slope
818 Drumlinized Onaway fine sandy loam, fsal 19.8 Cunard-Nahma
ground moraine 6-18% slopes Association
819 Drumlinized Emmet fine sandy loam, fsal 20.1 Kalkaska-Ishpeming-
ground moraine 6-18% slopes Rock Outcrop
Association
820 Bedrock-controlled Keewaydin-Dishno cfsal / cs1 18.6/ Sundog-Minocqua-
ground moraine complex, 6-18% slopes, 18.3 Charming Association
rocky, bouldery
821 Bedrock-controlled Keewaydin-Dishno cfsal / csl 18.6 / Sundog-Minocqua-
ground moraine complex, 6-18% slopes, 18.3 Charming Association
rocky, bouldery
822 Bedrock-controlled Keewaydin- cfsal 18.6 / Kalkaska-Carbondale-
ground moraine Michigamme-Rock 18.3 Deford Association
outcrop complex, 6-25%
slopes, very bouldery
823 Bedrock-controlled Keewaydin- cfsal 18.6 / Kalkaska-Carbondale-
ground moraine Michigamme-Rock 18.3 Deford Association

outcrop complex, 6-25%
slopes, very bouldery

 

114

Table A.2. continued.

 

 

Stlagd Landfomr Soil type Rifle 3150 Soil series
824 Bedrock-controlled Keewaydin—Dishno cfsal / cs1 18.6 / Sundog-Minocqua-
ground moraine complex, 6-18% slopes, 18.3 Charming Association
rocky, bouldery
825 Dissected moraine Garlic-Alcona—Voelker fsa / 1vfsa 18.9 / Kalkaska-Ishpeming—
complex, 15-70% slopes / fsa 18.6/ Rock Outcrop
1 8.6 Association
826 Dissected moraine Garlic-Alcona—Voelker fsa / 1vfsa 18.9 / Kalkaska-lshpeming—
complex, 15—70% slopes / fsa 186/ Rock Outcrop
1 8.6 Association
829 Dissected moraine Garlic-Alcona-Voelker fsa / 1vfsa 18.9 / Kalkaska-Ishpeming—
complex, 8-35% slopes / fsa 186/ Rock Outcrop
. 1 8.6 Association
832 Drumlinized Onaway fine sandy loam, fsal 19.8 Cunard-Nahma
ground moraine 6-l8% slopes Association
834 Disintegration Amasa very fine sandy vfsal 18.6 Cunard-Nahma
moraine loam, 1.6% slopes Association
837 Ground moraine Petticoat-Wabeno silt S] 20.1 / Petticoat-Wabeno
Ioams, 6-8% slopes, very 20.4
stony
842 Till-floored lake Munising-Yalmer fsal / 19.2 / Zeba-Jacobsville
plain complex, 6-18% slopes fsa 18.6 Association
844 Till-floored lake Munising-Yalmer fsal / 19.2 / Zeba-Jacobsville
plain complex, l-6% slopes fsa 18.6 Association
847 Ground moraine Shoepac-Trenary silt S] 19.8 / Shoepac-Carbondale
Ioams, 1-6% slopes 18.6 Association
851 Ground moraine Emmet-Pemene fine fsal 20.1 / Emmet-Pemene
sandy Ioams, 6-18\% 18.3
slope
852 Ground moraine Emmet fine sandy loam, fsal 20.1 Emmet
6-18% slope
1042 Bedrock-controlled Wabeno silt loam, 1-6% 51 20.4 Wabeno
ground moraine slopes, very stony
1043 Bedrock-controlled Wabeno silt loam, 1-6% 51 20.4 Wabeno
ground moraine slopes, very stony
1044 Bedrock-controlled Wabeno silt loam, 1-6% sl 20.4 Wabeno
ground moraine slopes, very stony
1081 Disintegration Goodman-Wabeno- s1/ 21.0 / Goodman-Wabeno-
moraine Sundog sandy substratum sl / 20.4 / Sundog
complex, 6-18% slopes, vfsal NA
stony
1133 Bedrock-controlled Emmet fine sandy loam, fsal 20.1 Emmet
ground moraine 6-18% slope
1272 Drumlinized Emmet fine sandy loam, fsal 20.1 Emmet
ground moraine 0-6% slope
1603 Ground moraine Petticoat-Wabeno silt s1 20.1 / Petticoat-Wabeno
Ioams, l-6% slopes, very 20.4
stony
4024 Drumlinized Onaway fine sandy loam, fsal 19.8 Onaway
ground moraine 3-9% slopes
4042 Drumlinized Onaway fine sandy loam, fsal 19.8 Onaway
ground moraine 3-9% slopes

 

115

Table A.2. continued.

 

Stand

Soil

 

ID Landforrn Soil type texture 3150 Soil series
4052 Drumlinized Onaway fine sandy loam, fsal 19.8 Onaway
ground moraine 12-35% slopes
4082 Bedrock-controlled Emmet-Rock outcrop fsal 20.1 Emmet
ground moraine complex, 6-18% slope
4102 Drumlinized Emmet fine sandy loam, fsal 20.1 Emmet
ground moraine 0-6% slope
4162 Ground moraine Emmet fine sandy loam, fsal 20.1 Emmet
0-6% slope
4172 Ground moraine Emmet fine sandy loam, fsal 20.1 Emmet
0-6% slope
4232 Drumlinized Emmet-Escanaba fsal / lfsa 20.1 / Rubicon-Sayner
ground moraine complex, 1-6% slopes 18.3 Association
4324 Disintegration Keweenaw-Kalkaska lsa/ 18.6/ Rubicon-Sayner
moraine complex, 6-18% slopes sa 19.5 Association
4332 Disintegration Keweenaw-Kalkaska lsa/ 18.6/ Rubicon-Sayner
moraine complex, 6-18% slopes sa 19.5 Association
4342 Disintegration Keweenaw loamy sand, lsa 18.6 Zeba-Jacobsville
moraine 6-18% slopes Association
4422 Drumlinized Onaway fine sandy loam, fsal 19.8 Onaway
ground moraine 3-9% slopes
4423 Drumlinized Onaway fine sandy loam, fsal 19.8 Onaway
ground moraine 3-9% slopes
4429 Drumlinized Onaway fine sandy loam, fsal 19.8 Onaway
ground moraine 3-9% slopes
4524 Ground moraine Longrie fine sandy loam, fsal 18.6 Longrie
0—6% slope
5032 Bedrock-controlled Emmet fine sandy loam, fsal 20.1 Emmet
ground moraine 6-18% slope
5052 Ground moraine Emmet fine sandy loam, fsal 20.1 Emmet
0-6% slope
5104 Disintegration Emmet fine sandy loam, fsal 20.1 Kalkaska—Ishpeming-
moraine 6-18% slopes Rock Outcrop
Association
5292 Bedrock-controlled Petticoat-Wabeno silt S] 20.1 / Petticoat-Wabeno
ground moraine Ioams, 6-8% slopes, very 20.4
stony
5314 Disintegration Stambaugh silt loam, 2- s1 18.6 Stambaugh
moraine 6% slopes, stony
5342 Drumlinized Emmet fine sandy loam, fsal 20.1 Emmet

gound moraine

0-6% slope

 

aJerome, Dwight S. and Upper Peninsula, Michigan Soil Survey Staff. 2006 drafi. Landforms of the Upper
Peninsula of Michigan. USDA NRCS.

b Linsemier, Lyle H. 1997. Soil survey of Iron County, Michigan. USDA NRCS and USPS. Washington,
DC; Schwenner, Chalres. 1989. Soil survey of Menominee County, Michigan. USDA NRCS.
Washington, DC; Schwenner, Chalres. 2007. Soil survey of Marquette County, Michigan. USDA NRCS.
Washington, DC; USDA. 1989. Soil survey of Dickinson County, Michigan. USDA NRCS. Washintgon,
D.C.

cAbbreviations: cfsa1= cobbly fine sandy loam; csl= cobbly silt loam; fsa= fine sand; lfsa= loamy fine sand;
fsa1= fine sandy loam; lsa= loamy sand; lvfsa= loamy very fine sand; sa= sand; sal= sandy loam; sl= silt
loam; vfsal= very fine sandy loam

116

APPENDIX B: Additional Details on Methodology

8.]. Pellet count layout and conversion of pellet counts to deer density
Pellet groups were double counted within ten 50 m by 4 m transects oriented in an

hour glass shape (LeBouton et a1 2005) (Figure B.1).

N
50 m x 4 m A
l] —>
hisseits / 5 m from plot
center to
/fiansect start

50 m T
between—p /
81 m

transects
<-—— between
transects

Figure 81 Layout of pellet count transects.

Average pellet group counts from all ten transects were used to calculate deer
density (equation 1). The pellet deposition rate estimated for the Upper Peninsula is 13.4
times per deer in a 24-h period (Hill 2001) and the period of deposition was the number

of days from November I (assumed date of leaf-off) to the date of the pellet surveys in

April and May (180 to 200 days). A constant was required to convert to deer / kmz.

average pellet group count x 5000

 

Estimated deer density / km2 = (1)

period of deposition x13.4

117

 

B.2. Predicting diameter at breast height from stump basal diameter and
measurement height

Stump basal diameter was converted to diameter at breast height using the
equation below developed by Demaerschalk and Omule 1978 for merchantable tree
species in British Columbia:

DBH = StDiam + b x StDiam x In W

StDiam is the stump diameter outside bark in centimeters, Sth is the ht of the diameter
measurement in meters and b is a species specific constant. Values for the constant “b”
from the southern central forest inventory zones in British Columbia were instead of the
higher latitude or coastal regions because this region is more similar to the Western
Upper Peninsula in terms of climatic conditions. Constant values were available for
balsam fir, quaking aspen, birch species and spruce species. Values for western hemlock
(Tsuga heterophylla) were used for eastern hemlock (Tsuga canadensis), western red-
cedar (Thuja plicata) for northern white-cedar (Thuja occidentalis), and bigleaf maple

(A cer macrophyllum) for sugar maple and all other hardwood species (Table B. 1).

Table B. 1. Species-specific parameter values for “b” in equation estimating dbh from
stump diameter (equation 2) and standard error of estimated dbh
(Demaerschalk and Omule 1978). “Species” refers to the species from this
study and “Equivalent species” refers to the species from Demaerschalk
and Omule 1978.

 

 

Species Equivalent species b Standard erroiicm)
Balsam fir 0.301495 2.26
Birch species 0.297892 2.07
Spruce species 0.415792 3.65
Quaking aspen 0.281240 1.68
Northern white-cedar Western red-cedar 0.486935 4.53
Eastern hemlock Western hemlock 0.264130 3.02
Other hardwoods Bigleaf maple 0.285892 0.91

 

118

 

B.3. FVS diameter growth simulations to “grow back” diameter of mature seed
trees to their diameter at the time of harvest

I wanted to predict seed production potential at the time of harvest because some
saplings present in harvest gaps may have come from seeds produced by trees that were
removed by the harvest. USDA Forest Vegetation Simulator Lake States Variant (Bush
and Brand 1993) diameter growth equations were implemented in NetLogo (Wilensky

1999) to “grow back” the diameter at breast height (dbh) of living trees to dbh at harvest

(dbhHarvest). Observed stand conditions were simulated and annual diameter growth for

each measured tree of a given dbh and species were estimated. For each time step: I)
predicted diameter growth from the current year was subtracted from the diameter, 2)
stand basal area and quadratic mean diameter were re-calibrated, and 3) annual diameter
growth was recalculated. This process was reiterated for the number of years since
harvest.

Species-specific FVS diameter growth equations require Site Index (obtained from
soil survey data) and stand basal area (BA) and quadratic mean diameter (QMD) of trees
with dbh > 2.54 cm. However, I only measured trees with dbh >= 20 cm. In order to
predict the basal area of trees with dbh < 20 cm (BAu) and the quadratic mean diameter
of trees with dbh > 2.54 cm (QMDt) in Netlogo, linear regression models were developed
with covariates that could produced from my data set (i.e. variables involving trees with
dbh > 20 cm). Models were calibrated with stand BA and dbh data from 55 managed
northern hardwood stands in the study area collected using prism point-sampling
(Millington et al 2009 in review). Twelve of these stands coincide with those sampled for

gap regeneration.

119

Variables included in regression models were BA (BAo), quadratic mean diameter
(QMDo), the percentage of trees per hectare not sugar maple (PerNotSM) and the species
richness (SR) for trees with dbh > 20 cm. Time since harvest (TSHarvest) was also
included, but only known for 23 stands. AIC values were compared in order to determine
the model that best predicted BAu and QMDt (i.e. lower AIC values by 2 indicate better
model performance). BAu and species richness (SR) were log transformed to comply
with normality assumptions. Two stands where no saplings < 20 cm were sampled were
removed because they were outliers. The best model for both variables (i.e. lowest AIC
of models tested) included BAo and PerNotSM (Tables B2 and B3). TSHarvest and SR
were not significant predictors.

I decided that incorporating stochasticity into the BAu and QMDt predictions was

unnecessary. This was determined based on twenty 15 year simulations for one randomly

selected stand for which the average standard deviation of the predicted dbhHarvest for all

183 trees at this stand was monitored. For each simulation, [El-values for the BAu and
QMDt regression models were randomly selected from a uniform distribution bounded
by the B point-estimate plus or minus one standard error. The average standard deviation

stabilized after 10 simulations at around 0052-0059 cm. Diameter estimates in the field

were measured to the nearest tenth of a centimeter, so a difference of ~ 1 cm in dbhHawest

estimates was insignificant.

Sensitivity analysis were also run to determine the effect of BAu and QMDt on one
year diameter growth estimates for trees with the average dbh for sugar maple (30.1 cm),
basswood (34.0 cm) and ironwood (24.0 cm). BAo and site index equal were set to the

average from the 59 sites (19.5 mz/ha and 19.7 m respectively). For baseline diameter

120

Table B.2.

Parameter estimates ([3), standard errors (SE), AIC and adjusted R2 values
from linear regression models to predict log transformed basal area (BAu)

 

 

 

 

 

 

of trees with dbh < 20 cm.
Log(BAu)
Model predictors [3 (SE) AIC Adj R2
Models using all stand (n=53)
Intercept 3.634 (0.291)**"‘ 100.3 0.484
BAo -0.017 (0.004)***
PerNotSM 1.360 (0.327)***
Intercept 4.237 (0.291)*** 114.0 0.319
BAo -0.021 (0.004)***
Intercept 4.053 (0.347)**"' 115.1 0.319
BAo -0.021 (0.004)**"'
Log(SR) 0.155 (0.160) NS
Intercept (null) 2.851 (0.114)*** 133.5
Models using stands with known time since harvest (n=25)
Intercept 3.511 (0.474)*** 45.8 0.409
BAo -0.017 (0.006)"I
PerNotAcesac 1.229 (0.643) NS
Intercept 3.686 (0.433)*** 46.0 0.403
BAD -0.022 (0.006)"
TSHarvest 0.067 (0.036) NS
Intercept 4.022 (0.415)*** 47.6 0.334
BAo -0.022 (0.006)"
Intercept (null) 2.655 (0.160)*** 56.1
Intercept 2.328 (O.312)*** . 56.5 0.021
TSHarvest 0.057 (0.047) NS

 

stup<0000| **p<0,01*p<0.05 NSp>0.05

BAo= basal area of trees with dbh > 20 cm, PerNotSM= percentage of trees with dbh > 20
not sugar maple, SR= species richness of trees with dbh > 20 cm, TSHarvest=time since

harvest.

growth estimates, BAu was set to 2.62 mz/ha and QMDt to 22.8 cm (the estimates

predicted from the average 8A0 and PerNot_SM (17.8%)). BAu was varied from the

average estimate to the average estimate plus the first quartile residual from the

121

Table B.3. Parameter estimates ([3), standard errors (SE), AIC and adjusted R2 values
from linear regression models to predict quadratic mean diameter (QMDt)

 

 

 

 

 

 

of trees with dbh > 2.54 cm.

QMDt

Model predictors (3 (s13) AIC Adj R2
Models using all stand (n=53)
Intercept ' 5.296 (0.948)*** 225.3 0.363
BAo 0.050 (0.012)***
PerNotSM -2.908 (1.064)"
Intercept 4.008 (0.873)*** 230.7 0.282
BAo 0.058 (0.013)*“
Intercept 4.753 (1.035)*** 230.9 0.293
BAo 0.057 (0.013)***
Log(SR) -0.628 (0.477) NS
Intercept -0.163 (0.962) NS 243.9 0.080
QMDo 0.683 (0.290)*
Intercept (null) 7.837 (O.333)*** 247.3
Models using stands with known time since harvest (n=25)

Intercept 6.156 (1.442)*** 97.0 0.383
BAo 0.042 (0.019)*
PerNotSM -4.327 (1.959)*
Intercept 4.357 (1.295)" 100.0 0.269
BAo 0.059 (0.019)M
Intercept 4.153 (1.459)" 101.9 0.237
BAo 0.058 (0.020)"
TSHarvest 0.041 (0.123)
Intercept (null) 8.066 (0.477)*** 106.3
Intercept 3.930 (6.500) NS 107.9 -0.028
QMDo 0.342 (0.535) NS
Intercept 7.670 (0.958)*** 108.1 -0.036
TSHarvest 0.068 (0.143) NS

 

*“p < 0.001, "p < 0.01, *p < 0.05, NS p > 0.05

BAo= basal area of trees with dbh > 20 cm, PerNotSM= percentage of trees with dbh > 20
not sugar maple, SR= species richness of trees with dbh > 20 cm, QMDo=quadratic mean
diameter of trees with dbh > 20 cm, TSHarvest=time since harvest.

122

predictive equations (-1.02 mz/ha), plus the third quartile residual (2.38 mz/ha) and plus

the maximum observed residual (9.33 mZ/ha). Subtracting the third quartile residual and

maximum residual produced unrealistic negative estimates of BAu. QMDt was varied
from the average estimate to the average plus the first quartile residual (-2.5 cm), and
plus or minus the third quartile (3.6 cm) and maximum residual (9.3 cm).

The sensitivity index (the percent change in diameter growth divided by the

percent change in BAu or QMDt) was on average 0.055 for BAu and 0.279 for QMDt,

 

well under the 1.5 threshold indicating significant sensitivity to a variable (Table B4). A
356% increase in BAu (adding the maximum residual) only decreased diameter grth
by 17.8 - 19.6 % and a 41% increase in QMDt (adding the maximum residual) only
decreased diameter growth by 3.8 - 21.1% for the three species. The maximum
differences in growth rate was a 0.36 cm decrease in annual sugar maple growth rate with
a 41% increase in QMDt and a 0.30 cm decrease in annual basswood growth rate with a

356% increase in BAu (adding the maximum residual). Over fifteen years, this would

translate into at most a 5.3 and 4.6 cm bias in dbhHawest. However, this represents a worst
case scenario in which the stand-level BAu and QMDt are most poorly predicted. Fifty
percent of the residuals (the interquartile range) lie within a much narrower range for
BAu (-1.02 and +2.38 m2 / ha) and QMD (-2.45 and +3.55 cm), and over- or under-

predictions within this range only results in a maximum change of 0.11 cm diameter

growth/yr, or a biased dbhHaWest estimation of 1.65 cm over 15 years. This error is

acceptable given the error already inherent to F VS growth estimates (Canavan and Ramm

2000; Pokharel and Froese 2008).

123

Table 8.4.

Sensitivity analysis: percent change in diameter growth (DG) with given
percent change in BAu or QMDt and sensitivity index (absolute percent
change in diameter grth divided by absolute percent change in BAu or

 

 

 

QMDt).
Species % Change DG % Change BAu Sensitivity Index
7 2
BAu average (2.62 m lha) + max residual (9.33 m lha)
Basswood -17.77 356.13 0.050
Sugar maple -19.64 356.13 0.055
Ironwood -18.97 356.13 0.053
2 2
BAu average (2.62 111 /ha) + 3rd quartile residual (2.38 m /ha)
Basswood -4.62 90.65 0.051
Sugar maple -5.24 90.65 0.058
Ironwood -S. l 8 90.65 0.057
2 2
BAu average (2.62 m /ha) + 1st quartile residual (-1.02 m /ha)
Basswood 2.03 -39.07 0.052
Sugar maple 2.34 -39.07 0.060
Ironwood 2.33 -39.07 0.060
Average Sensitivity Index BAu 0.055

 

% Change DG % Change QMDt

Sensitivity Index

 

 

Species
QMDt average (22.8 cm) + max residual (9.3 cm)
Basswood -21.1 1 40.90 0.516
Sugar maple -3.78 40.90 0.092
Ironwood -15.02 40.90 0.367
QMDt average (22.8 cm) - max residual (9.3 cm)
Basswood 3.74 -40.9 0.091
Sugar maple -1 1.59 -40.9 0.283
Ironwood 15.85 -40.9 0.387
QMDt average (22.8 cm) + 3rd quartile residual (3.6 cm)
Basswood -6.37 15.57 0.409
Sugar maple 0.45 15.57 0.029
Ironwood -6.01 15.57 0.386
QMDt average (22.8 cm) - 3rd quartile residual (3.6 cm)
Basswood 3.53 -15.57 0.226
Sugar maple -2.80 -15.57 0.180
Ironwood 6.22 -15.57 0.400
QMDt average (22.8 cm) + 1st quartile residual (-2.5 cm)
Basswood 2.74 -10.74 0.255
Sugar maple -l .70 -10.74 0.158
Ironwood 4.29 -10.74 0.399
0.279

 

Average Sensitivity Index QMDt

124

B.4. Determining the sensitivity of the random stand-level effect standard deviation
to prior specification

Multilevel model estimates with Bayesian inference can be sensitive to the prior
specification on the random-effect variance component because of its impact on fixed-
effect coefficients and random effect estimates (N atarajan and Kass 2000; Spiegelhalter
et al 2004; Browne and Draper 2006). To test the robustness of the posterior distribution
of the random effect variance to different priors, I ran the full generalized linear
multilevel model for sugar maple 1-2 m sapling abundance (Chapter 1) and the Model 4
logistic multilevel model for the gap-level percentage of hemlock seedlings browsed

(Chapter 2) with three commonly recommended priors: inverse gamma for the precision
(1/ 02) (Browne and Draper 2006), half-normal for the standard deviation (0)
(Spiegelhalter et a1 2004; Gelman 2006), and uniform for 0 on a wide range (Clayton
1996; Spiegelhalter et a1 2004; Gelman 2006). The standard deviation for the stand-level
random intercept was insensitive to prior specification (i.e. all 95% credible intervals

overlapped) for the sapling abundance model (Figure 8.2) and percentage of seedlings

browsed model (Figure 8.3).

125

 

16 (0.01, 0.01) :

16 (0.001, 0.001) 4 :

 

Half Normal -

"l '

Pfior

1)

Uniform (0, 50) _

 

Uniform (0, 25) JV

 

 

I j I

2.0 2.5 3.0 3.5 4.0
Standard Deviation of Random Stand-level Intercept

Figure B.2. Effect of the prior on the standard deviation (0) of the random stand-level
intercept in the sugar maple sapling abundance generalized linear
multilevel model: mean, 90% credible interval (thick black line), and 95%
credible interval (thin black line) of the sterior distribution. IG (inverse

gamma) prior was on the precision (1/o ) and half normal and uniform

density priors were on o. The half normal distribution was ~N(0, 10,000),
bounded above 0.

 

16 (0.01, 0.01) -

16 (0.001, 0.001) -

Pfior

—_———
—__——-
Half Normal - w
Uniform (o, 100) - “—-
—————

 

Uniform (0, 50) “

 

1.0 1.5 2.0 2.5 3.0
Standard Deviation of Random Stand-level Intercept

Figure B.3. Effect of the prior on the standard deviation (0) of the random stand-level
intercept in the gap-level percentage browsed hemlock seedlings logistic
multilevel model: mean, 90% credible interval (thick black line), and 95%
credible interval (thin black line) of the posterior distribution. IG (inverse
gamma) prior was on the precision (1/02) and half normal and uniform
density priors were on o. The half normal distribution was ~N(0, 10,000),
bounded above 0.

126

B.S. Creating image of deer density across the study area using DNR pellet count
data collected from 1996-2000 and Universal Kriging
The Wildlife Division within the Michigan Department of Natural Resources
collected deer pellet count data across the Upper Peninsula of Michigan from 1981 to
2000 (Doepker et a1 1994; Hill 2001). I developed an interpolated map from this dataset
for comparison with deer density estimates at my northern hardwood study stands to
validate the spatial pattern observed in deer densities.
I used a subset of the DNR data collected from 1996-2000 within 113 Public Land
Survey sections within the four study area counties (Iron, Dickinson, Marquette and 1

Menominee). The DNR randomly selected sections from three strata and surveyed deer

pellets within five randomly selected 1/50-acre (81 m2) rectangular plots (3.7 m by 22.1

m) (Hill 2001). The average pellet count from these five plots was converted into deer

density (deer / kmz) (Appendix B, Section B.1). The years 1996 to 2000 were selected

because pellet count surveys were consistently made in more sections during these years
than they were during other years. Only sections where pellet surveys were performed
more than four years during the five year timeframe were included in analysis.

An empirical semivariogram was developed to describe the spatial correlation

between deer densities estimates from sections. Deer density estimates were loglo +1

transformed to comply with normality assumptions. The spatial coordinates for the deer
density estimates were the UTM Northing and Easting for the center of the section. Deer
density generally decreased with increasing UTM Northing (Figure B.4) so a
semivariogram was developed using residuals from the relationship between transformed

deer density and UTM Northing. Generalized Least Squared (GLS) estimation was used

127

 

2.0 - 0 °

Log 10 + 1 Transformed

 

Average Deer Density (deerl kmz) 1996-2000

 

 

 

I I I
5,050,000 5,1 00,000 5,150,000

PLS Section UTM Northing

Figure 8.4. Average estimated deer density (deer/kmz) from DNR pellet count surveys
1996-2000 within Iron, Dickinson, Marquette and Menominee counties
versus UTM Northing for the Public Land Survey section. Log
transformed deer density = 23.8518 — 4.0e-06 >< UTM Northing based on
GLS-estimation. AIC is 107.7 compared to 85.7 for the null model.

to take into account the correlation between observations. Residuals from the Ordinary

Least Squared estimated linear regression between log transformed deer density and

UTM Northing were first used to estimate the nugget, partial sill and range for the

semivariogram. The semivariogram values were used to develop the correlation matrix

necessary for GLS-estimation. GLS residuals were used to develop the final

semivariogram which was best approximated with an exponential model in the form

128

y(h) = Co + c(1 — exp(— 27-) ) where co (the nugget) is 0.04, c (the partial sill) is 0.10, a
a

(the range) is 20,000 m and h (the lag distance) is 2,000 m (Figure B.5). There were very
few points separated by a distance less than 1,600 m (Table B.5), but the sample size was

adequate for the other size bins.

1 l 1

 

0.15 -

0.10 —

Semivariance v(2,000 m)

0.05 -

 

 

y(2,000) = 0.04 + 0.10 X (1 — exp[

_ 3 x 2,000 )
20,000

 

 

 

 

 

 

l I 1

20,000 40,000 80,000
Distance (meters)

Figure B.5. Semivariogram of the residuals from the GLS-estimated relationship
between estimated average deer density (deer/kmz) and UTM Northing. A
lag distance of 2,000 m and an exponential model form were used.

129

The semivariogram was used to perform universal kriging (Goovaerts 1997) and
develop an interpolated map of back-transformed deer density across the study area
counties using geostatistical analyst in ERSI ArcMap v 9.2. The average standard error
between the observed deer density and predicted from cross validation (Figure 3.6) was

0.35 and the root-mean-square error was 0.32.

Table B.5. Number of section (N. sections) centers separated by different distances.
The average of the semivariance values from the points in the different
distance bins was used to develop the semivariogram (Figure B.3).

 

 

N. sections Distance (m) N. sections Distance (In)

6 1610 134 40967
24 3258 140 42999
25 5108 142 45057
28 6854 146 46933
52 9049 145 48994
64 11096 154 50940
92 13104 134 53045
92 15068 128 55027
97 16895 145 56908
101 19012 120 58948
94 21079 116 60972
98 23027 108 62994
120 25008 117 64984
131 27059 155 66887
125 29070 147 69137
134 30938 117 70963
131 32907 110 72946
157 35056 128 74979
129 37014 11 76084
149 38918

 

130

 

 

 

 

 

1.6 -
O O
O
. .
_ o e o
E. 1.4 . . O 0 .
3 ° . 0.. \ o 0 0 °. 0 o
O 0
é o 0 .. 0~ o o g . .0 0
£3 1.21 '1; o. O... O
.0 g e ’6 o 0’ o ‘
3 '13 o 0 °
.2 5 ° 0 . .0
E E 1 0 — o o . : o
O
a ° 0
2 g 0.8 - .
3 _l 0 . 0 o
> .0
an ° .
§
0 0.6 _ Q 0
O
0.4 _ °
I I I I
0.5 1.0 1.5 2.0
Average Estimated Deer Density (deer! kmz) 19915-2000
Log 10 + 1 Transformed
Figure B.6. Observed deer density estimates versus values predicted from cross

validation.

131

 

APPENDIX C: Sample R2WinBUGS code and WinBUGS models
Models were run with WinBUGS v.1.4.3 (Spiegelhalter et a1 2003) through R
v.2.9.1 with the package R2WinBUGS (Sturtz et a1 2005). Below are WinBUGS models
for Chapter 1 seedling and sapling abundance generalized linear multilevel models
(GLMM) (Section C.1), log transformed sapling stocking levels and square root
transformed Shannon Diversity Index linear multilevel models (LMM) (Section C.2), log
transformed sapling grth rates linear models (LM) (Section C3) and Chapter 2
percentage of seedlings browsed multilevel logistic regression models (Section C.4).
WinBUGS model code was derived from examples in Ntzoufras 2009 and Gelman and
Hill 2007. Code is also provided for running the sugar maple seedling model with
R2WinBUGS, with R code modified from Qian and Shen 2007 (Section C.5).
Section C]. WinBUGS code: Seedling and sapling GLMM
model{
for (i in 1:n){
y[i] ~ dnegbin(p[i], r2)
p[i]<- r2/(r2+mu[i])
Iog(mu[i])<—a[site[i]] + b.cv*cv[i] + b.co*co[i] + b.ss*ss[i]
}
b.cv~dnorm(0,0.0001) #prior for cv (competing vegetation) B
b.co~dnorm(0,0.0001) #co = canopy openness
b.ss~dnorm(0,0.0001) #ss= seed production potential
r2~dgamma(0.1,0. 1)

for (i in 1:J){

132

 

a[j] ~ dnorm(a.hat[j], tau.a)
a.hat[j] <- g.0 + g.atdhp*atdhp[j] + g.atd*atd[j] + g.atm*atmfi] + g.tmc*tmc[i] +
g.d*d[j] + g.tsh*tshfi]
}
g.0 ~ dnorm(O, 0.0001) #g0= overall intercept (AOCa is reference Habitat Type)
g.atdhp ~ dnorm(O, 0.0001) #adthp= ATD-Hp Habitat Type
g.atd ~ dnorm(O, 0.0001) #atd= ATD
g.atm ~ dnorm(O, 0.0001) #atm= ATM
g.tmc ~ dnorm(0, 0.0001) #tmc= TMC
g.d ~ dnorm(O, 0.0001) #d= deer density
g.tsh ~ dnorm(O, 0.0001) #tsh= time since harvest
tau.a <- pow(sigma.a, -2) #random stand-level intercept precision
sigma.a ~ dunif (0, 50) #prior on random stand-level intercept stdev
}
Section C.2. WinBUGS code: Log transformed sapling stocking level and square
root transformed Shannon Diversity Index LMM
model{
for (i in 1:n){
y[i] ~ dnorm(mu[i], tau.b)
mu[i]<— a[site[i]] + b.cv*cv[i] + b.co*co[i] + b.hm*hm[i]
}
b.cv~dnorm(0,0.0001 )

b.co~dnorm(0,0.0001)

133

b.hm~dnorm(0,0.0001) #hm= overstory diversity

tau.b<-pow(sigma.b, -2)

sigma.b~dunif (0, 50)

for (j in l:.1){
a[j] ~ dnorm(a.hat[j], tau.a)
a.hat[j] <- g.0 + g.atdhp*atdhp[j] + g.atd*atd[j] + g.atm*atm[i] + g.tmc*tmc[j] +
g.d*d[j] + g.tsh*tsh[j]

}

g.0 ~ dnorm(O, 0.0001)

g.atdhp ~ dnorm(O, 0.0001)

g.atd ~ dnorm(O, 0.0001)

g.atm ~ dnorm(O, 0.0001)

g.tmc ~ dnorm(O, 0.0001)

g.d ~ dnorm(O, 0.0001)

g.tsh ~ dnorm(O, 0.0001)

tau.a <- pow(sigma.a, -2)

sigma.a ~ dunif (0, 50)
}
Section C.3. WinBUGS code: Log transformed sapling growth rates LM
model{

for (i in 1:n){

y[i] ~ dnorm(mu[i], tau.b)

I34

 

mu[i]<— b.0 + b.ht*ht[i] + b.ar*ar[i] + b.cv*cv[i] + b.co*co[i] + b.atdhp*atdhp[i]
+ b.atd*atd[i] + b.atm*atm[i] + b.tmc*tmc[i] + b.d*d[i] + b.tsh*tsh[i]
}
b.0 ~ dnorm(O, 0.0001)
b.ht~dnorm(0,0.0001) #ht=sapling height
b.ar~dnorm(0,0.0001) #ar= advanced regeneration
b.cv~dnorm(0,0.0001 )
b.co~dnorm(0,0.0001)
b.atdhp ~ dnorm(O, 0.0001)
b.atd ~ dnorm(O, 0.0001)
b.atm ~ dnorm(O, 0.0001)
0th ~ dnorm(0, 0.0001)
b.d ~ dnorm(0, 0.0001)
b.tsh ~ dnorm(O, 0.0001)
tau.b~dgamma(0.001 ,0.001 ) #prior on sapling-level precision

sigma.b<- 1/sqrt(tau.b) #sapling-level standard deviation

Section C.4. WinBUGS code: Gap-level percentage of seedlings browsed or severe
browsed multilevel logistic regression model
model{
for(i in 1:n.gaps){
r[i]~ dbin(p.bound[i],n[i]) #r[i] = predicted number browsed of n[i] seedlings

p.bound[i] <- max(0, min(1, p[i])) #constrains probability estimates btwn 0 and 1

135

logit(p[i]) <- a.stand[stand[i]] + b.pel*pel[i] + b.dlc*dlc[i]
}
b.pel~ dnorm(O, 0.0001) #pel= pellet counts / gap transect
b.dlc~ dnorm(O, 0.0001) #dlc= gap distance to lowland conifer

for(k in 1:n.stands){

a.stand[k]~dnorm(a.hat[k], tau.a)

a.hat[k]<- g.0 + g.dsd*dsd[k]
}
g.0~ dnorm(O, 0.0001) #g.0= overall intercept
g.dsw~ dnorm(O, 0.0001) #dsw= deep snow weeks
tau.a<-pow(sigma.a, -2)
sigma.a~dunif(0, l 00)

}

 

Section C.5: Code for running sugar maple seedling code with R2WinBUGS

AcesacGLMMl <- function(infile=AS, infi162=ASsite){
y <- infile$AcesacSeedTot

n <- length(y)

cv<-infile$CompetingS

co<- infile$CO

ss<-infile$RibbensAcesac

tsh<-infi1e2$TSHarvest

d<-infile2$DeerDen

atdhp=infile2$ATDHp

136

atd=infile2$ATD

atm=infile2$ATM

tmc=infile2$TMC

J<- length(tsh)

site <- as.numeric(ordered(infile$Site))

bugs.dat <- list(n=n, y=y, J=J, cv=cv, co=co, ss=ss, tsh=tsh, d=d,
atdhp=atdhp, atd=atd, atm=atm, tmc=tmc, site=site)

## input data needed to run the WinBUGS program

initsl <- 1ist(a=morm(length(infi162$Site)), g.atdhp=rnorm(1), g.atd=morm(1), 1
g.atm=morm(1), g.tmc=rnorm(1), g.tsh=rnorm(l), g.d=morm(1),
g.0=morm(l), sigma.a=runif( 1), r2=1, b.cv=morm(1), b.co=rnorm(1),
b.ss=morm(1))

inits2 <- Iist(a=morm(length(infileZ$Site)), g.atdhp=rnorm(1), g.atd=morm(1),
g.atm=morm(1), g.tmc=morm(l), g.tsh=morm(1), g.d=morm(1),
g.0=rnorm(1), sigma.a=runif(1), r2=2, b.cv=morm(1), b.co=morm(l),
b.ss=morm(1))

inits3 <- list(a=morm(length(infile2$Site)), g.atdhp= morm(1), g.atd=morm(1), g.atm=
morm(1), g.tmc=rnorm(1), g.tsh=morm(1), g.d=morm(1), g.0=rnorm(1),
sigma.a=runif(1), r2=5, b.cv=morm(1), b.co=morm(l), b.ss=morm(1))

inits <- list (initsl, inits2, inits3)

## variables to be monitored

137

parameters <- c("g.0", "g.atdhp", "g.atd", "g.atm", "g.tmc", "g.tsh", "g.d", "sigma.a",
"b.cv", "b.co", "b.ss", "r2", "mu") return(list(para=parameters, data=bugs.dat,
inits=inits))

}

input.to.bugs <- AcesacGLMM1()

AcesacGLMM] <- bugs(input.to.bugs$data, input.to.bugs$inits,
input.to.bugsSpara, model.file="Seedling_GLMM.txt", n.chains=3,
n.iter=70000, n.burnin=10000, n.thin=5, DIC=T, bugs.directory="

C:/Program Files/WinBUGS14", debug=TRUE)

138

APPENDIX D: Assessing Convergence

Seedling and sapling abundance generalized linear multilevel models (GLMM),
log transformed stocking-level linear multilevel models (LMM), square root transformed
Shannon Diversity Index linear multilevel models (LMM), log transformed sapling
growth rate linear models (LM) (Chapter 1) and percentage of seedlings browsed
multilevel logistic models (Chapter 2) were fit using Markov Chain Monte Carlo
sampling in WinBugs. Three parallel chains were run with dispersed randomly selected
starting values. Model convergence was assessed using the R package coda (Plummer et
a1 2009). If the Rafiery-Lewis diagnostic for a model indicated that an insufficient
number of iterations had been run to have a 95% probability of estimating the 2.5%
quantile with a tolerance of +/-0.5%, it was rerun with an increased number of iterations
and an increased thinning factor. I is the proportional increase in N attributable to
autocorrelation and indicates poor mixing of the MCMC chains if values are > 5.

Based on both diagnostics, the chains converged for all variables across all
analyses. All Rhat values (potential scale reduction factor) from the Gelman-Rubin
diagnostic were 1.0 (Gelman and Hill 2007) and the number of iterations run were greater

than N from the Raftery-Lewis diagnostic and I was less 5 (Table D.1 - D. 10)

139

Table D.1. Convergence in models of spatial trends in seedling and sapling
abundance and log transformed stocking-level: sample size (N obs.), bum—
in length (N .bumin), iterations afier burn-in (N .iter), thinning factor
(T.factor) and maximum Raftery-Lewis diagnostic N (number of iterations
required to have a 95% probability of estimating the 2.5% quantile with a

tolerance of +/-0.5%) and I (Dependence Factor) from three parallel

 

 

MCMC chains.

Species Sugar maple Ironwood Other species
Seedling abundance GLMM
Nobs. gaps /stands 337 /59 337 /59 337 /59
N .bumin 10,000 15,000 10,000
N.iter 60,000 85,000 60,000
T.factor 5 5 5
Variable Raftery-Lewis Diagnostics (N / I )
Intercept 4001 / 1.07 3945/ 1.05 4140/ 1.11
Northing 3940/ 1.05 4338/ 1.16 3972/ 1.06
Easting 3787/ 1.01 4290/ 1.15 3839/ 1.02
Stand Stdev 4344/ 1.16 10550/ 2.82 4487/ 1.20
Dispersion parameter 3881 / 1.04 4150/ 1.11 3903 / 1.04
Sapling abundance GLMM
N obs. gaps/stands 331 /58 331 /58 331 /58
N.burnin 10,000 10,000 10,000
N.iter 60,000 60,000 60,000
T.factor 5 5 5
Variable Rafiery-Lewis Diagnostics (N / I )
Intercept 9786/ 2.61 4001 / 1.07 4151 / 1.11
Northing 4275 / 1.14 3892 / 1.04 3865/ 1.03
Easting 3865 / 1.03 3865 / 1.03 3972/ 1.06
Stand Stdev 9364/ 2.50 4384/ 1.17 4340/ 1.16
Dispersion parameter 3881 / 1.04 3839/ 1.02 4067/ 1.09
Sapling stocking LMM
Nobs. gaps/stands 164 I45 184 /49 192 /48
N.burnin 5,000 5,000 5,000
N.iter 15,000 15,000 15,000
T.factor 1 1 1
Variable Rafiery-Lewis Diagnostics (N / I )
Intercept 4845 / 1.29 4307/ 1.15 4290/ 1,15
Northing 4713 / 1.26 4394/ 1.17 4503/ 1.20
Easting 4346/ 1.16 4299/1.15 4361 / 1.16
Stand Stdev 11360 / 3.03 9518 / 2.54 9024 / 2_41
Gap Stdev 5007/ 1.34 5257/ 1.4 5045/ 135

 

Stand Stdev= standard deviation for the random stand-level intercept.

140

Table D2. Convergence in LMM of spatial trends in square root transformed
Shannon Diversity Index (H’): Sample size (N obs.), burn-in length

(N .bumin), iterations after burn-in (N .iter), thinning factor (T.factor) and
maximum Raftery-Lewis diagnostic N (number of iterations required to
have a 95% probability of estimating the 2.5% quantile with a tolerance of

+/-0.5%) and I (Dependence Factor) from three parallel MCMC chains.

 

 

Species H’ Seedlings H’ Saplings 1-2 m H’ Saplings 1-7 m
N obs. gaps/stands 182 / 54 180 / 42 201 /46
N.burnin 10,000 5,000 1 0,000
N.iter 60,000 15,000 60,000
T.factor I 1 1
Variable Raftery-Lewis Diagnostics ( N / I )

Intercept 3918/ 1.05 4583 / 1.22 3918/ 1.05
Northing 3839/ 1.02 4541 / 1.21 3978/ 1.06
Easting 3892/ 1.04 4338/ 1.16 3761 / 1.00
Stand Stdev 6609/ 1.76 11274 / 3.01 4527/ 1.21
Gap Stdev 3952 / 1.05 4939/ 1.32 3836/ 1.02

 

Stand Stdev= standard deviation for the random stand-level intercept.

Table D.3. Convergence in LM of spatial trends in log transformed sapling growth
rate after harvest: Sample size (N obs.), burn-in length (N .bumin),
iterations after burn-in (N .iter), thinning factor (T.factor) and maximtun
Rafiery-Lewis diagnostic N (number of iterations required to have a 95%
probability of estimating the 2.5% quantile with a tolerance of +/—O.5%)

and 1 (Dependence Factor) from three parallel MCMC chains.

 

 

SPCCieS Sugar maple Ironwood

N obs. gaps/stands 127 / 26 98 / 30
N.bumin 5,000 5,000
N.iter 15,000 15,000
T.factor 1 1
Variable Rafiery-Lewis Diagnostics (N /1 )
Intercept 3803 / 1.02 3879/ 1.04
Northing 3908/ 1.04 3908 / 1.04
53511118 3802/ 1.01 3823 / 1.02
Sapling Stdev 3846/ 1.03 3946/ 1.05

 

141

 

Table D4. Convergence in GLMM of seedling abundance with gap- and stand-level
covariates: Sample size (N obs.), burn-in length (N .bumin), iterations after
burn-in (N .iter), thinning factor (T.factor) and maximum Rafiery-Lewis
diagnostic N (number of iterations required to have a 95% probability of
estimating the 2.5% quantile with a tolerance of +/-0.5%) and 1
(Dependence Factor) from three parallel MCMC chains.

 

 

Species Sugar maple Ironwood Other spp
N obs. gaps/stands 337 /59 337 / 59 337 / 59
N.bumin 10,000 15,000 15,000
N.iter 60,000 84,000 1 ,600,000
T.factor 5 7 20
Variable Raftery-Lewis Diagnostics ( N / I )
Full model
Canopy openness 7615 /2.03 5593 / 1,49 7217/ 1.93
Competing veg 6796/ 1.81 5705 / 1,52 6852/ 1.83
SP2003 6138/ 1.64 4886/ 1.30 6273 / 1.67
Intercept 3933 / 1.05 4067 / 109 4404/ 1.18
ATD-Hp 3945 / 1.05 4023 / 1,03 3918/ 1.05
ATD 3761 / 1.00 3242 /2,20 3918/ 1.05
ATM 3945 / 1.05 4234/ 1,14 3945 / 1.05
TMC 3918 / 1.05 10430 / 2,30 3865 / 1.03
Deer density 3945 / 1.05 4001 / 1,07 3945 / 1.05
TSHarvest 3839/ 1.02 4547/ 121 3918/ 1.05
Dispersion parameter 3972 / 1.06 4151 / 1_11 4011 / 1.07
Stand Stdev 4226/ 1.13 11494 / 307 9284 / 2.48
Null model

Intercept 3865 / 1.03 4272 / 1.14 4056/ 1.08
Dispersion parameter 3956/ 1.06 3888/ 1.04 3945 / 1.05
Stand Stdev 4167 / 1.11 8066/ 2.15 4665 / 1.25

 

SP2003= seedling production potential (EdIameter/distancez) in 2008, TSHarvest=time
since harvest, Stand Stdev= standard deviation for the random stand-level intercept.

142

Table D5. Convergence in GLMM of sapling (1-2 in) abundance with gap- and
stand-level covariates: Sample size (N obs.), burn-in length (N .bumin),
iterations after burn-in (N .iter), thinning factor (T.factor) and maimum
Raftery-Lewis diagnostic N (number of iterations required to have a 95%
probability of estimating the 2.5% quantile with a tolerance of +/-0.5%)

and 1 (Dependence Factor) from three parallel MCMC chains.

 

 

Species Sugar maple Ironwood Other spp
N obs. gaps/stands 331 /58 331 /58 331 /58
N.bumin 15,000 15,000 10,000
N.iter 84,000 84,000 60,000
T.factor 5 5 5
Variable Raftery-Lewis Diagnostics ( N/ 1)
Full model
Canopy openness 9812 / 2.62 6500/ 1.74 6393 / 1-84
Competing veg 9216 / 2.46 9008 / 2.40 9369 / 2.50
SPHar-vest 4768/ 1.27 5130/ 1.37 5341 / 1-43
111191091)t 9940/ 2.65 3972/ 1.06 4140/ 1-11
ATD-Hp NA 3685 /0,93 3839 / 1.02
ATD 4151/1.11 3735/100 3972/1.06
ATM 3967/ 1.06 3721 /0.99 3813 / 1.02
TMC 4028/ 1.08 3787 / 1.01 3813 / 1.02
Deer density 4496/ 1.20 3787/ 1.01 3929/ 1-05
TSHarvest 3813 / 1.02 3761 /1.00 3839/ 1.02
Dispersion parameter 4021 / 1,07 3319/ 102 3972/ 1.06
Stand Stdev 9612 / 2.57 4094/ 1.09 4190/ 1.12
Null model

Intercept 8576 / 2.29 4098/ 1.09 3945 / 1.05
Dispersion parameter 3892/ 1.04 3929/ 1.05 4011 / 1.07
Stand Stdev 9858/ 2.63 4295 / 1.15 4255 / 1.14

 

SPHarvest= seedling production potential (£diameter/distance2) at time of harvest,
TSHarvest=time since harvest, Stand Stdev= standard deviation for the random stand-

level intercept.

143

Table D6. Convergence in LMM of log transformed sapling 1-7 m tall stocking with
gap- and stand-level covariates: Sample size (N obs.), burn-in length
(N.burnin), iterations after burn-in (N .iter), thinning factor (T.factor) and
maximum Raftery-Lewis diagnostic N (number of iterations required to
have a 95% probability of estimating the 2.5% quantile with a tolerance of

+/-0.5%) and 1 (Dependence Factor) from three parallel MCMC chains.

 

 

Species Sugar maple Ironwood Other Spp

N obs. gaps/stands 164 I45 184 /49 192 /48
N.bumin 5,000 5,000 5,000
N.iter 15,000 15,000 15,000
T.factor 1 I 1
Variable Rafiery-Lewis Diagnostics ( N/ I )

Full model
Canopy openness 5094/ 1.36 9618 /2.57 7610 / 2,03
Competing veg 5835 / 1.56 9182 /2.45 8618 / 2,30
SPHarvest 4267 / 1.14 4661 / 1.24 4643 / 1.24
Intercept 7652 / 2.04 4508/ 1.20 4677/ 1.25
ATD-Hp 4479/ 1.20 4558/ 1.22 4583 / 1,22
ATD 4410/1.18 4508/1.20 4290/1,15
ATM 4385 / 1.17 4387/ 1.17 4338/1,16
TMC 5152/ 1.38 8032 /2.I4 4533/ 121
Deer density 4267/1.14 4252 / 1.14 4410/1,18
TSHarvest 4299/1.15 4338 / 1.16 4129/ 1.10
Stand Stdev 10438 / 2.79 10440 / 2.79 10152 / 2,71
Gap Stdev 5094/ 1.36 5343 / 1.43 5044/ 135
Null model '

Intercept 4402/ 1.18 4418/ 1.18 4229/ 1.13
Stand Stdev 5825 / 1.55 9988/ 2.67 5770/ 1.54
Gap Stdev 5021 / 1.34 5021 / 1.34 5094 / 1.36

 

SPHarvest= seedling production potential (Ediameter/distancez) at time of harvest,
TSHarvest=time since harvest, Stand Stdev= standard deviation for the random stand-

level intercept.

144

Table D7. Convergence in LMM of square root transformed of Shannon Diversity
Index (H’) with gap- and stand-level covariates: Sample size (N obs.),
burn-in length (N .bumin), iterations afier burn-in (N .iter), thinning factor
(T.factor) and Rafiery-Lewis diagnostic N (number of iterations required
to have a 95% probability of estimating the 2.5% quantile with a tolerance
of +/-0.5%) and 1 (Dependence Factor) from three parallel MCMC chains.

 

 

 

Species H’ Seedlings H’ Saplings 1-2 m H’ Saplings 1-7 in
N obs. gaps/stands 182 / 54 I80 / 42 201 /46
N.bumin 20,000 5,000 5,000
N.iter 200,000 15,000 15,000
T.factor 10 1 1
Variable Raftery-Lewis Diagnostics ( N / l )

Full model
Canopy openness 3740/ 1.00 5010/ 1.34 9214/ 2.46
Competing veg 3834/ 1.02 9784 / 2.61 8990 /2.40
11’ memory] 3872 / 1.03 9120 / 2.43 5038/ 1.34
Intercept 3818/ 1.02 4800/ 1.28 8762 / 2.34
AT D-Hp 3740/ 1.00 4687/ 1.25 4765 / 1.27
ATD 3787/ 1.01 4299/ 1.15 8800/2.35
ATM 3818 / 1.02 4669/ 1.25 4982 / 1.33
TMC 3802/ 1.01 7864/2.10 8752/2.34
Deer density 3802/ 1.01 4660/ 1.24 4765/ 1.27
TSHarvest 3802/ 1.01 4483 / 1.20 4661 / 1.24
Stand Stdev 9650 / 2.58 14559 / 3.89 12726 / 3.40
Gap Stdev 3845 / 1.03 5103 / 1.36 4997/ 1.33

Null model
Intercept 3853 / 1.03 4385 / 1.17 4687/ 1.25
Stand Stdev 4095/ 1.09 9358 / 2.50 12498 / 3.34
Gap Stdev 3978/ 1.06 4970/ 1.33 5048/ 1.35
TSHarvest=time since harvest, Stand Stdev= standard deviation for the random stand-level
intercept.

l
H’Overstory for H’Seedlings is Shannon Diversity of overstory trees, H’Overstory for
H’Saplings 1-2 and 1-7 m tall is Shannon Diversity of overstory trees and stumps

145

 

Table D8.

Convergence in LM of log transformed sapling growth rates with gap- and
stand-level covariates: Sample size (N obs.), burn-in length (N .bumin),
iterations after burn-in (N .iter), thinning factor (T.factor) and Raftery-
Lewis diagnostic N (number of iterations required to have a 95%
probability of estimating the 2.5% quantile with a tolerance of +/-0.5%)
and I (Dependence Factor) from three parallel MCMC chains.

 

 

Species Sugar maple Ironwood
N obs. gaps/stands 127 / 26 98 / 30
N.bumin 5,000 5,000
N.iter 15,000 15,000
T.factor 1 1
Variable Rafiery-Lewis Diagnostics ( N/ I )
Full model
Intercept 4043 / 1.08 3913 / 1.04
Advanced regen 3832/ 1.02 3802/ 1.01
Sapling height 3740/ 1.00 3802/ 1.01
Canopy openness 3844/ 1.03 3823 / 1.02
Competing veg 3782/ 1.01 3908/ 1.04
ATD-Hp 3844/ 1.03 3811 / 1.02
ATD 3938/ 1.05 3782/ 1.01
ATM 3799/ 1.01 3782 / 1.01
TMC 3802/1.01 3761 /1.00
Deer density 3938/ 1.05 3865 / 1.03
TSHarvest 3746/ 1.00 3770/ 1.01
Sapling Stdev 3860/ 1.03 3959/ 1.06
Null model
Intercept 3841 / 1.03 3834/ 1.02
Sapling Stdev 3841 / 1.03 3782/ 1.01

 

Advanced regen= bivariate variable indicating if a variable is advanced regeneration,
TSHarvest=time since harvest.

146

 

Table D9.

Convergence in multilevel logistic models of gap-level percentage of
seedlings browsed within gaps with gap- and stand-level covariates:
Sample size (N obs.), burn-in length (N .bumin), iterations after burn-in
(N .iter), thinning factor (T.factor) and Rafiery—Lewis diagnostic N
(number of iterations required to have a 95% probability of estimating the
2.5% quantile with a tolerance of +/-0.5%) and 1 (Dependence Factor)
from three parallel MCMC chains.

 

 

Species Hemlock White pine

N obs. gaps/stands 338 / 34 338 / 34

N.bumin 10,000 10,000

N.iter 30,000 30,000

T.factor 5 5

Variable Raftery-Lewis Diagnostics ( N/ l )
Model 1 (null)

Intercept 3887/ 1.04 3761 / 1.00
Stand Stdev 3962 / 1.06 3900/ 1.04
Model 2 (Deer + DSW)

Intercept 4084/ 1.09 3710/ 0.99
Deer density 4028/ 1.08 3973 / 1.06
DSW 3973 / 1.06 3866/ 1.03
Stand Stdev 4028/ 1.08 4028/ 1.08
Model 3 (PGC + DSW)

Intercept 3973 / 1.06 3710/ 0.99
Centered PGC 3973 / 1.06 3710 / 0.99
DSW 3866/ 1.03 3761 / 1.00
Stand Stdev 3973 / 1.06 4084/ 1.09

 

2
Deer= winter deer density in 10’s of deer / km , DSW=deep snow weeks (number of
weeks from November 2007 to April 2008 with snow depth > 40 cm), PGC= pellet

group count (the number of deer fecal pellets found within a 40 in area in a planted
harvest gap), Stand Stdev= standard deviation for the random stand-level intercept.

147

 

Table D.10. Convergence in multilevel logistic models of gap-level percentage of
seedlings browsed within gaps with gap- and stand-level covariates for
stands where DLC could be determined: Sample size (N obs.), burn-in
length (N .bumin), iterations after burn-in (N .iter), thinning factor
(T.factor) and Raftery-Lewis diagnostic N (number of iterations required
to have a 95% probability of estimating the 2.5% quantile with a tolerance
of +/-0.5%) and I (Dependence Factor) from three parallel MCMC chains.

 

 

Species Hemlock White pine

N obs. gaps/stands 126/ 13 126/ 13

N.bumin 10,000 10,000

N.iter 30,000 30,000

T.factor 5 5

Variable Raflery-Lewis Diagnostics (N / l )
Model 1 (null)

Intercept 3973 / 1.06 3866/ 1.03
Stand Stdev 4084/ 1.09 3973 / 1.06
Model 4 (DLC + PGC + DSW)

Intercept 4050/ 1.08 3973 / 1.06
Centered DLC 4558/ 1.22 5485/ 1.46
Centered PGC 5408/ 1.44 5705 / 1.52
DSW 3866/ 1.03 3866/ 1.03
Stand Stdev 3973 / 1.06 3866/ 1.03
Model 5 (DLC)

Intercept 4084/ 1.09 3761 / 1.00
Centered DLC 5095/ 1.36 4623 / 1.23
Stand Stdev 4028/ 1.08 4028 / 1.08

 

 

DSW=deep snow weeks (number of weeks fi'om November 2007 to April 2008 with
snow depth > 40 cm), PGC= pellet group count (the number of deer fecal pellets

found within a 40 m area in a planted harvest gap), DLC= distance to lowland

conifer in 105 of m, Stand Stdev= standard deviation for the random stand-level
intercept.

 

148

APPENDIX E: Model Parameters

Bayesian inference was used with Markov chain Monte Carlo (MCMC) sampling

to estimate parameter values and obtain credible intervals for generalized linear mixed

models (GLMM), linear mixed models (LMM) and linear models (LM).

 

 

 

Table E.1. Effects of gap and stand-level variables on seedling abundance <1 m tall:
mean parameter estimates and 95% credible interval from posterior
distributions. Positive values indicate increases in seedling abundance
with a one standard deviation increase in that covariate on the log scale.
AOCa is the reference Habitat Type. A decrease in DIC of 5 or more
indicates improved performance for full model over the null model.

Variable Sugar maple Ironwood Other Spp
Mean CI Mean CI Mean CI
GLMM full model
Gap level
Canopy openness -0.08 -0.31 — 0.14 -0.13 -0.52 — 0.26 0.21 -0.07 — 0.48
Competing veg -0.23** -0.42 — -0.05 -0.34"' -0.71 — 0.02 -0.43** -0.70 — -0.l6
3pm,; 0.12 -0.06 — 0.30 0.37" 0.08 — 0.74 -0.03 -0.25 — 0.21
Stand level
Intercept 1.40" 0.80 — 2.00 -0.37 -I .23 - 0.43 -0.90** -1.82 — -0.08
ATD-Hp 0.18 -0.70 — 1.08 -0.02 -l.27 - 1.21 1.49“ 0.30 —- 2.75
ATD 1.50“ 0.55 — 2.48 -2.1 I" -3.81 — -0.61 0.83 -0.49 — 2.19
ATM 1.15" 0.16 — 2.15 -1.86** -3.53 - -0.29 0.94 -0.44 — 2.35
TMC 1.74“ -0.13 — 3.65 -2.91 -7.14 — 0.55 3.68" 1.18 - 6.31
Deer density -0.24 -0.61 - 0.12 0.11 -0.43 - 0.68 -0.08 -0.61 - 0.42
Time since harvest 0.24 -0.10 - 0.58 -0.53* -l.l l - 0.04 0.21 -0.27 -— 0.70
Variance components
r 0.94 0.76 — 1.16 0.50 0.32 — 0.75 0.78 0.56 — 1.06
Stand-level Stdev 1.14 0.88 - 1.45 1.45 0.94— 2.10 1.52 1.11 — 2.06
DIC 2172.61 701.20 1132.04
GLMM null model
Intercept 2.06“ 1.69 - 2.42 -0.99** -1 .58 - -0.47 -0.02* -0.57 — 0.48
r 0.90 0.73 — 1.09 0.45 0.29 — 0.68 0.77 0.56 — 1.05
Stand-level Stdev 1.32 1.05 — 1.66 1.54 1.06 - 2.15 1.74 1.33 — 2.28
DIC 2184.88 714.41 1132.46
N. gaps / stands 337 / 59 337 / 59 337 / 59

 

"indicates 95% credible interval does not overlap 0, * indicate 90% CI does not overlap 0

r= negative binomial gap-level dispersion parameter, SP2003= seedling production potential

2
(EdIameter/distance ) for 2008, Stand-level Stdev= standard deviation for the random stand-level intercept.

149

 

 

Table B.2. Effects of gap and stand-level variables on sapling abundance 1-2 m tall:
mean parameter estimates and 95% credible interval from posterior
distributions. Positive values indicate increases in sapling abundance with
a one standard deviation increase in that covariate on the log scale. AOCa
is the reference Habitat Type for ironwood and other species. ATD-Hp
and AOCa were combined as the reference Habitat Types to improve
convergence of the sugar maple model. A decrease in DIC of 5 or more
indicates improved performance for full model over the null model.

 

 

 

Variable Sugar maple Ironwood Other Spp

Mean CI Mean CI Mean CI

GLMM full model

Gap level
Canopy openness 0.29“ 0.03 — 0.56 0.05 -0.30 - 0.40 0.59" 0.33 — 0.86
Competing veg -0.49** -0.72 — -0.26 -O.15 -0.48 — 0.19 -0.42** -0.70 — -0.I4
SPHarvest -0.07 -0.23 — 0.09 -0.04 -O.23 — 0.21 0.07 -0.10 — 0.27
Stand level
Intercept -5.32** -7.70 — -3.53 0.38 -l.30 — 1.91 -0.67 -2.18 — 0.71
ATD-Hp NA NA -0.48 -2.80 — 1.86 0.76 -l.30 — 2.85
ATD 7.63" 5.23 — 10.58 -1.48 —4.06 — 1.08 1.64 -0.51 — 3.89
ATM 7.40" 4.90 — 10.46 -0.41 -2.96 — 2.23 1.47 -0.78 - 3.78
TMC 8.34” 3.59 — 13.80 -I.75 -6.93 — 3.19 3.35 -0.80 — 7.71
Deer density -1.24** -2.47 — -0.08 0.42 -0.55 — 1.37 -0.49 -I .34 - 0.36

Time since harvest 0.75 -0.17 — 1.73 -0.02 -0.90 —- 0.88 0.39 -0.40 — 1.19
Variance components

r 1.49 1.06 — 2.01 0.64 0.48 — 0.82 0.84 0.64 — 1.07
Stand-level Stdev 2.93 2.10 — 4.11 3.03 2.25 — 4.08 2.68 2.05 — 3.50
DIC 1 197.10 1299.96 1468.80
GLMM null model
Intercept -2.05** -4.10 — -0.41 -0.08 -0.95 — 0.70 0.44 -0.35 - 1.18
r 1.28 0.94 — 1.72 0.64 0.49 — 0.82 0.74 0.58 - 0.94
Stand-level Stdev 5.46 3.99 — 7.55 2.78 2.12 — 3.65 2.64 2.05 — 3.41
DIC 1218.17 1296.83 1489.85
N. gaps / stands 331 /58 331 /58 331 /58

 

"indicates 95% credible interval does not overlap 0

F negative binomial gap-level dispersion parameter, SPHarvest= seedling production potential

2
(£diameter/distance ) at the time of harvest, Stand-level Stdev= standard deviation for the random stand—
level intercept.

150

Table E.3. Effects of gap and stand-level variables on log transformed sapling 1-7m
tall: mean parameter estimates and 95% credible interval from posterior
distributions. Only gaps with stocking-levels > 0 by species were used in
these analyses. Positive values indicate increases in stocking-levels on the

log scale with a one standard deviation increase in that covariate. AOCa is

the reference Habitat Type. A decrease in DIC of 5 or more indicates
improved performance for full model over the null model.

 

 

 

Variable Sugar maple Ironwood Other Spp
Mean CI Mean CI Mean CI

LMM full model
Gap level
Canopy openness -0.05 -O.30 - 0.21 0.09 -020 - 0.38 03 1 .1... 0,02 _ 0.61
Competing veg -0.36** -0.61 — -O.l I -0.38** -O.65 _ -0_ 10 -O.36** -0.64 _ -003
SPHarvest 0-02 '0.” - 0-22 0-00 -0.16 —- 0.16 0.02 -0.17 — 0.23
Stand level
Intercept -2.10** -2.99 - -1.24 -2.30** -304 _ -1 _57 -343" -4.26 _ -2.63
ATD-Hp -0-09 -1.40 - 1.22 -0. 14 -1.21 — 0.93 0.63 -0.61 - 1.87
ATD 139“ 0.18 — 2.60 0.00 -1.22 — 1.21 1.21" 0.01 — 2.42
ATM 0.69 -0.53 - 1.93 -0.16 -1.37 — 1.04 0.61 -0.65 — 1.89
TMC 0-56 -1-59 - 2.70 0-83 -2.05 — 3.70 1.28 -1.02 — 3.59
Deer density -O.26 -0.74 — 0.22 0.20 -029 _ 0,70 -023 .074 _ 023
Time since harvest 0.02 -0.40 — 0.43 0.18 -020 .. 0,57 0,24 -023 _ 0,72
Variance components
Stand-level Stdev 1.20 0.86- 1.62 1.20 0,33- 1.60 1,30 094. 1,73
Gap-level Stdev 1.12 0.99— 1.28 1.18 1,04- 1,33 1,25 1,12- 1,41
DIC 540.62 622.61 674.85

LMM null model
Intercept -1.50** -1.96 — -1.05 -2.34** -2.72 — -l.96 -2.73** -3.17 - -2.30
Stand-level Stdev 1.35 1.03 — 1.75 1.14 0.84 — 1.50 1.34 1.02 — 1.74
Gap-level Stdev 1.14 1.01 — 1.30 1.20 1.07 — 1.36 1.27 1.13 — 1.42
DIC 544.32 625.81 675.80
N. gaps / stands 164 / 45 184 / 49 I93 / 48

 

"indicates 95% credible interval does not overlap 0

2
SPHarvest= seedling production potential (Ediameter/distance ) at the time of harvest, Stand-level Stdev=
standard deviation for the random stand-level intercept.

151

 

Table E.4. Effects of gap and stand-level variables on square root transformed
Shannon Diversity Index (H’) for seedlings, saplings 1-2 m and saplings
1-7 m tall: mean parameter estimates and 95% credible interval from
posterior distributions. Only gaps with H’ > 0 by size class were used in
these analyses. Positive values indicate increases in H’ on the square root
scale with a one standard deviation increase in that covariate. AOCa is the
reference Habitat Type. A decrease in DIC of 5 or more indicates
improved performance for full model over the null model.

 

 

 

Variable H’ Seedlings H’ Saplings 1-2 m H’ Saplings 1-7 m

Mean Cl Mean Cl Mean CI
LMM full model

Gap level

Canopy openness 0.02 -0.02 — 0.05 0.00 -0.04 — 0.05 0.00 -0.04 — 0.04

Competing veg -0.01 -0.04 — 0.03 0.03 -0.01 — 0.07 0.02 -0.02 — 0.06

H’ Overstoryl 0.07" 0.04 — 0.10 0.02 -0.02 — 0.07 0.02 -0.02 — 0.06

Stand level

Intercept 0.75" 0.68 — 0.82 0.78" 0.66 — 0.90 0.80“ 0.69 -— 0.90

ATD-Hp 0.04 -0.06 — 0.14 -0.1 1 -0.28 — 0.07 -0.09 -0.24 — 0.06

ATD -0.10* -0.21 — 0.01 -0.17** -0.34 - -0.005 -0.14* -0.29 — 0.01

ATM -0.10* -0.22 — 0.01 -0.13 -0.30 — 0.04 -0.10 -0.25 — 0.05

TMC 0.00 -0.21 — 0.21 0.10 -0.21 — 0.40 0.06 -0.21 — 0.33

Deer density 0.01 -0.04 — 0.05 -0.04 -0.10 — 0.03 -0.03 -0.09 - 0.03

Time since harvest 0.02 -0.02 — 0.06 0.00 -0.06 — 0.06 0.00 -0.05 - 0.05
Variance components
Stand-level Stdev 0.10 0.06 — 0.14 0.16 0.10 —- 0.22 0.14 0.09 — 0.19

Gap-level Stdev 0.16 0.14 — 0.18 0.19 0.17 — 0.22 0.19 0.17 — 0.22
DIC -112.27 -46.81 -54.46
LMM null model

Intercept 0.73* 0.69 — 0.77 0.69" 0.63 — 0.75 0.72* 0.67 — 0.77
Stand-level Stdev 0.12 0.08 — 0.17 0.16 0.12 — 0.21 0.14 0.10 — 0.18
Gap-level Stdev 0.17 0.15 — 0.19 0.19 0.17 — 0.22 0.19 0.17 — 0.22
DIC -100.26 -50.06 -57.47

N. gaps / stands 182 / 54 180 / 42 201 / 46

 

"indicates 95% credible interval does not overlap 0, * indicates 90% CI does not overlap 0
Stand-level Stdev= standard deviation for the random stand-level intercept

l
H’Overstory for H’Seedlings is Shannon Diversity of overstory trees, H’Overstory for H’Saplings 1-2 and
1-7 m tall is Shannon Diversity of overstory trees and stumps

152

Table E.5. Effects of gap and stand-level variables on log transformed sapling after
harvest growth rates (m/yr): mean parameter estimates and 95% credible
interval from posterior distributions. Positive values indicate increases in
growth rate on the log scale with a one standard deviation increase in that
covariate. AOCa is the reference Habitat Type. A decrease in DIC of 5 or
more indicates improved performance for full model over the null model.

 

 

 

Variable Sugar maple Ironwood
Mean CI Mean CI
LM full model

Intercept -l.28** -I.40—-l.16 -1.26** -l.42—-l.10
Advanced regen -0.40** -0.53 — -0.28 -0.23** -0.38 — -0.08
Sapling height 0.32“ 0.26 - 0.38 0.34” 0.26 - 0.41
Canopy openness -0.03 -0.08 — 0.02 0.06“ -0.01 — 0.13
Competing veg 0.00 -0.04 — 0.05 0.02 -0.03 - 0.07
ATD-Hp -0.14 -0.48 — 0.21 -0.09 -0.24 - 0.06
ATD -0.09 -0.21 — 0.04 0.06 -0.10— 0.23
ATM -0.06 -0.21 — 0.08 -0.07 -0.24 — 0.09
TMC -0.12 -0.33 — 0.08 -0.03 -0.28 — 0.23
Deer density -0.02 -0.09 — 0.05 0.03 -0.03 — 0.08

Time since harvest -0.26** -O.31 — -0.21 -0.10** -0.17 - -0.02
Variance components

Sapling Stdev 0.23 0.21 — 0.27 0.24 0.21 - 0.28
DIC 2.34 13.64

LM null model
Intercept 41.64" -I .70 — -l .57 -1.46* -1.53 — -I .38
Sapling Stdev 0.38 0.34 — 0.43 0.37 0.32 — 0.43
DIC 1 14.86 84.03
N. gaps / stands 127 / 26 98 / 30

 

"indicates 95% credible interval does not overlap 0, * indicates 90% CI
does not overlap 0
Advanced regen= bivariate variable indicating if a variable is advanced regeneration.

153

Literature Cited

154

LITERATURE CITED

Ahlgren, CE. 1976. Regeneration of red pine and white pine following wildfire and
logging in northeastern Minnesota. Journal of Forestry 74 (3): 135-140.

Albert, DA. 1995. Regional landscape ecosystems of Michigan, Minnesota and
Wisconsin: a working map and classification. USDA Forest Service, North
Central Forest Experiment Station General Technical Report 178, St. Paul,
Minnesota.

Aldrich, P.R, G.R. Parker, J. Romero-Severson, and CH. Michler. 2005. Confirmation of
oak recruitment failure in Indiana old-growth forest: 75 years of data. Forest
Science 51 (5): 406-416.

Arbogast, C., Jr. 1957. Marking guides for northern hardwoods under the selection
system. USDA Forest Service, Lake States Forest Experimental Station Paper 56,
St. Paul, Minnesota.

Augustine, DJ. and SJ. McNaughton. 1998. Ungulate effects on the functional species
composition of plant communities: Herbivore selectivity and plant tolerance. The
Journal of Wildlife Management 62 (4): 1165-1183.

Bailey, RE. and R]. Putman. 1981. Estimation of fallow deer (Dama dama) populations
from faecal accumulation. Journal of Applied Ecology 18 (3): 697—702.

Barrett, A. 2003. National Operational Hydrologic Remote Sensing Center Snow Data
Assimilation System (SNODAS) Products at NSIDC. National Snow and Ice Data
Center Special Report 11, Boulder, Colorado.

Bates, D. and M. Maechler. 2009. lme4: Linear mixed-effects models using S4 classes. R
package version 0.9993 75-32. Online: http://cran.r-
project.org/web/packages/lme4/index.htrnl. Accessed: 4 June 2009.

Belden, AC. and S.G. Pallardy. 2009. Successional trends and apparent Acer saccharum
regeneration failure in an oak-hickory forest in central Missouri, USA. Plant
Ecology 204 (2): 305-322.

Bigelow, SW. and CD. Canham. 2002. Community organization of tree species along
soil gradients in a north-eastem USA forest. Journal of Ecology 90 (1): 188-200.

Blouch, RI. 1984. Northern Great Lake states and Ontario forests. pp 391-410. In: L.K.

Halls (ed.), White-Tailed Deer: Ecology and Management. Stackpole Books,
Harrisburg, Pennsylvania.

155

Bolker, B.M., M.E. Brooks, C.J. Clark, S.W. Geange, J .R. Poulsen, M.H.H. Stevens and
J .S. White. 2009. Generalized linear mixed models: A practical guide for ecology
and evolution. Trends in Ecology and Evolution 24 (3): 127-135.

Brokaw, N. and RT. Busing. 2000. Niche versus chance and tree diversity in forest gaps.
Trends in Ecology and Evolution 15 (5): 183-188.

Browne, W.J. and D. Draper. 2006. A comparison of Bayesian and likelihood-based
methods for fitting multilevel models. Bayesian Analysis 1 (3): 473-514.

Burger, T.L. and J. Kotar. 2003. A Guide to Forest Communities and Habitat Types of
Michigan. The Department of Forest Ecology and Management, University of
Wisconsin-Madison.

Bush, RR. and G. Brand. 1993. Lake States (LS) Variant Overview: Forest Vegetation
Simulator. USDA Forest Service, Forest Management Service Center, Fort
Collins, Colorado.

Canavan, SJ. and CW. Ramm. 2000. Accuracy and precision of 10 year predictions for
Forest Vegetation Simulator-Lake States. Northern Journal of Applied Forestry 17
(2), 62-70.

Canham, CD. 1988. Growth and canopy architecture of shade-tolerant trees: Response to
canopy gaps. Ecology 69 (3): 786-795.

Canham, C.D., J.S. Denslow, W.J. Platt, J .R. Runkle, T.A. Spies and RS. White. 1990.
Light regimes beneath closed canopies and tree-fall gaps in temperate and tr0pical
forests. Canadian Journal of Forest Research 20 (5): 620-631.

Canham, CD, and PL. Marks. 1985. The response of woody plants to disturbance:
patterns of establishment and growth. pp. 197-216. In: S.T.A. Pickett and RS.
White (eds.). The Ecology of Natural Disturbance and Patch Dynamics.
Academic Press, Inc, Orland, Florida.

Carleton, T.J., P.F. Maycock, R. Amup and AM. Gordon. 1996. In situ regeneration of
Pinus strobus and P. resinosa in the Great Lakes forest communities of Canada.
Journal of Vegetation Science 7(3): 431-444.

Carmean, W.H. 1972. Site index curves for upland oaks in the central states. Forest
Science 18 (2): 109-120.

Caspersen, J .P. and R.K. Kobe. 2001. Interspecific variation in sapling mortality in
relation to growth and soil moisture. Oikos 92 (1): 160-168.

Clark, J .S. 2007. Models for Ecological Data: An Introduction. Princeton University
Press, Princeton, New Jersey.

156

Clark, J.S., E. Macklin and L. Wood. 1998. Stages and spatial scales of recruitment

limitation in southern Appalachian forests. Ecological Monographs 68 (2): 213-
235.

Clayton, D.G. 1996. Generalized linear mixed models. pp. 276-301. In: Gilks, W.R., S.
Richardson and DJ. Spiegelhalter (eds.). Markov Chain Monte Carlo in Practice.
Chapman & Hall: New York, New York.

Coates, K.D. 2000. Conifer seedling response to northern temperate forest gaps. Forest
Ecology and Management 127: 249-269.

Cole, W.G. and CG. Lorimer. 2005. Probabilities of small-gap capture by sugar maple
saplings based on height and crown growth data from felled trees. Canadian
Journal of Forest Research 35 (3): 643-655.

Collins, B.S., K.P. Dunne and S.T.A. Pickett. 1985. Responses of forest herbs to canopy
gaps. pp. 217-234. In: S.T.A. Pickett and RS. White (eds.). The Ecology of
Natural Disturbance and Patch Dynamics. Academic Press, Inc: Orland, Florida.

Congdon, P. 2003. Applied Bayesian Modeling. John Wiley & Songs, Ltd: New York,
New York.

Cook, B. 2008. Deer depredation on the forests of Michigan. Michigan Society of
American Foresters.

Cote, S.D., T.P. Rooney, J. Tremblay, C. Dussault, and D.M. Waller. 2004. Ecological
impacts of deer overabundance. Annual Review of Ecology, Evolution and
Systematics 35: 113-147.

Crow, T.R., A. Haney, and D.M. Waller. 1994. Report on the scientific roundtable on
biological diversity convened by the Chequamegon and Nicolet National Forests.
USDA Forest Service, North Central Forest Experiment Station General
Technical Report 166, St. Paul, Minnesota.

DeGraaf, R.M., J .B. Hestbeck and M. Yamasaki. 1998. Associations between breeding
bird abundance and stand structure in the White Mountains, New Hampshire and
Maine, USA. Forest Ecology and Management 103: 217-233.

Demaerschalk, J .P. and S.A.Y. Omule. 1978. Stump and breast height diameter tables for
the British Columbia merchantable tree species. Faculty of Forestry Report,
University of British Columbia.

De Stevens, D. 1991. Experiments on mechanisms of tree establishment in old-field
succession: Seedling survival and growth. Ecology 72 (3): 1076-1088.

157

Didier, K.A. and W.F. Porter. 2003. Relating spatial patterns of sugar maple reproductive
success and relative deer density in northern New York State. Forest Ecology and
Management 181 (1-2): 253-266.

Dietze, M. and J.S. Clark. In prep. The role of individual variability and dieback in
sapling response to large forest gaps.

Doepker, R., D.E. Beyer, and M. Donovan. 1994. Deer population trends in Michigan’s
Upper Peninsula. Michigan Department of Natural Resources, Wildlife Division
Report 3254.

Domke, G.M., J .P. Caspersen and TA. Jones. 2007. Light attenuation following selection

harvesting in northern hardwood forests. Forest Ecology and Management: 239:
182-190.

Donovan, G. 2005. Chronic regeneration failure in northern hardwood stands: A liability
to certified forest landowners. Proceedings from Forests & Whitetails: Striving
for Balance. Michigan Society of American Foresters. 9-10 June 2005. St. Ignace,
Little Bear Conference Center.

Donovan, M. 2005. Michigan GAP analysis project. pp 47-49. In: Maxwell et al. (eds).
Gap Analysis Bulletin No. 13 USGS/BRD/Gap Analysis Program. Moscow,
Idaho, USA.

Dovciak, M., PB. Reich and LE. Frelich. 2003. Seed rain, safe sites, competing
vegetation, and soil resources spatially structure white pine regeneration and
recruitment. Canadian Journal of Forest Research 33 (10): 1892-1904.

Dumont, A., J. Oullet, M. Créte, and J. Huot. 2005. Winter foraging strategy of white-
tailed deer at the northern limit of its range. Ecoscience 12 (4): 476-484.

Enders, OK. and D. Tofighi. 2007. Centering predictor variables in cross-sectional
multilevel models: A new look at an old issue. Psychological Methods 12 (2):
121-138.

Eschtruth, AK. and J .J . Battles. 2008. Deer herbivory alters forest response to canopy
decline caused by an exotic insect pest. Ecological Applications 18 (2): 360-376.

Eyre, PH. and W.M. Zillgitt. 1953. Partial cuttings in northern hardwoods of the Lake
States: Twenty-year experimental results. USDA Forest Service, Lake States
Forest Experiment Station Technical Bulletin 1076, St. Paul, Minnesota.

F ei, S. and KC. Steiner. 2008. Relationships between advance oak regeneration and
biotic and abiotic factors. Tree Physiology 28: 1111-1119.

158

 

Felix, A.B., H. Campa, III, K.F. Millenbah, S.R. Winterstein and WE. Moritz. 2004.
Development of landscape-scale habitat-potential models for forest wildlife
planning and management. Wildlife Society Bulltein 32 (3): 795-806.

Filip, Stanley M. 1973. Cutting and cultural methods for managing northern hardwoods
in the Northeastern United States. USDA Forest Service, Northeastern Forest
Experiment Station General Technical Report 5. Upper Darby, Pennsylvania.

Finzi, AC. and CD. Canham. 2000. Sapling growth in. response to light and nitrogen
availability in a southern New England forest. Forest Ecology and Management
131:153-165.

Forsyth, D.M., R.J. Barker, G. Morriss, M.P. Scroggie. 2007. Modeling the relationship
between fecal pellet indices and deer density. The Journal of Wildlife
Management 71 (3): 964-970.

F rawley, B.J. 2005. 2004 quality deer management (QDM) survey: Deer management
units in the Upper Peninsula. Michigan Department of Natural Resources Wildlife
Division Report No. 3432.

Frazer, G.W., C.D. Canham, and KP. Lertzman. 1999. Gap Light Analyzer (GLA),
Version 2.0: Imaging Software to Extract Canopy Structure and Gap Light
Transmission Indices from True-colour Fisheye Photographs, Users Manual and
Program Documentation. Simon Fraser University, Burnaby, British Columbia
and the Institute of Ecosystem Studies, Millbrook, New York.

Frelich, LE. and CG. Lorimer. 1985. Current and predicted long-term effects of deer
browsing in hemlock forests in Michigan, USA. Biological Conservation 34: 99-
120.

Fuller, T.K. 1991. Do pellet counts index white-tailed deer numbers and population
change? Journal of Wildlife Management 55 (3): 393-396.

Garrett, P.W. and RE. Graber. 1995. Sugar maple seed production in northern New
Hampshire. USDA Forest Service, Northeastern Forest Experiment Station
Research Paper 697, Radnor, Pennsylvania.

Gelman, A. 2006. Prior distributions for variance parameters in hierarchical models.
Bayesian Analysis: 1(3): 515-533.

Gelman, A. and J. Hill. 2007. Data Analysis Using Regression and
Multilevel/Hierarchical Models. Cambridge University Press, New York, New
York.

George, LC. and FA. Bazzaz. 1999. The fern understory as an ecological filter:
Emergence and establishment of canopy-tree seedlings. Ecology 80 (3): 83 3-845.

159

 

Goldstein, H. 1995. Multilevel Statistical Models. Edward Arnold, London.

Goovaerts, P. 1997. Geostatistics for Natural Resources Evaluation. Oxfort University
Press: New York, New York.

Graber, R.E., and W.B. Leak. 1992. Seed fall in an old-growth northern hardwood forest.
USDA Forest Service, Northeastern Forest Experiment Station Research Paper
663, Radnor, Pennsylvania.

Graham, SA. 1954. Changes in northern Michigan forests from browsing by deer.
Transactions of the Nineteenth North American Wildlife Conference 19: 526-533.

Gray, A.N., T.A. Spies and M.J. Easter. 2002. Microclimatic and soil moisture responses
to gap formation in coastal Douglas-fir forests. Canadian Journal of Forest
Research 32 (2): 332- 343.

Green, J .C. Ecological features of white pine stands for wildlife. 1992. In: R.A. Stine and
M.J. Baughman (eds.). Proceedings from The White Pine Symposium. 16-18
September 1992. Duluth, Minnesota.

Grisez, T.J. 1960. Slash helps protect seedlings from deer browsing. Journal of Forestry
58 (5): 385-387.

Guldin, J .M. 1996. The role of uneven-aged silviculture in the context of ecosystem
management. Western Journal of Applied Forestry 11 (1): 4-12.

Heezen, KL. and J .R. Tester. 1967. Evaluation of radio-tracking by triangulation with

special reference to deer movements. Journal of Wildlife Management 31 (1):
124-141.

Henne, P.D., F.S. Hu and D.T. Cleland. 2007. Lake-effect snow as the dominant control
of mesic-forest distribution in Michigan, USA. Journal of Ecology 95 (3): 517-
529.

Herrnan, K., M. Joseph, T. Oliver, D. Wagner, H.W. Scullon, J. Ferris and D. Kuhr.
2004. A process for implementing mesic conifer restoration on State land,
Western Upper Peninsula, Michigan. Michigan Department of Natural Resources,
Wildlife Division, Western Upper Peninsula Management Unit.

Hill, HR. 2001. The 2001 deer pellet group surveys. Michigan Department of Natural
Resources Wildlife Report No. 3349. Lansing, Michigan.

Horsley, SB. and DA. Marquis. 1983. Interference by weeds and deer with Allegheny
hardwood reproduction. Canadian Journal of Forest Research 13: 61-69.

160

 

Horsley, S.B., S.L. Stout, and D.S. deCalesta. 2003. White-tailed deer impact on the
vegetation dynamics of a northern hardwood forest. Ecological Applications 13
(1): 98-118.

Hough, AF. 1949. Deer and rabbit browsing and available winter forage in Allegheny
hardwood forests. Journal of Wildlife Management 13 (1): 135-141.

Houle, G. 1992. Spatial relationship between seed and seedling abundance and mortality

in a deciduous forest of northeastern North America. Journal of Ecology 80: 99-
108.

Hughes, J .W. and T.J. F ahey. 1988. Seed dispersal and colonization in a disturbed
northern hardwood forest. Bulletin of the Torrey Botanical Club 115 (2): 89-99.

Jenkins, J. 1997. Hardwood regeneration failure in the Adirondacks: Preliminary studies

of incidence and severity. The Wildlife Conservation Society, Working Paper No.
9.

Jenkins, K.J. and B.F.J. Manly. 2008. A double-observer method for reducing bias in
faecal pellet surveys of forest ungulates. Journal of Applied Ecology 45 (5): 1339-
1348.

Jerome, D.S. 2006 draft. Landforms of the Upper Peninsula, Michigan. USDA Natural
Resource Conservation Service.

Juice, S.M., Timothy J .F ., Thomas G.S., C.T. Driscoll, E.G. Denny, C. Eagar, N.L.
Cleavitt, R. Minocha and AD. Richardson. 2006. Response of sugar maple to
calcium addition to northern hardwood forest. Ecology 87 (5): 1267-1280.

Kelsall, J. P. 1969. Structural adaptations of moose and deer for snow. Journal of
Mammology 50 (2): 302-310.

Kendeigh, SC. 1945. Community selection by birds on the Helderberg Plateau of New
York. The Auk 62 (3): 418-436.

Kittredge, DB. and P.M.S. Ashton. 1995. Impact of deer browsing on regeneration in
mixed stands in southern New England. Northern Journal of Applied Forestry 12
(3): 115-120.

Kobe, R.K. 2006. Sapling growth as a function of light and landscape-level variation in
soil water and foliar nitrogen in northern Michigan. Oecologia 147: 119-133.

Kobe, R.K., G.E. Likens and C. Eagar. 2002. Tree seedling growth and mortality

responses to manipulations of calcium and aluminum in a northern hardwood
forest. Canadian Journal of Forest Research 32 (6): 954-966.

161

Kohn, BE. and J .J . Mooty. 1971. Summer habitat of white-tailed deer in north-central
Minnesota. The Journal of Wildlife Management 35 (3): 476-487.

Kraft, L.S., T.R. Crow, D.S. Buckley, E.A. Nauertz and J .C. Zasada. 2004. Effects of
harvesting and deer browsing on attributes of understory plants in northern

hardwood forests, Upper Michigan, USA. Forest Ecology and Management 199:
219-230.

Krasny, ME. and MC. Whitrnore. 1992. Gradual and sudden forest canopy gaps in

Allegheny northern hardwood forests. Canadian Journal of Forest Research 22
(2): 139-143.

Krueger, L.M. and C]. Peterson. 2006. Effects of white-tailed deer on Tsuga canadensis

regeneration: Evidence of microsites as refugia from browsing. American
Midland Naturalist 156 (2): 353-362.

Krueger, J .A. and K.J. Puettmann. 2004. Growth and injury patterns of eastern white pine
(Pinus strobus L.) seedlings as affected by hardwood overstory density and
weeding treatments. Northern Journal of Applied Forestry 21 (2): 61-68.

Langdon, CA. 2001. A comparison of white-tailed deer population estimation methods
in Western Virginia. Master’s thesis. Division of Forestry, West Virginia
University, Morgantown, West Virginia.

Laurent, B.J., H.J. Shi, D. Gatziolis, J .P. LeBouton, M.B. Walters, and J .G. Liu. 2005.
Using the spatial and spectral precision of satellite imagery to predict wildlife
occurrence patterns. Remote Sensing of Environment 97 (2): 249-262.

Leak, W.B. 1999. Species composition and structure of a northern hardwood stand after

61 years of group/patch selection. Northern Journal of Applied Forestry 16:151—
153.

LeBouton, J .P., B.J. Laurent, M.B. Walters, J .G. Liu. 2005. Forests for dinner: Exploring
a model of how deer affect advance regeneration at stand and landscape scales.
Proceedings from Forests & Whitetails: Striving for Balance. Michigan Society of
American Foresters. 9-10 June 2005. St. Ignace, Little Bear Conference Center.

Liang, S.Y. and SW. Seagle. 2002. Browsing and microhabitat effects on riparian forest
woody seedling demography. Ecology 83 (1): 212-227.

Linden, D.W. 2006. Modeling current and historic habitat for Canada lynx (Lynx
canadensis) in the Upper Peninsula of Michigan. Master’s thesis. Department of
Fisheries and Wildlife, Michigan State University, East Lansing, Michigan.

Marquis, DA. 1975. Seed storage and germination under northern hardwood forests.
Canadian Journal of Forest Research 5: 478-484.

162

 

Marquis, DA. and R. Brenneman. 1981. The impact of deer on forest vegetation in
Pennsylvania. USDA Forest Service, Northeastern Forest Experiment Station
General Technical Report 65, Broomall, Pennsylvania.

Marx, L.M. 2005. Substrate Limitations to T suga canadensis and Betula alleghaniensis
seedling establishment. PhD thesis. Department of Forestry, Michigan State
University, East Lansing, Michigan.

Marx, L.M. and MB. Waters. 2008. Survival of tree seedlings on different species of
decaying wood maintains tree distribution in Michigan hemlock-hardwood
forests. Journal of Ecology 96: 505-513.

McClure, J .W., T.D. Lee and W.B. Leak. 2000. Gap capture in northern hardwoods:
Patterns of establishment and height growth in four species. Forest Ecology and
Management 127: 181-189.

MDNR (Michigan Department of Natural Resources). 2006. Draft 2006 State forest
management plan. Forest, Mineral and Fire Management Division, IC 4600.

MDNR (Michigan Department of Natural Resources). 2008. Weekly winter severity
index report. Wildlife Division.

Metzger, FT. and CH. Tubbs. 1971. The influence of cutting method on regeneration of
second-growth northern hardwoods. Journal of Forestry 69 (9): 559-564.

Miller, R.O. 2004. Regeneration in a heavily browsed northern hardwood stand twelve
years after scarification and fencing. Michigan State University, Upper Peninsula
Tree Improvement Center Research Report.

Millington, J .D.A., M.B. Walters, M.S. Matonis and J .G. Liu. 2009. Effects of local and
regional landscape characteristics on wildlife distribution across managed forests.
In review Forest Ecology and Management.

Mladenoff, DJ. and F. Stearns. 1993. Eastern hemlock regeneration and deer browsing in
the northern Great Lakes region: A re-examination and model simulation.
Conservation Biology 7: 889-900.

N'atarajan, R. and RE. Kass. 2000. Reference Bayesian methods for generalized linear
mixed models. Journal of the American Statistical Association 95 (449): 227-23 7.

National Climatic Data Center. 2009. US. Department of Commerce. Online:
http://www.ncdc.noaa.gov/oa/ncdc.html. Accessed: 11 June 2009.

Neff, DJ. 1968. The pellet-group count technique for big game trend, census, and
distribution: A review. Journal of Wildlife Management 32 (3): 597-614.

163

 

Nobis, M. 2005. SideLook 1.1- Imaging software for the analysis of vegetation structure
with true-colour photographs. Online: http://www.appleco.ch. Accessed: 5
January 2009.

NOHRSC (National Operational Hydrologic Remote Sensing Center). 2004. Snow Data
Assimilation System (SNODAS) Data Products at NSIDC. Boulder, Colorado
USA: National Snow and Ice Data Center.

Ntzoufras, I. 2009. Bayesian Modeling Using WinBUGS. John Wiley & Sons, Inc:
Hoboken, New Jersey.

Nyland, RD. 1998. Selection system in northern hardwoods. Journal of Forestry 96 (7):
18-21.

OMNR. 1998. A silvicultural guide for the tolerant hardwood forest in Ontario. Ontario
Ministry of Natural Resources. Queen’s Printer for Ontario: Toronto, Ontario.

OMNR. 2003. Old growth policy for Ontario’s Crown forests, v1. Ontario Ministry of
Natural Resources: Toronto, Ontario.

Ozoga, J .J . 1968. Variations in microclimate in a conifer swamp deeryard in northern
Michigan. Journal of Wildlife Management 32 (3): 574-585.

Ozoga, J .J and L.W. Gysel. 1972. Response of white-tailed deer to winter weather.
Journal of Wildlife Management 36 (3): 892-896.

Ozoga, J .J . and L.J. Verme. 1970. Winter feeding patterns of penned white-tailed deer.
Journal of Wildlife Management 34 (2): 431-439.

Patmos, W. 1995. New Hampshire’s native trees, shrubs and vines with wildlife value.
University of New Hampshire, Cooperative Extension Report 449.

Pauley, G.R., J .M. Peek, P. Zager. 1993. Predicting white-tailed deer habitat use in
northern Idaho. Journal of Wildlife Management 57: 904-913.

Plummer, M., N. Best, K. Cowles and K. Vines. 2009. Coda: Output analysis and
diagnostics for MCMC. R package version 0.13-4. Online: http://cran.r-
project.org/web/packages/coda/index.htrnl. Accessed: 10 August 2009.

Pokharel, B. and RE. Froese, 2008. Evaluating alternative implementations of the Lake

States FVS diameter increment model. Forest Ecology and Management 255 (5-
6): 1759-1771.

Poole, K.G. and G. Mowat. 2005. Winter habitat relationships of deer and elk in the

temperate interior mountains of British Columbia. Wildlife Society Bulletin 33
(4): 1288-1302.

164

 

Potter-Witter, K. and J .T. Lacksen. 1993. The status of the maple-birch forest type in
Michigan. Michigan State University Department of Forestry Research Report
533. East Lansing, Michigan.

Qian, SS. and Z. Sheri. 2007. Ecological applications of multilevel analysis of variance.
Ecology 88 (10): 2489—2495.

R Development Core Team. 2009. R: a language and environment for statistical
computing. R Foundation for Statistical Computing, Vienna, Austria. Online:
http://www.R-project.org. Accessed: 10 January 2009.

Racevskis, LA. and F. Lupi. 2006. Comparing urban and rural perceptions of and
familiarity with the management of forest ecosystems. Society and Natural
Resources 19(6): 479-495.

Randall, J .A. 2007. Deer and sedge effects on tree seedling dynamics in northern
temperate forests. PhD Dissertation. Department of Forestry, Michigan State
University, East Lansing, Michigan.

Ray, D.G., R.D. Nyland, and RD. Yanai. 1999. Patterns of early cohort development
following Shelterwood cutting in three Adirondack northern hardwood stands.
Forest Ecology and Management 119 (1-3): 1-11.

Ribbens, E., J.A. Silander, Jr., S.W. Pacala. 1994. Seedling recruitment in forests:

Calibrating models to predict patterns of tree seedling dispersion. Ecology 75 (6):
1 794-1 806.

Romagosa, MA. and DJ. Robison. 2003. Biological constraints on the growth of

hardwood regeneration in upland Piedmont forests. Forest Ecology and
Management 175: 545-561.

Rooney, T.P. 2001. Deer impacts on forest ecosystems: A North American perspective.
Forestry 74 (4): 201-208

Rooney, T.P, R.J. McCormick, S.L. Solheim, and D.M. Waller. 2000. Regional variation
in recruitment of hemlock seedlings and saplings in the upper Great Lakes, USA.
Ecological Application 10: 1119-1 132.

Rooney, T.P., S.L. Solheim, D.M. Waller. 2002. Factors affecting the regeneration of
northern white cedar in lowland forests of the Upper Great Lakes region, USA.
Forest Ecology and Management 163 (1): 119-130.

Rooney, T.P. and D.M. Waller. 2003. Direct and indirect effects of white-tailed deer in
forest ecosystems. Forest Ecology and Management 181 (1-2): 165-176.

165

Runkle, J .R. 1981. Gap regeneration in some old-growth forests of the eastern United
States. Ecology 62 (4): 1041-1051.

Runkle, J .R. 1982. Patterns of disturbance in some old-growth mesic forests of eastern
North America. Ecology 63 (5): 1533-1546.

Runkle, J .R. 1984. Development of woody vegetation in treefall gaps in a beech-sugar
maple forest. Holarctic Ecology 7: 157-164.

Runkle, J .R. 1985. Disturbance regimes in temperate forests. pp 17-33. In: S.T.A. Pickett
and RS. White (eds.). The ecology of natural disturbance and patch dynamics.
Academic Press, New York, New York.

Sage, R.W. Jr., W.F. Porter, and H. B. Underwood. 2003. Windows of opportunity:
White-tailed deer and the dynamics of northern hardwood forests of the
northeastern US. Journal for Nature Conservation 10 (4): 213-220.

Saunders, MR. and K.J. Puettmann. 1999. Use of vegetational characteristics and
browsing patterns to predict deer damage in eastern white pine (Pinus strobus)
plantations. Northern Journal of Applied Forestry 16 (2): 96-102.

Schreeg, L.A., R.K. Kobe and MB. Walters. 2005. Tree seedling growth, survival and
morphology in response to landscape-level variation in soil resource availability
in northern Michigan. Canadian Journal of Forest Research 35 (2): 263-273.

Schumann, M.E., A.S. White, J .W. Witham. 2003. The effects of harvest-created gaps on
plant species diversity, composition and abundance in a Maine oak-pine forest.
Forest Ecology and Management 176: 543-561.

Shafer, E.L., T.J. Grisez and E. Sowa. 1961. Results of deer exclosure studies in
northeastern Pennsylvania. USDA Forest Service Northeastern Forest Experiment
Station Research Note 121, Upper Darby, Pennsylvania.

Shi, H., E.J. Laurent, J. LeBouton, L. Racevskis, K.R. Hall, M. Donovan, R.V. Doepker,
M.B. Walters, F. Lupi, J .G. Liu. 2006. Local spatial modeling of white-tailed deer
distribution. Ecological Modelling 190: 171-189.

Shibata, M. and T. Nakashizuka. 1995. Seed and seedling demography of four co-
occurring Carpinus species in a temperate deciduous forest. Ecology 76 (4):
1099-1108.

Shields, J .M. and CR. Webster. 2007. Ground-layer response to group selection with
legacy-tree retention in a managed northern hardwood forest. Canadian Journal of
Forest Resources 37 (10): 1797-1807.

166

Shields, J .M., C.R. Webster, and L.M. Nagel. 2007. Factors influencing tree species

diversity and Betula alleghaniensis establishment in silvicultural openings.
Forestry 80 (3): 293-307.

Smith, D.M. 1951. The influence of seedbed conditions on the regeneration of eastern
white pine. Connecticut Agricultural Experiment Station Bulletin 545, New
Haven, Connecticut.

Smith, D.M., B.C. Larson, M.J. Kelty, P.M.S. Ashton. 1997. The Practice of
Silvilculture: Applied Forest Ecology 9th ed. John Wiley & Sons, Inc: New York,
New York.

Smith, W.B., P.D. Miles, J .S. Vissage and SA. Pugh. 2004. Forest Resources of the
United States, 2002. USDA Forest Service, North Central Research Station
General Technical Report 241, St Paul, Minnesota.

Spiegelhalter, D.J., N.G. Best, B.P. Carlin and A. Linde. 2002. Bayesian measures of
model complexity and fit. Journal of the Royal Statistical Society, Series B
(Statistical Methodology) 64(4): 583-639.

Spiegelhalter, D.J., A. Thomas, N.G. Best, W.R. Gilks and D. Lunn. 2003. BUGS:
Bayesian inference using Gibbs sampling. MRC Biostatistics Unit, Cambridge,
England. Online: www.mrc-bsu.cam.ac.uk/bugs. Accessed: 15 July 2009.

Spiegelhalter, D.J., K.R. Abrams and J .P. Myles. 2004. Bayesian approaches to clinical
trails and health-care evaluation. John Wiley & Sons, Ltd, West Sussex, England.

Stearns, F.W. 1997. History of the Lake States forests: Natural and human impacts. pp 8-
29.1n: J .M. Vasievich and H.H. Webster (eds.). Lake States Regional Forest
Resources Assessment: Technical Papers. USDA Forest Service, North Central
Forest Experiment Station General Technical Report 189, St. Paul, Minnesota.

Stoeckeler, J .H., R.O. Strothmann, and L.W. Krefting. 1957. Effect of deer browsing on
reproduction in the northern hardwood-hemlock type in northeastern Wisconsin.
Journal of Wildlife Management 21 (1): 75-80.

Sturtz, S., U. Ligges, and A. Gelman. 2005. R2WinBUGS: A package for running
WinBUGS from R. Journal of Statistical Software, 12(3), 1-16.

Swifi, E. 1949. Wisconsin’s deer damage to forest reproduction survey: Final report.
Wisconsin Conservation Department Publication 347.

Swihart, R.K. and RH. Yahner. 1983. Browse preferences of jackrabbits and cottontails
for species used in shelterbelt plantings. Journal of Forestry 81 (2): 92-94.

167

Telfer, ES. 1972. Browse selection by deer and hares. Journal of Wildlife Management
36(4): 1344-1349.

Thome, S.G. 1992. White pine: Challenges for the future. In: R.A. Stine and M.J.
Baughman (eds.). Proceedings from the White Pine Symposium. 16-18 September
1992. Duluth, Minnesota.

Tierson, W.C., G.F. Mattfeld, R.W. Sage, Jr., D.F. Behrend. 1985. Seasonal movements
and home ranges of white-tailed deer in the Adirondacks. Journal of Wildlife
Management 49 (3): 760-769.

Todd, J .B. 1927. Winter food of cottontail rabbits. Journal of Mammalogy 8 (3): 222-
228.

Tubbs, CH. 1968. The influence of residual stand densities on regeneration in sugar
maple stands. USDA Forest Service, North Central Forest Experiment Station
Research Note 47, St. Paul, Minnesota.

Tubbs, CH. 1977. Natural regeneration of northern hardwoods in the northern Great
Lakes region. USDA Forest Service, North Central Forest Experiment Station
Research Paper 150, St. Paul, Minnesota.

VanDeelen, T.R., H. Campa, III, M. Hamady, J .B. Haufler. 1998. The Journal of Wildlife
Management 62 (l ): 205-213.

Verme, L.J. 1965. Swamp conifer deeryards in northern Michigan: Their ecology and
management. Journal of Forestry 63 (7): 523-529.

Verme, L.J. 1973. Movements of white-tailed deer in upper Michigan. Journal of
Wildlife Management 37 (4): 545-552.

Walters, MB. and RB. Reich. 1997. Growth of Acer saccharum seedlings in deeply
shaded understories of northern Wisconsin: Effects of nitrogen and water
availability. Canadian Journal of Forest Research 27 (2): 237-247.

Ward, J .S. and GR. Stephens. 1995. Protection of tree seedlings from deer browsing. In:
K.W. Gottschalk and S.L.C. Fosbroke (eds.). Tenth Central Hardwood Forest
Conference Proceedings. 5-8 March 1995. Morgantown, West Virginia. USDA
Forest Service, Northeastern Forest Experiment Station General Technical Report
197, Radnor, Pennsylvania.

Webb, W.L., D.F. Behrend, and B. Saisom.1977. Effect of logging on songbird
populations in a northern hardwood forest. Wildlife Monographs 55: 1-35.

168

Webster, CR. and CG. Lorimer. 2002. Single-tree versus group selection in hemlock-
hardwood forests: Are smaller opening less productive? Canadian Journal of
Forest Research 32 (4): 591-604.

Webster, CR. and CG. Lorimer. 2005. Minimum opening sizes for canopy recruitment
of midtolerant tree species: A retrospective approach. Ecological Applications 15
(4): 1245-1262.

Weisberg, P.J. and H. Bugmann. 2003. Forest dynamics and ungulate herbivory: From
leaf to landscape. Forest Ecology and Management 181: 1-12.

Whitney, G.G. 1987. An ecological history of the Great Lakes forest of Michigan.
Journal of Ecology 75: 667-684.

Wilensky, U. 1999. NetLogo. Center for Connected Learning and Computer-Based
Modeling, Northwestern University. Evanston, Illinois. Online:
http://ccl.northwestem.edu/netlogo/. Accessed: 16 September 2009.

Williamson, SJ. and DH. Hirth. 1985. An evaluation of edge use by white-tailed deer.
Wildlife Society Bulletin 13 (3): 252-257.

Yamamoto, S. 2000. Forest Gap dynamics and tree regeneration. J oumal of Forest
Research 5: 223-229.

Yawyney, H.W. and CM. Carl, Jr. 1970. A sugar maple planting study in Vermont.
USDA Forest Service, Northeastern Forest Experiment Station Research Paper
175, Upper Darby, Pennsylvania.

Zhao, Y., J. Staudenmayer, B.A. Coull and MP. Wand. 2006. General design Bayesian
generalized linear mixed models. Statistical Science 21 (1): 35-51.

169

 

M11111111111111111111111111111111111“