EVALUATING HOW STAND LEVEL FOREST DYNAMICS INFLUENCE
CARIBOU WINTER RANGE

by
Erin Colleen Hall
BSc., University of Northern British Columbia, 1999

THESIS SUBMITTED IN PARTIAL FULFILLMENT OF
THE REQUIREMENTS FOR THE DEGREE OF
MASTER OF SCIENCE
IN
NATURAL RESOURCES AND ENVIRONMENTAL STUDIES

UNIVERSITY OF NORTHERN BRITISH COLUMBIA
August 2023
© Erin Hall, 2023

Abstract
Forest ecosystems and the ecological services they provide are complex and dynamic.
Disturbance and successional processes impact habitat for animals such as caribou (Rangifer
tarandus caribou Gmelin 1788, Banfield 1961, 1974) and moose (Alces alces L.). The
forested winter habitat of the Tweedsmuir-Entiako caribou herd has been significantly
affected by recent disturbance. In order to evaluate the mid- and long-term consequences for
habitat, we require an understanding of how habitat attributes change as stand structure
changes over time.
I developed a framework to link a stand dynamics model (SORTIE-ND), via linker functions,
to three ecosystem elements relevant to caribou populations: terrestrial forage lichen, moose
forage, and vertical cover. I first used SORTIE-ND to model stand structure dynamics
following mountain pine beetle (Dendroctonus ponderosae Hopkins) attack and clearcut
harvesting. I developed linker functions and refined them with empirical data. These
functions were applied to stand structure output from SORTIE-ND to project the response of
habitat to the disturbance agents. I addressed three questions: 1. What is the influence of
edaphic site on the provision of elements of caribou habitat over time? 2. What is the
influence of stand type on the provision of elements of caribou habitat over time? 3. How do
recovery trajectories of caribou habitat elements differ in response to disturbances such as
mountain pine beetle and clearcut harvest?
I found that site and stand type interacted to influence the suitability of terrestrial lichen,
moose forage and vertical cover. Poor productivity sites were most favourable for terrestrial
forage lichen. Higher productivity stands and broadleaf stands provided the greatest
suitability for moose shrubs. Both mountain pine beetle attack and clearcut harvest improved
stand suitability for lichen and moose browse though for different sites and stand types.
ii

Interestingly, I found that severe mountain pine beetle attack resulted in conditions that
favoured forage lichen for longer than clearcuts. Clearcuts resulted in conditions that were
more suitable than mountain pine beetle attack for winter forage for moose. When operating
on landscapes with caribou winter habitat, forest managers should consider disturbance
history, the distribution of sites and stand types, silviculture systems and the future desired
state of ecosystem services. The SORTIE-ND modelling framework could be used to support
the evaluation of alternative management decisions in caribou winter habitat including
alternative silviculture strategies. This approach could also provide stand- and site-specific
trajectories of the recovery of caribou habitat that could be applied to habitat supply models.

iii

Table of Contents

Abstract

ii

Table of Contents

iv

List of Tables

vi

List of Figures

vii

Acknowledgements

x

Chapter 1

1

Introduction

1

1.1 Background

1

1.1.1 Forest Ecosystem Dynamics

1

1.1.2 Caribou Forested Winter Habitat

3

1.1.3 Effects of Mountain Pine Beetle and Timber Harvesting on Caribou Winter Habitat
8
1.2 Thesis Objectives and Structure

10

Chapter Two

13

Developing Linker Functions to Enable Projection of Elements of Caribou Winter
Habitat

13

2.1 Introduction

13

2.2 Methods

17

2.2.1 Study Area and Caribou Habitat Context

17

2.2.2 SORTIE-ND

20

2.2.3 Linker Functions

21

2.2.4 Model Description

22

2.2.5 Field Data Collection

29

2.2.6 Linker Function Fitting and Calibration

35

2.3 Results

36

2.3.1 Lichen

36

2.3.2 Moose Shrubs

39

2.3.3 Moose Cover

43

2.4 Discussion

43

Chapter Three

50

Modelling the Effects of Disturbance on Caribou Winter Habitat

50

3.1 Introduction

50

3.2 Methods

54

3.2.1 Study Area and Caribou Habitat Context

54
iv

3.2.2 Stand Modelling
3.3 Results

55
60

3.3.1 Component Linker Functions

60

3.3.2 Terrestrial Lichen Dynamics

64

3.3.3 Moose Shrub Dynamics

66

3.3.4 Vertical Cover Dynamics

68

3.4 Discussion

70

3.5 Conclusions

78

Chapter Four

79

Conclusions

79

4.1 Main Findings and Management Implications

80

4.2 Contributions

82

4.3 Next Steps and Future Considerations

83

4.4 Concluding Remarks

86

Literature Cited

87

Personal Communication

114

Appendix 1

115

Appendix 2

119

Appendix 3

122

Appendix 4

124

Appendix 5

125

Appendix 6

137

Appendix 7

145

v

List of Tables
Table 1. Structure and evaluation of linker functions developed to project the relative value
of a forest stand to support forage lichens.
36
Table 2. Structure and evaluation of linker functions developed to project the relative value
of a forest stand to support forage shrubs used by moose.
37
Table 3. Parameter values for linker functions developed to project the relative value of a
forest stand to support forage lichen in west-central BC
39
Table 4. Parameter values for linker functions developed to project the relative value of a
forest stand to support forage shrubs used by moose in west-central BC
41
Table 5. Structure of stands used to initiate SORTIE-ND model scenarios for mature stands
across gradients of edaphic site and stand type.
57
Table 6. Stand establishment characteristics for sub-boreal clearcut scenarios modelled with
SORTIE-ND including planting densities (in stems/ha) and natural regeneration (‘seed’). 58

vi

List of Figures
Figure 1. The range of the Tweedsmuir-Entiako caribou population (red outline) located in
west-central British Columbia, Canada.
5
Figure 2. Distribution of empirical plots (black triangles) established within the Sub-Boreal
Spruce moist cold Babine variant (SBSmc2) (purple) during the summer of 2019 in westcentral British Columbia, Canada.
18
Figure 3. Linker function components that represent how the suitability of a stand to support
forage lichen varies as a function of stand state
38
Figure 4. Observed lichen (a) and moose shrub (b) suitability estimated from empirical plots
(n = 76) vs. predicted suitability derived from developed linker functions.
40
Figure 5. Linker function components that represent how the suitability of a stand to support
moose shrubs for forage varies as a function of stand state
41
Figure 6. Distribution of species of moose shrubs across soil moisture and nutrient gradients.
42
Figure 7. The linker function that represents the effect of conifer basal area (trees taller than
10 m) on the suitability of a stand to provide vertical cover in winter.
43
Figure 8. Output of the component linker functions for forage lichen on subxeric sites.

61

Figure 9. Output of the component linker functions for moose browse on subhygric sites.

62

Figure 10. The projected response of lichen suitability over time in undisturbed mature forest
(1st column), after MPB disturbance (2nd column) and after clearcut (3rd column) across four
edaphic sites (rows) and four stand types (colours).
65
Figure 11. The projected response of moose shrub suitability over time in undisturbed mature
forest (1st column), after MPB disturbance (2nd column) and after clearcut (3rd column)
across four edaphic sites (rows) and four stand types (colours).
67
Figure 12. The projected response of suitability to provide vertical cover over time in
undisturbed mature forest (1st column), after MPB disturbance (2nd column) and after
clearcut (3rd column) across four edaphic sites (rows) and four stand types (colours).

69

Figure 13. Comparison of the linker functions for lichen parameterized with the calibrated
parameters (blue) versus parameterized with the maximum likelihood estimates (black). 122
Figure 14. Comparison of the linker functions for moose shrubs used for forage
parameterized with the calibrated parameters (blue) versus parameterized with the maximum
likelihood estimates (black).
123
Figure 15. Distribution of the abundance of forage lichen, moose shrubs, aspen (0.5 - 3 m in
height), subalpine fir (0.5 - 3 m in height), and conifer basal area in the empirical dataset (n =
76).
124

Figure 16. Output of the component linker functions for forage lichen in three different stand
types (colours) on subxeric sites.
125
Figure 17. Output of the component linker functions for forage lichen in four different stand
types (colours) on submesic sites.
126
Figure 18. Output of the component linker functions for forage lichen in four different stand
types (colours) on mesic sites.
127
Figure 19. Output of the component linker functions for forage lichen in four different stand
types (colours) on subhygric sites.
128
Figure 20. Output of the component linker functions for moose browse in three different
stand types (colours) on subxeric sites.
129
Figure 21. Output of the component linker functions for moose browse in four different stand
types (colours) on submesic sites.
131
Figure 22. Output of the component linker functions for moose browse in four different stand
types (colours) on mesic sites.
133
Figure 23. Output of the component linker functions for moose browse in four different stand
types (colours) on subhygric sites.
135
Figure 24. The projected response of lichen suitability over time in undisturbed mature forest
(1st column), after MPB disturbance (2nd column) and after clearcut (3rd column) across
four edaphic sites (rows) and four stand types (colours).
137
Figure 25. The projected response of moose shrub suitability over time in undisturbed mature
forest (1st column), after MPB disturbance (2nd column) and after clearcut (3rd column)
across four edaphic sites (rows) and four stand types (colours).
138
Figure 26. The projected response of moose browse suitability (includes aspen and subalpine
fir 0.5 – 3 m in height and moose shrubs) over time in undisturbed mature forest (1st
column), after MPB disturbance (2nd column) and after clearcut (3rd column) across four
edaphic sites (rows) and four stand types (colours).
139
Figure 27. The projected response of the suitability of a stand to provide vertical cover over
time in undisturbed mature forest (1st column), after MPB disturbance (2nd column) and
after clearcut (3rd column) across four edaphic sites (rows) and four stand types (colours).140
Figure 28. The projected response of lichen suitability over time in undisturbed mature forest
(1st column), after MPB disturbance (2nd column) and after clearcut (3rd column) across
four edaphic sites (colours) and four stand types (rows).
141
Figure 29. The projected response of moose shrub suitability over time in undisturbed mature
forest (1st column), after MPB disturbance (2nd column) and after clearcut (3rd column)
across four edaphic sites (colours) and four stand types (rows).
142
Figure 30. The projected response of moose browse suitability (includes aspen and subalpine
fir 0.5 – 3 m in height and moose shrubs) over time in undisturbed mature forest (1st
viii

column), after MPB disturbance (2nd column) and after clearcut (3rd column) across four
edaphic sites (colours) and four stand types (rows).
143
Figure 31. The projected response of the suitability of a stand to provide vertical cover over
time in undisturbed mature forest (1st column), after MPB disturbance (2nd column) and
after clearcut (3rd column) across four edaphic sites (colours) and four stand types (rows).144
Figure 32. Comparison of the variability in lichen suitability for four stand types (colours) on
four edaphic sites (rows) as represented by the difference between the 95th suitability
quantile and the 5th suitability quantile at each timestep.
145
Figure 33. Comparison of the variability in lichen suitability between four edaphic sites
(colours) for four stand types (rows) as represented by the difference between the 95th
suitability quantile and the 5th suitability quantile at each timestep.

146

Figure 34. Comparison of the variability in moose shrub suitability for four stand types
(colours) on four edaphic sites (rows) as represented by the difference between the 95th
suitability quantile and the 5th suitability quantile at each timestep.

147

Figure 35. Comparison of the variability in moose shrub suitability between four edaphic
sites (colours) for four stand types (rows) as represented by the difference between the 95th
suitability quantile and the 5th suitability quantile at each timestep.
148
Figure 36. Comparison of the variability in suitability of a stand to provide vertical cover for
four stand types (colours) on four edaphic sites (rows) as represented by the difference
between the 95th suitability quantile and the 5th suitability quantile at each timestep.
149
Figure 37. Comparison of the variability in suitability of a stand to provide vertical cover
between four edaphic sites (colours) for four stand types (rows) as represented by the
difference between the 95th suitability quantile and the 5th suitability quantile at each
timestep.
150

ix

Acknowledgements
There are many people who have been involved and encouraged me along this journey, all of
whom enriched my experience and the final product. I sincerely thank each and every one of
you.
To my supervisory committee, Ché Elkin, Dr. Chris Johnson, and Erica Lilles, thank you for
your support throughout this project. I sincerely appreciate all that I’ve learned from each of
you as well as your time, patience, mentoring and encouragement. Additionally, Ché, thank
you for your enthusiastic support and many contributions that turned an idea and vision for
this project into a reality.
I am very appreciative of the funding I received from the Province of British Columbia’s
Caribou Recovery Program and the Bulkley Valley Research Centre. Further support from
Jocelyn Campbell, Anne-Marie Roberts, and the British Columbia Ministry of Forests was
key to initiating and carrying out this project. The Cheslatta Carrier Nation generously
provided me with barge access to my research sites on the south side of Ootsa Reservoir.
So many people have contributed to this project. Working with Dave Coates and Allen
Banner studying stand dynamics and ecosystem recovery played a large role in inspiring this
project and provided a strong base to build upon. I’ve appreciated the camaraderie, support,
and perspectives provided by members of the Mixed Wood Ecology lab and our informal
Vertebrate Ecology lab group. Weekend sessions with fellow MSc students in Smithers were
invaluable to completing this project and made working on the weekends much more
enjoyable – thank you Nadia Nowak for instigating and Emily Mason and Bryce O’Connor
for continuing to join me. I appreciate the knowledgeable and thoughtful insight that Deb
Cichowski and Karine Pigeon provided in their reviews of my thesis. I have enjoyed the

many insightful discussions I’ve had with numerus people throughout this project including:
Alana Clason, Anne-Marie Roberts, Dave Coates, Sybille Haeussler, Benita Kaytor,
Genevieve Perkins, Jodi Axelson, John Hagen, Dave Wilford, Charlie Canham, Paula
Bartemucci, Roy Rea, Darwyn Coxson, Shannon Crowley, and Dexter Hodder. I appreciated
and enjoyed time spent with Paula Bartemucci honing my lichen and moss identification
skills. Thank you all for readily agreeing to take the time out of your busy schedules to
participate in these discussions and provide valuable perspective, suggestions and insight.
I have great memories from the many weeks of field work across the breadth of the SBSmc2.
Many thanks are due to my field assistants Bronwyn Hobson and Jonathan VanElslander as
well as several others who filled in over the summer: Erin Holtzman, Marie-Lou LeFrançois,
Brandy Hughes, Joanna Lee and Jeremy Roscoe. Ben Steele, Lauren Runge, Alex Ritz, and
Erin Holtzman volunteered to assist with data entry – their contributions sped analysis along
considerably.
Finally, I am forever grateful to my partner in life, Jeremy Roscoe, for his unwavering
support throughout my MSc.

xi

Chapter 1
Introduction
1.1 Background
1.1.1 Forest Ecosystem Dynamics
Ecosystems are dynamic and change over time, driven by shifting environmental conditions,
abiotic and biotic disturbances, competition dynamics, successional processes, and the
interplay between these factors. As ecosystems change, the type and provisioning level of
ecosystem services are expected to shift. For example, biodiversity and habitat quality are
spatially and temporally dynamic ecosystem services that reflect changes in ecosystem
structure and disturbance history (Carroll and Bliss 1982, Halpern et al. 2012, Lilles et al.
2018, Greuel et al. 2021).
Vegetation dynamics often underpin terrestrial ecosystem dynamics. At the landscape level
this is reflected in the distribution and dynamics of vegetation biomes (grasslands, savanna,
forests etc.). At a finer scale, such as a forest stand, vegetation dynamics also drive changes
in the provisioning of ecosystem services through time. For example, the vegetation
community at the stand level determines the amount of forage available to ungulates and, by
extension, the habitat quality of the stand (Crête and Courtois 1997, Keim et al. 2017). Some
ecosystem components, such as lichen communities used by caribou (Rangifer tarandus
Lin.), may require decades or centuries to develop (Maikawa and Kershaw 1976, Boudreault
et al. 2009). Others, such as shrub communities favoured by moose (Alces alces L.) develop
more rapidly following disturbance (Eastman 1977, Julianus et al. 2019).
Forest management, implemented through different silviculture systems, influences forest
composition and structure, the successional trajectory of stands, and the associated provision
1

of ecosystem services (Marzluff et al. 2002, Maleki et al. 2021). Trees play a pivotal role in
forested ecosystems, modifying ecosystems through influencing the availability and
distribution of resources such as light and nutrients and therefore driving understory
vegetation dynamics (Bartemucci et al. 2006, Haughian and Burton 2015). As a result, forest
management activities directly influence the spatial and temporal provisioning of a broad
variety of ecosystem services in forested systems (Marzluff et al. 2002, Mönkkönen et al.
2014, Lilles et al. 2018). To manage for these multiple ecosystem services, we need to
understand the temporal dynamics of forests and ecosystem services, and how recovery
trajectories of these services are affected by management activities and natural disturbance
(Marzluff et al. 2002). Improving our understanding of how stand structure and stand
dynamics influence forest-based ecosystems services will improve our ability to robustly
manage for a range of economic, ecological, and social values (Mönkkönen et al. 2014).
In recent years there has been considerable improvement in the methods for modelling the
dynamics of some ecosystem services. Many of these methods use existing models of forest
dynamics that have a fifty-year history of development (Bugmann and Seidl 2022). These
models enable the projection of stand dynamics including stand structure and tree species
composition over time (Bugmann 2001, Bugmann and Seidl 2022). However, most models
do not explicitly and dynamically project other ecosystem services provided by the forest.
Recently, modellers and ecologists have begun to link ecosystem services to forest dynamics
models via models called linker functions (Elkin et al. 2013, Cordonnier et al. 2014,
Mönkkönen et al. 2014, Lafond et al. 2017). Linker functions enable the projection of
ecosystem service dynamics from stand structure output. Thus, linker functions extend our
ability to dynamically model not only stand structure attributes but also ecosystem services
such as habitat, carbon, and biodiversity (Elkin et al. 2013, Cordonnier et al. 2014,
Mönkkönen et al. 2014, Lafond et al. 2017).
2

The overarching aim of this thesis is to improve our understanding of how forest structure
and composition affect key habitat elements for caribou, an ungulate of conservation concern,
and how these habitat elements recover following mountain pine beetle (MPB)
(Dendroctonus ponderosae Hopkins) outbreaks and clearcut harvesting. In the following
section I provide an overview of the forest habitat elements that are important for caribou
winter habitat. I then describe the agents of habitat disturbance evaluated in this thesis and
their impact on forest structure, and by extension caribou habitat. I conclude Chapter 1 by
describing the thesis chapters and how they relate to the thesis objectives.
1.1.2 Caribou Forested Winter Habitat
Caribou are a species of conservation concern with an uncertain future in Canada (FestaBianchet et al. 2011, Johnson et al. 2015). As of 2011, most caribou populations in Canada
were listed under the Species at Risk Act as either Endangered, Threatened or Special
Concern. Predation and, for some populations, overharvest are proximate causes of decline
(Festa-Bianchet et al. 2011). However, habitat change is the primary cause of decline for
many populations of caribou (Festa-Bianchet et al. 2011). Activities such as forestry and oil
and gas exploration and development create early seral forests that provide habitat for the
apparent competitors of caribou, namely moose and deer (Odocoileus spp.). These apparent
competitors support larger and more abundant populations of wolves (Canis lupus L.) and
bears (Ursus spp.) that ultimately lead to unsustainable predation for caribou (Wittmer et al.
2007, Festa-Bianchet et al. 2011).
Although only one species of caribou is recognized throughout its global distribution, caribou
exhibit a wide range of migratory and movement behaviours in addition to morphological and
genetic variability (COSEWIC 2011, Festa-Bianchet et al. 2011). In North America, four
extant subspecies were described by Banfield: woodland (Rangifer tarandus caribou Gmelin

3

1788, Banfield 1961, 1974), Grant’s caribou (R. t. granti Allen, 1902), barren-ground caribou
(R. t. groendlandicus L. 1767), and Peary caribou (R. t. pearyi Allen, 1902) (COSEWIC
2011). Due to challenges with this classification and availability of improved information,
COSEWIC proposed twelve Designatable Units to classify caribou populations within
Canada (COSEWIC 2011). These Designatable Units were developed using available
information regarding caribou phylogenetics; genetic diversity and structure; morphology;
movements, behaviour and life history strategies; and distribution (COSEWIC 2011).
Designatable Units are intended to recognize discrete and evolutionarily significant
populations of species in Canada that are considered irreplaceable (COSEWIC 2009, 2011).
The Tweedsmuir-Entiako caribou herd (TEC) is the focus of this project. These caribou are
classified as the northern mountain ecotype of woodland caribou and are included in the
Northern Mountain Designatable Unit (COSEWIC 2002, 2011, 2014). The herd’s home
range is located in west-central British Columbia (BC), east of the Coast Mountains, and
southwest of Prince George (Figure 1).

4

Figure 1. The range of the Tweedsmuir-Entiako caribou population (red outline) located in westcentral British Columbia, Canada.

As with many northern mountain caribou, TEC migrate annually in spring and fall. From
December to the onset of spring migration in early May, the majority of winter is spent in the
east in low-elevation lodgepole pine (Pinus contorta var. latifolia Douglas ex Loudon) forests
in the Entiako Park area south of Tetachuk Lake. However, during some years, caribou may
also move to alpine and subalpine areas in the Fawnie Mountains or reside in the Cheslaslie
River area (northwest of Entiako Park and the core of the winter range) (Cichowski 1993,
Steventon 1996). In late winter the caribou congregate near the Entiako River where they
inhabit mature lodgepole pine and mature mixed pine-interior spruce (Picea glauca (Moench)
Voss) forests on mesic sites. By early May the caribou cross Tetachuk Lake and migrate
west to the mountains in the northern portion of Tweedsmuir Park and mountains further to

5

the west and northwest. By late November almost all of the caribou return to the winter range
in the Entiako Park area (Cichowski 1989, 1993).
The TEC herd, and specifically its forested winter habitat, was selected as the focus of this
study. Forested winter habitat was selected for two reasons. First, the herd is in decline, bluelisted (designated as threatened by the Province of British Columbia), listed as Threatened
under the federal Species at Risk Act (COSEWIC 2014), and habitat management in the
herd’s range is a provincial priority. Second, the herd’s forested winter habitat has been
subject to significant levels of recent disturbance from wildfire, mountain pine beetle and
timber harvest. For simplicity, I will refer to forested winter habitat as winter habitat for the
remainder of the thesis. Winter habitat is important because forage availability, diversity and
quality are most limited in winter. In summer, the diet of northern mountain caribou reflects
that of a generalist herbivore consuming deciduous shrubs, forbs, lichens, and mushrooms
(Denryter et al. 2017). In the winter, caribou transition to a specialized diet consisting
primarily of lichen (Cichowski 1993, Johnson et al. 2000). However, forage comprises only
one component of habitat. Habitat is defined as “the resources and conditions present in an
area that produce occupancy, which may include survival and reproduction by a given
organism” (Hall et al. 1997 p. 175).
TEC winter range is characteristic of northern mountain caribou and reflects prevailing
topography, vegetation, and weather of their range (COSEWIC 2011, 2014). Winter ranges
of northern mountain caribou are generally gentle with rolling mountains interspersed with
lower elevation plateaus (COSEWIC 2011). Winters are colder and the snow pack generally
lower than conditions typical of neighbouring ranges of central and southern mountain
caribou to the east and southeast (COSEWIC 2011). In winter, northern mountain caribou
feed primarily on terrestrial lichens found either in low-elevation mature conifer forests

6

(lodgepole pine and black spruce (Picea mariana (Mill.) Britton, Sterns & Poggenb.)) or on
high windswept alpine slopes (COSEWIC 2011). Arboreal lichens are also consumed during
winter especially where they are more abundant such as in mixed stands of pine and spruce
and along the fringes of wetlands. These lichen are often sourced later in the winter when
snow pack is deeper or less amenable to cratering (Cichowski 1993, 1993, Johnson et al.
2001, 2003). In addition to the availability of lichen, an important characteristic of caribou
winter habitat is the canopy condition that influences the availability of terrestrial lichen. A
conifer canopy intercepts snow reducing the depth of the snowpack; it may also protect the
snow from developing as thick a crust as may form in open conditions (Cichowski 2010, Seip
and Jones 2010, Kivinen et al. 2010).
In addition to stand types that provide caribou forage, another important attribute of northern
mountain caribou habitat relates to their primary defense strategy. Similar to all woodland
caribou, northern mountain caribou avoid predation through the use of habitats with a
relatively low density of deer and moose and their shared predators including wolves and
bears (Rettie and Messier 2000, Wittmer et al. 2007, Venier et al. 2014). Landscapes with
habitat that supports abundant prey species have the potential to support larger wolf
populations; this heightens mortality risk for caribou. This relationship, where the population
of one species indirectly results in the decline of another prey species via a common predator,
is termed apparent competition (Holt 1977, Wittmer et al. 2007, Anderson et al. 2018).
Therefore, an important metric of caribou habitat is the presence and quality of moose
habitat. Elements of moose habitat include forage, vertical and security cover, and
abundance of predators (Peek 2007).
In this study, I focus on caribou habitat and stand dynamics at the stand scale. Oliver and
Larson (1990 p. 1) define a stand as “a spatially continuous group of trees and associated

7

vegetation having similar structures and growing under similar soil and climatic conditions”.
Habitat elements that can be assessed at a stand scale include abundance of individual forage
species or basal area of a particular tree species. In contrast, landscape-scale elements of
habitat include patterns of disturbances or networks of roads and linear features. Studies of
caribou habitat are generally situated at the landscape scale (Johnson et al. 2002, DeCesare et
al. 2012, Stockdale et al. 2019, Leblond et al. 2022). However, stand-scale studies and
process-based modelling of stand and habitat dynamics can inform our understanding of
landscape-scale patterns (Holt et al. 1995, Leblond et al. 2022). In addition, process-based
models at the stand-scale allow us to explore the implications of different severities of
disturbance or forest management decisions such as the application of alternative silviculture
systems (Maleki et al. 2021).
1.1.3 Effects of Mountain Pine Beetle and Timber Harvesting on Caribou Winter Habitat
There are a number of interacting reasons that habitat disturbance, either natural or industrial,
is reported as the ultimate cause of caribou population decline. In the previous section I
outlined two important aspects of winter habitat for northern mountain caribou: forage
availability and co-occurrence of moose habitat. Disturbance affects these as well as other
elements of caribou habitat (Festa-Bianchet et al. 2011). Severe disturbances that remove old
forest result in reduced lichen abundance (Cichowski et al. 2022). The resulting younger
stand types provide a greater abundance of winter forage for moose and deer and thereby
increase predation risk via apparent competition (Holt 1977, Wittmer et al. 2007, FestaBianchet et al. 2011, Barber et al. 2018). Development of roads and linear corridors also
increases predation risk by providing wolves with access to caribou habitat and increased
hunting efficiency (Whittington et al. 2011, Dickie et al. 2017, Mumma et al. 2018).

8

Over the last two decades disturbance within TEC winter habitat has been significant with
several large wildfires and the largest MPB epidemic documented in Canadian history (Dhar
et al. 2016, Cichowski and Haeussler 2020, Province of British Columbia 2023). Extensive
salvage harvest of standing dead trees has occurred following the wildfire and MPB events.
MPB attack affects winter habitat through killing pine trees and thus altering the structure and
species composition of the overstory. The degree to which MPB attack affects the provision
of caribou habitat will depend upon the intensity of the attack and the proportion of pine trees
killed (Nobert et al. 2020, Cichowski et al. 2022). Where pine are killed, microclimate and
light levels change gradually over several years as the dead trees slowly lose needles (needle
loss occurs for up to five years following attack) and then fall (Cichowski and Haeussler
2020). Monitoring plots within the TEC winter range have documented that tree fall occurs
gradually. Few trees fell within the first ten years following MPB attack, but after an
additional five years, 26% of MPB attacked trees had fallen and fall rates had increased
(Cichowski and Haeussler 2020).
The most immediate ecological impact of MPB attack is that, with the death of pine trees, the
availability of water and nutrients increases for other plants. A study in the TEC winter range
found a significant increase in kinnikinnick cover (Arctostaphylos uva-ursi (L.) Spreng.) for
the first several years following MPB attack (Cichowski and Haeussler 2020). A similar
response in dwarf Vaccinium species was observed in central BC following MPB attack (Seip
and Jones 2010). In contrast, abundance of terrestrial lichen declined significantly. After
about ten years, these trends reversed on the poorest productivity sites; kinnikinnick
abundance declined and lichen decline stopped (in the Tweedsmuir-Entiako lichen
monitoring plots, Cichowski and Haeussler 2020). When considering these dynamics it is

9

important to keep in mind that lichen are slow growing and poor competitors in comparison
to other ground vegetation and mosses (Rowe 1984).
Timber harvest can have significantly different effects on caribou winter habitat, depending
on the silviculture system and the timing of harvest. The effects of silviculture systems could
be considered along an intensity gradient, from low intensity single-tree selection to
gradually higher intensity group selection, shelterwoods, clearcuts with dispersed retention or
reserves, to high-intensity clearcuts. Clearcuts completely remove the mature trees and
associated vegetation structure that caribou often select (Joly et al. 2003, Fortin et al. 2008).
Further, caribou avoid young and immature cutblocks (Fisher and Wilkinson 2005,
Metsaranta and Mallory 2007, Courtois et al. 2008, Courbin et al. 2009). Partial harvest
systems provide the option to extract revenue from timber yet also retain some stand
structure. Although there is extensive information regarding caribou interactions with
clearcuts, there is limited information regarding caribou interactions with partial harvest
systems (Fisher and Wilkinson 2005, Courtois et al. 2007, Waterhouse et al. 2011). Some
partial harvest systems are being tested and evaluated to assess the degree to which they
retain elements of caribou habitat (Waterhouse et al. 2011, Nadeau Fortin et al. 2016).
Harvesting methods also play a role in the degree to which winter habitat is affected.
Harvesting on a thick snowpack can minimize the physical disturbance to terrestrial lichen
associated with the machinery necessary to cut and remove trees (Webb 1998, Coxson and
Marsh 2001). Harvest methods that limit the distribution of logging debris can also minimize
the direct impacts to lichen (Kivinen et al. 2010, Cichowski et al. 2022)
1.2 Thesis Objectives and Structure
In this thesis I evaluated how elements of caribou winter habitat changed through time in
response to stand dynamics following disturbance events. This study addresses the need for a
10

greater understanding of the recovery and successional dynamics of winter habitat for
northern mountain caribou. Researchers have developed the capability to model the dynamics
of stand structure from a tree perspective. However, we lack the corresponding ability to
model and project how habitat elements within forests respond and change as the stands that
they are embedded within regenerate, mature, and respond to disturbance over time. It is
important to be able to project the influence of disturbance events and proposed treatments on
both present and future habitat values to manage for caribou habitat and prioritize treatments
that favour the recovery of winter habitat (Holt et al. 1995, Mäkelä et al. 2012).
In Chapter 2, I develop a framework that uses the output of the stand dynamics model
SORTIE-ND (Coates et al. 2003) to project the temporal dynamics of caribou habitat. To
develop the framework, I developed linker functions based on literature review and
interviews with subject experts. I then conducted a field study across a broad range of stand
types and collected data representative of the occurrence and abundance of caribou habitat
and associated site and stand structure. I used the empirical data to calibrate and evaluate the
linker functions.
In Chapter 3, I used the framework developed in Chapter 2 to assess the effect of disturbance
on the trajectory of forest stands and caribou habitat. I used SORTIE-ND to model mountain
pine beetle attack and clearcut harvest. I used the resulting model projections to compare the
temporal recovery of caribou winter habitat in disturbed stands to undisturbed mature stands.
Based on those comparisons, I examined three specific questions:
1.

What is the influence of edaphic site on the provision of elements of caribou
habitat over time?

2.

What is the influence of stand type on the provision of elements of caribou habitat
over time?
11

3.

How do recovery trajectories of caribou habitat elements differ in response to
disturbances such as mountain pine beetle and clearcut harvest?

In Chapter 4, I synthesize the results of the thesis and discuss future work and management
implications.

12

Chapter Two
Developing Linker Functions to Enable Projection of Elements of Caribou Winter
Habitat
2.1 Introduction
“The complexity of a forest ecosystem makes difficult any attempt to
synthesize knowledge about forest dynamics or to perceive the implications
of information and assumptions regarding forest growth.” (Botkin et al.
1972a p. 849)
Forest ecosystems are complex and dynamic with interactions across spatial and temporal
scales (Botkin et al. 1972a, Filotas et al. 2014). Many processes influence ecosystem
dynamics including disturbance, competition dynamics, and seasonal change in resource
availability (Angelstam and Kuuluvainen 2004, Filotas et al. 2014). Models are useful for
testing our understanding of complex ecosystem dynamics and for providing a framework
within which to set existing knowledge (Botkin et al. 1972b, Filotas et al. 2014). Models also
enable us to visualize and project the effects of ecosystem dynamics and management actions
on ecosystem elements of interest (Mönkkönen et al. 2014, Leblond et al. 2022). These
elements may include indices of economic output, carbon sequestration or habitat for species
of conservation concern (Marzluff et al. 2002, Spies et al. 2007, Mönkkönen et al. 2014).
Caribou are a species of conservation concern in Canada and habitat management is a priority
in many jurisdictions (COSEWIC 2014, Leblond et al. 2014). As presented in Chapter 1, the
ultimate cause of decline is habitat disturbance (Festa-Bianchet et al. 2011). Timber
harvesting, wildfire, and insect outbreaks can change the quality or result in the loss of
habitat for caribou across Canada. As introduced in Chapter 1, winter habitat is an important
component of caribou habitat and is sensitive to disturbance. Two components of the winter
13

habitat for northern mountain caribou that can be assessed at the stand scale are terrestrial
lichen for winter forage and moose habitat that indirectly leads to greater rates of predation
for caribou. Moose habitat at the stand scale consists of forage as well as thermal and
security cover (Peek 2007).
Although there has been much research focused on understanding stand conditions that
provide high-value caribou habitat, the response of these habitat types to the dynamics of the
stands that they are embedded within is not well understood (Cichowski 1993, Johnson et al.
2000, Cichowski and Haeussler 2020, Nobert et al. 2020). Dynamic, stand-scale models of
caribou habitat would support our understanding of the longer-term implications of
disturbance on caribou habitat. Those models could be used to identify current forest types
and management actions that best support development of high-value habitat in the mid- to
long-term. Trade-offs exist between management strategies that retain high-value habitat in
the short-term versus those that promote habitat supply in the long-term (Coxson 2015,
Nobert et al. 2020). Models could be used to support strategic application of management
approaches to balance the short-term and long-term provision of habitat on the landscape.
Habitat models are available for informing management of ungulate habitat. Species
distribution models have a long history of use and inform how landscape- and stand-scale
ecological conditions influence the space use of animals (Metsaranta et al. 2003, Johnson et
al. 2004, Leblond et al. 2014, Whitman et al. 2017). Some recent studies have used models
of vegetation dynamics to forecast habitat dynamics into the future (Barber et al. 2018,
Leblond et al. 2022). The effect of disturbance depends upon the severity of the disturbance.
Some events are severe and stand-replacing, such as a high-intensity wildfire or clearcutting.
However, other disturbances affect only portions of a stand resulting in complex stand
structures (e.g., partial-cut silviculture methods, light to moderate severity fire, and MPB

14

attack). Variability in disturbance types and the resulting stand structure can lead to different
seral stand types (species composition, density, size of trees, etc.) on similar edaphic sites
(Meidinger and Pojar 1991, Banner et al. 1993, Bergeron 2000). Incorporating the ability to
model the dynamics of habitat in response to complex changes in stand structure would add
further value to habitat models. This is especially relevant for habitat elements that develop
over relatively long time frames or that are sensitive to stand structure, such as terrestrial
lichen communities (Holt et al. 1995).
There are a number of forest dynamics models that represent the stand- and landscape-scale
processes related to the tree components of forest stands (Bugmann 2001, Bugmann and Seidl
2022). SORTIE-ND is one such model that has been developed for use in northwest BC
including forests that provide habitat for the TEC. SORTIE-ND is an individual-tree,
spatially explicit model of forest dynamics (Coates et al. 2003, Astrup et al. 2008a). SORTIEND has been developed to model complex stand dynamics including multi-layered and multispecies stands and partial harvest silviculture systems (Canham et al. 2004, Coates et al.
2009, Maleki et al. 2021) (Appendix 1).
To manage for the sustainable provision of wildlife habitat over time, the ability to project
future dynamics of habitat elements and their response to disturbance is important (Holt et al.
1995). In Europe, recent studies have developed frameworks that enable researchers to
extend the results of stand dynamics models to project stand-level ecosystem services
(Bugmann et al. 2017, Lafond et al. 2017, Cristal et al. 2019). Several of these studies have
labelled the functions used to project ecosystem services as “linker functions”. Bugmann et
al. (2017 p. 8) define linker functions as: “algorithms that provide an unbiased and accurate
quantification of ecosystem service provision from forest simulation model outputs”. Linker
functions have been used to relate broad suites of ecosystem services such as vegetation

15

biodiversity and deadwood to stand dynamics resulting from forest management (Lafond et
al. 2015, 2017, Bugmann et al. 2017, Blattert et al. 2017). Less common are the development
of linker functions that relate species-specific habitat elements to stand structure dynamics;
however, some have been developed for species such as the flying squirrel, the Eurasian tree
creeper, and numerous woodpeckers and grouse species (Mönkkönen et al. 2014, Bugmann et
al. 2017).
SORTIE-ND, as a spatially explicit individual-tree model, is well-suited to applying linker
functions to model ecosystem services. Although output from SORTIE-ND is limited to
variables that represent tree, light, and forest floor substrate (Pacala et al. 1996, Coates et al.
2003), linker functions can relate those outputs to ecosystem services such as wildlife habitat.
SORTIE-ND can be used to model the response of stand structure to disturbances that result
in complex stand structure such as natural disturbances (windthrow patches or speciesspecific pests such as mountain pine beetle), partial harvest silviculture systems and a variety
of silviculture treatments (e.g., thinning or broadleaf tree removal). Applying linker functions
allows researchers and forest managers to address questions regarding habitat dynamics
following disturbance events or management activities.
Linker functions extend forest dynamics models to ecosystem services by representing
relationships between stand structure and ecosystem elements (Elkin et al. 2013, Mönkkönen
et al. 2014, Bugmann et al. 2017). The overstory tree canopy plays an important role in
determining the structure and composition of the understory (Coxson and Marsh 2001,
Sulyma and Coxson 2001, Bartemucci et al. 2006, Haughian and Burton 2015). Thus, habitat
elements within forested stands are intrinsically linked to overstory dynamics. Developing
linker functions that represent these relationships enables modelling not only tree dynamics
but the corresponding dynamics of habitat as it responds to changes in the overstory. Most

16

studies develop linker functions from existing literature and expert knowledge (e.g., Elkin et
al. 2013, Mönkkönen et al. 2014, Bugmann et al. 2017, Cristal et al. 2019). However, Lafond
et al. (2015, 2017) provided instances where linker functions were derived from forest
inventory data.
I used literature and informal interviews with experts to develop linker functions that
represented the foraging habitat of TEC caribou and their apparent competitors, moose. I
then conducted a field study and used empirical data to fit and calibrate the linker functions.
These functions could be applied to a number of modelling platforms, but I developed them
for use with SORTIE-ND as it was the stand dynamics model best suited to modelling the
complex stand dynamics of the region.
The objectives of this chapter were to describe the method I used to develop linker functions
premised on empirical data. When applied to the output from SORTIE-ND, these linker
functions can be used to model the winter habitat of northern mountain caribou over space
and time at the scale of the forest stand. These linker functions can support research or
conservation planning designed to understand and maintain the habitat for northern mountain
caribou, a species of high conservation priority.
2.2 Methods
2.2.1 Study Area and Caribou Habitat Context
The study area is located in west-central BC, Canada. The area is classified as the SubBoreal Spruce moist cold Babine variant (SBSmc2) (Figure 2) and characterized as a
montane zone with a continental climate (Meidinger and Pojar 1991, Banner et al. 1993).
Mature upland coniferous forests are dominated by hybrid spruce (Picea glauce x
engelmannii) and subalpine fir (Abies lasiocarpa (Hook.) Nutt.) with lodgepole pine and

17

trembling aspen (Populus tremuloides Michx.) as common seral species (lodgepole pine is
also common in drier mature stands) (Meidinger and Pojar 1991). This area was selected
because it comprises a portion of TEC winter range that is subject to timber harvesting.
Wildlife managers are interested in developing timber management strategies that support the
provision of winter habitat for caribou.

Figure 2. Distribution of empirical plots (black triangles) established within the Sub-Boreal Spruce
moist cold Babine variant (SBSmc2) (purple) during the summer of 2019 in west-central British
Columbia, Canada.

TEC winter habitat consists predominantly of mature (greater than 80 years of age) pine
forests on poor productivity sites (Cichowski 1989, 1993). These forests support the primary
winter forage of Tweedsmuir-Entiako caribou: terrestrial lichens in the genus Cladonia (P.
18

Browne) (the Cladina morphotype (Cladonia subgenus Cladina) and the Cladonia
morphotype (squamulose species in Cladonia subgenus Cladonia)), Cetraria (Ach).,
Flavocetraria (Kärnefelt & Thell), and Stereocaulon (Hoffm.) (Cichowski 1989, Cichowski
and Haeussler 2013, Denryter et al. 2017). I will refer to this group of terrestrial lichen as
forage lichen or simply lichen for the remainder of the thesis. This is the response variable for
the first linker function.
There have been few studies of the structural attributes of pine stands within the TEC range
that support abundant forage lichens. A study by Lance and Mills (1996) assessed the
structure of stands used by the TEC in the Chelaslie migration corridor (this area is
occasionally also used as winter habitat (Steventon 1996)). These stands were noted as
having similar habitat characteristics to the winter habitat described by Cichowski (Lance and
Mills 1996). The heavily used stands in the Chelaslie corridor had less than 1600 stems/ha of
trees greater than 12.5 cm diameter at breast height (dbh) and stand density was rarely greater
than 2000 stems/ha (Lance and Mills 1996).
High-quality caribou habitat should have abundant forage, but also relatively little risk of
predation due to poor-quality moose habitat (Rettie and Messier 2000, Wittmer et al. 2007,
Venier et al. 2014; see apparent competition discussion in Chapter 1). This is consistent with
the known distribution of TEC as they avoid habitats that coincide with vegetation
communities that are used by moose. That includes more productive site types and mixed
broadleaf stands, cutblocks between 11 and 35 years of age, and spruce-dominated stands
(Cichowski 2015, DeMars and Serrouya 2020).
Peek (2007) stated that forage is the most important component of habitat for moose.
Although there are regional differences, the primary winter browse species for moose in
northern B.C. are trembling aspen, birch (Betula spp.), subalpine fir, willow (Salix spp.), red19

osier dogwood (Cornus stolonifera Michx.), saskatoon (Amelanchier alnifolia (Nutt.) Nutt.
Ex M. Roem.), highbush-cranberry (Viburnum edule (Michx.) Raf.), and mountain ashes
(Sorbus scopulina (Greene) and S. sitchensis (M. Roem.)) (Eastman 1974, Goulet 1985,
Westworth et al. 1989, Rea and Booth 2011). The suite of shrubs suitable as moose browse
will hereafter be referred to as moose shrubs. The majority of these moose shrubs are
associated with mesic to rich sites and younger stands (Eastman 1977, Beaudry et al. 1999).
Canopy cover, primarily from conifer trees, is also recognized as an important component of
moose winter habitat (Mysterud and Østbye 1999, Peek 2007). Canopy cover provides
thermal cover and snow interception that enables moose to conserve energy throughout the
winter months (Eastman 1977). Throughout this paper, the cover provided by a conifer
canopy is referred to as vertical cover. Moose shrub and vertical cover are response variables
for the second and third linker functions developed in this study.
2.2.2 SORTIE-ND
I used SORTIE-ND to simulate within-stand dynamics. SORTIE-ND was initially developed
from a small-scale disturbance model of transitional oak-northern hardwood forests in the
northeastern United States (Pacala et al. 1996). In 1995, research scientists at the British
Columbia Ministry of Forests in northwest BC began collaborating with the developers of
SORTIE and initiated a SORTIE research program to model the dynamics of two local forest
types (Canham et al. 2004, Astrup et al. 2008a). As a result, SORTIE-ND has now been
parameterized for sub-boreal spruce (SBS) and northern temperate interior cedar-hemlock
forests (ICHmc2) (Meidinger and Pojar 1991, Banner et al. 1993, Canham et al. 2004).
SORTIE-ND simulates changes in tree populations at the stand scale. The model uses a
combination of empirical and mechanistic behaviours. Empirical behaviours (sub-models)
were developed from field experiments that measured fine-scale and short-term interactions
20

among individual trees, to project forest dynamics (Canham et al. 2004, Coates and Hall
2005, Maleki et al. 2020). The model structure uses field experiments and testing of alternate
hypotheses to best parameterize the demographic processes and tree growth relationships
found in the model (Kobe and Coates 1997, Wright et al. 1998, 2000, Canham et al. 1999,
2004, LePage et al. 2000, Coates and Hall 2005, Coates et al. 2013).
SORTIE-ND is a spatially explicit, individual-tree model. The model tracks not only
individual tree size, regeneration, growth and mortality information but also the spatial
distribution of the trees and the simulated light environment in a grid across the entire stand
(Canham et al. 2004). A snag sub-model represents the dynamics of dead trees (Bose et al.
2015). Appendix 1 provides an in-depth description of SORTIE-ND following the protocol
presented by Grimm et al. (2006, 2010). I selected SORTIE-ND for this study because it is a
neighborhood dynamics (i.e. stand) model that has been parameterized to represent
interactions between all major tree species in the study area, and allows growth and
regeneration dynamics following mountain pine beetle attack and partial cutting silviculture
systems to be simulated.
2.2.3 Linker Functions
I developed functions that represented how stand structural attributes influenced target habitat
elements (see Mönkkönen et al. (2014) and Cordonnier et al. (2013) for the development of a
comparable framework). I used a four-step process to develop the linker functions. First,
mathematical functions were chosen for forage lichen, moose shrub cover, and vertical
conifer cover based on a literature review and discussions with regional experts (see
Appendix 2). The predictor variables included: understory light, time since disturbance,
edaphic site condition, competing low vegetation, humus depth, conifer basal area and
broadleaf basal area. Second, these linker functions informed a field sampling effort where

21

corresponding stand structure, edaphic site condition, and vegetation community data were
collected. Third, the empirical data were used to fit the linker functions. Fourth, the functions
were calibrated and the predictor variables were limited to only those that could be output
from the stand dynamics model SORTIE-ND.
The linker functions were designed to output a “suitability index” that represented the
suitability of a stand to provide vertical cover or support the occurrence and growth of forage
lichen or moose shrubs. This is not equivalent to the conventional habitat suitability index
(HSI) developed by the United States Fish and Wildlife Service (1981) for the quantification
of wildlife habitat (Roloff and Kernohan 1999, Pecchi et al. 2019). My linker functions did
not directly model species abundance, cover or biomass nor did they evaluate habitat at the
scale represented by most HSI models (Roloff and Kernohan 1999, Pecchi et al. 2019). My
linker functions were also not full indices of habitat availability because the suite of predictor
variables was limited to those that can be derived from SORTIE-ND output. Also, I was not
able to test the relevance of all potential stand-scale variables (e.g., stand initiating
disturbance intensity). The linker functions provided a proxy measure of whether a certain
stand had the potential to support the occurrence and growth of lichen or moose shrubs, or
provide vertical cover. As a short form, throughout this thesis I refer to this potential for a
stand to support or provide a habitat element as simply “suitability” (e.g., lichen suitability or
suitability for lichen). Appendix 2 describes the linker functions following the protocol
presented by Grimm et al. (2006, 2010).
2.2.4 Model Description
2.2.4.1. Lichen
The suitability of a forest stand to support terrestrial forage lichen was represented as the
product of time (Equation 2) since the area was disturbed, light availability in the understory
(Equation 3), abundance of broadleaf trees (Equation 4), and the stand’s edaphic conditions
22

(Equation 5). Equation 1 illustrates how these four variables were combined multiplicatively
to define the full linker function for lichen suitability. Although competing low vegetation
and humus were identified as potential predictor variables, these are not presented here
because they were not selected for the final set of linker functions and are not simulated by
SORTIE-ND.
ℎ

=

∗

ℎ ∗

∗

ℎ

(1)

Forage lichen requires time to establish following a stand replacing disturbance that removes
or severely disturbs the forest floor. Other lichen may establish more quickly following
disturbance but this linker function represents only those lichen species selected by caribou as
forage lichen (Maikawa and Kershaw 1976). Previous studies have found that the cover and
biomass of lichen is much less following disturbance but increases as lichens re-colonize the
site (Maikawa and Kershaw 1976, Johnson 1981, Rowe 1984, Coxson and Marsh 2001).
Although there are regional and ecosystem differences, studies have shown that, in general,
approximately 30–50 years are required before species of forage lichen establish in quantities
similar to pre-disturbance levels (Maikawa and Kershaw 1976, Johnson 1981, Carroll and
Bliss 1982, Rowe 1984, Coxson and Marsh 2001). I represented this post-disturbance
recovery period using a Weibull function that assumed that suitability for forage lichens was
initially low and increased to an asymptote of 1 by approximately year 50 (Equation 2, Figure
3a).
∗
=1−
where tdisturbance represented time since stand replacing disturbance

(2)

Lichen are symbiotic organisms that consist of fungi and associated photobionts
(cyanobacteria or algae) (Brodo et al. 2001, Spribille 2018). Although lichen may be present
in stands with low light levels, a minimum amount of light is required to support

23

photosynthesis. Numerous studies have found that abundance of terrestrial lichen is
positively associated with higher light levels typically found under open conifer canopies
(Carroll and Bliss 1982, Sulyma and Coxson 2001, Haughian and Burton 2015, Lewis et al.
2019).
I used a logistic function to represent the positive relationship between light (GLI) and site
suitability for lichen (Equation 3, Figure 3b). The inflexion point in the function represented
the significant increase in lichen suitability, evident in the empirical data, at light levels
between 25 to 40% (Figure 3b). Lichen suitability approached full suitability above 50% of
full light.

ℎ =

+

(3)
1+

( – ∗

)

Terrestrial lichen is generally limited in abundance in broadleaf stands (Longton 1992, Crites
and Dale 1998, Russell and Johnson 2019). The abundant and repeated litterfall of deciduous
broadleaf trees and vascular plants limits the abundance of slow growing lichen (Canters et
al. 1991, Cornelissen et al. 2001). However, terrestrial lichen can be found in low abundance
in some stands with a broadleaf component (Figure 3c) (Johnson 1981). The negative
logistic function provided flexibility in the rate at which lichen suitability declined in
response to increasing broadleaf abundance.

= 1−

(4)

1
1 + e(

–

∗

)

where BAdecid represented the average basal area per hectare of all
broadleaf stems
Terrestrial forage lichens grow slowly and are poor competitors in comparison to other
ground vegetation and moss (Johnson 1981, Rowe 1984). These lichen survive and persist in
conditions unfavourable for other vegetation (Rowe 1984, Haughian and Burton 2015).
Primary determinants of vegetation abundance are soil nutrient level and water holding
24

capacity (Maikawa and Kershaw 1976, Meidinger and Pojar 1991, Banner et al. 1993,
Coxson and Marsh 2001, Sulyma and Coxson 2001). In British Columbia, relative soil
moisture regime (SMR) and relative soil nutrient regime (SNR) are indices for soil water
holding capacity and nutrient availability (respectively) that are commonly used to describe
edaphic site characteristics (Meidinger and Pojar 1991). These indices correlate well with
field measurements of moisture and nutrient availability (Kranabetter et al. 2007, Kranabetter
and Simard 2008). SMR and SNR are highly correlated in the study region and were
therefore combined into one edaphic variable (Equation 5).
I used a negative logistic equation to represent the relationship between lichen suitability and
edaphic condition (Equation 5, Figure 3d). Suitability is greatest on low productivity edaphic
sites (low moisture or nutrients) and declines as site productivity increases. The inflexion
point provided flexibility to reflect how steeply suitability declines as edaphic site
productivity increases.

ℎ

= 1−

a4
1 + e(

–

(5)
∗

)

where smnr =

2.2.4.2. Moose
Browse for moose during winter consists of both shrubs and low trees (Eastman 1974, Goulet
1985, Westworth et al. 1989, Rea and Booth 2011). Browse availability from trees can be
derived directly from SORTIE-ND (individual tree data). However, linker functions were still
required to project the relative suitability of the tree component within the stand (on a scale of
0 to 1). Linker functions were also developed to project the suitability of the stand to support
understory shrubs browsed by moose. The moose shrub (Equation 7) linker function was
combined with the aspen (Equation 11) and subalpine fir (Equation 12) linker functions to
project the overall suitability of the stand to provide winter browse for moose (Equation 6).
25

ℎ

=

i.

+

(6)

+
3

Moose Shrubs

The suitability of a forest stand to support shrub species used by moose for winter browse
(Equation 7) was represented as a function of time since disturbance (Equation 8), conifer
basal area (Equation 9) and edaphic site condition (Equation 10) (Figure 5).
ℎ

=

∗

∗

(7)

ℎ

As with lichen, shrubs require time to recover or establish following stand replacing
disturbances. Abundance of moose forage has been shown to peak approximately 25 years
post-disturbance (Julianus et al. 2019); a lognormal function enabled a peak in suitability
(Equation 8, Figure 5a) (Canham et al. 2004).
(8)
=

. ∗

where tdisturbance represented time since disturbance
Conifer trees outgrow moose shrubs in height and size and over time limit the light, water
and nutrients available for shrubs in the understory. I used a negative logistic function to
represent the competitive effect of conifers on moose shrubs, with conifer basal area
representing conifer abundance (Equation 9, Figure 5b). The shape of this function reflected
the decline in suitability for moose shrubs as conifer basal area increases.

=1−

(9)
1+

(

–

∗

)

where conBA represents total basal area of conifers greater than 12.5 cm dbh
Understory plants require moisture and nutrients to grow. Following bark beetle outbreaks,
competition from overstory trees is rapidly reduced resulting in an increase in available soil
26

moisture and nutrients. A corresponding increase in understory vegetation abundance has
been documented (Pec et al. 2015, Cichowski and Haeussler 2020). Therefore, I
hypothesized that in general, as site productivity increases, the biomass of understory
vegetation will increase. I used a simple linear function to represent the positive relationship
between edaphic condition and suitability to support moose shrubs (Equation 10, Figure 5c).

= 5

ℎ

ii.

(10)

+ 6
2

Aspen

Young aspen trees often comprise a component of moose winter diet (Eastman 1977, Rea
2014). Aspen is available as browse until fresh shoots are out of reach of moose (Cumming
1987, Timmermann and McNicol 1988). Vegetation between 50 cm and 3 m in height is
considered accessible to moose in the winter (Peek et al. 1976, Cumming 1987, Timmermann
and McNicol 1988). I used a linear function to reflect the increase in suitability as stem
density of young aspen increases (Equation 11, Figure 5d). I set a threshold, above which,
additional stems did not provide a further increase in browse suitability.

=

(

, 1)

(11)

where at represented the density of aspen 0.5 - 3 m in height and
maxat represented the density threshold for aspen
iii.

Subalpine fir

Advanced regeneration of subalpine fir constitutes a significant component of moose diet in
some areas and some winters (Eastman 1977, Rea 2014). For example, Hodder et al. (2013)
found that approximately 45% of moose diets consisted of subalpine fir during a deep-snow
winter. I used a linear function to reflect the increase in browse suitability with increased
availability of subalpine fir stems (Equation 12, Figure 5e). Similar to aspen, I applied a
27

maximum density threshold and only trees 0.5 - 3 m in height contributed to browse
suitability.

=

(

, 1)

(12)

where bl represented the density of subalpine fir 0.5 - 3 m in
height and maxbl represented the density threshold for
subalpine fir
iv.

Vertical cover

Canopy cover has been shown to provide significant value for moose winter habitat through
snow interception and thermal protection (Peek et al. 1976, Eastman 1977). I selected a
positive linear function to represent the increase in suitability with increased conifer
presence. Conifers greater than 10 m in height are considered to provide suitable vertical
cover (Allen et al. 1987, Timmermann and McNicol 1988). Because the function was linear,
a maximum threshold was required. In discussions with local experts, fifty years of age was
regarded as the approximate age at which stands in the region provided sufficient cover for
moose (Crowley, Hodder, Johnson, pers. comm. 2020). This was supported by results from
Norway that found moose increased selection of stands greater than 40 years of age in winter
(Bjørneraas et al. 2011). Local growth and yield curves (Goudie 1984) were used to
determine that, at age 50, hybrid spruce and lodgepole pine stands achieved a basal area of 30
m2/ha. Therefore the maximum threshold for vertical cover was set at a conifer basal area of
30 m2/ha.

=

(

, 1)

(13)

where con represented basal area of conifers greater than 10 m in
height and maxcon represented the threshold for conifer
basal area
28

2.2.5 Field Data Collection
I established plots and collected data across the SBSmc2 subzone (Figure 2). Due to
logistical constraints, I did not restrict sampling to solely within the winter habitat of the
TEC. This improved the dataset as I was able to sample mature stands with live pine trees
that were located in the western portion of the SBSmc2. There was lower incidence of MPB
attack in this area. In addition, the larger sample area contained unlogged, mature stands
across a broader range of edaphic site types.
I used a purposeful sampling approach to acquire a balanced dataset across three gradients of
site and stand conditions. Purposeful sampling was employed because it provides an efficient
means to acquire a balanced data set that includes rare stand types and ecosystems (Lilles and
Astrup 2012). This site selection method was employed by other studies in the region that
used the same analysis methods (Lilles and Astrup 2012, Coates et al. 2013, 2018).
Gradients of stand age, trees species composition, and edaphic condition were selected to
guide sampling with the intent to acquire a representation of the full range of stand and site
variables.
I developed maps that stratified the study area according to: (1) stand age (four age classes:
0–40, 41–80, 81–140, > 140 years of age), (2) tree species composition (aspen leading stands,
pine leading stands, spruce leading stands, subalpine fir leading stands), and, where available,
(3) edaphic condition (from poor and subxeric to rich and subhygric). I used the Province of
BC’s Vegetation Resource Inventory mapping (British Columbia Data Catalogue n.d.) and
available ecosystem mapping to develop these maps and identify potential sample sites. To
develop a sample plan, I identified stands that were mapped as containing rare combinations
of stand conditions missing from my sample set (e.g., dry sites, aspen stands, mature spruce
stands on mesic or productive sites). Other suitable stands, encountered in the vicinity, would

29

be opportunistically sampled. I limited sampling to stands that were within approximately
2000 m of a road or a reasonable walking distance.
I sampled stands if they were uniform in stand structure, uniform in edaphic condition (soil
moisture and nutrient regime), and at least 50 m from a disturbance edge. To be considered
uniform in stand structure, the stand needed to have a similar species composition and similar
canopy structure throughout (e.g., either uniformly dense or uniformly open). This was to
control for a relatively uniform light level within the plot. Once a stand was determined to be
suitable, plot center was established by moving 10 paces along a randomly selected azimuth.
Each sample site consisted of two nested plots. Large plots were 7.98 m in radius (200 m2)
and were used to describe the site, measure large trees (all live and dead trees > 12.5 cm dbh),
perform a vegetation cover assessment, and survey for signs of moose browsing and pellet
piles. I also assessed the percent cover of ground cover in the large plots (low growing
vegetation, terrestrial cryptogams, litter, wood, and rock).Within a smaller 3.99-m radius plot
(50 m2), I measured trees less than 12.5 cm dbh.
Site description
I followed the Province of British Columbia’s “Field Manual for Describing Terrestrial
Ecosystems” (hereafter referred to as the Field Manual) to record sample site information and
describe the ecological conditions at each sampled site (2010). Site information included:
GPS location, elevation, slope, aspect, mesoslope position, surface shape, and stand structural
stage. For naturally regenerated stands, I determined stand age by taking a tree core from one
of the dominant canopy trees (at 1.3 m height, corrected for the time to grow to 1.3 m using
height-age site index curves (Thrower et al. 1994)). For plantations, I used the age of the
planted trees recorded in the Province of BC’s silviculture database. In mountain pine beetle
attacked stands, I estimated the time since mountain pine beetle attack by noting the years

30

since release when counting the age of a tree. I cross-referenced this estimate by counting the
number of years of release exhibited in the growth whorls of understory advance
regeneration. I dug a soil pit in each plot to a minimum depth of 30 cm or until a restricting
soil layer was reached. I recorded the following soil information as per the Field Manual:
depth of humus layer (including litter, fermented, and humic organic horizons (LFH)), humus
form identification (mor, mormoder, moder or mull), seepage depth, rooting depth, root
restricting layer, A and B horizon identification, A horizon depth, dominant soil texture (of
the rooting zone), and percentage of coarse fragments (in the rooting zone) (Province of
British Columbia 2010). I used the soil and site description information to determine the soil
moisture regime and soil nutrient regime of each site. Vegetation was used to cross-reference
the site determination and site series (Banner et al. 1993).
I took hemispherical photos at five locations within the plot: one at plot center and one in
each ordinal direction, 3.99 m from plot center. The photos were taken 1.3 m above the
ground with a Nikon Coolpix 5000 digital camera and Nikon FC-E8 fisheye lens converter. I
analysed the images using Gap Light Analyzer software (GLA 2.0) to obtain the gap light
index (GLI), a measure of the percentage of full sunlight (sum of direct and diffuse light) that
a point receives throughout the growing season (Frazer et al. 1999). GLA models the seasonal
and daily solar radiation transmitted through the canopy (Canham 1988, Wright et al. 1998). I
used the mean GLI of the five photos to provide a metric of light available in the understory
of the stand.
Ground Cover
I assessed lichen cover and other classes of ground cover in each of the four quadrants of the
7.98-m plot. I used an ocular method to estimate the percent cover of 12 classes of ground
cover: forage lichen, non-forage lichen (all other species of lichen), feathermoss (Pleurozium
schreberi, Ptilium cristae-cristensis, Hylocomium splendens), non-feathermoss moss, dwarf
31

shrubs, grass, herbs, conifer litter (woody debris < 5 cm in diameter), deciduous litter, wood
(> 5 cm in diameter), rock, and bare ground (Province of British Columbia 2010). Lichen
cover was recorded by species and was subsequently averaged across all quadrants to provide
the estimate of percent cover by species for the plot.
I used percent cover of forage lichen as a surrogate for forage lichen biomass. Bergerud
(1971) and Thomas et al. (1996) demonstrated a strong correlation between lichen cover and
biomass and recommend the use of cover where efficient measurement is required. Johnson
et al. (2000) noted that caribou browsing would reduce lichen biomass; browsing pressure
could be inconsistent and percent cover was considered to be less sensitive to bias. Also,
lichen cover is measured much more quickly than volume. That efficiency was important
considering that I collected an extensive amount of data within in each plot. Lichen cover was
normalized to produce a relative index of cover that ranged from 0 to 1.
Vegetation and Tree Mensuration
I recorded vegetation cover for each species located in the 7.98-m radius plot. I assessed
vegetation by layer: tree (A), shrub (B), herb (C), and moss/lichen/seedling layer (D)
(Province of British Columbia 2010). I assessed cover of tree species by sub-layer with cover
of each species (any woody species > 10 m in height) differentiated into cover of the sublayer: dominant trees (A1), main canopy trees (A2), sub-canopy trees (A3). I differentiated
the cover of shrub species (woody plants including trees and shrubs) into a tall shrub layer
(2–10 m in height) and a low shrub layer (˂ 2 m in height except low-growing woody or
trailing plants that are usually < 15 cm in height) (Province of British Columbia 2010). I
evaluated the cover of each species of terrestrial moss, lichen, liverwort, and tree seedling < 2
years old (D layer). I followed the methods presented in the Field Manual for a full
vegetation plot with the only differences being the plot size and shape. Unknown species
were collected for later identification.
32

I measured the diameter (dbh) of all standing trees ≥ 12.5 cm dbh in the 7.98-m radius plot.
In addition to dbh, I recorded tree species, live/dead status, canopy position (dominant, codominant, intermediate or suppressed) and cause of death (if possible). I tallied the number
of fallen trees killed by mountain pine beetle and estimated the dbh of each fallen tree. I
measured the height of a tree that represented the average canopy height of the stand.
I assessed smaller trees in the 3.99-m radius plot. I measured the diameters of trees with 7.5–
12.4 cm dbh. I tallied trees ˂ 7.5 cm dbh by the following size classes: 0–10 cm in height,
10–30 cm in height, 30–130 cm in height, 0–2.4 cm dbh, 2.5–4.9 cm dbh, and 5–7.4 cm dbh.
I recorded tree species and live/dead status for trees of all size classes.
Browse on Moose Shrubs
I recorded the species and percent cover of each shrub and shrub-sized tree within the 7.98-m
plot. For the purposes of quantifying available moose browse, I collected further data on
individuals that had branches available for browse (< 3 m height). For these individuals, I
recorded percent cover and modal height (measured to the nearest 10 cm). I also assessed the
trees and shrubs for signs of moose browse. Moose and deer browse were difficult to
distinguish. The areas sampled were not considered high-quality winter habitat for deer
therefore browsing was attributed to moose. Signs of moose browse was differentiated from
rabbit browse because moose browsing results in a torn, ragged edge; rabbit browsing leaves
a sharp edge at a 45° angle.
I adapted the “Procedures for Environmental Monitoring in Range and Wildlife Habitat
Management Version 5.0” method to assess the degree of browse (Habitat Monitoring
Committee 1996). For each species, I estimated the percent of each plant browsed instead of
assigning an ordinal browse class. I also assessed the hedging on any species that was
browsed. I used the utilization form classes outlined in the B.C. Ministry of Environment

33

manual “Describing Ecosystems in the Field” (Luttmerding et al. 1990). I condensed the sixcode system into two simple codes: code 1 indicated little or no hedging and code 2 indicated
moderate to severe hedging. The number of piles of moose pellets in each plot was also
recorded.
I used the moose browse data in combination with regional literature sources (Eastman 1974,
Goulet 1985, Westworth et al. 1989, Baker 1990, Rea and Booth 2011) to determine which
plant species should be considered moose forage in the study area. The majority of browse in
our field plots occurred on species that were previously documented as winter browse:
trembling aspen, subalpine fir, willow, red-osier dogwood, saskatoon, highbush-cranberry,
and mountain ash. However, additional species were also browsed in some field plots. I
compared the amount that each species was browsed to the abundance of the species in the
plots. Red elderberry (Sambucus racemosa L.) and Douglas maple (Acer glabrum Torr.)
occurred in low abundance but moderate to heavy browse was evident when the shrubs were
present; therefore, I included these species as moose shrubs. Black cottonwood (Populus
balsamifera ssp. Trichocarpa (T. & G.) was also browsed when present but was not included
because it is a tree. Twinberry (Lonicera involucrata (Richardson) Banks ex Spreng.), prickly
rose (Rosa acicularis Lindl.), and thimbleberry (Rubus parviflorus Nutt.) were only
occasionally browsed when the species were present. Hedging was not evident, therefore, I
did not include these species as moose shrubs.
I derived the moose shrub variable from the total cover of all moose shrub species in a plot.
Total cover was normalized across all plots to provide an index of cover between 0 and 1 for
each plot.

34

2.2.6 Linker Function Fitting and Calibration
I used the empirical data (76 sample points) to fit and calibrate the lichen and moose shrub
linker functions. During calibration, predictor variables were limited to only those available
as SORTIE-ND output.
First, I used likelihood methods presented in Canham and Uriarte (2006) to parameterize and
select the most parsimonious models. The likelihood method follows a 4-step process: “(1)
model specification, (2) parameter estimation using maximum likelihood methods, (3) model
comparison, and (4) model evaluation” (Canham and Uriarte 2006 p. 64). I compared a
nested suite of models to determine the degree of support for including each of the
hypothesized predictor variables (Table 1 and Table 2; see Lilles and Astrup (2012) for a
similar application of the likelihood method)). I used the Cran R ‘likelihood’ package (v1.7,
(Murphy 2015)) in R (version 4.0; R Core Team 2021) for parameterization and model
comparison. This package uses simulated annealing to produce maximum likelihood
estimates for each parameter in each model and calculates support intervals for parameter
estimates (Edwards 1992, Goffe et al. 1994). Parameters within each model were estimated
simultaneously (Canham and Uriarte 2006, Murphy 2015). The likelihood estimation used in
this package requires a probability model; here I used a gaussian distribution for the
probability density function. In this package R2 is calculated as the sum of squared expected
divided by the sum of squared total (Murphy 2015). Akaike information criterion corrected
for small sample size (AICc) was used to assess model fit and evaluate the suitability of
functional forms (Burnham and Anderson 2002).
After parameterizing the functions and selecting the model with the most support from the
data, I evaluated the best fit parameters and model forms. For some functions, calibration was
required. Due to the small sample size, there was the potential that the likelihood method

35

Table 1. Structure and evaluation of linker functions developed to project the relative value of a
forest stand to support forage lichens. Linker functions were developed using data collected in subboreal forests of west-central British Columbia during the summer of 2019. Goodness of fit and
parsimony of each model was assessed using ΔAICc and R2.

Model
Null
Light
Time
Broadleaf
SMR2
SNR3
SMNR4
Full Model
(C.Veg5)
Full Model
(SMNR)
Full Model +
Humus
1

Variables

light
time
BAdecid1
smr
snr
smnr
time * light * decid *
veg.effect
time * light * decid *
smnr
time * light * decid *
smnr * humus

Number of
parameters
1
4
3
3
4
4
4
8

AICc

ΔAICc

R2

13929.3
605.7
616.9
613.5
578.2
574.8
574
582.7

13380.5
56.9
68.1
64.7
29.4
26.0
25.2
33.9

0
0.17
0
0.05
0.42
0.44
0.45
0.46

8

562.4

13.6

0.58

9

548.8

0

0.66

BAdecid = basal area of broadleaf trees (m2/ha) 2 SMR = relative soil moisture regime; 3 SNR =

relative soil nutrient regime; 4 SMNR = (SMR + SNR) / 2; 5 C.Veg = competing low vegetation.

overfit the data. Some best fit statistical models produced inappropriate values that did not
enable suitability on sites where the empirical data showed habitat elements did occur. For
these situations I calibrated the linker function to ensure a minimum level of suitability was
possible (Table 3, Table 4, Appendix 3).
Metrics of the tree species that are browsed by moose (i.e., aspen and subalpine fir stem
density) and conifer cover (BA) are directly output by SORTIE-ND models, and therefore did
not require calibration.
2.3 Results
2.3.1 Lichen
In my field plots, the maximum cover of terrestrial forage lichen was 61% (mean cover was
5.3%, median cover was 0.04%) (Figure 15 in Appendix 4). I used 61% to define the

36

Table 2. Structure and evaluation of linker functions developed to project the relative value of a
forest stand to support forage shrubs used by moose. Linker functions were developed using data
collected in sub-boreal forests of west-central British Columbia during the summer of 2019. Goodness
of fit and parsimony of each model was assessed using ΔAICc and R2.

Model
Null
SMR1
SNR2
Light
Time
Conifer BA3
SMR, SNR, Light
SMNR4, Light
SMNR, Light, Time
SMNR, Light, Time,
Conifer BA
SMNR, Light, Time,
Conifer BA 12.55
SMNR, Light, Time,
Conifer BA 17.56
SMNR, Time,
Conifer BA 12.5

Number of
parameters
1
2
2
3
4
4
3
3
5
7

AICc

ΔAICc

R2

479.9
324.8
332
334.7
325.5
328.3
324.8
325.6
321
320.4

162
6.9
14.1
16.8
7.6
10.4
6.9
7.7
3.1
2.5

0
0.11
0.02
0.01
0.15
0.12
0.13
0.12
0.22
0.27

smnr * time * conBA12.5 * light

7

320.3

2.4

0.28

smnr * time * conBA17.5 * light

7

318.9

1

0.29

smnr * time * conBA12.5

6

317.9

0

0.28

Variables

smr
snr
light
time
conBA
SNR * SMR * Light
SMNR * light
smnr * time * light
smnr * time * conBA * light

1

SMR = relative soil moisture regime; 2 SNR = relative soil nutrient regime; 3 BA = basal area
(m2/ha) of all conifer trees; 4 SMNR = (SMR + SNR) / 2; 5 BA 12.5 = basal area of all conifer trees
with diameter at breast height (dbh) ≥ 12.5 cm; 6 BA 17.5 = basal area of all conifer trees with dbh ≥
17.5 cm

maximum Lichen Suitability Index of 1.0 (Figure 3). For lichen, the non-linear model with
time, light, broadleaf basal area, edaphic site and humus was the model that best fit the data
(Full Model + Humus, Table 1). However, humus data can not be derived from SORTIE-ND
output so I selected the second-best model for the lichen suitability linker function (Full
Model (SMNR)).
Lichen suitability approached maximum suitability within 50 years following a standreplacing disturbance (Figure 3a). Lichen suitability, and observed lichen cover, increased
with increasing light availability (Figure 3b). For the light component of the lichen linker
function, I set the minimum suitability to 0.2 (Table 3, d = 0.2) because forage lichen
37

Figure 3. Linker function components that represent how the suitability of a stand to support forage
lichen varies as a function of stand state: (a) years since stand replacing disturbance, (b) light level
represented by Gap Light Index (GLI), (c) basal area of broadleaf trees in the stand and (d) soil
nutrients and soil moisture (normalized gradients: SNR/5 and SMR/6 respectively). Linker functions
are represented by the blue lines in (a) - (c) and by shading in (d). Empirical plot data (n = 76)
collected during the summer of 2019 in sub-boreal forests of west-central British Columbia are
represented by black dots.

occasionally occurred when understory light was very low (Figure 3b). The suitability of a
stand to support lichen decreased as broadleaf tree basal area increased (Figure 3c). I
calibrated the logistic function to enable stands with a small broadleaf component to support
lichen suitability since lichen occurred in plots with up to a 10% broadleaf component. I set
the inflection point at a low broadleaf basal area of 5 m2/ha (Figure 3c). Lichen suitability,
and observed lichen cover, was greatest in dry stands with low soil nutrients (Figure 3d). The
edaphic component of the lichen linker function was calibrated so that the poorest
productivity edaphic sites in the empirical data (that had the greatest amount of lichen) had
38

Table 3. Parameter values for linker functions developed to project the relative value of a forest stand
to support forage lichen in west-central BC: maximum likelihood estimates (with support intervals)
and calibrated parameter values (where calibration occurred) are shown for each of the parameters.

Lichen
Maximum
Linker
Likelihood Estimate
Function
Parameter
z
0.09
a
b
c
d
b1
b2
a4
b4
c4

1 (not estimated)
4.98
0.15
0 (not estimated)
n/a1
n/a1
0 (not estimated)
6.15
17.70

Support Interval

0.04 – 5.00 (upper bound set
to 5.0)
n/a
4.33 – 5.68
0.14 – 0.18
n/a

5.72 – 6.53
16.63 – 18.87

Calibrated
Parameter
Value

0.8
0.2
2.5
0.5
0.9
4
8

1

b1 and b2 did not have a maximum likelihood estimate because a Weibull model was initially
included in the model set. It was not included in the final linker function because the Weibull form did
not allow for lichen occurrence at low levels of overstory broadleaf abundance.

suitability approaching 1.0 (Figure 3d) (see also Haughian and Burton 2015). The minimum
suitability of 0.1 (a = 0.9) was set to reflect the occasional sparse presence of lichen on some
higher moisture sites (Figure 3d, (Ahti and Oksanen 1990)).
Compared to the empirical data, the lichen linker function predicted lichen suitability well for
plots with > 0.2 suitability (Figure 4). The linker function over-predicted lichen suitability for
many plots with low observed suitability (< 0.2). The root mean squared error was 0.2045.
2.3.2 Moose Shrubs
The maximum horizontal cover of moose browse in the empirical data was 10%; this was
used to define the maximum Moose Shrub Suitability Index of 1 (mean cover was 1.5%,
median cover was 0.8%) (Figure 15 in Appendix 4). I selected the model with the best fit to
the data as the linker function for suitability for moose shrubs: time, conifer basal area and
edaphic site (SMNR, Time, Conifer BA 12.5, Table 2).

39

Figure 4. Observed lichen (a) and moose shrub (b) suitability estimated from empirical plots (n = 76)
vs. predicted suitability derived from developed linker functions. Empirical data collected in subboreal stands in west-central British Columbia during the summer of 2019.

Abundance of moose shrubs was greatest in stands 20–30 years post-disturbance (Figure 5a),
and subsequently decreased as stands aged. The potential for elevated moose shrub
suitability was maintained until the basal area of large conifers (> 12.5 cm dbh) reached ~20
m2/ha (Figure 5b); above this basal area, moose shrub cover declined. I set a lower threshold
for moose suitability at 0.25 for stands with conifer basal areas greater than 40 m2/ha (Table
4). This reflected the presence of moose shrubs in some stands with greater conifer basal
area (Figure 5b). In contrast to the response of lichen cover to edaphic conditions, the
occurrence of moose shrubs was greatest in stands with high soil moisture and nutrient levels
(Figure 5c).
Shrubs eaten by moose grow across a broad range of edaphic sites (Figure 6). For example,
willow (Salix spp.) was found across all edaphic sites whereas saskatoon was generally found
on poorer productivity sites and red-osier dogwood on higher productivity sites.

40

Figure 5. Linker function components that represent how the suitability of a stand to support moose
shrubs for forage varies as a function of stand state: (a) years since stand replacing disturbance, (b)
basal area of all conifers greater than 12.5 cm dbh, c) soil nutrients and soil moisture (normalized
gradients: SNR/5 and SMR/6 respectively), and (d) and (e) density of aspen and subalpine fir trees
(respectively) 0.5 - 3 m in height. Linker functions are represented by shading in (c) and blue lines in
the other panes. Empirical plot data (n = 76) collected during the summer of 2019 in sub-boreal
forests of west-central British Columbia are represented by the black dots.
Table 4. Parameter values for linker functions developed to project the relative value of a forest stand
to support forage shrubs used by moose in west-central BC: maximum likelihood estimates (with
support intervals) and calibrated parameter values (where calibration occurred) are shown for each of
the parameters.

Moose Shrub
Linker Function
Parameter
Xo
Xb
aaa1
bbb
ccc
maxat
maxbl
maxcon

Maximum
Likelihood
Estimate
19.58
1.64
1 (not estimated)
7.28
0.22

Support Interval
10.77 – 40.88
1.15 – 3.03
n/a
5.61 – 9.95
0.16 – 0.29

Calibrated
Parameter

0.75
5000
3000
30

1

aaa did not have a maximum likelihood estimate because it was set to 1 during model fitting and
therefore was not estimated (aaa = 1).

41

Figure 6. Distribution of species of moose shrubs across soil moisture and nutrient gradients. Data
points are scaled in size to represent the percent cover of each observation in the empirical dataset (n
= 76). Data points are displaced slightly to display individual points that would otherwise entirely
overlap at the same moisture/nutrient junction. Species codes: ACERGLA – Douglas maple,
AMELALN – saskatoon, CORNSTO – red-osier dogwood, SORBSCO – western mountain-ash,
SORBSIT – Sitka mountain-ash, VIBUEDU – highbush-cranberry, SAMBRAC – red elderberry,
SALIX – willow. Data collected during the summer of 2019 in sub-boreal forests of west-central
British Columbia.

Density of aspen and subalpine fir between 50 cm and 3 m in height are direct outputs of the
stand simulation model. The suitability of a stand to support moose browse increased with
increasing abundance of these trees until an asymptote was reached at 5000 and 3000
stems/ha (Figure 5d and Figure 5e, respectively). The distribution of aspen and subalpine fir
saplings varied significantly within the empirical data (Appendix 4).
Compared to the empirical data, the moose shrub linker function predicted shrub suitability
better for plots with > 0.4 suitability than for plots with < 0.4 suitability (Figure 4). The linker

42

function over-predicted moose shrub suitability for many plots. The root mean squared error
was 0.2887.
2.3.3 Moose Cover
Conifer trees provide vertical cover, contributing to the suitability of a stand as moose
habitat. Vertical cover suitability increased as conifer basal area increased until an asymptote
was reached at a conifer basal area of 30 m2 (conifers > 12.5 cm dbh, Figure 7). Conifer
basal area varied within the empirical dataset between 0 and 58 m2/ha (Appendix 4).

Figure 7. The linker function that represents the effect of conifer basal area (trees taller than 10 m) on
the suitability of a stand to provide vertical cover in winter.

2.4 Discussion
Past studies have developed linker functions primarily through extensive literature review,
related studies, and application of expert knowledge (Cordonnier et al. 2014, Mönkkönen et
al. 2014, Blattert et al. 2017, Cristal et al. 2019). I built on that approach by also
incorporating empirical data I collected during a focused field study. I used the empirical data
43

to improve the predictive ability of these linker functions for the target study area and local
caribou population. The empirical data enabled me to select predictor variables that best
explained patterns in the data and to evaluate the fit of the linker functions.
I used simulated annealing to test and simultaneously parameterize the linker functions. This
resulted in linker functions that isolated the effect of each individual predictor variable on
habitat suitability. For example, the linker function for time reflects solely the time required
for lichen to re-establish on a site following stand-replacing disturbance. The function is not
initially concave as is often observed in stands where forage lichen becomes established
following fire, declines as conifer trees establish, and then increases again as stands mature
(Thomas et al. 1996). However, when the conifer basal area linker function is combined with
the time linker function, the expected concave pattern becomes evident for stands that follow
this successional pattern. The separation of these two linker functions provided the flexibility
to represent the pattern observed in other stand types that transition from lichen-dominated to
moss-dominated understories due to increasing overstory density and biomass (Sulyma and
Coxson 2001, Horstkotte and Moen 2019). Fitting these functions simultaneously enabled the
linker functions to model this interaction of time and overstory structure.
At low suitability, predicted values exceed observed values because predicted values only
reflect the influence of stand structure and edaphic site on lichen and moose shrub
development (Figure 4). Other ecological variables also influence stand suitability and, if
incorporated, would dampen the suitability of sites to support the elements. These variables
include past disturbance type and severity, competing vegetation and bryophytes, and edaphic
differences between moose shrub species. In addition, landscape and microsite scale factors
are also known to influence the provision of caribou habitat (Sulyma and Coxson 2001,
Wittmer et al. 2007, Haughian and Burton 2015). Although the suite of linker functions may

44

not represent a complete ecological model, the linker functions allow for comparisons of the
relative effect of changes in stand structure on forage lichen and moose shrub suitability.
Disturbance type and severity are not reflected in the lichen and moose shrub linker functions
though they play a significant role in the regeneration success of these communities (Viereck
and Schandelmeier 1980, Cichowski et al. 2022). Disturbances such as wildfire, timber
harvest and insect attack differ in terms of the degree of impact that they have on vegetation
communities and the potential for persistence or regeneration for individual species.
Disturbances impact vegetation either directly through physical damage or indirectly through
the provision of substrate suitable for regeneration or abundance of competing vegetation
(most significant for lichen) (Viereck and Schandelmeier 1980, Haughian and Burton 2015,
Lewis et al. 2019).
Linker functions did not directly reflect disturbance type and severity because historical
disturbance data were not available. Even if available, expanding sampling to extend across
disturbance type and severity gradients would not have been feasible. Balanced sampling of
each disturbance type across gradients of intensity, would have required more resources and
time than available for this study. In addition, modelling the effects of disturbance type and
severity are not fully available in SORTIE-ND.
Humus depth is affected by past disturbance severity and type: low intensity surface fires,
high intensity crown fire, and season of burn will influence the degree to which the fire burns
the humus layer (Viereck and Schandelmeier 1980, MacCracken and Viereck 1990). Lichen
is associated with microhabitats that have a shallow humus layers (Haughian and Burton
2015). The lichen model that included the depth of the forest floor (Table 1, model
“Fullmodel + Humus”) had substantially more support from the data than the comparable
model without humus thickness (Δ AICc = 13.6). Humus depth could be considered a
45

surrogate metric that reflects intensity of historic stand-level disturbance. However, many
other factors also influence humus depth including dominant tree species, soil type, humus
form, and associated microbial activity, etc. Humus layer dynamics are not modelled in
SORTIE-ND so it was not included in the final linker function set.
Understory vegetation influences lichen abundance through litter fall in addition to
competition for light, space, and resources (Cornelissen et al. 2001). Increases in vegetation
cover from dwarf shrub species such as kinnikinnick and dwarf Vaccinium species have been
correlated with reductions in terrestrial lichen communities following mountain pine beetle
attack (Seip and Jones 2010, Cichowski and Haeussler 2020). Understory vegetation other
than tree regeneration is not currently modelled in SORTIE-ND. I tested a linker function
that included a function representing the suitability of a site to support low-growing (e.g.,
dwarf shrubs) vegetation (see Table 1, model “Full Model (C.Veg)”). However, limiting the
suite of linker functions to incorporating SMR and SNR directly had greater support in the
data. Therefore, I selected the more parsimonious suite of linker functions that did not
incorporate a linker function representing the presence of low-growing vegetation.
It was challenging to develop one index to represent the total amount of browse for moose.
This required combining a diverse array of species into one response variable; the linker
function needed to represent species that differed in edaphic niche and regeneration strategies
(Figure 6, Beaudry et al. 1999). I did not have a large enough dataset to develop linker
functions for the individual species that were browsed by moose. I also did not identify
willow to the species-level due to differentiating features such as catkins being unavailable at
the time of sampling. In general, however, greater diversity and abundance of moose shrubs
was associated with higher productivity sites. Schrempp et al. (2019) also found higher
abundance of moose shrubs on higher productivity sites.

46

I expected light to be an important predictor variable of moose shrub suitability. However,
the likelihood analysis revealed that adding light to the linker function did not improve
support for a model that already incorporated conifer basal area. It should be noted that
incorporating conifer basal area rather than light results in a linker function that is not
sensitive to the effect of edge on understory species. There is a gradation of light conditions
between full light in the center of a large gap, partial light at the edge of gap adjacent to a
forest edge, and low light under a full forest canopy (Chen et al. 1995). This gradual change
in light and corresponding vegetation growth, is not reflected by simply incorporating conifer
basal area in the immediate vicinity. The suitability for moose shrubs under a forest canopy,
but in a location with higher light levels due to an adjacent gap, may be underestimated.
Alternatively, suitability may be overestimated for a location in a gap that is directly adjacent
to a forest canopy that provides shade. The presence and abundance of understory vegetation
species are influenced by proximity to edge as well as the orientation of the site relative to the
edge (i.e., south side vs. north side of a gap) (Canham 1988, Matlack 1994, Nelson and
Halpern 2005).
Moose often select mixedwood and broadleaf stands in winter (Osko et al. 2004, Street et al.
2015, Mumma et al. 2021) and mixedwood and broadleaf stands can provide a greater
quantity of moose forage (Crête and Courtois 1997). I did not explicitly include the broadleaf
component of the stand in the linker function because I was unable to find evidence in the
literature providing a direct mechanistic link between broadleaf tree presence and moose
shrubs. Interestingly, the support for conifer basal area rather than light, does indirectly
reflect the presence of the broadleaf component in a stand. A mixed stand with a broadleaf
component would have a reduced conifer basal area in comparison to a comparable pure
conifer stand with the same number of total stems.

47

Past studies of wildlife habitat often have focused on landscape-scale analyses and
applications (Whitman et al. 2017, Stockdale et al. 2019, Stewart et al. 2020, Leblond et al.
2022). This modelling framework, projecting habitat dynamics using a stand dynamics
model and linker functions, can be used to inform landscape-scale habitat models. Modelling
of the forest cover component of caribou habitat has historically been static or dynamic at the
landscape scale (e.g., resetting stand age and cover type due to wildfire or timber harvest or
tracking change in stand age through time but not representing the corresponding changes in
forest cover) (Johnson et al. 2015, Barber et al. 2018, Nagy-Reis et al. 2021, Leblond et al.
2022). Some disturbances affect only a portion of the stand or vary in terms of the intensity
of disturbance (e.g., different intensities of mountain pine beetle attack or partial harvest).
Employing process-based, individual-tree, stand dynamics models provides the ability to
evaluate the effect of these complex disturbances on stand structure and, via linker functions,
on habitat dynamics (Holt et al. 1995, Maleki et al. 2021). Scaling output from this modelling
framework to the landscape scale could enable modelling of habitat dynamics customized to
the dynamics of unique stands on the landscape (Barber et al. 2018, St-Laurent et al. 2022).
There are elements of caribou habitat that are beyond the scope of this study because they are
best modelled at a landscape scale. For example, linear disturbances play a key role in
predation risk (Dickie et al. 2017, DeMars and Boutin 2018, Mumma et al. 2018), but cannot
be effectively modelled at the stand scale by SORTIE-ND. Roads have been shown to
improve the hunting efficiency of wolves through improved access to caribou habitat and
predation refugia as well as increasing the efficiency at which they hunt (Dickie et al. 2017,
DeMars and Boutin 2018). The effect of linear features on habitat value could be modelled
with a buffer along the linear corridor representing the zone of avoidance of these features
(Johnson et al. 2015). However, at the scale at which stands are modelled in SORTIE-ND (a
300 x 300-m or 500 x 500-m polygon), applying a buffer of over 1 km is not feasible.
48

SORTIE-ND can be used to model the dynamics of stand revegetation along linear features
such as roads or skid trails and to model the effect of higher light levels on the stand
dynamics of the adjacent stands.
I developed a simple but robust modelling framework to relate stand structure to three key
ecological elements of the winter habitat of northern mountain caribou: terrestrial forage
lichen, moose shrubs, and vertical cover. I developed linker functions that can be used to
relate the output from forest simulation models to habitat suitability. As with all models there
are limitations, but these limitations should not shadow the contributions that these linker
functions make for providing insight into the effects of operational forest management
decisions on habitat for a species of conservation concern. Where disturbance events result in
stands with complex structure, informing landscape-scale habitat models with habitat
dynamics that are stand and site specific would further improve our understanding of the
overall dynamics of habitat across landscapes and longer time scales.

49

Chapter Three
Modelling the Effects of Disturbance on Caribou Winter Habitat
3.1 Introduction
Abiotic and biotic disturbances can significantly change forest ecosystems and substantially
alter the ecosystem services that forests provide (Buma 2015). The impact of disturbances on
forest ecosystems will vary with the type of disturbance (Thom and Seidl 2016, Viljur et al.
2022), the severity and magnitude of the disturbance (Seidl et al. 2016), and the ecosystem
services being evaluated (Thom and Seidl 2016, Sánchez et al. 2021). Disturbances will
often immediately impact the provisioning of forest ecosystem services, but can also
influence the recovery trajectory of the ecosystem and rate at which ecosystem services are
regained (Seidl et al. 2016, Bartels et al. 2016, Sutherland et al. 2016, Bergeron et al. 2017).
For example, a severe wildfire removes most vegetation from the landscape, significantly
reduces the depth of forest floor, and adds a significant pulse of nutrients. One of the
outcomes of that fire could be a reduction in habitat for wildlife that depend on older trees
and associated forest structure. As a corollary, numerous wildlife species will benefit from
the early successional plant communities that establish following the fire (Viereck and
Schandelmeier 1980).
Managing for diverse ecosystem services is complex and interacting disturbances add further
challenges (Schwenk et al. 2012, Mönkkönen et al. 2014, Bowd et al. 2021). For example,
forest canopy removal is positive for some species but negative for others and the response
can differ depending on the pattern and intensity of canopy removal (Halpern et al. 2012). In
many regions, planning decisions are made for landscapes where not only timber harvest but
also natural disturbance (e.g., wildfire, insect infestations) influences the provision of

50

ecosystem services. Interactions between disturbances, management actions, and
successional forest dynamics add further complexity (Seidl et al. 2011, Elkin et al. 2013,
Buma 2015). We need the ability to project the influence of both management decisions and
natural disturbances to support the sustainable provision of these services (Seidl et al. 2011,
Whitman et al. 2017).
Numerous factors play a role in the response of ecosystems to a particular disturbance. At a
landscape scale, these include past disturbance history, current disturbance severity and
intensity, and landscape structure (Racey et al. 1996, Angelstam and Kuuluvainen 2004).
Type and intensity of disturbance influence the species composition and structure of the
regenerating stands following disturbance. For example, shade tolerant species are more
likely to survive or regenerate following a low intensity fire or MPB attack (Astrup et al.
2008b, Kielland and Brown 2015). At a stand scale, factors such as stand structure, species
composition, and edaphic site conditions influence the response and recovery of the
ecosystem and ecosystem services (Brulisauer et al. 1996, Clark et al. 2003, Lilles et al.
2018). Understanding how these factors interact and influence the response of ecosystem
features is required to manage and maintain these features following disturbance.
Caribou is an ungulate species of conservation concern in Canada with habitat that has been
affected by a variety of disturbance agents (COSEWIC 2014). Habitat change due to both
anthropogenic factors and natural disturbance are considered ultimate factors in the decline of
mountain caribou populations that reside throughout the interior of British Columbia, Canada
(COSEWIC 2014). The forests of the central interior of BC have experienced the largest
mountain pine beetle epidemic in Canadian history (Dhar et al. 2016) and numerous large
wildfires (Cichowski and Haeussler 2020, Province of British Columbia 2023). Extensive
timber salvage resulted in the removal of a large proportion of pine leading stands that serve

51

as winter habitat of caribou. Therefore, to ensure the sustainable provision of caribou winter
habitat, managers need to understand the effects of disturbance and the recovery trajectories
of critical habitat elements.
Forest managers require tools that support tactical and strategic land-use decisions for timber
harvest within caribou habitat. Models allow us to project what we know across longer time
frames and integrate the interactions of multiple factors. As such, models provide useful tools
for exploring the impacts of disturbance on ecosystem services (Seidl et al. 2011). Dynamic
forest models have a long history of development and use for improving our understanding of
stand dynamics and projecting changes in tree populations and stand structure over time
(Bugmann 2001, Bugmann and Seidl 2022). These models were developed to quantitatively
represent stand dynamics and are typically focused on the demographics of individual trees or
stands of trees (Bugmann and Seidl 2022). Less common are models that project the
response of non-tree ecosystem elements to changes in stand structure (e.g., the response of
specific ungulate habitat attributes such as terrestrial lichen to changes in stand structure due
to differing disturbance events such as clearcutting, partial harvest or MPB attack). A tool
that models habitat dynamics at the stand scale provides the ability to evaluate the effects of
alternative management strategies over time.
SORTIE-ND is a complex stand dynamics model that can be used to project the response of
stands to harvest or disturbance over a timeframe of hundreds of years. As introduced in
Chapter 2, SORTIE-ND enables the modelling of mixed-species stands and stands with
complex vertical structure including snags. Linker functions provide a link between stand
models and ecosystem features that enables modelling ecosystem features that are responsive
to changes in stand structure. By providing the linkage from stand structure to ecosystem
features, stand dynamics models such as SORTIE-ND can be used to project not only stand

52

dynamics following disturbance but also the dynamics of ecosystem features, such as
important habitat elements. Projecting the response of habitat elements over time can inform
management decisions on a disturbed landscape.
Forests in the central interior are valued for timber production as well as providing habitat for
wildlife such as caribou (Cichowski 1996). The ability of forests to grow trees and provide
wildlife habitat is expected to depend on forest management decisions and the legacies of
disturbances such as MPB and wildfire. In order to develop improved forest management
strategies we need to improve our understanding of how forest composition and structure
influences caribou habitat, and how stand dynamics will alter caribou habitat value through
time. I applied linker functions to SORTIE-ND model output to evaluate the influences of
different disturbance events on the dynamics of elements of caribou habitat. I examined three
specific questions:
1.

What is the influence of edaphic site on the provision of elements of
caribou habitat over time?

2.

What is the influence of stand type on the provision of elements of caribou
habitat over time?

3.

How do recovery trajectories of caribou and moose habitat elements differ
in response to disturbances such as mountain pine beetle and clearcut
harvest?

I evaluated these three questions by applying them to three important elements of caribou
winter habitat: terrestrial forage lichen, moose shrubs, and vertical cover.

53

3.2 Methods
3.2.1 Study Area and Caribou Habitat Context
The focus of this study is the winter habitat of the Tweedsmuir-Entiako caribou herd (TEC).
This herd of northern mountain caribou is in decline and is blue-listed or designated as
threatened by the Province of British Columbia (BC) (COSEWIC 2014, Min. of Env. and
Climate Change 2018). The herd’s home range is located in west-central BC southwest of
Prince George (Figure 1). The herd’s winter habitat has been heavily affected by MPB,
wildfire, and timber harvesting. Maintaining both healthy caribou herds and sustainable
timber supplies are provincial and federal priorities.
The study area is defined by the Sub-Boreal Spruce moist cold Babine variant (SBSmc2), a
biogeoclimatic subzone and variant described by the local biogeoclimatic ecosystem
classification system (BEC) (Figure 1, Banner et al. 1993). Chapter 2 includes a description
of the geography and forests that are characteristic of this subzone. Although the TEC winter
habitat extends across several BEC subzones, the SBSmc2 was selected for this study
because it comprises a portion of the winter habitat that is subject to industrial forest
management and caribou winter habitat management.
Three elements of TEC winter habitat are the focus of this paper: terrestrial lichens favoured
as caribou forage (henceforth referred to as forage lichen or lichen), shrubs favoured by
moose for winter forage (moose shrubs), and vertical canopy cover. These elements are
described in detail in the methods section of Chapter 2. Other stand-scale habitat elements
are important for caribou during winter and throughout other seasons, however this study is
focused solely on the three elements of forested winter habitat described above.

54

3.2.2 Stand Modelling
SORTIE-ND simulates changes in tree populations across small stands (generally less than
nine hectares) over time. As a spatially explicit model, output includes not only individual
tree size, regeneration, growth and mortality information but also the spatial distribution of
the trees and the resulting light environment in a grid across the entire stand (Canham et al.
2004). A snag sub-model enables modeling of snag dynamics (Bose et al. 2015). SORTIEND was parameterized to model four edaphic sites of the SBSmc2 by applying growth
relationships derived from previous studies (Lilles and Astrup 2012, Coates et al. 2013).
SORTIE-ND plots can be conceptualized as a torus (donut) where the trees on the apparent
‘edge’ of a plot are neighbours of trees on the opposite ‘edge’ of the plot. Therefore, trees on
the western edge of the plot are modelled as growing next to the trees located on the eastern
edge (at the same latitude); the same relationship holds for trees on the north and south edges
of plots. Thus, 9ha clearcut plots model the interior condition of a stand with no forest edge.
I modelled four stand types crossed with four edaphic sites on a 9-ha plot. The stand types
were categorized as: lodgepole pine dominant, mixed conifer and broadleaf (lodgepole pine,
hybrid spruce, subalpine fir and trembling aspen), mixed conifer (lodgepole pine, hybrid
spruce and subalpine fir), and hybrid spruce dominant. Throughout this paper, these stand
types will be referenced as pine, conifer/broadleaf, conifer mix, and spruce (respectively) and
individual tree species as spruce, pine, subalpine fir and aspen. The four edaphic sites are
situated along a productivity gradient from a poor subxeric site to a rich subhygric site and
referenced as: subxeric, submesic, mesic, and rich subhygric. All four stand types were
modelled on submesic, mesic, and rich subhygric sites. Spruce-dominated stands aren’t
found on subxeric sites and therefore were not modelled on this site.

55

Edaphic conditions at a given forest site can give rise to a range of stand conditions (i.e.
species composition, age class distribution, DBH distribution) as a result of interactions
between stand age, disturbance history and the forest’s regeneration pathway (Meidinger and
Pojar 1991, Banner et al. 1993, Coates et al. 1994). I evaluated the impact of stand type by
simulating a range of forest compositions that naturally occur under a given edaphic
condition. Initialization of each stand type was based on tree lists from the empirical plot
data described in Chapter 2. I selected one field plot that was representative of each
stand/edaphic site combination (Table 5). I used the following criteria to select stands: 1) the
trees in the dominant canopy were between 75 and 125 years of age; 2) all stands within each
site type had similar basal areas; 3) for pine and spruce stands, the pine or spruce component
(respectively) comprised at least 60 percent of the total basal area; 4) for conifer/broadleaf
stands, a significant component of aspen was required but the stand could not be dominated
by aspen (I targeted 30–60% aspen by basal area); and 5) conifer mix stands had no broadleaf
component and a relatively even representation (by basal area) of all three conifer species.
Some stand/edaphic site combinations did not have suitable field plots so the tree list from a
stand for a similar edaphic site was used (the tree list for the mesic spruce stand was used to
initiate the tree composition of the submesic spruce stand, the submesic conifer/broadleaf
tree list was used for mesic conifer/broadleaf stand, and the mesic pine tree list was used for
the subhygric pine stand). Stands were initiated as pre-mountain pine beetle attack stands so
all pine found in the field plots (live and dead from MPB attack) were initiated as live pine.
Three scenarios were modelled: no disturbance, mountain pine beetle attack (MPB), and
clearcut harvest. All scenarios were simulated for 200 timesteps (years). The scenarios were
initiated with the stand conditions described in the previous paragraph. For the MPB
scenario I applied 100% mortality to all pine greater than 15 cm dbh and 50% mortality to all

56

Table 5. Structure of stands used to initiate SORTIE-ND model scenarios for mature stands across
gradients of edaphic site and stand type. This stand structure data was drawn from data collected
during the summer of 2019 in sub-boreal forests of west-central British Columbia, Canada.
Site Type

Site
Series

Stand Type

Species
composition
(% basal area)

Basal
Area
(m2/ha)

Conifer
Basal
Area
(m2/ha)

Broadleaf
Basal
Area
(m2/ha)

Density
(> 5 cm
dbh)

Light
(%,
GLI)

Age
(yrs)

subxeric

02

pine

Pl 100

28.2

28.2

0

1000

42

88

subxeric

02

conifer/broadleaf

23.8

16.6

7.2

700

35

77

subxeric

02

conifer mix

13.2

13.2

0

750

48

77

submesic

01c

pine

28.0

28.0

0

800

26

76

submesic

02-01

conifer/broadleaf

22.4

8.7

13.7

750

33

95

submesic

01c

conifer mix

47.7

47.7

0

1800

22

88

submesic

01

spruce

38.1

38.1

0

3400

22

87

mesic

01

pine

Pl 37 At 31
Bl 20 Sx 12
Bl 44 Sx 32
Pl 24
Pl 79 Sx 18
Bl 3
At 60 Pl 30
Bl 10
Bl 50 Pl 41
Sx 9
Sx 61 Pl 33
Bl 6
Pl 98 Sx 2

37.9

37.9

0

1900

34

121

mesic

02-01

conifer/broadleaf

At 60 Pl 30
Bl 10

22.4

8.7

13.7

750

33

95

mesic

01

conifer mix

17.9

17.9

0

600

22

76

mesic

01

spruce

38.1

38.1

0

3400

22

87

subhygric

01

pine

Bl 44 Pl 29
Sx 27
Sx 61 Pl 33
Bl 6
Pl 98 Sx 2

37.9

37.9

0

1900

34

121

subhygric

05

conifer/broadleaf

Sx 69 At 31

43.3

30.5

12.8

650

15

84

subhygric

05-06

conifer mix

38.4

38.4

0

3750

10

108

subhygric

05

spruce

Bl 52 Sx 37
Pl 11
Sx 100

63.7

63.7

0

2700

15

92

trees between 10 and 15 cm dbh ( SORTIE-ND randomly selected 50% of the trees in the 10
to 15 cm tree population). This differential level of attack across pine diameter classes
mimicked the intensity of MPB attack observed during the recent MPB epidemic in BC
(Safranyik 2004, Dhar et al. 2015). The MPB attacked trees became snags, thus, decay and
57

fall-down rates were simulated. All trees were removed in the clearcut scenario and the
establishment of an even-aged plantation was simulated. Seedlings were planted at a density
of 1800 stems/ha (average height of seedlings was 21 cm; see Table 6). Natural regeneration
of all conifer species (seed rain) was simulated for all stand types in each timestep. Natural
regeneration of aspen was included for the conifer/broadleaf stand type only.
I tested the amount of variation in output that resulted from stochastic elements within
SORTIE-ND to determine if multiple simulations should be run for each
stand/site/disturbance combination. Stochastic elements included the location of individual
trees, the assignment of tree diameters (within a diameter class) to trees of each species, the
probability of mortality, the spatial distribution of seed rain and the probability of seed
establishment. To do this, I conducted 10 replicate runs of a clearcut harvest scenario for each
stand type and edaphic site combination. The variance between the 10 replicated simulations
was minor. Therefore, I determined one simulation for each treatment/edaphic site/stand type
combination was sufficient because limited additional information on variability was gained
from 10 simulations.
Table 6. Stand establishment characteristics for sub-boreal clearcut scenarios modelled with
SORTIE-ND including planting densities (in stems/ha) and natural regeneration (‘seed’).

Stand type

Hybrid Spruce

Lodgepole
Pine

Subalpine Fir

Trembling
Aspen

Pine

seed

1800 + seed

seed

0

Conifer/Broadleaf 700 + seed

700 + seed

seed

400 + seed

Conifer Mix

400 + seed

1000 + seed

400 + seed

0

Spruce

1400 + seed

400 + seed

seed

0

58

Output from all simulations consisted of tree maps and light maps at each timestep. A tree
map is a list of all trees (live and dead) in the simulated 9-ha stand with the following data for
each tree: tree species, live/dead indicator, dbh, diameter at 10 cm height, height, and location
(x and y coordinate). The light map provided the Gap Light Index (GLI, see Chapter 2
Methods for a description) estimated at 1.3 m height for each 1×1-m grid cell in the simulated
stand.
The calibrated linker functions presented in Chapter 2 were applied to the SORTIE-ND
output to project the dynamics of the three habitat elements. The linker functions related
stand metrics (light, basal areas of conifer and broadleaf trees, densities of aspen and
subalpine fir saplings, edaphic site condition, and time since disturbance) to the suitability of
a site to support terrestrial lichen, moose shrub communities, and vertical cover. Each
timestep of SORTIE-ND output was rasterized into 10 x 10-m rasters using the raster
package (v3.4-13; (Hijmans 2021) in R (v. 4.1.2; R Core Team 2021). All seven stand
metrics were calculated for each raster. The linker functions were used to calculate the
suitability indices for each raster using the applicable stand metric value. The median
suitability index of the stand and the 5th and 95th quantiles were reported. The 5th and 95th
quantiles provide a measure of the variability of the data around the median for data that may
not be normally distributed.
The terrestrial lichen and moose shrub linker functions were comprised of component linker
functions. Each component linker function represented the influence of one stand variable on
the suitability of the stand to support the habitat element. The component linker functions
were combined multiplicatively to project the relative suitability for terrestrial lichen and
moose shrubs across the plot (scaled from 0 to 1) (Equation 1 and 7). The vertical cover

59

linker function was a single function projecting relative canopy cover based on conifer basal
area (Equation 13).
The linker functions projected a relative measure of suitability of a stand. They did not
project the expected abundance of each habitat element in a stand. Other factors are known
to influence the abundance of terrestrial lichen and moose shrubs (e.g., disturbance history,
abundance of competing vegetation). Because I developed the linker functions with the
specific intent of applying them to output from the stand model SORTIE-ND, the predictor
variables were constrained to those available from SORTIE-ND. The linker functions were
not predictive, ecological models of abundance of terrestrial lichen or moose shrubs.
3.3 Results
3.3.1 Component Linker Functions
Applying the component linker functions to SORTIE-ND output allowed me to parse out the
effects and interactions of individual stand characteristics on lichen and moose shrub
suitability over time. Lichen dynamics were best illustrated by focusing on subxeric sites
(Figure 8, Appendix 5) and moose habitat dynamics were best illustrated on the subhygric
sites (Figure 9, Appendix 5).
The full lichen linker function was comprised of four component functions representing the
combined effects of stand-level determinants on the suitability of a stand to support terrestrial
lichen: edaphic site, time since stand-replacing disturbance, understory light, and overstory
broadleaf abundance. Edaphic effect was static and did not change over time (Figure 8).
Time since stand-replacing disturbance also appeared static for mature and MPB-attacked
stands, reflecting the length of time since a severe stand-replacing disturbance (Figure 8).
Applied to the clearcut scenario, the time since disturbance linker function reflected the

60

Figure 8. Output of the component linker functions for forage lichen on subxeric sites. The top four
rows present the output of the component linker functions, with the full suitability index for forage
lichen shown in row five. Coloured lines represent different stand types. Where only the conifer mix
stand suitability appears on the figure, the other stand types have the same suitability (e.g., the
broadleaf effect for the mature pine-dominated stand is the same as the conifer mix stand).

physical damage to lichen communities from timber harvest and the time required for lichen
communities to re-populate and regenerate on the site. The light effect was dynamic
reflecting changes in stand structure including relative tree species abundance and growth.
The light effect decreased with increasing canopy closure for conifer mix and
conifer/broadleaf stands. The pine stand exhibited an increase in light effect over time in
response to tree mortality due to senescence. The broadleaf effect suggested that once
broadleaf species were well-established in a mature or MPB attacked stand, the stands were

61

unsuitable for lichen. Following clearcutting, the suitability for lichen was reduced during
mid-succession for stands with a low-density broadleaf component. Considering the
combination of all these variables, the best subxeric stands for lichen (> 0.35 suitability) were
mature conifer mix stands before the slow growing trees matured further, pine stands
experiencing mortality due to senescence, and clearcut stands of all types 20–45 years old.

Figure 9. Output of the component linker functions for moose browse on subhygric sites. The top
three rows present the output of the component linker functions for moose shrubs, with the full
suitability index for moose shrubs shown in row four. The suitability of a stand for subalpine fir and
aspen browse are presented in rows five and six and the total moose browse suitability of a stand in
row seven. The bottom row illustrates the suitability of a stand to provide vertical cover. Coloured
lines represent different stand types. Where only the spruce stand suitability appears on the figure,
the other stand types have the same suitability (e.g., vertical cover, bottom row).

62

The full moose shrub linker function was comprised of three component functions (top 4
panels of Figure 9). Similarly to terrestrial lichen, the edaphic effect was static over time.
The time since severe disturbance curve showed a slow decline in the suitability of mature
and MPB stands to support early-seral moose shrubs as they aged whereas clearcuts reached
high suitability by year 20 as shrubs regenerated following harvest. The competitive effect of
conifer basal area meant that mature stands declined in suitability as their basal area
increased, MPB attacked pine stands had high suitability for about 50 years after attack (until
their basal area recovered), and clearcut stands had high suitability that started to decline
around 50 years. Overall for moose shrubs, the combined effects of the component functions
suggested that the best subhygric stands were 5 to 50-year-old clearcuts (suitability > 0.35).
The full moose browse response combined moose shrub, aspen sapling and subalpine fir
sapling suitabilities (Figure 9: 4th to 7th panels). Aspen and subalpine fir browse was above
zero only at the very beginning of the scenario for mature and MPB stands and immediately
following clearcutting. Disturbance resulted in good conditions for natural regeneration and
release of advanced regeneration. However, the regenerating trees quickly grew beyond a
height where they provided forage opportunities for moose. The suitability function for full
moose browse showed very little difference between treatments, stand type or site over a 200year timeframe because the shrub response was weakened by averaging it with the weak
responses of fir and aspen. Therefore, the full moose browse suitability function will not be
further evaluated in this thesis (figures are available in Appendix 6).
The vertical cover linker function reflected the dynamics of conifer basal area (Figure 9).
Thus, mature stands provided full suitability for cover. MPB attacked pine stands exhibited a
reduction in vertical cover until the basal area recovered after 50 years. Clearcut stands, on

63

these high productivity sites, also required approximately 50 years to develop suitable
vertical cover for moose.
For the remainder of the results section I will evaluate the dynamics of habitat suitability by
focussing on the response of the full lichen, full moose shrub, and vertical cover linker
functions.
3.3.2 Terrestrial Lichen Dynamics
Effect of Site: Poorer productivity sites provided much greater lichen suitability than richer
and wetter stands (Figure 10). However, an interaction between stand type and site was
evident where conifer/broadleaf stands exhibited poorer suitability regardless of edaphic site.
Improvements in lichen suitability following both MPB and clearcutting were greater and
extended for longer as stands moved along the gradient from higher to lower productivity
sites (subhygric to subxeric).
Effect of Stand Type: Differences in lichen suitability related to stand type became evident
as site productivity declined (Figure 10). The greatest difference was on subxeric sites,
however on submesic sites, differentiation between stand types was also evident. In general,
subxeric and submesic pine stands provided the highest suitability for lichen. However the
subxeric pine stand initially had lower suitability than the conifer/broadleaf and conifer mix
stands. This reflected the greater conifer basal area of the subxeric pine stand compared to
the other stand types. Following mountain pine beetle, pine stands provided greater
suitability for lichen because the larger proportion of pine in the overstory (compared to other
stand types with less pine) resulted in a greater proportion of the overstory being killed. The
resulting increase in understory light was beneficial for lichen. There was limited effect of
stand type following clearcutting for stands managed primarily for timber production. These
stands were initiated with a relatively high density of conifers and only 400 stems per hectare
64

of aspen. Following clearcutting, all stand types initially increased in lichen suitability but
the pulse was short lived. Although conifer/broadleaf clearcuts (with relatively low aspen
density) initially also provide improved lichen conditions, in the mid-term (50 to 125 years)
suitability is reduced in comparison to other stand types.

Figure 10. The projected response of lichen suitability over time in undisturbed mature forest (1st
column), after MPB disturbance (2nd column) and after clearcut (3rd column) across four edaphic
sites (rows) and four stand types (colours). The 4th column illustrates the difference in lichen
suitability between a MPB attacked stand and an undisturbed stand and the 5th column shows the
difference between a clearcut stand and a MPB attacked stand.

65

Effect of Disturbance: MPB and clearcutting improved lichen suitability in certain stand
types (Figure 8 and Figure 10) due to the reduction of the live overstory. However, MPB only
improved suitability significantly in pine stands (compared to mature stands) because the
proportion of susceptible host (pine) was limited in the other stand types; thus, the effect of
MPB was also limited in comparison to mature stands. Interestingly, terrestrial lichen
suitability was greater for drier MPB affected pine stands than clearcuts (subxeric and
submesic sites) and suitability extended for a longer time frame following MPB attack.
Clearcutting, through complete canopy removal, elevated the suitability of all stand types to
support lichen in the short term, though only until canopy closure occurred (Figure 10).
Within-stand variability in lichen suitability differed across site, stand, and disturbance types
(Appendix 7). The greatest variability was evident in subxeric and submesic pine stands
(mature and MPB-attacked) and submesic conifer/broadleaf stands. In general, variability
within clearcuts and spruce stands was low (except clearcut conifer/broadleaf stands)
however, a pulse in variability did occur between 25 and 50 years following disturbance in
conjunction with crown closure.
3.3.3 Moose Shrub Dynamics
Effect of Site: Across all disturbance scenarios, moose shrub habitat was least suitable on the
poorer sites, regardless of stand type (Figure 11, Appendix 6, Figure 29). Following
clearcutting, moose shrub suitability declines as site productivity declines. Time interacted
with site effect following disturbance: as site productivity increased, the period of time that
clearcuts and MPB stands supported greater suitability for moose shrubs declined (though
subhygric and mesic stands were similar). In general, for mature stands (timesteps 0 to 50
which equated to approximately 80- to 130-year-old stands), subxeric and subhygric sites
exhibited poorer moose shrub suitability than submesic and mesic sites.

66

Figure 11. The projected response of moose shrub suitability over time in undisturbed mature forest
(1st column), after MPB disturbance (2nd column) and after clearcut (3rd column) across four edaphic
sites (rows) and four stand types (colours). The 4th column illustrates the difference in moose shrub
suitability between a MPB attacked stand and an undisturbed stand and the 5th column shows the
difference between a clearcut stand and a MPB attacked stand.

Effect of Stand Type: All mature conifer stand types without a broadleaf component
provided lower suitability for moose shrubs than mature conifer/broadleaf stands (Figure 11).
Although all mature stand types exhibited a decline in moose shrub suitability over time,
conifer/broadleaf stands retained elevated suitability for longer than the other stand types.
For the first 75 years following MPB attack, pine stands provided equivalent or higher moose
shrub suitability than conifer/broadleaf stands (except on mesic sites). MPB attacked pine
67

stands provided elevated moose shrub suitability in comparison to other MPB attacked
conifer stands on all sites except mixed conifer stands on mesic sites. In mature stands there
was not a consistent differentiation in suitability between the conifer stand types (e.g., mature
mixed conifer stands provided the poorest suitability on submesic sites but greater suitability
on mesic, subxeric, and subhygric sites).
Effect of Disturbance: Disturbance (MPB and clearcutting) increased moose shrub
suitability; the complete removal of overstory canopy from clearcutting resulted in a greater
relative effect compared to MPB (Figure 11). MPB attack had the greatest influence on
moose shrubs in pine stands due to the larger proportion of pine in the stands.
There were no strong trends in within-stand variability for moose shrubs in mature and MPBaffected stands (Appendix 7). However, for clearcut sites within-stand variability decreased
as sites decreased in productivity. For the first 100 years, the greatest within-stand variability
for moose shrubs was evident in mature and MPB-attacked conifer/broadleaf stands on
subhygric sites.
3.3.4 Vertical Cover Dynamics
Effect of Site: Following disturbance, stands on poorer sites developed vertical conifer cover
more slowly (Figure 12). This reflected the slower growth of trees on those poor productivity
sites. There was an interaction with stand type where, on higher productivity sites, there was
little differentiation between stand type. However, as site productivity declined, the influence
of stand type on the recovery trajectory of vertical cover grew stronger. On the poorest sites,
conifer/broadleaf stands developed vertical cover most slowly.
Effect of Stand Type: Across most stand types and sites, vertical cover suitability was least
for conifer/broadleaf stands (Figure 12). However, variability was also greatest for this stand
type, so portions of these stands did have the potential to provide greater cover in addition to
68

less cover in other areas (Appendix 7). On the richest sites, although conifer/broadleaf stands
provided full cover suitability (median at 1.0), there remains very high variability within the
stand such that portions still provided low cover. Spruce stands consistently provided the
highest cover of all stand types and the least within-stand variability.

Figure 12. The projected response of suitability to provide vertical cover over time in undisturbed
mature forest (1st column), after MPB disturbance (2nd column) and after clearcut (3rd column) across
four edaphic sites (rows) and four stand types (colours). The 4th column illustrates the difference in
suitability to provide vertical cover between a MPB attacked stand and an undisturbed stand and the
5th column shows the difference between a clearcut stand and a MPB attacked stand.

Effect of Disturbance: Not surprisingly, both forms of disturbance reduced vertical cover.
MPB had an effect proportional to the amount of pine in the stand. Therefore, pine stands
69

had the greatest reduction in cover following MPB attack. Conifer/broadleaf stands exhibited
a reduction in vertical conifer cover following MPB attack that extended for a longer period
of time than the other stand types (except on xeric sites where pine stands had a longer period
of time with reduced cover). Clearcutting reduced cover below MPB-attacked cover values
for a brief period following MPB attack before planted stands regained cover value.
However, in some poorer productivity stand types (drier pine (subxeric and submesic) and
submesic and mesic conifer/broadleaf stands), MPB-attacked stands sustained reduced cover
for longer than similar clearcut stands (for 25 to 125 years following disturbance).
3.4 Discussion
Model results show that forest disturbances differ in their effects on habitat elements and
their recovery trajectories. Results also show that site and stand structure provide important
context when evaluating the recovery of habitat elements. Managers should incorporate site,
stand structure, and disturbance type when determining the best approach to managing for the
recovery of habitat elements following disturbance events. On a mountain pine beetle
affected landscape where recovery of caribou winter habitat is a priority, these model results
indicate that certain forest management decisions, such as reserving MPB attacked subxeric
stands from harvest, will promote the provision and recovery of caribou winter habitat more
effectively than others.
Model results indicated that low productivity MPB attacked pine stands provided the greatest
lichen suitability post-disturbance and suitability remained elevated for the longest period of
time. MPB attacked stands were also less suitable for moose shrubs than clearcuts. In
addition, low-productivity MPB attacked pine stands were projected to be slower to provide
vertical cover when compared to clearcut stands. Comparing the effects of MPB versus
clearcutting on all these habitat elements suggests that leaving low productivity MPB
70

attacked pine stands to recover naturally has greater value for caribou winter habitat than
clearcutting the stands.
Our understanding of lichen ecology supports the finding that MPB attacked stands on
subxeric sites have the potential to support suitable conditions for terrestrial lichen.
Terrestrial lichens are well adapted to dry poor sites as they are not dependent on soil to
provide moisture and nutrients (Rowe 1984). Therefore, lichen have a competitive advantage
on low productivity sites where resources aren’t available for competing vegetation (Rowe
1984, Payette and Delwaide 2018). MPB disturbance does not physically damage lichen
whereas lichen may be damaged by machinery or smothered by debris that falls to the ground
during tree harvest and processing (Coxson and Marsh 2001, Waterhouse et al. 2011). In
MPB attacked stands, the shade and competition from surviving trees and germination
substrate limitations suppress natural regeneration (Astrup et al. 2008b) so stands may not
achieve crown closure as quickly as fast-growing plantations following clearcutting.
Although these considerations indicate that poor productivity MPB attacked stands should
provide suitable conditions to support terrestrial lichen, long-term studies monitoring the
effect of MPB on lichen have found lichen declined in the first 10 years following MPB
attack (Cichowski and Haeussler 2020, Nobert et al. 2020). This decline was associated with
a positive response of understory vegetation following MPB attack (Cichowski and Haeussler
2020). After 10 years, the lichen decline stopped and abundance appeared to stabilize on
subxeric sites but not submesic or mesic sites (Cichowski and Haeussler 2020). The
discrepancy between the model projection and these empirical results indicate that the linker
functions may not represent this competitive effect. More time is required to determine if the
initial reduced lichen suitability is sustained or if, following a lag, suitability increases as
projected over the 50–100-year time frame.

71

Studies monitoring lichen abundance following MPB attack have recommended protecting
poor productivity stands from harvest to minimize the impact of the MPB on low-elevation
caribou winter habitat (Seip and Jones 2010, Nobert et al. 2020, Cichowski et al. 2022).
These studies found that, although lichen declines post-MPB attack, lichen abundance
remains greater than in corresponding clearcuts (Nobert et al. 2020, Cichowski et al. 2022).
These studies determined that clearcuts created conditions suitable for lichen recovery and
expansion, however, the recovery trajectory was slow and uncertain. Cichowski and
Haeussler (2020) found that forage lichen communities in clearcut stands had peaked and
were declining 16–17 years post-harvest. They attributed the decline in lichen to growth of
regenerating conifers; a similar finding to Lewis et al. (2019). Interestingly, numerous studies
have found that although lichen is positively correlated with greater understory light levels in
mature stands, as the percent of trees attacked by mountain pine beetle increases, there is a
greater corresponding decline in terrestrial lichen abundance (Cichowski and Haeussler 2020,
Nobert et al. 2020, Cichowski et al. 2022). This could also be attributed to competitive
interactions with understory vascular plant and bryophyte communities. Despite reduced
abundance of forage lichen in stands, caribou continued to forage in these stands for up to six
years post-MPB attack (Cichowski 2010, Seip and Jones 2010).
These results could inform harvest practices: where harvest is planned within caribou winter
habitat, employing harvest patterns that mimic mountain pine beetle conditions could
promote the maintenance of caribou habitat (e.g., uniform partial harvest). Model results
showed that stand thinning due to MPB caused tree mortality suppressed moose shrubs and
provided conditions more suitable for forage lichen than those in a clearcut or closed canopy
mature forest. Numerous studies have found higher forage lichen abundance associated with
more open stand conditions and thus recommend stand thinning (citing reduced basal area
and/or higher light levels) (Horstkotte and Moen 2019, Nobert et al. 2020). Thinning is also
72

recommended to maintain conditions that benefit lichen and reduce the environmental
conditions that favour bryophytes in maturing stands (Coxson 2015, Vitt et al. 2019,
Horstkotte and Moen 2019). However, thinning too heavily risks creating conditions suitable
for vascular plants that compete with lichen and produce litter that suppresses lichen
(Cornelissen et al. 2001, Coxson 2015, Cichowski and Haeussler 2020). Where harvest
occurs, methods should be employed that minimize mechanical disturbance of the lichen
(e.g., employ winter harvest or dedicated use of skid trails) and that limit the spread of debris
over lichen (Coxson 2015, Horstkotte and Moen 2019, Cichowski et al. 2022). A stand
thinning study in Alberta found that lichen cover remained either stable or slightly increased
and moss cover significantly declined 19 years following treatment (Vitt et al. 2019).
Although lichen abundance didn’t increase with thinning in that study area, it didn’t decrease
as has occurred following clearcutting and mountain pine beetle (Webb 1998, Cichowski and
Haeussler 2020, Cichowski et al. 2022). When considering partial harvest, the potential for
increasing the density and network of active roads and the associated mortality risk must be
taken into consideration (Dickie et al. 2017, Mumma et al. 2018).
Model results indicated that clearcutting should be avoided on mesic and high productivity
sites within caribou winter habitat. Simulated clearcutting resulted in heightened moose
shrub suitability on these sites in comparison to mature and MPB attacked stands. Although
clearcutting on submesic sites was not projected to increase moose shrub suitability to as
great an extent, it was projected to increase suitability above MPB attacked and mature stand
levels and remained greater for a longer duration than on more productive sites. Although
more productive sites are generally not suitable for terrestrial lichen, these results suggest that
predation risk for caribou increases where subxeric lichen sites are located adjacent to sites
with greater moose shrub suitability.

73

Numerous studies have found that clearcuts enhance moose forage and that moose select
cutblocks over other habitat types (Carleton and MacLellan 1994, Fisher and Wilkinson
2005, Mumma et al. 2021). Clearcuts on mesic and higher productivity sites were projected
to provide greater suitability for moose shrubs. This is supported by the work of Chen et al.
(2004) where understory woody species was positively associated with site index.
Understory shrub species on most productive sites included preferred moose shrubs: Cornus
stolonifera and Viburnum edule (Chen et al. 2004). Submesic clearcuts were projected to
sustain moose shrub suitability above that found in MPB stands for longer than on richer sites
due to the lower productivity of the sites. Growth of conifer trees is slower on these sites
therefore more time elapses before crown closure limits light in the understory shrub layer.
On poorer productivity sites (subxeric and submesic), clearcuts were projected to provide
conditions suitable for terrestrial lichen within 25 years; however, these conditions weren’t
sustained beyond 50 years (Figure 10). This mimics the common pattern of succession
following wildfire where lichen communities initially establish but then diminish as tree
crown closure occurs (Maikawa and Kershaw 1976, Coxson and Marsh 2001, Lewis et al.
2019). Where harvest methods protect the lichen community (e.g., winter harvest with a deep
snow pack), lichen cover can be sustained or recover quickly (Coxson and Marsh 2001, Seip
and Jones 2010). Although total lichen cover may be relatively high following harvest or
fire, the Cladonia morphotypes often dominate the initial lichen community (Coxson and
Marsh 2001). Within 40 to 60 years the lichen community transitions to Cladina
morphotypes, providing greater forage biomass for caribou (Coxson and Marsh 2001).
Although lichen may be abundant, the characteristics of the snowpack in clearcuts can render
terrestrial lichen inaccessible or limit detectability (Johnson et al. 2001); without a mature
canopy to intercept snow, the snowpack is deeper and a hard crust may form (Kivinen et al.
2010). Caribou have been found to stop foraging on terrestrial lichen when the snow pack
74

exceeds a threshold height or hardness (Johnson et al. 2001, Seip and Jones 2010). Seip and
Jones (2010) found that post-MPB attack, the dead snags continued to provide an
ameliorative effect such that snowpack conditions did not limit cratering.
Model results suggest that where timber harvest is planned, broadleaf tree regeneration
should be minimized to maximize caribou winter habitat. Conifer/broadleaf stand types
support greater forage for moose in the form of both regenerating broadleaf saplings and
greater moose shrub communities (Eastman 1977, Rea 2014). In addition, stands with
significant broadleaf components are associated with greater moose forage (Eastman 1977,
Crête and Courtois 1997, Rea 2014). Also, lichens were projected to be least abundant under
conifer/broadleaf stands regardless of site or disturbance. Other studies have found a
decreased likelihood of lichen presence under these stand types (Russell and Johnson 2019,
Hillman and Nielsen 2020). Numerous studies provide evidence that deciduous litter, in
addition to shade from vascular plants, suppresses lichen communities (Canters et al. 1991,
Cornelissen et al. 2001). Studies of habitat selection have found that caribou avoid
mixedwood and broadleaf stands (Metsaranta et al. 2003, Metsaranta and Mallory 2007,
Fortin et al. 2008).
Interestingly, the mature subxeric conifer mix was projected to support greater lichen
suitability than the pine stand type. This was due to the lower initial basal area of the stand;
there was very little pine regeneration and reduced regeneration success and growth of
subalpine fir and spruce on this very poor productivity site. As these trees grow and natural
ingress occurs, tree cover will increase and lichen suitability will decline (Horstkotte and
Moen 2019). However, in the short term, there was greater light and space available to
support lichen than in a pine stand on a similar site where pine is best adapted to subxeric
conditions.

75

Limitations and Qualifications
I modelled three habitat elements considered significant for caribou winter habitat however,
other elements do play a role. Deadfall in stands limits caribou movement and hinders
caribou’s ability to crater for lichen (Metsaranta et al. 2003). Also, dead trees cover and
eliminate lichen (Cichowski and Haeussler 2020). Horizontal cover affects the ability of
caribou to see and escape from predators (Metsaranta et al. 2003). Dense young regenerating
stands weren’t used by caribou in southeast Manitoba despite available lichen; Schaeffer and
Pruitt (1991) suggest this was due to the density of the trees limiting visibility. As discussed
previously, the degree of vertical cover provided by a stand influences the snowpack
characteristics and thus accessibility and detectability of terrestrial lichen. Although these
additional factors influence the suitability of a stand for caribou winter habitat, I was not able
to develop linker functions for these factors or model them with SORTIE-ND.
Planning of forest management activities should consider the provision of caribou habitat at a
landscape scale. My work has focused on how stand level conditions influence key
determinants of caribou and moose habitat. These findings can be used to inform both stand
level management decisions, as well as landscape level planning. By combining the output of
my work into a spatially explicit representation of a real forest landscape mosaic the
landscape implications of forest disturbances and management decisions on caribou
population dynamics could be better projected (Holt et al. 1995). For example, in Labrador
SORTIE was linked to two other models, LANDIS-II (landscape scale disturbance and
succession model) and Patchworks (spatial modelling tool for harvest planning) using SELES
(a tool that links models) (Sturtevant et al. 2007). This interactive model system was
developed to support sustainable management of a landbase with multiple resource values.
SORTIE provided the succession trajectories and yield curves for complex stand types that

76

were needed to model alternative silviculture treatments throughout the landscape (Sturtevant
et al. 2007). Scaling up to the landscape scale with LANDIS-II enabled modelling the
interaction of alternative forest management options with natural disturbance across the
region (Sturtevant et al. 2007).
Consideration of the density and location of roads and the size, distribution, and adjacency of
stand types is important for caribou population management and requires a landscape scale
approach. Roads provide wolves with efficient pathways to travel and access to areas that
previously acted as refugia from predation. Both of those outcomes can increase predation or
risk-related displacement (Whittington et al. 2011, Dickie et al. 2017, DeMars and Boutin
2018, Mumma et al. 2018).
Minimizing fragmentation and maintaining large areas of older forest are important
landscape-scale management strategies for caribou (Racey et al. 1991). Longer rotations and
maintaining MPB attacked stands can sustain arboreal lichen that provide forage when snow
conditions limit availability of terrestrial lichen (Johnson et al. 2001, Seip and Jones 2010,
Kivinen et al. 2010). Minimizing fragmentation is important for maintaining a landscape
where caribou may forage more efficiently and also minimizes the co-occurrence of good
lichen foraging sites adjacent to stands with moose forage (Kivinen et al. 2010). Among the
first documented occurrences of caribou foraging in cutblocks, Racey et al. (1996) suggested
the caribou used the young cutblocks, in part due to the close proximity to an adjacent
extensive peatland complex that was undisturbed and had high early winter habitat value
(Racey et al. 1991). Minimizing the provision of moose habitat should also be considered at
the landscape scale to minimize the effects of apparent competition. The landscape primarily
used by TEC for winter habitat is characterized as low-lying, flat to gently rolling terrain;
subxeric pine stands with abundant forage lichen often occur in close proximity to sites with

77

high potential for moose shrubs when disturbed (mesic and subhygric stands) (Cichowski
1993, Cichowski and Haeussler 2020). Therefore, planned management activities on a
landscape should consider proximity to high-value caribou habitat, density and location of
roads, and the distribution of moose habitat following disturbance.
3.5 Conclusions
I applied a framework that uses linker functions to quantitatively relate habitat elements to
stand dynamics. That modelling framework allowed me to evaluate the post-disturbance
recovery dynamics of caribou winter habitat. I projected the effects of mountain pine beetle
attack and clearcut harvest on four different stand types across four different sites. My results
suggest that managers must consider site, stand type, and disturbance type and history when
planning forest harvest that will sustain winter range for caribou and reduce habitat for
moose. In general, disturbance improves conditions for either lichen or moose shrubs
depending on the site and disturbance type: MPB sustains greater terrestrial lichen suitability
than clearcuts on dry, poor sites and there is little improvement in moose shrub suitability.
Moose shrubs are enhanced on richer sites and following clearcut harvest. This framework
can be used to inform landscape-scale habitat modelling of the supply of caribou winter
habitat in areas subject to multiple disturbances of differing intensities.

78

Chapter Four
Conclusions
Forest ecosystems are dynamic, complex systems with many components interacting on
different scales (Messier et al. 2013). Disturbances influence forest structure and
composition, while also influencing the trajectory of succession. Forest ecosystems provide
habitats for many species. Species also interact with forest ecosystems at a diversity of
scales. Some organisms, such as lichen, interact with a stand at a small, local scale, others,
such as northern mountain caribou, migrate over one hundred kilometres across broad
landscapes and diverse forest ecosystems (Cichowski 1989, Haughian and Burton 2015).
Each species relies on different components of the ecosystem. We are becoming increasingly
aware of the importance of sustaining the diversity of components and processes within
ecosystems to support resilient and robust systems (Messier et al. 2015).
Extensive forest management has far reaching impacts and affects a multitude of ecosystem
components and services provided by forested communities (Mönkkönen et al. 2014, Suzuki
and Parker 2016). Although forest management acknowledges the importance of the many
ecosystem services, it is challenging to manage explicitly for ecosystem services beyond
economic values (Schwenk et al. 2012). We have process-based models that represent the
dynamics of trees within forest ecosystems, however, models of other ecosystem services
such as wildlife habitat often operate at the landscape-scale and don’t reflect stand-scale
dynamics (Holt et al. 1995, Marzluff et al. 2002).
The objective of this study was to build upon an existing stand dynamics model to enable
modelling habitat elements at the stand scale. I developed and applied a modelling
framework to explore the habitat dynamics of a species of conservation concern, northern

79

mountain caribou. This chapter presents a synthesis of my findings, the contributions of this
study, and future considerations.
4.1 Main Findings and Management Implications
Linking habitat elements to a stand dynamics model has allowed me to explore the
implications of disturbance, stand structure, and edaphic site on habitat. Model results
indicate that stand type and edaphic site must be taken into consideration when considering
habitat suitability, whether in response to disturbance events or, in the absence of disturbance,
in the context of ongoing stand structural development. For northern mountain caribou, these
results can inform management decisions to support the provision of forested winter habitat
over mid- to long-term timeframes. The results highlight sites and stand types where certain
management actions would best support the provision of winter habitat.
Conifer stands on the poorest productivity sites provided the highest quality winter habitat for
caribou of all the mature stands and sites. These stand types were the most suitable for
terrestrial forage lichen and the least suitable for elements of moose habitat (moose shrubs
and vertical cover). Maintaining these stands as mature stands should be prioritized to
maintain existing caribou winter forage and limit predation risk. Pine stands on subxeric sites
are projected to increase in suitability to support lichen as they mature (these stands were
initiated at approximately 80 years of age). This is especially important for herds such as the
TEC with heavily disturbed winter habitats. This is reflected in most planning for caribou
habitat and recent calls for protection of critical caribou habitat (McNay et al. 2008, FestaBianchet et al. 2011, Environment Canada 2014, Johnson et al. 2015).
Disturbance on the poorest productivity pine stands was projected to generate the most
suitable lichen habitat for caribou (for 50–75 years post-disturbance). For these same sites,

80

moose shrub suitability remained low and vertical cover was slow to recover from
disturbance. These results suggested that leaving MPB attacked stands to recover naturally
would maintain greater forage opportunities for caribou than clearcutting and limit the
creation of higher value moose habitat. Recent studies evaluating the response of terrestrial
lichen to MPB attack also recommend reserving poor productivity MPB attacked stands from
harvest in order to support declining caribou populations (Cichowski and Haeussler 2020,
Nobert et al. 2020, Cichowski et al. 2022). Although these studies have documented an
initial decline in forage lichen following MPB attack, unharvested MPB-attacked stands still
support more lichen than clearcuts (on comparable sites) for up to 14 years post-MPB attack
(Cichowski et al. 2022). The decline post-MPB may also be short lived on the poorest
productivity sites, as lichen did not decline between 10 and 16 years post-MPB (Cichowski
and Haeussler 2020).
Stands with a substantial component of broadleaves are not beneficial for caribou winter
habitat. These stands are projected to have very low suitability for lichen regardless of site.
In addition, these stands negatively influence caribou habitat by supporting the growth of
shrubs and trees that serve as moose forage. Numerous habitat selection studies have found
that caribou select against mixed and broadleaf stand types (Cichowski 1989, Metsaranta and
Mallory 2007, Fortin et al. 2008). Therefore, when managing for high-value caribou winter
habitat, minimizing the occurrence of broadleaf trees on the landscape is desirable. This can
be accomplished either through promoting and maintaining mature and old stands (Maleki et
al. 2020) or minimizing the regeneration of broadleaf trees in regenerating stands.
Clearcut harvest of mesic and subhygric stand types are also projected to negatively affect
caribou winter habitat. Clearcutting these sites can result in an increase of moose shrubs for
up to 70 years. Although clearcutting reduces the availability of vertical cover for moose,

81

cover is reduced for a limited period of time in most stands. Caribou have been found to
select against young stand types and harvested cutblocks (Cichowski 1989, Rettie and
Messier 2000, Courtois et al. 2008). Thus, my research supports other work and management
strategies that suggest limiting clearcut harvest especially in areas that promote habitat for
apparent competitors of caribou, including moose and deer (Serrouya et al. 2011, Johnson et
al. 2015, St-Laurent et al. 2022, Leblond et al. 2022).
4.2 Contributions
This thesis was initiated in part by my desire to explore the ability to use an existing stand
dynamics model to evaluate the effects of stand management decisions on habitat elements
beyond the tree components of forest ecosystems. The framework that has resulted
contributes to caribou habitat management by providing resource managers with a means to
evaluate the effects of management decisions on caribou winter habitat and could be used to
support habitat supply modelling.
This work developed a new framework for modelling the influence of a variety of forest
disturbances and silviculture approaches on ecosystem services. Developing the framework
involved identifying key habitat elements, creating linker functions, and then collecting
empirical data that were used to fit and calibrate linker functions. This framework allows one
to evaluate the provision of ecosystem services following a variety of disturbances and
alternative silviculture systems. Although I only evaluated the results of three treatments (no
treatment, clearcutting and severe MPB attack), SORTIE-ND is very flexible and a variety of
silviculture treatments and disturbances of varying intensities and spatial patterns can be
modelled. For example, the effects of spruce budworm, partial harvest with gap cuts, partial
harvest with individual tree removal, thinning treatments, and shelterwoods can be modelled
(Ameztegui et al. 2017, Cristal et al. 2019, Maleki et al. 2020, 2021). This is an important
82

step towards evaluating the effects of forest management decisions on not only the tree
component of the ecosystem but also other ecosystem services. This framework may be most
useful for ecosystem elements that require longer time frames to develop such as terrestrial
forage lichens (Holt et al. 1995).
This approach could also be scaled up to landscape-scale models of ecosystem services such
as habitat supply. This could be highly valuable on a landscape heavily affected by a variety
of disturbances such the TEC winter range. Enabling the dynamic (and site- and standspecific) projection of the recovery of winter habitat elements allows the evaluation of the
effects of alternate stand management decisions.
4.3 Next Steps and Future Considerations
Ecological models are powerful, useful, and provide a framework to test our understanding of
ecological systems and explore the implications of perturbations to those systems (Filotas et
al. 2014). Implicit in testing our understanding, is the recognition that we will discover areas
that require further study. In addition, as models are applied and ecological problems are
explored, additional potential uses become evident. Here I explore future considerations for
extending the application of this framework to the TEC winter range as well as identify areas
that would benefit from further study.
A valuable future step would be developing an optimization framework for evaluating
alternative management strategies to address the sustained provision of caribou winter habitat
for the TEC. Linking dynamic habitat modelling to caribou habitat supply models would be
the first step. Optimization could be used to determine the type and distribution of treatments
that ensured continual provision of winter habitat. In conjunction with evaluating habitat
supply, the provision of additional ecosystem services such as timber value and carbon

83

sequestration could be simultaneously evaluated (Mönkkönen et al. 2014, Triviño et al. 2017,
Lafond et al. 2017). Optimization requires assigning weights to each ecosystem service or
habitat element; this can be challenging as weights significantly influence the selection of
optimal outcomes (Schwenk et al. 2012, Lafond et al. 2017). Results of optimization studies
indicate that the application of a diversity of silvicultural approaches can supply a broad
range of ecosystem services and resolve conflicts between competing objectives (Schwenk et
al. 2012, Mönkkönen et al. 2014, Triviño et al. 2017).
Another logical extension of this project is to evaluate the effect of partial harvesting on
caribou habitat. SORTIE-ND is ideally suited to exploring partial harvest scenarios as it was
designed to model the complex dynamics of stands with diverse mixtures of species and tree
size classes. Partial harvest strategies have been suggested as a means to maintain lichen in
mature stands where displacement of lichens by moss is a concern (Coxson 2015, Stevenson
and Coxson 2015, Haughian and Burton 2015, Horstkotte and Moen 2019). Partial harvest
strategies could also reduce the emergence of moose shrubs that follow clearcut harvest
(Eastman 1974, Carleton and MacLellan 1994, Mumma et al. 2021). The results for MPB
attacked pine stands suggest that it would be informative to explore partial cutting scenarios
that mimic stand structure following MPB attack such as uniform partial cutting (e.g.,
enhanced forage lichen suitability and limited moose shrub response compared to clearcut
harvest and mature forests).
Validation of the linker functions would add value to this study and lend support to model
projections. Although I did not validate the linker functions in this thesis, the results from
Chapter 3 are consistent with theoretical expectations and thus provide confidence that the
linker functions adequately project lichen, moose, and cover suitability. Validation by
withholding a portion of the dataset was not possible due to the small sample size.

84

It would be desirable to expand the dataset to provide greater insight into several ecological
relationships. The existing empirical data illustrate that humus depth was a strong predictor
variable for lichen suitability and I hypothesize that this could be related to disturbance
history. Further data collection across disturbance types and severities could inform both the
lichen and moose shrub linker functions. As discussed in Chapter 2, modelling all moose
shrubs as one ecological unit is also problematic. Adding to the dataset to provide enough
data to differentiate between some groups of moose shrubs could enable refinement of the
moose shrub linker functions and improve projections of shrub suitability. Further work to
investigate the interactions between terrestrial lichen, moss, and low growing vegetation
(especially ericaceous shrubs) is important. Numerous studies have found lichen abundance
correlated with changes in low vegetation or moss abundance (Canters et al. 1991, Pharo and
Vitt 2000, Horstkotte and Moen 2019, Cichowski and Haeussler 2020). Developing a linker
function that reflects the relationship between these groups of vegetation may improve
modelling of low-intensity disturbances. Insights from indigenous knowledge were not
included in this study but could enrich and inform our understanding and modelling of stand
and habitat dynamics. This study, and the TEC herd range, falls within the traditional
territories of several nations who could provide valuable perspective and knowledge of
habitat and wildlife dynamics.
Climate change is an important factor when considering mid- to long-term modelling (Barber
et al. 2018, Leblond et al. 2022). Further work is necessary to determine how best to
consider the potential effects of climate change for these elements of winter habitat. The
scale at which climate change is best modeled will need to be explored (e.g., the interaction
between habitat elements and site, at the stand scale, or via disturbance dynamics at the
landscape scale). The effects of climate change can be implemented in SORTIE-ND via tree
growth and mortality parameters (Morán-Ordóñez et al. 2020, Moran et al. 2021).
85

Several recent studies have modelled the influence of climate change on the distribution of
caribou habitat (Barber et al. 2018, Leblond et al. 2022, Neilson et al. 2022). Barber et al.
(2018) reported significant shifts in upland habitat types as a result of the influence of climate
change on dominant vegetation cover and wildfire. The other two studies found that habitat
change via harvesting or wildfire played a greater role in affecting habitat suitability than
climate induced changes to caribou habitat suitability itself (Leblond et al. 2022, Neilson et
al. 2022). Although we need to be mindful that the drivers of ecosystem dynamics will shift
as climate changes, the learnings in this study are important for projecting how systems will
change in the future and what that means for caribou habitat.
4.4 Concluding Remarks
Here I have demonstrated that process-based forest models can be used to describe how forest
structure changes and how these changes can then be extrapolated to other important
ecosystem components such as habitat. Modelling the effects of disturbance on caribou
winter habitat revealed that stand type, edaphic site, and disturbance type all influenced the
suitability of forested habitat. This modelling approach can be used to evaluate silviculture
decisions and habitat management options. Results can inform decisions regarding
management of critical winter habitat for caribou at the stand scale and inform landscapescale habitat supply modelling. This approach, though focused on caribou winter, can also be
applied to questions of moose habitat where managing for moose is a priority.

86

Literature Cited
Ahti, T., and J. Oksanen. 1990. Epigeic lichen communities of taiga and tundra regions.
Vegetatio 86:39–70.
Allen, A. W., P. A. Jordan, and J. W. Terrell. 1987. Habitat suitability index models: moose,
Lake Superior region. Biol. Rep., U.S. Fish Wildl. Serv. 47 pp.
Ameztegui, A., A. Cabon, M. De Cáceres, and L. Coll. 2017. Managing stand density to
enhance the adaptability of Scots pine stands to climate change: A modelling
approach. Ecological Modelling 356:141–150.
Anderson, M., B. N. McLellan, and R. Serrouya. 2018. Moose response to high-elevation
forestry: Implications for apparent competition with endangered caribou. The Journal
of Wildlife Management 82:299–309.
Angelstam, P., and T. Kuuluvainen. 2004. Boreal forest disturbance regimes, successional
dynamics and landscape structures - a European perspective. Ecological
Bulletins:117–136.
Astrup, R., K. D. Coates, and E. Hall. 2008a. Finding the appropriate level of complexity for
a simulation model: an example with a forest growth model. Forest Ecology and
Management 256:1659–1665.
Astrup, R., K. D. Coates, and E. Hall. 2008b. Recruitment limitation in forests: lessons from
an unprecedented mountain pine beetle epidemic. Forest Ecology and Management
256:1743–1750.
Baker, B. G. 1990. Winter habitat selection and use by moose in the West-Chilcotin region of
British Columbia. M.Sc. Thesis, University of British Columbia, Vancouver, B.C.
Banner, A., W. MacKenzie, S. Haeussler, S. Thomson, J. Pojar, and R. Trowbridge. 1993. A
field guide to site identification and interpretation for the Prince Rupert Forest
Region. Province of British Columbia, Victoria, B.C.
87

Barber, Q. E., M.-A. Parisien, E. Whitman, D. Stralberg, C. J. Johnson, M.-H. St‐Laurent, E.
R. DeLancey, D. T. Price, D. Arseneault, X. Wang, and M. D. Flannigan. 2018.
Potential impacts of climate change on the habitat of boreal woodland caribou.
Ecosphere 9:e02472.
Bartels, S. F., H. Y. H. Chen, M. A. Wulder, and J. C. White. 2016. Trends in postdisturbance recovery rates of Canada’s forests following wildfire and harvest. Forest
Ecology and Management 361:194–207.
Bartemucci, P., C. Messier, and C. D. Canham. 2006. Overstory influences on light
attenuation patterns and understory plant community diversity and composition in
southern boreal forests of Quebec. Canadian Journal of Forest Research 36:2065–
2079.
Beaudry, L., R. Coupé, C. Delong, and J. Pojar. 1999. Plant indicator guide for northern
British Columbia: Boreal, sub-boreal, and subalpine biogeoclimatic zones. Forestry
Division Services Branch. Province of B.C., Victoria, B.C.
Bergeron, J. A. C., J. Pinzon, S. Odsen, S. Bartels, S. E. Macdonald, and J. R. Spence. 2017.
Ecosystem memory of wildfires affects resilience of boreal mixedwood biodiversity
after retention harvest. Oikos 126:1738–1747.
Bergeron, Y. 2000. Species and Stand Dynamics in the Mixed Woods of Quebec’s Southern
Boreal Forest. Ecology 81:1500–1516.
Bergerud, A. T. 1971. Abundance of forage on the winter range of Newfoundland caribou.
The Canadian Field-Naturalist 85:39–52.
Bjørneraas, K., E. J. Solberg, I. Herfindal, B. V. Moorter, C. M. Rolandsen, J.-P. Tremblay,
C. Skarpe, B.-E. Sæther, R. Eriksen, and R. Astrup. 2011. Moose Alces alces habitat
use at multiple temporal scales in a human-altered landscape. Wildlife Biology 17:44–
54.

88

Blattert, C., R. Lemm, O. Thees, M. J. Lexer, and M. Hanewinkel. 2017. Management of
ecosystem services in mountain forests: Review of indicators and value functions for
model based multi-criteria decision analysis. Ecological Indicators 79:391–409.
Bose, A. K., B. D. Harvey, K. D. Coates, S. Brais, and Y. Bergeron. 2015. Modelling stand
development after partial harvesting in boreal mixedwoods of eastern Canada.
Ecological Modelling 300:123–136.
Botkin, D. B., J. F. Janak, and J. R. Wallis. 1972a. Some ecological consequences of a
computer model of forest growth. Journal of Ecology 60:849–872.
Botkin, D. B., J. F. Janak, and J. R. Wallis. 1972b. Rationale, limitations, and assumptions of
a northeastern forest growth simulator. IBM Journal of Research and Development
16:101–116.
Boudreault, C., Y. Bergeron, and D. Coxson. 2009. Factors controlling epiphytic lichen
biomass during postfire succession in black spruce boreal forests. Canadian Journal of
Forest Research 39:2168–2179.
Bowd, E., W. Blanchard, L. McBurney, and D. Lindenmayer. 2021. Direct and indirect
disturbance impacts on forest biodiversity. Ecosphere 12:e03823.
British Columbia Data Catalogue. n.d. VRI - 2022 - Forest Vegetation Composite Polygons.
https://catalogue.data.gov.bc.ca/dataset/vri-2022-forest-vegetation-compositepolygons/resource/85724845-3973-4a26-b715-100306768bbc.
Brodo, I. M., S. D. Sharnoff, and S. Sharnoff. 2001. Lichens of North America. Yale
University Press, New Haven and London.
Brulisauer, A. R., G. E. Bradfield, and J. Maze. 1996. Quantifying organizational change
after fire in lodgepole pine forest understorey. Canadian Journal of Botany 74:1773–
1782.
Bugmann, H. 2001. A review of forest gap models. Climatic Change 51:259–305.

89

Bugmann, H., T. Cordonnier, H. Truhetz, and M. J. Lexer. 2017. Impacts of business-asusual management on ecosystem services in European mountain ranges under climate
change. Regional Environmental Change 17:3–16.
Bugmann, H., and R. Seidl. 2022. The evolution, complexity and diversity of models of longterm forest dynamics. Journal of Ecology 110:2288–2307.
Buma, B. 2015. Disturbance interactions: characterization, prediction, and the potential for
cascading effects. Ecosphere 6:art70.
Burnham, K. P., and D. R. Anderson. 2002. Model selection and multimodel inference: a
practical information-theoretic approach. 2nd ed. Springer, New York.
Canham, C. D. 1988. An index for understory light levels in and around canopy gaps.
Ecology 69:1634–1638.
Canham, C. D., K. D. Coates, P. Bartemucci, and S. Quaglia. 1999. Measurement and
modeling of spatially explicit variation in light transmission through interior cedar–
hemlock forests of British Columbia 29:1775–1783.
Canham, C. D., P. T. LePage, and K. D. Coates. 2004. A neighborhood analysis of canopy
tree competition: effects of shading versus crowding. Canadian Journal of Forest
Research 34:778–787.
Canham, C. D., and L. E. Murphy. (n.d.). SORTIE-ND Documentation. http://www.sortiend.org/help/manuals/index.html.
Canham, C. D., and M. Uriarte. 2006. Analysis of neighborhood dynamics of forest
ecosystems using likelihood methods and modeling. Ecological Applications 16:62–
73.
Canters, K. J., H. Schöllerj, S. Ott, and H. M. Jahns. 1991. Microclimatic influences on lichen
distribution and community development. The Lichenologist 23:237–252.

90

Carleton, T. J., and P. MacLellan. 1994. Woody vegetation responses to fire versus clearcutting logging: a comparative survey in the central Canadian boreal forest.
Écoscience 1:141–152.
Carroll, S. B., and L. C. Bliss. 1982. Jack pine – lichen woodland on sandy soils in northern
Saskatchewan and northeastern Alberta. Canadian Journal of Botany 60:2270–2282.
Chen, H. Y. H., S. Légaré, and Y. Bergeron. 2004. Variation of the understory composition
and diversity along a gradient of productivity in Populus tremuloides stands of
northern British Columbia, Canada. Canadian Journal of Botany 82:1314–1323.
Chen, J., J. F. Franklin, and T. A. Spies. 1995. Growing-season microclimatic gradients from
clearcut edges into old-growth Douglas-fir forests. Ecological Applications 5:74–86.
Cichowski, D. 2010. Tweedsmuir-Entiako caribou project: Effects of a mountain pine beetle
epidemic on northern caribou habitat use: final report. 77 pp.
Cichowski, D. 2015. Tweedsmuir-Entiako caribou population status and background
information summary. BC Ministry of Forests, Lands and Natural Resource
Operations, Smithers, B.C. 145 pp.
Cichowski, D. B. 1989. Seasonal movements, habitat use, and winter feeding ecology of
woodland caribou in West-Central British Columbia. M.Sc. Thesis, University of
British Columbia, Vancouver, B.C.
Cichowski, D. B. 1993. Seasonal movements, habitat use, and winter feeding ecology of
woodland caribou in west-central British Columbia. Research Branch, Ministry of
Forests, Victoria, B.C. 65 pp.
Cichowski, D. B. 1996. Managing woodland caribou in west-central British Columbia.
Rangifer 16:119–126.
Cichowski, D., and S. Haeussler. 2013. The response of caribou terrestrial forage lichens to
mountain pine beetles and forest harvesting in the East Ootsa and Entiako areas.

91

Annual report – 2012/13 – year 11. Prepared for Bulkley Valley Centre for Natural
Resources Research and Management, Habitat Conservation Trust Foundation, B.C.
Min. FLNRO, Smithers, B.C. 55 pp.
Cichowski, D., and S. Haeussler. 2020. The response of caribou terrestrial forage lichens to
mountain pine beetles, forest harvesting and fire in the East Ootsa and Entiako areas:
2016/17 (year 16) field season report. Bulkley Valley Centre for Natural Resources
Research and Management, Habitat Conservation Trust Foundation, BC Ministry of
Forests, Lands, Natural Resource Operations and Rural Development, BC Parks, and
Canadian Forest Products Ltd., Smithers, B.C. 97 pp.
Cichowski, D., G. D. Sutherland, R. S. McNay, and R. Sulyma. 2022. Direct and indirect
effects of habitat disturbances on caribou terrestrial forage lichens in montane forests
of British Columbia. Forests 13:251.
Clark, D. F., J. A. Antos, and G. E. Bradfield. 2003. Succession in sub-boreal forests of westcentral British Columbia. Journal of Vegetation Science 14:721–732.
Coates, K. D., C. D. Canham, M. Beaudet, D. L. Sachs, and C. Messier. 2003. Use of a
spatially explicit individual-tree model (SORTIE/BC) to explore the implications of
patchiness in structurally complex forests. Forest Ecology and Management 186:297–
310.
Coates, K. D., C. D. Canham, and P. T. LePage. 2009. Above- versus below-ground
competitive effects and responses of a guild of temperate tree species. Journal of
Ecology 97:118–130.
Coates, K. D., S. Haeussler, S. Lindeburgh, R. Pojar, and A. J. Stock. 1994. Ecology and
silviculture of interior spruce in British Columbia. FRDA Report, Canada/British
Columbia Partnership Agreement on Forest Resource Development, Victoria, B.C.
190 pp.

92

Coates, K. D., and E. C. Hall. 2005. Implications of alternate silvicultural strategies in
mountain pine beetle damaged stands: Technical report for forest science program
project y051161. Bulkley Valley Centre for Natural Resources Research and
Management, Smithers, B.C. 31 pp.
Coates, K. D., E. C. Hall, and C. D. Canham. 2018. Susceptibility of trees to windthrow
storm damage in partially harvested complex-structured multi-species forests. Forests
9:199.
Coates, K. D., E. B. Lilles, and R. Astrup. 2013. Competitive interactions across a soil
fertility gradient in a multispecies forest. Journal of Ecology 101:806–818.
Cordonnier, T., F. Berger, C. Elkin, T. Lämås, and M. Martinez. 2013. Models and linker
functions (indicators) for ecosystem services. ARANGE Deliverable Report,
ARANGE Advanced Multifunctional Management of European Mountain Forests
Project. 94 pp.
Cordonnier, T., F. Berger, C. Elkin, T. Lämås, and M. Martinez. 2014. ARANGE deliverable
D2.2: Models and linker functions (indicators) for ecosystem services. ARANGE
Deliverable Report, ARANGE Advanced Multifunctional Management of European
Mountain Forests Project.
Cornelissen, J. H. C., T. V. Callaghan, J. M. Alatalo, A. Michelsen, E. Graglia, A. E. Hartley,
D. S. Hik, S. E. Hobbie, M. C. Press, C. H. Robinson, G. H. R. Henry, G. R. Shaver,
G. K. Phoenix, D. Gwynn Jones, S. Jonasson, F. S. Chapin III, U. Molau, C. Neill, J.
A. Lee, J. M. Melillo, B. Sveinbjörnsson, and R. Aerts. 2001. Global change and
arctic ecosystems: is lichen decline a function of increases in vascular plant biomass?
Journal of Ecology 89:984–994.

93

COSEWIC. 2002. COSEWIC assessment and update status report on the Woodland caribou
Rangifer tarandus caribou in Canada. Committee on the Status of Endangered
Wildlife in Canada. xi + 98 pp.
COSEWIC. 2009. COSEWIC / COSEPAC - Guidelines for Recognizing Designatable Units.
https://www.cosewic.ca/index.php/en-ca/reports/preparing-status-reports/guidelinesrecognizing-designatable-units.html.
COSEWIC. 2011. Designatable units for caribou (Rangifer tarandus) in Canada. Committee
on the Status of Endangered Wildlife in Canada, Ottawa. 88 pp.
COSEWIC. 2014. COSEWIC assessment and status report on the caribou Rangifer tarandus:
Northern Mountain population, Central Mountain population, Southern Mountain
population in Canada. Committee on the Status of Endangered Wildlife in Canada,
Ottawa, Canada. 135 pp.
Courbin, N., D. Fortin, C. Dussault, and R. Courtois. 2009. Landscape management for
woodland caribou: the protection of forest blocks influences wolf-caribou cooccurrence. Landscape Ecology 24:1375–1388.
Courtois, R., A. Gingras, D. Fortin, A. Sebbane, B. Rochette, and L. Breton. 2008.
Demographic and behavioural response of woodland caribou to forest harvesting.
Canadian Journal of Forest Research 38:2837–2849.
Courtois, R., J.-P. Ouellet, L. Breton, A. Gingras, and C. Dussault. 2007. Effects of forest
disturbance on density, space use, and mortality of woodland caribou1. Écoscience;
Sainte-Foy 14:491–498.
Coxson, D. S. 2015. Using partial-cut harvesting to conserve terrestrial lichens in managed
landscapes. Canadian Wildlife Biology & Management 4:150–162.

94

Coxson, D. S., and J. Marsh. 2001. Lichen chronosequences (postfire and postharvest) in
lodgepole pine (Pinus contorta) forests of northern interior British Columbia.
Canadian Journal of Botany 79:1449–1464.
Crête, M., and R. Courtois. 1997. Limiting factors might obscure population regulation of
moose (Cervidae: Alces alces) in unproductive boreal forests. Journal of Zoology
242:765–781.
Cristal, I., A. Ameztegui, J. R. González-Olabarria, and J. Garcia-Gonzalo. 2019. A decision
support tool for assessing the impact of climate change on multiple ecosystem
services. Forests 10:440.
Crites, S., and M. R. Dale. 1998. Diversity and abundance of bryophytes, lichens, and fungi
in relation to woody substrate and successional stage in aspen mixedwood boreal
forests 76:14.
Cumming, H. G. 1987. Sixteen years of moose browse surveys in Ontario. Alces 23:125–155.
DeCesare, N. J., M. Hebblewhite, F. Schmiegelow, D. Hervieux, G. J. McDermid, L.
Neufeld, M. Bradley, J. Whittington, K. G. Smith, L. E. Morgantini, M. Wheatley,
and M. Musiani. 2012. Transcending scale dependence in identifying habitat with
resource selection functions. Ecological Applications 22:1068–1083.
DeMars, C. A., and S. Boutin. 2018. Nowhere to hide: Effects of linear features on predatorprey dynamics in a large mammal system. Journal of Animal Ecology 87:274–284.
DeMars, C., and R. Serrouya. 2020. Tweedsmuir-Entiako primary prey strategy: phase 1.
Caribou Monitoring Unit, Alberta Biodiversity Monitoring Institute. 41 pp.
Denryter, K. A., R. C. Cook, J. G. Cook, and K. L. Parker. 2017. Straight from the caribou’s
(Rangifer tarandus) mouth: detailed observations of tame caribou reveal new insights
into summer–autumn diets. Canadian Journal of Zoology 95:81–94.

95

Dhar, A., N. Balliet, K. Runzer, and C. Hawkins. 2015. Impact of a mountain pine beetle
outbreak on young lodgepole pine stands in central British Columbia. Forests 6:3483–
3500.
Dhar, A., L. Parrott, C. Hawkins, A. Dhar, L. Parrott, and C. D. B. Hawkins. 2016. Aftermath
of mountain pine beetle outbreak in British Columbia: Stand dynamics, management
response and ecosystem resilience. Forests 7:171.
Dickie, M., R. Serrouya, R. S. McNay, and S. Boutin. 2017. Faster and farther: wolf
movement on linear features and implications for hunting behaviour. Journal of
Applied Ecology 54:253–263.
Eastman, D. S. 1974. Habitat use by moose of burns, cutovers and forests in north-central
British Columbia. Tenth North American Moose Conference and Workshop. Duluth,
Minnesota. 19 pp.
Eastman, D. S. 1977. Habitat selection and use in winter by moose in sub-boreal forests of
north-central British Columbia, and relationships to forestry. PhD Thesis, University
of British Columbia, Vancouver, B.C.
Edwards, A. W. F. 1992. Likelihood. Page Likelihood—expanded edition. Johns Hopkins
University Press, Baltimore, Maryland, USA.
Elkin, C., A. G. Gutiérrez, S. Leuzinger, C. Manusch, C. Temperli, L. Rasche, and H.
Bugmann. 2013. A 2 °C warmer world is not safe for ecosystem services in the
European Alps. Global Change Biology 19:1827–1840.
Environment Canada. 2014. Recovery strategy for the woodland caribou, Southern Mountain
population (Rangifer tarandus caribou) in Canada. Species at Risk Act Recovery
Strategy Series., Environment Canada, Ottawa, Ontario. viii + 103 pp.

96

Festa-Bianchet, M., J. C. Ray, S. Boutin, S. D. Côté, and A. Gunn. 2011. Conservation of
caribou (Rangifer tarandus) in Canada: an uncertain future. Canadian Journal of
Zoology 89:419–434.
Filotas, E., L. Parrott, P. J. Burton, R. L. Chazdon, K. D. Coates, L. Coll, S. Haeussler, K.
Martin, S. Nocentini, K. J. Puettmann, F. E. Putz, S. W. Simard, and C. Messier.
2014. Viewing forests through the lens of complex systems science. Ecosphere 5:art1.
Fisher, J. T., and L. Wilkinson. 2005. The response of mammals to forest fire and timber
harvest in the North American boreal forest. Mammal Review 35:51–81.
Fortin, D., R. Courtois, P. Etcheverry, C. Dussault, and A. Gingras. 2008. Winter selection of
landscapes by woodland caribou: behavioural response to geographical gradients in
habitat attributes. Journal of Applied Ecology 45:1392–1400.
Frazer, G. W., C. D. Canham, and K. P. 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.
Goffe, W. L., G. D. Ferrier, and J. Rogers. 1994. Global optimization of statistical functions
with simulated annealing. Journal of Econometrics 60:65–99.
Goudie, J. W. 1984. Height growth and site index curves for lodgepole pine and white spruce
and interim managed stand yield tables for lodgepole pine in British Columbia.
Unpubl. Rep., B.C. Min. For., Res. Branch., F.P.D.S. Section. 75 pp.
Goulet, L. A. 1985. Winter habitat selection by moose in northern British Columbia. Alces
21:103–25.
Greuel, R. J., G. É. Degré-Timmons, J. L. Baltzer, J. F. Johnstone, E. J. B. McIntire, N. J.
Day, S. J. Hart, P. D. McLoughlin, F. K. A. Schmiegelow, M. R. Turetsky, A.

97

Truchon-Savard, M. D. van Telgen, and S. G. Cumming. 2021. Predicting patterns of
terrestrial lichen biomass recovery following boreal wildfires. Ecosphere 12:e03481.
Grimm, V., U. Berger, F. Bastiansen, S. Eliassen, V. Ginot, J. Giske, J. Goss-Custard, T.
Grand, S. K. Heinz, G. Huse, A. Huth, J. U. Jepsen, C. Jørgensen, W. M. Mooij, B.
Müller, G. Pe’er, C. Piou, S. F. Railsback, A. M. Robbins, M. M. Robbins, E.
Rossmanith, N. Rüger, E. Strand, S. Souissi, R. A. Stillman, R. Vabø, U. Visser, and
D. L. DeAngelis. 2006. A standard protocol for describing individual-based and
agent-based models. Ecological Modelling 198:115–126.
Grimm, V., U. Berger, D. L. DeAngelis, J. G. Polhill, J. Giske, and S. F. Railsback. 2010.
The ODD protocol: A review and first update. Ecological Modelling 221:2760–2768.
Habitat Monitoring Committee. 1996. Procedures for environmental monitoring in range and
wildlife habitat management Version 5.0. British Columbia Ministry of Environment,
Lands and Parks and British Columbia Ministry of Forests, Victoria, B.C. 225 pp.
Hall, L. S., P. R. Krausman, and M. L. Morrison. 1997. The habitat concept and a plea for
standard terminology. Wildlife Society Bulletin (1973-2006) 25:173–182.
Halpern, C. B., J. Halaj, S. A. Evans, and M. Dovčiak. 2012. Level and pattern of overstory
retention interact to shape long-term responses of understories to timber harvest.
Ecological Applications 22:2049–2064.
Haughian, S. R., and P. J. Burton. 2015. Microhabitat associations of lichens, feathermosses,
and vascular plants in a caribou winter range, and their implications for understory
development. Botany 93:221–231.
Hijmans, R. J. 2021. Raster: geographic data analysis and modeling.
Hillman, A. C., and S. E. Nielsen. 2020. Quantification of lichen cover and biomass using
field data, airborne laser scanning and high spatial resolution optical data: A case
study from a Canadian boreal pine forest. Forests 11:682.

98

Hodder, D. P., R. V. Rea, and S. M. Crowley. 2013. Diet content and overlap of sympatric
mule deer, moose, and elk during a deep snow winter in north-central British
Columbia, Canada. Canadian Wildlife Biology & Management 2:43–50.
Holt, R. D. 1977. Predation, apparent competition, and the structure of prey communities.
Theoretical Population Biology 12:197–229.
Holt, R. D., S. W. Pacala, T. W. Smith, and J. Liu. 1995. Linking Contemporary Vegetation
Models with Spatially Explicit Animal Population Models. Ecological Applications
5:20–27.
Horstkotte, T., and J. Moen. 2019. Successional pathways of terrestrial lichens in changing
Swedish boreal forests. Forest Ecology and Management 453:117572.
Johnson, C. J., N. D. Alexander, R. D. Wheate, and K. L. Parker. 2003. Characterizing
woodland caribou habitat in sub-boreal and boreal forests. Forest Ecology and
Management 180:241–248.
Johnson, C. J., L. P. W. Ehlers, and D. R. Seip. 2015. Witnessing extinction – Cumulative
impacts across landscapes and the future loss of an evolutionarily significant unit of
woodland caribou in Canada. Biological Conservation 186:176–186.
Johnson, C. J., K. L. Parker, and D. C. Heard. 2000. Feeding site selection by woodland
caribou in north-central British Columbia. Rangifer:158–172.
Johnson, C. J., K. L. Parker, and D. C. Heard. 2001. Foraging across a variable landscape:
behavioral decisions made by woodland caribou at multiple spatial scales. Oecologia
127:590–602.
Johnson, C. J., K. L. Parker, D. C. Heard, and M. P. Gillingham. 2002. A multiscale
behavioral approach to understanding the movements of woodland caribou.
Ecological Applications 12:1840–1860.

99

Johnson, C. J., D. R. Seip, and M. S. Boyce. 2004. A quantitative approach to conservation
planning: using resource selection functions to map the distribution of mountain
caribou at multiple spatial scales. Journal of Applied Ecology 41:238–251.
Johnson, E. A. 1981. Vegetation organization and dynamics of lichen woodland communities
in the Northwest Territories, Canada. Ecology 62:200–215.
Joly, K., B. W. Dale, W. B. Collins, and L. G. Adams. 2003. Winter habitat use by female
caribou in relation to wildland fires in interior Alaska. Canadian Journal of Zoology
81:1192.
Julianus, E., T. N. Hollingsworth, A. D. McGuire, and K. Kielland. 2019. Availability and
use of moose browse in response to post-fire succession on Kanuti National Wildlife
Refuge, Alaska. Alces 55:24.
Keim, J. L., P. D. DeWitt, J. J. Fitzpatrick, and N. S. Jenni. 2017. Estimating plant abundance
using inflated beta distributions: Applied learnings from a lichen-caribou ecosystem.
Ecology and Evolution 7:486–493.
Kielland, K., and C. Brown. 2015. Understanding the effects of wildfire severity on moose
habitat characteristics and use in Interior, Alaska. University of Alaska Fairbanks,
Fairbanks, Alaska.
Kivinen, S., J. Moen, A. Berg, and Å. Eriksson. 2010. Effects of modern forest management
on winter grazing resources for reindeer in Sweden. AMBIO 39:269–278.
Kobe, R. K., and K. D. Coates. 1997. Models of sapling mortality as a function of growth to
characterize interspecific variation in shade tolerance of eight tree species of
northwestern British Columbia 27:10.
Kranabetter, J. M., C. R. Dawson, and D. E. Dunn. 2007. Indices of dissolved organic
nitrogen, ammonium and nitrate across productivity gradients of boreal forests. Soil
Biology and Biochemistry 39:3147–3158.

100

Kranabetter, J. M., and S. W. Simard. 2008. Inverse relationship between understory light
and foliar nitrogen along productivity gradients of boreal forests. Canadian Journal of
Forest Research 38:2487–2496.
Lafond, V., T. Cordonnier, and B. Courbaud. 2015. Reconciling biodiversity conservation
and timber production in mixed uneven-aged mountain forests: Identification of
ecological intensification pathways. Environmental Management 56:1118–1133.
Lafond, V., T. Cordonnier, Z. Mao, and B. Courbaud. 2017. Trade-offs and synergies
between ecosystem services in uneven-aged mountain forests: evidences using Pareto
fronts. European Journal of Forest Research 136:997–1012.
Lance, A. N., and B. Mills. 1996. Attributes of woodland caribou migration habitat in westcentral British Columbia. Rangifer 16:355–364.
Leblond, M., Y. Boulanger, J. Pascual Puigdevall, and M.-H. St-Laurent. 2022. There is still
time to reconcile forest management with climate-driven declines in habitat suitability
for boreal caribou. Global Ecology and Conservation 39:e02294.
Leblond, M., C. Dussault, and M.-H. St-Laurent. 2014. Development and validation of an
expert-based habitat suitability model to support boreal caribou conservation.
Biological Conservation 177:100–108.
LePage, P. T., C. D. Canham, K. D. Coates, and P. Bartemucci. 2000. Seed abundance versus
substrate limitation of seedling recruitment in northern temperate forests of British
Columbia 30:415–427.
Lewis, K. J., C. J. Johnson, and M. N. Karim. 2019. Fire and lichen dynamics in the Taiga
Shield of the Northwest Territories and implications for barren-ground caribou winter
forage. Journal of Vegetation Science 0:13.

101

Lilles, E. B., and R. Astrup. 2012. Multiple resource limitation and ontogeny combined: a
growth rate comparison of three co-occurring conifers. Canadian Journal of Forest
Research.
Lilles, E., A. Dhar, K. D. Coates, and S. Haeussler. 2018. Retention level affects dynamics of
understory plant community recovery in northern temperate hemlock-cedar forests.
Forest Ecology and Management 421:3–15.
Longton, R. E. 1992. The role of bryophytes and lichens in terrestrial ecosystems. Pages 32–
76 in J. W. Bates and A. M. Farmer, editors. Bryophytes and lichens in a changing
environment. Oxford University Press.
Luttmerding, H. A., D. A. Demarchi, E. C. Lea, D. V. Meidinger, and T. Vold, editors. 1990.
Describing ecosystems in the field. 2nd ed. Ministry of Environment in cooperation
with Ministry of Forests, Victoria, B.C.
MacCracken, J. G., and L. A. Viereck. 1990. Browse regrowth and use by moose after fire in
interior Alaska. Northwest Science 64:11–18.
Maikawa, E., and K. A. Kershaw. 1976. Studies on lichen-dominated systems. XIX. The
postfire recovery sequence of black spruce – lichen woodland in the Abitau Lake
Region, N.W.T. Canadian Journal of Botany 54:2679–2687.
Mäkelä, A., M. del Río, J. Hynynen, M. J. Hawkins, C. Reyer, P. Soares, M. van Oijen, and
M. Tomé. 2012. Using stand-scale forest models for estimating indicators of
sustainable forest management. Forest Ecology and Management 285:164–178.
Maleki, K., M. A. Gueye, B. Lafleur, A. Leduc, and Y. Bergeron. 2020. Modelling postdisturbance successional dynamics of the Canadian boreal mixedwoods. Forests
11:28.

102

Maleki, K., B. Lafleur, A. Leduc, and Y. Bergeron. 2021. Modelling the influence of
different harvesting methods on forest dynamics in the boreal mixedwoods of western
Quebec, Canada. Forest Ecology and Management 479:118545.
Marzluff, J. M., J. J. Millspaugh, K. R. Ceder, C. D. Oliver, J. Withey, J. B. McCarter, C. L.
Mason, and J. Comnick. 2002. Modeling changes in wildlife habitat and timber
revenues in response to forest management. Forest Science 48:191–202.
Matlack, G. R. 1994. Vegetation dynamics of the forest edge -- trends in space and
successional time. Journal of Ecology 82:113–123.
McNay, R. S., D. Heard, R. Sulyma, and R. Ellis. 2008. A recovery action plan for northern
caribou herds in north-central British Columbia. FORREX Forest Research Extension
Partnership, Kamloops, B.C. 105 pp.
Meidinger, D. V., and J. Pojar, editors. 1991. Ecosystems of British Columbia. Research
Branch, Ministry of Forests, Victoria, B.C.
Messier, C. C., K. J. Puettmann, and K. D. Coates. 2013. Managing forests as complex
adaptive systems: Building resilience to the challenge of global change. Routledge,
New York.
Messier, C., K. Puettmann, R. Chazdon, K. P. Andersson, V. A. Angers, L. Brotons, E.
Filotas, R. Tittler, L. Parrott, and S. A. Levin. 2015. From management to
stewardship: Viewing forests as complex adaptive systems in an uncertain world.
Conservation Letters 8:368–377.
Metsaranta, J. M., and F. F. Mallory. 2007. Ecology and habitat selection of a woodland
caribou population in west-central Manitoba, Canada. Northeastern Naturalist
14:571–588.

103

Metsaranta, J. M., F. F. Mallory, and D. W. Cross. 2003. Vegetation characteristics of forest
stands used by woodland caribou and those disturbed by fire or logging in Manitoba.
Rangifer 23:255–266.
Min. of Env. and Climate Change. 2018. Caribou herd locations for BC.
https://catalogue.data.gov.bc.ca/dataset/2b217585-f48d-4d9f-b7ba-746909ac35ca.
Mönkkönen, M., A. Juutinen, A. Mazziotta, K. Miettinen, D. Podkopaev, P. Reunanen, H.
Salminen, and O.-P. Tikkanen. 2014. Spatially dynamic forest management to sustain
biodiversity and economic returns. Journal of Environmental Management 134:80–89.
Moran, E. V., N. Vannest, and M. Aubry-Kientz. 2021. Modeling the forest dynamics of the
Sierra Nevada under climate change using SORTIE-ND. Annals of Forest Science
78:75.
Morán-Ordóñez, A., A. Ameztegui, M. De Cáceres, S. de-Miguel, F. Lefèvre, L. Brotons, and
L. Coll. 2020. Future trade-offs and synergies among ecosystem services in
Mediterranean forests under global change scenarios. Ecosystem Services 45:101174.
Mumma, M. A., M. P. Gillingham, S. Marshall, C. Procter, A. R. Bevington, and M.
Scheideman. 2021. Regional moose (Alces alces) responses to forestry cutblocks are
driven by landscape-scale patterns of vegetation composition and regrowth. Forest
Ecology and Management 481:118763.
Mumma, M. A., M. P. Gillingham, K. L. Parker, C. J. Johnson, and M. Watters. 2018.
Predation risk for boreal woodland caribou in human-modified landscapes: Evidence
of wolf spatial responses independent of apparent competition. Biological
Conservation 228:215–223.
Murphy, L. 2015. Likelihood package R documentation. https://cran.rproject.org/web/packages/likelihood/likelihood.pdf

104

Mysterud, A., and E. Østbye. 1999. Cover as a habitat element for temperate ungulates:
effects on habitat selection and demography. Wildlife Society Bulletin (1973-2006)
27:11.
Nadeau Fortin, M.-A., L. Sirois, and M.-H. St-Laurent. 2016. Extensive forest management
contributes to maintain suitable habitat characteristics for the endangered AtlanticGaspésie caribou. Canadian Journal of Forest Research 46:933–942.
Nagy-Reis, M., M. Dickie, A. M. Calvert, M. Hebblewhite, D. Hervieux, D. R. Seip, S. L.
Gilbert, O. Venter, C. DeMars, S. Boutin, and R. Serrouya. 2021. Habitat loss
accelerates for the endangered woodland caribou in western Canada. Conservation
Science and Practice 3:e437.
Neilson, E. W., C. Castillo-Ayala, J. F. Beckers, C. A. Johnson, M. H. St-Laurent, N.
Mansuy, D. Price, A. Kelly, and M. A. Parisien. 2022. The direct and habitatmediated influence of climate on the biogeography of boreal caribou in Canada.
Climate Change Ecology 3:100052.
Nelson, C. R., and C. B. Halpern. 2005. Edge-related responses of understory plants to
aggregated retention harvest in the Pacific Northwest. Ecological Applications
15:196–209.
Nobert, B. R., T. A. Larsen, K. E. Pigeon, and L. Finnegan. 2020. Caribou in the cross-fire?
Considering terrestrial lichen forage in the face of mountain pine beetle
(Dendroctonus ponderosae) expansion. PLOS ONE 15:e0232248.
Oliver, C. D., and B. C. Larson. 1990. Forest stand dynamics. McGraw-Hill, Inc.
Osko, T. J., M. N. Hiltz, R. J. Hudson, and S. M. Wasel. 2004. Moose habitat preferences in
response to changing availability. The Journal of Wildlife Management 68:576–584.

105

Pacala, S. W., C. D. Canham, J. Saponara, J. A. Silander, R. K. Kobe, and E. Ribbens. 1996.
Forest models defined by field measurements: estimation, error analysis and
dynamics. Ecological Monographs 66:1–43.
Payette, S., and A. Delwaide. 2018. Tamm review: The North-American lichen woodland.
Forest Ecology and Management 417:167–183.
Pec, G. J., J. Karst, A. N. Sywenky, P. W. Cigan, N. Erbilgin, S. W. Simard, and J. F. C. Jr.
2015. Rapid increases in forest understory diversity and productivity following a
mountain pine beetle (Dendroctonus ponderosae) outbreak in pine forests. PLOS
ONE 10:e0124691.
Pecchi, M., M. Marchi, V. Burton, F. Giannetti, M. Moriondo, I. Bernetti, M. Bindi, and G.
Chirici. 2019. Species distribution modelling to support forest management. A
literature review. Ecological Modelling 411:108817.
Peek, J. M. 2007. Habitat Relationships. Pages 351–375 in A. W. Franzmann and C. C.
Schwartz, editors. Ecology and management of the North American moose. 2nd
edition. University Press of Colorado, Boulder, Colorado, USA.
Peek, J. M., D. L. Urich, and R. J. Mackie. 1976. Moose Habitat Selection and Relationships
to Forest Management in Northeastern Minnesota. Wildlife Monographs:3–65.
Pharo, E. J., and D. H. Vitt. 2000. Local Variation in Bryophyte and Macro-lichen Cover and
Diversity in Montane Forests of Western Canada. The Bryologist 103:455–466.
Province of British Columbia. 2010. Field manual for describing terrestrial ecosystems. 2nd
ed. B.C. Ministry of Forests and Range and B.C. Ministry of Environment, Victoria,
B.C.
Province of British Columbia. 2023. Wildfire Season Summary.
https://www2.gov.bc.ca/gov/content/safety/wildfire-status/about-bcws/wildfirehistory/wildfire-season-summary.

106

R Core Team. 2021. R: a language and environment for statistical computing. R Foundation
for Statistical Computing, Vienna, Austria.
Racey, G. D., K. Abraham, W. R. Darby, H. R. Timmermann, and Q. Day. 1991. Can
woodland caribou and the forest industry coexist: the Ontario scene. Rangifer 11:108–
115.
Racey, G. D., A. Harris, and R. F. Foster. 1996. Caribou winter habitat in the new forest:
lessons from Lucy Lake. Ontario Ministry of Natural Resources, Thunder Bay,
Ontario. 16 pp.
Rea, R. V. 2014. A preliminary assessment of moose (Alces alces) winter diets in the Aleza
Lake Research Forest in north-central British Columbia. Wildlife Afield 11:50–53.
Rea, R. V., and A. L. Booth. 2011. Use of trembling aspen bark by moose in a browseabundant habitat during winter. Wildlife Afield 8:104–107.
Rettie, W. J., and F. Messier. 2000. Hierarchical habitat selection by woodland caribou: its
relationship to limiting factors. Ecography 23:466–478.
Roloff, G. J., and B. J. Kernohan. 1999. Evaluating reliability of habitat suitability index
models. Wildlife Society Bulletin (1973-2006) 27:973–985.
Rowe, J. S. 1984. Lichen woodland in northern Canada. Pages 225–236 in R. Olson, R.
Hastings, and F. Geddes, editors. Northern ecology and resource management:
Memorial essays honouring Don Gill. The University of Alberta Press, Edmonton,
Alberta.
Russell, K. L. M., and C. J. Johnson. 2019. Post-fire dynamics of terrestrial lichens:
Implications for the recovery of woodland caribou winter range. Forest Ecology and
Management 434:1–17.

107

Safranyik, L. 2004. Mountain pine beetle epidemiology in lodgepole pine. Mountain pine
beetle symposium: challenges and solutions. Natural Resources Canada, Canadian
Forest Service, Pacific Forestry Centre, Victoria, B.C. 287 pp.
Sánchez, J. J., R. Marcos-Martinez, L. Srivastava, and N. Soonsawad. 2021. Valuing the
impacts of forest disturbances on ecosystem services: An examination of recreation
and climate regulation services in U.S. national forests. Trees, Forests and People
5:100123.
Schaefer, J. A., and W. O. Pruitt. 1991. Fire and woodland caribou in southeastern Manitoba.
Wildlife Monographs:3–39.
Schrempp, T. V., J. L. Rachlow, T. R. Johnson, L. A. Shipley, R. A. Long, J. L. Aycrigg, and
M. A. Hurley. 2019. Linking forest management to moose population trends: The role
of the nutritional landscape. PLOS ONE 14:e0219128.
Schwenk, W. S., T. M. Donovan, W. S. Keeton, and J. S. Nunery. 2012. Carbon storage,
timber production, and biodiversity: comparing ecosystem services with multi-criteria
decision analysis. Ecological Applications 22:1612–1627.
Seidl, R., P. M. Fernandes, T. F. Fonseca, F. Gillet, A. M. Jönsson, K. Merganičová, S.
Netherer, A. Arpaci, J.-D. Bontemps, H. Bugmann, J. R. González-Olabarria, P.
Lasch, C. Meredieu, F. Moreira, M.-J. Schelhaas, and F. Mohren. 2011. Modelling
natural disturbances in forest ecosystems: a review. Ecological Modelling 222:903–
924.
Seidl, R., T. A. Spies, D. L. Peterson, S. L. Stephens, and J. A. Hicke. 2016. REVIEW:
Searching for resilience: addressing the impacts of changing disturbance regimes on
forest ecosystem services. Journal of Applied Ecology 53:120–129.

108

Seip, D., and E. Jones. 2010. Response of woodland caribou to partial retention logging of
winter ranges attacked by mountain pine beetle: 2010. FIA/FSP Project # Y102010
2009-10, British Columbia Ministry of Forests and Range, Prince George, B.C. 33 pp.
Serrouya, R., B. N. McLellan, S. Boutin, D. R. Seip, and S. E. Nielsen. 2011. Developing a
population target for an overabundant ungulate for ecosystem restoration. Journal of
Applied Ecology 48:935–942.
Spies, T. A., B. C. McComb, R. S. H. Kennedy, M. T. McGrath, K. Olsen, and R. J. Pabst.
2007. Potential effects of forest policies on terrestrial biodiversity in a multiownership province. Ecological Applications 17:48–65.
Spribille, T. 2018. Relative symbiont input and the lichen symbiotic outcome. Current
Opinion in Plant Biology 44:57–63.
Stevenson, S. K., and D. S. Coxson. 2015. Can partial-cut harvesting be used to manage
terrestrial lichen habitat? A review of recent evidence. Rangifer. Special Issue 35:11–
26.
Steventon, J. D. 1996. Caribou habitat use in the Chelaslie River migration corridor and
recommendations for management. Ministry of Forests, Research Program, Victoria,
B.C. 17 pp.
Stewart, F. E. C., J. J. Nowak, T. Micheletti, E. J. B. McIntire, F. K. A. Schmiegelow, and S.
G. Cumming. 2020. Boreal caribou can coexist with natural but not industrial
disturbances. The Journal of Wildlife Management 84:1435–1444.
St-Laurent, M.-H., Y. Boulanger, D. Cyr, F. Manka, P. Drapeau, and S. Gauthier. 2022.
Lowering the rate of timber harvesting to mitigate impacts of climate change on
boreal caribou habitat quality in eastern Canada. Science of The Total Environment
838:156244.

109

Stockdale, C., Q. Barber, A. Saxena, and M.-A. Parisien. 2019. Examining management
scenarios to mitigate wildfire hazard to caribou conservation projects using burn
probability modeling. Journal of Environmental Management 233:238–248.
Street, G. M., A. R. Rodgers, T. Avgar, and J. M. Fryxell. 2015. Characterizing demographic
parameters across environmental gradients: a case study with Ontario moose (Alces
alces). Ecosphere 6:art138.
Sturtevant, B. R., A. Fall, D. D. Kneeshaw, N. P. P. Simon, M. J. Papaik, K. Berninger, F.
Doyon, D. G. Morgan, and C. Messier. 2007. A toolkit modeling approach for
sustainable forest management planning: achieving balance between science and local
needs. Ecology and Society 12.
Sulyma, R., and D. S. Coxson. 2001. Microsite displacement of terrestrial lichens by feather
moss mats in late seral pine-lichen woodlands of north-central British Columbia. The
Bryologist 104:505–516.
Sutherland, I. J., E. M. Bennett, and S. E. Gergel. 2016. Recovery trends for multiple
ecosystem services reveal non-linear responses and long-term tradeoffs from
temperate forest harvesting. Forest Ecology and Management 374:61–70.
Suzuki, N., and K. L. Parker. 2016. Potential conflict between future development of natural
resources and high-value wildlife habitats in boreal landscapes. Biodiversity and
Conservation 25:3043–3073.
Thom, D., and R. Seidl. 2016. Natural disturbance impacts on ecosystem services and
biodiversity in temperate and boreal forests. Biological Reviews 91:760–781.
Thomas, D. C., S. J. Barry, and G. Alaie. 1996. Fire - caribou - winter range relationships in
northern Canada. Rangifer 16:57–67.

110

Thrower, J. S., A. F. Nussbaum, and C. M. Di Lucca. 1994. Site index curves and tables for
British Columbia - interior species. 2nd edition. Land Management Handbook Field
Guide Insert 6, B.C. Min. For. Res. Branch, Victoria, B.C. 50 pp.
Timmermann, H. R., and J. G. McNicol. 1988. Moose habitat needs. The Forestry Chronicle
64:238–245.
Triviño, M., T. Pohjanmies, A. Mazziotta, A. Juutinen, D. Podkopaev, E. Le Tortorec, and M.
Mönkkönen. 2017. Optimizing management to enhance multifunctionality in a boreal
forest landscape. Journal of Applied Ecology 54:61–70.
United States Fish and Wildlife Service. 1981. Standards for the development of habitat
suitability index models. Release No. 1-81, 103 ESM. Division of Ecological
Services, U.S. Fish and Wildlife Service, Department of the Interior.
Venier, L. A., I. D. Thompson, R. Fleming, J. Malcolm, I. Aubin, J. A. Trofymow, D.
Langor, R. Sturrock, C. Patry, R. O. Outerbridge, S. B. Holmes, S. Haeussler, L. De
Grandpré, H. Y. H. Chen, E. Bayne, A. Arsenault, and J. P. Brandt. 2014. Effects of
natural resource development on the terrestrial biodiversity of Canadian boreal
forests. Environmental Reviews 22:457–490.
Viereck, L. A., and L. A. Schandelmeier. 1980. Effects of fire in Alaska and adjacent
Canada: a literature review. BLM-Technical Report, USDA Forest Service,
Fairbanks, Alaska. 135 pp.
Viljur, M.-L., S. R. Abella, M. Adámek, J. B. R. Alencar, N. A. Barber, B. Beudert, L. A.
Burkle, L. Cagnolo, B. R. Campos, A. Chao, B. Chergui, C.-Y. Choi, D. F. R. Cleary,
T. S. Davis, Y. A. Dechnik-Vázquez, W. M. Downing, A. Fuentes-Ramirez, K. J. K.
Gandhi, C. Gehring, K. B. Georgiev, M. Gimbutas, K. B. Gongalsky, A. Y.
Gorbunova, C. H. Greenberg, K. Hylander, E. S. Jules, D. I. Korobushkin, K. Köster,
V. Kurth, J. D. Lanham, M. Lazarina, A. B. Leverkus, D. Lindenmayer, D. M. Marra,

111

P. Martín-Pinto, J. A. Meave, M. Moretti, H.-Y. Nam, M. K. Obrist, T. Petanidou, P.
Pons, S. G. Potts, I. B. Rapoport, P. R. Rhoades, C. Richter, R. A. Saifutdinov, N. J.
Sanders, X. Santos, Z. Steel, J. Tavella, C. Wendenburg, B. Wermelinger, A. S.
Zaitsev, and S. Thorn. 2022. The effect of natural disturbances on forest biodiversity:
an ecological synthesis. Biological Reviews 97:1930–1947.
Vitt, D. H., L. Finnegan, and M. House. 2019. Terrestrial bryophyte and lichen responses to
canopy opening in pine-moss-lichen forests. Forests 10:233.
Waterhouse, M. J., H. M. Armleder, and A. F. L. Nemec. 2011. Terrestrial lichen response to
partial cutting in lodgepole pine forests on caribou winter range in west-central British
Columbia. Rangifer 31:119–134.
Webb, E. T. 1998. Survival, persistence, and regeneration of the reindeer lichens, Cladina
stellaris, C. rangiferina, and C. mitis following clearcut logging and forest fire in
northwestern Ontario. Rangifer 18:41–47.
Westworth, D., L. Brusnyk, J. Roberts, and H. Veldhuzien. 1989. Winter habitat use by
moose in the vicinity of an open-pit copper mine in north-central British Columbia.
Alces 25:156–166.
Whitman, E., M.-A. Parisien, D. T. Price, M.-H. St‐Laurent, C. J. Johnson, E. R. DeLancey,
D. Arseneault, and M. D. Flannigan. 2017. A framework for modeling habitat quality
in disturbance-prone areas demonstrated with woodland caribou and wildfire.
Ecosphere 8:e01787.
Whittington, J., M. Hebblewhite, N. J. DeCesare, L. Neufeld, M. Bradley, J. Wilmshurst, and
M. Musiani. 2011. Caribou encounters with wolves increase near roads and trails: a
time-to-event approach: Wolf-caribou encounter rates. Journal of Applied Ecology
48:1535–1542.

112

Wittmer, H. U., B. N. Mclellan, R. Serrouya, and C. D. Apps. 2007. Changes in landscape
composition influence the decline of a threatened woodland caribou population.
Journal of Animal Ecology 76:568–579.
Wright, E. F., C. D. Canham, and K. D. Coates. 2000. Effects of suppression and release on
sapling growth for 11 tree species of northern, interior British Columbia. Canadian
Journal of Forest Research 30:10.
Wright, E. F., K. D. Coates, C. D. Canham, and P. Bartemucci. 1998. Species variability in
growth response to light across climatic regions in northwestern British Columbia.
Canadian Journal of Forest Research 28:16.

113

Personal Communication
Crowley, Shannon. February 6, 2019 and July 9, 2020. Ecological Monitoring Coordinator.
John Prince Research Forest. Fort St. James, B.C.
Hodder, Dexter. February 6, 2019 and July 9, 2020. Director of Research and Education.
John Prince Research Forest. Fort St. James, B.C.
Johnson, Chris. 2020. Professor. University of Northern British Columbia. Prince George,
B.C.
Attended half day meeting to workshop linker functions (March 27, 2019):
Cichowski, Debbie. March 27, 2019. Biologist. Caribou Ecological Consulting.
Smithers, B.C.
Clason, Alana. March 27, 2019. Research Associate and Ecologist. Bulkley Valley
Research Centre. Smithers, B.C.
Coates, Dave. March 27, 2019. Retired Research Silviculturist (British Columbia
Ministry of Forests). Smithers, B.C.
Elkin, Ché . March 27, 2019. Associate Professor. FRBC/Slocan Mixedwood
Ecology Chair. University of Northern British Columbia. Prince George, B.C.
Lilles, Erica. March 27, 2019. Research Generalist. British Columbia Ministry of
Forests. Smithers, B.C.

114

Appendix 1
SORTIE-ND: Model Description
http://sortie-nd.org/help/manuals/index.html
SORTIE-ND is an individual-tree, spatially explicit model. This appendix describes the
model using the Overview, Design concepts, Details (ODD) protocol recommended by
Grimm et al. (2006, 2010). This description of SORTIE-ND is sourced from SORTIE-ND
online documentation (Canham and Murphy n.d.). SORTIE-ND can be parameterized to
simulate and output the dynamics of a broad array of tree and stand variables.
Purpose: SORTIE-ND was developed to study neighbourhood processes and forest dynamics.
It was initially developed as a small-scale disturbance model to study transitional oaknorthern hardwood forests in the northeastern United States (Pacala et al. 1996). It has since
been parameterized for Interior Cedar Hemlock moist cool (ICHmc2) and Sub-boreal Spruce
moist cool Babine variant (SBSmc2) forests of northwest BC (Banner et al. 1993, Canham et
al. 2004, Astrup et al. 2008a).
Entities, State Variables, and Scales: There are four hierarchical levels of the model, the
individual (tree or grid cell), the neighbourhood (neighbouring trees within a specified radius
of a target location), the subplot (a specified area within the plot), and the stand (the entire
plot). There are three state entities: individual trees, grid cells (for substrate and light) and the
plot. State variables for trees include: species, diameter at 10 cm above the ground, diameter
at breast height (DBH, measured at 1.35 m above the ground), height, x and y coordinates
(location), live or dead status, crown radius and depth, and years since death (if dead). State
variables for substrate are cell size, substrate type, and time since substrate type was
assigned. State variables for the plot include: length (in X and Y direction) and latitude.
Numerous auxiliary variables can be calculated from state variables including measures of
115

tree volume, basal area, density, carbon, and light. These may be calculated at various scales
(tree, grid cell, subplot and stand level). This description of variables pertains to the model
set-up that has been established for northwest British Columbia. The specific variables that
are available for a model scenario are dependent on the model behaviours that are selected
(i.e., the ecological and biological processes selected for incorporation into the model).
Basic plot and scenario set-up options include timestep, timeframe and plot and grid cell size.
The typical timestep and minimum timeframe is 1 year; there is no maximum timeframe. The
minimum plot size is 100 m by 100 m. There is no maximum plot size however plot size
typically reflects the computational intensity of the simulation, computing capacity, and the
research question. 200 x 200-m or 300 x 300-m are typical plot sizes. The default setting for
grid cell size (pixel) is typically set at 8 x 8- m. A plot can be conceptualized as a torus (or
donut) where a tree on the west ‘edge’ of the plot is modelled as the neighbour of a tree on
the east ‘edge’ of the plot at the same latitude. Thus a 9-ha plot set up as a clearcut (9-ha
clearcut) models the interior condition of a stand with no neighbouring mature trees or edge
effect (i.e., all trees in the plot are surrounded by trees regenerating in the same clearcut
conditions).
Process overview and scheduling: A SORTIE-ND simulation flows in the following process
order: 1. harvest and disturbance activity, 2. light calculation, 3. growth, 4. mortality, 5.
substrate dynamics, 6. snag dynamics, 7. seed dispersal, 8. seedling establishment, 9. planting
activity, 10. output generation.
Design concepts: SORTIE-ND was initially designed to model neighbourhood interactions,
that is, the interactions of trees with their nearest neighbours to study local ecosystem
dynamics. SORTIE-ND extrapolates from measurable fine-scale and short-term interactions
among individual trees to large-scale and long-term dynamics of forest communities (Coates
116

and Hall 2005). SORTIE-ND was initially designed to study ecological processes but has
since been developed to also address questions regarding forest management, natural
disturbance, climate change and carbon dynamics.
Each tree is modelled as an individual with an explicit location; the establishment, growth,
mortality, and regeneration dynamics of each individual is determined by its neighbourhood
and the ecological processes included in the simulation. These processes are modelled in
SORTIE-ND by ‘behaviours’. Researchers can parameterize existing behaviours or develop
new ones specific to their forest type and research question. Different behaviours can be
applied to different tree species and tree life stages. The intent with the design of SORTIEND was the development of a flexible model with a simple structure that provides the user
with the opportunity to control almost every aspect of the model. The description of state
variables and behaviours of SORTIE-ND in this appendix describe the model as
parameterized for the SBS.
The design of SORTIE-ND results in emergent properties. The interaction of behaviours
specific to size and species of individual trees results in changes in species composition of
multi-species stands over time that reflects observed successional processes (Maleki et al.
2020). Light dynamics occurring within stands over time are also emergent properties of the
model and reflect the interaction of species-specific tree establishment, growth and mortality
behaviours.
Behaviours can incorporate stochasticity; currently mortality, snag fall, seed dispersal and
seedling establishment are behaviours with stochastic elements.
Initialization: Initialization of the model involves establishing the state variables for the plot
and initializing the tree population. Trees can be initialized in three ways: 1) a stem map that
identifies the location, species and DBH of each tree, 2) specifying initial densities of trees by
117

species and either height class (for seedlings) or DBH class (for trees with a DBH) or 3) the
tree planting behaviour. Height and DBH classes are established by the user. Trees are
initialized at random locations unless locations are specified with a stem map. Grid cell pixel
size can be specified for substrate and light. Substrate features of grid cells are initialized
from the input parameters.
Input data: Tree maps can be input for starting conditions. The current version of SORTIEND in use in northwest B.C. does not incorporate other external inputs.
Submodels: Grimm et al. (2010) state that it is important to describe sub-models in detail
with respect to model parameters, parameter dimensions, decisions guiding model form
selection, and parameter testing and validation. Behaviours can be viewed as submodels.
Extensive work has been conducted to develop each behaviour and parameterize those that
are applied to forests of northwest B.C. The online SORTIE-ND manual provides detailed
information regarding each behaviour (Canham and Murphy n.d.).

118

Appendix 2
Linker Functions: Model Description
Linker functions have been used to extend the ability of stand dynamics models to predict not
only tree related variables but also habitat elements (e.g. Cordonnier et al. 2013, Mönkkönen
et al. 2014, Bugmann et al. 2017). The linker functions developed in this study described, in a
mathematical form, the relationship between specific stand-level variables and habitat
elements. For example, one linker function described the relationship between understory
light and suitability for forage lichen. Mönkkönen et al. (2014) and Cordonnier et al. (2013)
provide examples of the development of linker functions to describe the suitability of a site to
provide habitat for a specific species or guild of species.
The Grimm et al. (2006, 2010) ODD protocol was developed to describe complex models but
can also be applied to simple models such as the linker functions in this study. In this context,
the purpose of linker functions was to link stand variables (e.g., light, tree, or site variables)
to habitat variables (e.g., relative suitability of a stand to support forage lichen or moose
shrub communities). Variables included understory light, time since disturbance, edaphic site
condition (indices of soil moisture and nutrients), conifer basal area and broadleaf basal area.
These variables can be acquired from SORTIE-ND output. The process overview and
scheduling can be described by describing the nested structure of the model. The full lichen
or moose shrub linker function represents the interactions and relationship between all the
identified stand variables and the habitat variable. The full linker function is comprised of
component linker functions. A component linker function is developed specific to each stand
variable and defines the relationship between the stand variable and the habitat variable.
Component linker functions are combined into a full linker function. Scheduling requires the
following order of calculation: first the component linker function and second the linker

119

function. The design concept behind the development of linker functions is that they are
developed through a combination of available literature, data and expert knowledge. They
can further be refined through field studies targeting the areas where data is least available or
analysis indicates further study is required. The mathematical form of each function was
selected on the principle that they should be simple but also flexible enough to represent the
hypothesized potential relationship.
The first step of the project was an extensive literature review. Informal meetings were then
organized with local experts in the fields of ecology, moose and caribou forage, and lichen
ecology to discuss and draft the linker and component functions. A half-day meeting was
held with Deb Cichowski, Dr. Ché Elkin, Dr. Alana Clason, Erica Lilles, and Dr. Dave
Coates and a shorter meeting occurred with Dexter Hodder, Shannon Crowley and Dr. Ché
Elkin. I also met with Dr. Darwyn Coxson. At these meetings we identified appropriate
predictor and response variables and hypothesized the appropriate functional forms for the
component and full linker functions.
Terrestrial lichen and moose shrub species were determined to be the primary forage groups
that determined the value of a stand for caribou winter habitat. Terrestrial lichens make up the
majority of northern caribou diet during the winter (Miller 1976, Thomas and Hervieux 1986,
Cichowski 1989, Thompson et al. 2015). Woody shrubs comprise the majority of moose
winter diet (Eastman 1977, Rea and Booth 2011, Hodder et al. 2013).
The linker functions were developed following a process similar to the habitat suitability
functions described in Mönkkönen et al. (2014). That is, for each forage type (lichen or
moose shrub) a linker function was developed that represented the relationship of each stand
variable to the index of habitat suitability (i.e., for terrestrial lichens, stand variables include
light availability and proportion of broadleaf trees in the stand). The interactions between
120

component linker functions were either multiplicative (if each component was fundamentally
non-replaceable) or additive (if each item could be considered an alternative to the other) (see
Whitman et al. 2017).
In developing component linker functions, it was important to keep in mind that each
component linker function represented solely the relationship between one stand variable and
the resultant response variable. For example for lichen, the component linker function for
understory light represented only the influence of understory light on the relative suitability
for lichen. It did not reflect any differences in soil moisture or vegetation community that
could occur at different light levels.
Personal Communication
Cichowski, Debbie. March 27, 2019. Biologist. Caribou Ecological Consulting. Smithers,
BC
Clason, Alana. March 27, 2019. Research Associate and Ecologist. Bulkley Valley Research
Centre. Smithers, B.C.
Coates, Dave. March 27, 2019. Retired Research Silviculturist (British Columbia Ministry of
Forests). Smithers, B.C.
Coxson, Darwin. March 25, 2019. Professor. University of Northern British Columbia.
Prince George, B.C.
Crowley, Shannon. February 6, 2019 and July 9, 2020. Ecological Monitoring Coordinator.
John Prince Research Forest. Fort St. James, B.C.
Elkin, Ché . March 27, 2019. Associate Professor. FRBC/Slocan Mixedwood Ecology Chair.
University of Northern British Columbia. Prince George, B.C.
Hodder, Dexter. February 6, 2019 and July 9, 2020. Director of Research and Education.
John Prince Research Forest. Fort St. James, B.C.
Lilles, Erica. March 27, 2019. Research Generalist. British Columbia Ministry of Forests.
Smithers, B.C.

121

Appendix 3
Linker Function Calibration Figures

Figure 13. Comparison of the linker functions for lichen parameterized with the calibrated parameters
(blue) versus parameterized with the maximum likelihood estimates (black). BA = basal area (m2/ha).
Where only a blue line is present, the maximum likelihood estimate was used in the final function.
Empirical data are represented by black dots (n=76) and were collected during the summer of 2019 in
sub-boreal forests of west-central British Columbia.

122

Figure 14. Comparison of the linker functions for moose shrubs used for forage parameterized with
the calibrated parameters (blue) versus parameterized with the maximum likelihood estimates (black).
Where only a blue line is present, the maximum likelihood estimate was used in the final function.
Empirical data are represented by black dots (n=76) and were collected during the summer of 2019 in
sub-boreal forests of west-central British Columbia.

123

Appendix 4
Empirical Data Distribution Figures

Figure 15. Distribution of the abundance of forage lichen, moose shrubs, aspen (0.5 - 3 m in height),
subalpine fir (0.5 - 3 m in height), and conifer basal area in the empirical dataset (n = 76). Data were
collected during the summer of 2019 in sub-boreal forests of west-central British Columbia.

124

.

Appendix 5
Linker Component Function Figures

Figure 16. Output of the component linker functions for forage lichen in three different stand types
(colours) on subxeric sites. The top four rows present the output of the component linker functions,
with the full suitability index for forage lichen shown in row five. Coloured lines represent different
stand types. Where only the conifer mix stand suitability appears on the figure, the other stand types
have the same suitability (e.g., the broadleaf effect for the mature pine stand is the same as the conifer
mix stand). The 5th and 95th quantiles (thin lines) provide measures of the variability in suitability
within the stand. The 4th column illustrates the difference in suitability between a MPB attacked stand
and an undisturbed stand and the 5th column shows the difference between a clearcut stand and a
MPB attacked stand.

125

Figure 17. Output of the component linker functions for forage lichen in four different stand types
(colours) on submesic sites. The top four rows present the output of the component linker functions,
with the full suitability index for forage lichen shown in row five. Where only the spruce stand
suitability appears on the figure, the other stand types have the same suitability (e.g., the broadleaf
effect for the mature pine stand is the same as the spruce stand). The 5th and 95th quantiles (thin lines)
provide measures of the variability in suitability within the stand. The 4th column illustrates the
difference in suitability between a MPB attacked stand and an undisturbed stand and the 5th column
shows the difference between a clearcut stand and a MPB attacked stand.

126

Figure 18. Output of the component linker functions for forage lichen in four different stand types
(colours) on mesic sites. The top four rows present the output of the component linker functions, with
the full suitability index for forage lichen shown in row five. Where only the spruce stand suitability
appears on the figure, the other stand types have the same suitability (e.g., the broadleaf effect for the
mature pine stand is the same as the spruce stand). The 5th and 95th quantiles (thin lines) provide
measures of the variability in suitability within the stand. The 4th column illustrates the difference in
suitability between a MPB attacked stand and an undisturbed stand and the 5th column shows the
difference between a clearcut stand and a MPB attacked stand.

127

Figure 19. Output of the component linker functions for forage lichen in four different stand types
(colours) on subhygric sites. The top four rows present the output of the component linker functions,
with the full suitability index for forage lichen shown in row five. Where only the spruce stand
suitability appears on the figure, the other stand types have the same suitability (e.g., the broadleaf
effect for the mature pine stand is the same as the spruce stand). The 5th and 95th quantiles (thin lines)
provide measures of the variability in suitability within the stand. The 4th column illustrates the
difference in suitability between a MPB attacked stand and an undisturbed stand and the 5th column
shows the difference between a clearcut stand and a MPB attacked stand.

128

Figure 20. Output of the component linker functions for moose browse in three different stand types
(colours) on subxeric sites. The top three rows present the output of the component linker functions
for moose shrubs, with the full suitability index for moose shrubs shown in row four. The suitability
of a stand for subalpine fir and aspen browse are presented in rows five and six and the total moose
browse suitability of a stand in row seven. The bottom row illustrates the suitability of a stand to

129

provide vertical cover. Coloured lines represent different stand types. Where only the conifer mix
stand suitability appears on the figure, the other stand types have the same suitability (e.g., aspen
browse for the mature pine stand is the same as the conifer mix stand). The 5th and 95th quantiles (thin
lines) provide measures of the variability in suitability within the stand. The 4th column illustrates the
difference in suitability between a MPB attacked stand and an undisturbed stand and the 5th column
shows the difference between a clearcut stand and a MPB attacked stand.

130

Figure 21. Output of the component linker functions for moose browse in four different stand types
(colours) on submesic sites. The top three rows present the output of the component linker functions
for moose shrubs, with the full suitability index for moose shrubs shown in row four. The suitability

131

of a stand for subalpine fir and aspen browse are presented in rows five and six and the total moose
browse suitability of a stand in row seven. The bottom row illustrates the suitability of a stand to
provide vertical cover. Coloured lines represent different stand types. Where only the spruce stand
suitability appears on the figure, the other stand types have the same suitability (e.g., aspen browse for
the mature pine stand is the same as the spruce stand). The 5th and 95th quantiles (thin lines) provide
measures of the variability in suitability within the stand. The 4th column illustrates the difference in
suitability between a MPB attacked stand and an undisturbed stand and the 5th column shows the
difference between a clearcut stand and a MPB attacked stand.

132

Figure 22. Output of the component linker functions for moose browse in four different stand types
(colours) on mesic sites. The top three rows present the output of the component linker functions for
moose shrubs, with the full suitability index for moose shrubs shown in row four. The suitability of a
stand for subalpine fir and aspen browse are presented in rows five and six and the total moose
browse suitability of a stand in row seven. The bottom row illustrates the suitability of a stand to

133

provide vertical cover. Coloured lines represent different stand types. Where only the spruce stand
suitability appears on the figure, the other stand types have the same suitability (e.g., aspen browse for
the mature pine stand is the same as the spruce stand). The 5th and 95th quantiles (thin lines) provide
measures of the variability in suitability within the stand. The 4th column illustrates the difference in
suitability between a MPB attacked stand and an undisturbed stand and the 5th column shows the
difference between a clearcut stand and a MPB attacked stand.

134

Figure 23. Output of the component linker functions for moose browse in four different stand types
(colours) on subhygric sites. The top three rows present the output of the component linker functions
for moose shrubs, with the full suitability index for moose shrubs shown in row four. The suitability
of a stand for subalpine fir and aspen browse are presented in rows five and six and the total moose
browse suitability of a stand in row seven. The bottom row illustrates the suitability of a stand to
provide vertical cover. Coloured lines represent different stand types. Where only the spruce stand
suitability appears on the figure, the other stand types have the same suitability (e.g., aspen browse for
the mature pine stand is the same as the spruce stand). The 5th and 95th quantiles (thin lines) provide

135

measures of the variability in suitability within the stand. The 4th column illustrates the difference in
suitability between a MPB attacked stand and an undisturbed stand and the 5th column shows the
difference between a clearcut stand and a MPB attacked stand.

136

Appendix 6
Figures of Suitability Indices for all Stand Types and all Sites

Figure 24. The projected response of lichen suitability over time in undisturbed mature forest (1st
column), after MPB disturbance (2nd column) and after clearcut (3rd column) across four edaphic
sites (rows) and four stand types (colours). The 4th column illustrates the difference in suitability
between a MPB attacked stand and an undisturbed stand and the 5th column shows the difference
between a clearcut stand and a MPB attacked stand. The 5th and 95th quantiles (thin lines) provide
measures of the variability in suitability within the stand.

137

Figure 25. The projected response of moose shrub suitability over time in undisturbed mature forest
(1st column), after MPB disturbance (2nd column) and after clearcut (3rd column) across four edaphic
sites (rows) and four stand types (colours). The 4th column illustrates the difference in suitability
between a MPB attacked stand and an undisturbed stand and the 5th column shows the difference
between a clearcut stand and a MPB attacked stand. The 5th and 95th quantiles (thin lines) provide
measures of the variability in suitability within the stand.

138

Figure 26. The projected response of moose browse suitability (includes aspen and subalpine fir 0.5 –
3 m in height and moose shrubs) over time in undisturbed mature forest (1st column), after MPB
disturbance (2nd column) and after clearcut (3rd column) across four edaphic sites (rows) and four
stand types (colours). The 4th column illustrates the difference in suitability between a MPB attacked
stand and an undisturbed stand and the 5th column shows the difference between a clearcut stand and
a MPB attacked stand. The 5th and 95th quantiles (thin lines) provide measures of the variability in
suitability within the stand.

139

Figure 27. The projected response of the suitability of a stand to provide vertical cover over time in
undisturbed mature forest (1st column), after MPB disturbance (2nd column) and after clearcut (3rd
column) across four edaphic sites (rows) and four stand types (colours). The 4th column illustrates
the difference in suitability between a MPB attacked stand and an undisturbed stand and the 5th
column shows the difference between a clearcut stand and a MPB attacked stand. The 5th and 95th
quantiles (thin lines) provide measures of the variability in suitability within the stand.

140

Figure 28. The projected response of lichen suitability over time in undisturbed mature forest (1st
column), after MPB disturbance (2nd column) and after clearcut (3rd column) across four edaphic
sites (colours) and four stand types (rows). The 4th column illustrates the difference in suitability
between a MPB attacked stand and an undisturbed stand and the 5th column shows the difference
between a clearcut stand and a MPB attacked stand. The 5th and 95th quantiles (thin lines) provide
measures of the variability in suitability within the stand.

141

Figure 29. The projected response of moose shrub suitability over time in undisturbed mature forest
(1st column), after MPB disturbance (2nd column) and after clearcut (3rd column) across four edaphic
sites (colours) and four stand types (rows). The 4th column illustrates the difference in suitability
between a MPB attacked stand and an undisturbed stand and the 5th column shows the difference
between a clearcut stand and a MPB attacked stand. The 5th and 95th quantiles (thin lines) provide
measures of the variability in suitability within the stand.

142

Figure 30. The projected response of moose browse suitability (includes aspen and subalpine fir 0.5 –
3 m in height and moose shrubs) over time in undisturbed mature forest (1st column), after MPB
disturbance (2nd column) and after clearcut (3rd column) across four edaphic sites (colours) and four
stand types (rows). The 4th column illustrates the difference in suitability between a MPB attacked
stand and an undisturbed stand and the 5th column shows the difference between a clearcut stand and
a MPB attacked stand. The 5th and 95th quantiles (thin lines) provide measures of the variability in
suitability within the stand.

143

Figure 31. The projected response of the suitability of a stand to provide vertical cover over time in
undisturbed mature forest (1st column), after MPB disturbance (2nd column) and after clearcut (3rd
column) across four edaphic sites (colours) and four stand types (rows). The 4th column illustrates
the difference in suitability between a MPB attacked stand and an undisturbed stand and the 5th
column shows the difference between a clearcut stand and a MPB attacked stand. The 5th and 95th
quantiles (thin lines) provide measures of the variability in suitability within the stand.

144

Appendix 7
Figures Comparing Variability in Suitability across Stand Types and Edaphic Sites

Figure 32. Comparison of the variability in lichen suitability for four stand types (colours) on four
edaphic sites (rows) as represented by the difference between the 95th suitability quantile and the 5th
suitability quantile at each timestep.

145

Figure 33. Comparison of the variability in lichen suitability between four edaphic sites (colours) for
four stand types (rows) as represented by the difference between the 95th suitability quantile and the
5th suitability quantile at each timestep.

146

Figure 34. Comparison of the variability in moose shrub suitability for four stand types (colours) on
four edaphic sites (rows) as represented by the difference between the 95th suitability quantile and the
5th suitability quantile at each timestep.

147

Figure 35. Comparison of the variability in moose shrub suitability between four edaphic sites
(colours) for four stand types (rows) as represented by the difference between the 95th suitability
quantile and the 5th suitability quantile at each timestep.

148

Figure 36. Comparison of the variability in suitability of a stand to provide vertical cover for four
stand types (colours) on four edaphic sites (rows) as represented by the difference between the 95th
suitability quantile and the 5th suitability quantile at each timestep.

149

Figure 37. Comparison of the variability in suitability of a stand to provide vertical cover between
four edaphic sites (colours) for four stand types (rows) as represented by the difference between the
95th suitability quantile and the 5th suitability quantile at each timestep.

150