THE DEVELOPMENT AND APPLICATION OF SPATIO-TEMPORAL METHODS TO
UNDERSTAND AND PREDICT BROAD-SCALE PATTERNS OF FOREST CHANGE
By
Malcolm S. Itter

A DISSERTATION
Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of
Forestry – Doctor of Philosophy
2018

ABSTRACT
THE DEVELOPMENT AND APPLICATION OF SPATIO-TEMPORAL METHODS TO
UNDERSTAND AND PREDICT BROAD-SCALE PATTERNS OF FOREST CHANGE
By
Malcolm S. Itter
The function, composition, and health of regional forest systems are driven by factors operating at
a range of spatio-temporal scales. Climate shapes regional species composition at centennial-tomillennial timescales, but may also contribute to more rapid forest change through the occurrence of
climate extremes. Disturbance events operate at scales ranging from individual trees to landscapelevel metacommunities impacting forest dynamics and resetting forest succession and development
over decadal-to-centennial time frames. At the local-scale, forest function, composition, and health
at a given time are determined by forest demographic processes including growth, mortality, and
regeneration. Understanding and predicting broad-scale patterns of forest change requires methods
to integrate these different factors synthesizing information across spatio-temporal scales. The
research presented here focuses on the development and application of spatio-temporal, Bayesian
hierarchical methods to advance understanding of the processes and factors driving large-scale
forest change. The methods seek to make inference about latent forest processes of interest based
on noisy observations of forest demographics, climate, and disturbance events.
The impacts of novel climatic conditions forecast to occur over the next century on forest
ecosystem function are difficult to predict given potential interactions between climate, disturbance
events, and forest characteristics such as species composition, density, and tree size/age distribution.
The first three chapters of the following dissertation focus on the development and application of
methods to advance understanding of such interactions. First, a dynamic Bayesian hierarchical
model is presented allowing forest growth responses to climate variables to vary over time in relation to past climate extremes, disturbance events, and forest dynamics. The model was applied
to tree-ring data from a range of sites within northeastern Minnesota. Results revealed significant
growth responses to soil water availability triggered by large climatic water deficits across multiple

seasons and years, forest tent caterpillar defoliation events, and high forest density following large
regeneration events. Building on these results, the interactive effects of past water deficit and insect
defoliation stress on forest growth were further explored using broad-scale tree-ring and defoliation data from two regions of the Canadian boreal forest with contrasting species compositions,
primary insect defoliators, and regional climates. A series of novel methods were developed to
quantify the ecological memory of boreal trees to antecedent water and insect defoliation stress.
Results highlighted the temporal persistence of drought and defoliation stress on boreal tree growth
dynamics and provided an empirical estimate of their interactive effects. Finally, a Bayesian state
space framework for the assimilation of tree-ring and forest inventory data with a forest growth
and yield model (Forest Vegetation Simulator) was developed to reconstruct forest dynamics with
explicit uncertainty. The framework allows for the use of tree-ring data to inform growth-climate
relationships and inventory data to inform estimates of past forest composition, density, and tree
size/age distribution. The unique inference afforded by the framework is demonstrated through its
application to red pine plantation data from northern Minnesota. The final chapter of the dissertation presents a Bayesian point process model for the reconstruction of past fire regimes using
sediment charcoal data. The framework was applied to a network of boreal forest lakes in interior
Alaska demonstrating a significant reduction in the uncertainty of past fire identification compared
to existing methodologies. Further, results highlighted shifts in the regional fire regime coincident
with changes in regional species composition over the past ∟10,000 years.
The methods developed herein and their application to a range of forest data types provide
increased understanding of the multi-scale factors contributing to changes in forest growth and
mortality over time and space. Still missing, however, is a process-based framework that integrates
the various spatio-temporal methods presented to gain mechanistic understanding of forest responses
to extreme climate and disturbance events. Future work is needed to develop such a framework
and apply it to extensive regional forest data sets to advance mechanistic understanding and predict
forest responses to the novel environmental conditions of the 21st century.

To my family and friends for their tireless love and support

iv

ACKNOWLEDGEMENTS

I would like to thank my advisor Dr. Andrew Finley and the rest of my graduate committee Dr.
Phoebe Zarnetske, Dr. Michael Walters, and Dr. Christopher Wikle for their mentorship, guidance,
support, and valuable feedback over the course of my doctoral program. My PhD was enriched by
a number of collaborators and colleagues on the Paleo-Ecological Observatory Network (PalEON)
and Forest Biodiversity and Climate at the Macro Scale projects. In particular, I would like to
thank Dr. Jason McLachlan, Dr. James Clark, Dr. Mevin Hooten, and Dr. Michael Dietze
for providing me with a number of exciting and challenging opportunities to learn and grow as
a graduate student and independent researcher. Further, I would like to thank forestry faculty
members Dr. David Rothstein, Dr. Scott Stark, and Dr. Richard Kobe for their guidance. Thanks
also to forestry administrative staff members Julie Mack, Katie James, and Sandra Dunnebacke
for their valuable assistance during my time in the Department of Forestry. I am grateful for the
support and camaraderie of fellow members (past and present) of the Geospatial Lab within the
Department of Forestry at Michigan State University: Chad Babcock, Gloria Desanker, Megan
Kress, and Neil Ver Planck. The following collaborators were integral to the analyses presented in
the following dissertation: Andria Dawson, Anthony D’Amato, Brian Palik, Daniel Kneeshaw, Jane
Foster, Jason McLachlan, Jennifer Marlon, John Bradford, Loïc D’Orangeville, Louis Duchesne,
Mevin Hooten, Philip Higuera, Ryan Kelly. In addition, I would like to thank and acknowledge R.
Ouimet for sharing Quebec tree-ring and defoliation data, M.-C. Lambert who generated weather
data using the BioSIM model, and C. Whitehouse for providing Alberta defoliation data used
in Chapter 3. The idea for the model presented in Chapter 4 came out of a working group at
the PalEON-sponsored, Climate Ecology and Tree Growth Workshop, Harvard Forest, September
26-29, 2016. I am indebted to Shannon Kay, Patrick Bartlein, and collaborators on the PalEON
project for their valuable contributions and feedback on the model presented in Chapter 5. Finally,
I would like to thank my family and friends without whose love and support this dissertation would
not have been possible.
v

TABLE OF CONTENTS

LIST OF TABLES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

ix

LIST OF FIGURES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

x

CHAPTER 1 INTRODUCTION . . . . . . . . . . . . . . . .
1.1 Factors Driving Forest Change . . . . . . . . . . . . .
1.2 Adaptive Forest Management . . . . . . . . . . . . . .
1.3 Ecological Process & Bayesian Hierarchical Models . .
1.4 Statistical Models for Forest Ecology and Management
1.5 Overview of Chapters . . . . . . . . . . . . . . . . . .

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

1
1
2
3
5
5

VARIABLE EFFECTS OF CLIMATE ON FOREST GROWTH IN RELATION TO CLIMATE EXTREMES, DISTURBANCE, AND FOREST
DYNAMICS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Modeling Approach . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.3.1 Fixed Climate Effects . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.3.2 Variable Climate Effects . . . . . . . . . . . . . . . . . . . . . . . . .
Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.4.1 Tree Growth Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.4.2 Forest Tent Caterpillar . . . . . . . . . . . . . . . . . . . . . . . . . .
2.4.3 Climate Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.5.1 Fixed Climate Effects (FCE) Model . . . . . . . . . . . . . . . . . . .
2.5.2 Variable Climate Effects (VCE) Model . . . . . . . . . . . . . . . . . .
Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.6.1 Fixed Climate Effects . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.6.2 Variable Climate Effects . . . . . . . . . . . . . . . . . . . . . . . . .

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

8
8
9
11
12
14
17
17
19
20
21
21
22
26
27
28

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

34
34
35
38
38
38
39
41
43

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

CHAPTER 2
2.1
2.2
2.3
2.4

2.5
2.6

CHAPTER 3
3.1
3.2
3.3

BOREAL TREE GROWTH EXHIBITS DECADAL-SCALE ECOSYSTEM MEMORY TO DROUGHT AND INSECT DEFOLIATION, BUT
NO NEGATIVE RESPONSE TO THEIR INTERACTION . . . . . . . .
Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.3.1 Study Regions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.3.1.1 West . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.3.1.2 East . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.3.2 Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.3.3 Model Specification . . . . . . . . . . . . . . . . . . . . . . . . . . . .

vi

3.4
3.5
3.6

Results . . . . . . . . . . .
3.4.1 Ecological Memory
3.4.2 Antecedent Effects
Discussion . . . . . . . . .
3.5.1 Ecological Memory
3.5.2 Antecedent Effects
Conclusions . . . . . . . .

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

47
47
48
51
51
54
56

ASSIMILATION OF TREE-RING AND FOREST INVENTORY DATA
TO MODEL INTERACTIONS BETWEEN CLIMATE AND FOREST
DYNAMICS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Data and Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.3.1 Forest Vegetation Simulator . . . . . . . . . . . . . . . . . . . . . . .
4.3.2 State Space Model . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.3.3 Birch Lake Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.5.1 Model Applications . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.5.2 Management Applications . . . . . . . . . . . . . . . . . . . . . . . .
4.5.3 Future Steps . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.

58
58
59
62
62
65
69
70
75
76
77
78
79

A MODEL-BASED APPROACH TO WILDLAND FIRE RECONSTRUCTION USING SEDIMENT CHARCOAL RECORDS . . . . . . . . . . .
Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Bayesian Point Process Model for Charcoal Deposition . . . . . . . . . . . . . .
5.3.1 Univariate Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.3.1.1 Probability of Fire . . . . . . . . . . . . . . . . . . . . . . . .
5.3.1.2 Mean Fire Return Interval . . . . . . . . . . . . . . . . . . . .
5.3.2 Multivariate Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.3.2.1 Regional Mean Fire Return Interval . . . . . . . . . . . . . . .
5.3.3 Bayesian Implementation . . . . . . . . . . . . . . . . . . . . . . . . . .
5.3.3.1 Univariate Model . . . . . . . . . . . . . . . . . . . . . . . . .
5.3.3.2 Multivariate Model . . . . . . . . . . . . . . . . . . . . . . . .
5.3.3.3 Identification of Optimal Penalties . . . . . . . . . . . . . . . .
5.3.3.4 Mean FRI . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Model Application . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.4.1 Simulation Study . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.4.2 Yukon Flats . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.4.2.1 Univariate Model . . . . . . . . . . . . . . . . . . . . . . . . .
5.4.2.2 Multivariate Model . . . . . . . . . . . . . . . . . . . . . . . .
Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

81
81
81
84
85
88
88
89
91
92
92
92
94
95
96
96
97
99
100
103

CHAPTER 4
4.1
4.2
4.3

4.4
4.5

4.6

CHAPTER 5
5.1
5.2
5.3

5.4

5.5

vii

CHAPTER 6 CONCLUSIONS
6.1 Research Synthesis . .
6.2 Future Research . . . .
6.3 Concluding Remarks .

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

108
108
110
112

APPENDICES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114
APPENDIX A
VARIABLE EFFECTS OF CLIMATE ON FOREST GROWTH
IN RELATION TO CLIMATE EXTREMES, DISTURBANCE,
AND FOREST DYNAMICS . . . . . . . . . . . . . . . . . . . . . 115
APPENDIX B
BOREAL TREE GROWTH EXHIBITS DECADAL-SCALE ECOSYSTEM MEMORY TO DROUGHT AND INSECT DEFOLIATION,
BUT NO NEGATIVE RESPONSE TO THEIR INTERACTION . . 130
APPENDIX C
ASSIMILATION OF TREE-RING AND FOREST INVENTORY
DATA TO MODEL INTERACTIONS BETWEEN CLIMATE
AND FOREST DYNAMICS . . . . . . . . . . . . . . . . . . . . . 138
APPENDIX D
A MODEL-BASED APPROACH TO WILDLAND FIRE RECONSTRUCTION USING SEDIMENT CHARCOAL RECORDS . 144
BIBLIOGRAPHY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 153

viii

LIST OF TABLES

Table 2.1: Summary of species sampled. . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
Table 2.2: Summary of seasonal aggregations of climate variables. Bullets indicate that
a seasonal aggregation is calculated for a given variable. . . . . . . . . . . . . . 20
Table 2.3: Summary of tree growth sensitivity to climate variables (reference for Figure 2.4). 23
Table 3.1: Posterior summary of stand-level, antecedent variable coefficients. Posterior
median coefficient values are provided for antecedent insect defoliation, antecedent climatic water deficit, and their interaction for each study region by
species (defoliator host vs. non-host) and size (large- vs. small-diameter)
categories. 95 percent credible intervals are given in parentheses. Coefficients
for which credible intervals do not include zero are bolded. . . . . . . . . . . . . 50
Table 4.1: Posterior summary of observation and process variance parameters. . . . . . . . 74
Table 5.1: Summary of univariate and multivariate Poisson process model results for each
lake in the Yukon Flats data set. Mean FRI is equal to the posterior mean fire
return interval with 95 percent credible intervals in parentheses. . . . . . . . . . 101
Table A.1: Posterior distribution summary for model variance terms. Posterior mean
values are provided with 95 percent credible intervals in parentheses. FCE
= fixed climate effects model, VCE = variable climate effects model. Note,
σy2 =

2
σpe
.
1−φ2

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129

Table B.1: Posterior summary of tree size effect coefficients and tree- and stand-level
variances for the East and West study regions. Posterior median values are
reported along with 95 percent credible intervals in parentheses. Note a single
linear term is used to control for tree size in the East, whereas host/non-host
specific linear and quadratic terms are used in the West. . . . . . . . . . . . . . . 133
Table C.1: Shape (a x ) and scale (b x ) hyperparameter values for inverse gamma prior
distributions for each observation and process variance parameter (x indicates
the variable: dtr , dfc , nfc , D, N). . . . . . . . . . . . . . . . . . . . . . . . . . . 140

ix

LIST OF FIGURES

Figure 2.1: Schematic of fixed climate effects and variable climate effects tree growth
models. Note: Fixed and variable climate effects sub-models are applied
separately; 0:T subscript indicates all time points from 0 to T. . . . . . . . . . . 15
Figure 2.2: Location of 105 forest plots (35 stands, 3 plots per stand) in relation to
Superior National Forest in northeastern Minnesota, USA. . . . . . . . . . . . . 18
Figure 2.3: Standardized coefficient values for fixed climate effects model. Note: Points
represent posterior median coefficient estimates; black lines indicate 95 percent credible interval bounds; a dashed grey line at a coefficient value of zero
is provided for reference. Variable abbreviations are as follows: fall deficit
(FAL-DEF), spring deficit (SPR-DEF), summer deficit (SUM-DEF), summer
deficit in the previous growing season (SUM-DEF-LAG), and mean annual
snow pack (SNOW). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
Figure 2.4: Evolution of climate coefficient values for each climate variable over the
study period (1897-2007). Solid black line and points indicate posterior
mean coefficient values. Dashed lines delineate 95 percent credible intervals.
Points are colored to indicate different response categories. Zero Response:
credible interval includes zero; Weak Response: credible interval does not
contain zero, but annual r 2 < 0.25; Threshold Response: strong response to
climate (credible interval does not contain zero and annual r 2 > 0.25) within
two years of threshold exceedance; Persistent Response: strong response to
climate in years immediately following a threshold exceedance; Disturbance
Response: strong response to climate in years where forest tent caterpillar is
present in study region; Unknown Response: strong response to climate not
attributable to threshold exceedance, persistent response, or disturbance. The
upper right panel indicates the number of study stands from which tree growth
data exists in relation to the study period. The grey shading in each panel
indicates the period during which tree growth data from fewer than 20 study
stand are available (1897-1929). Variable abbreviations are as follows: fall
deficit (FAL-DEF), spring deficit (SPR-DEF), summer deficit (SUM-DEF),
summer deficit in the previous growing season (SUM-DEF-LAG), and mean
annual snow pack (SNOW). . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24

x

Figure 2.5: Observed climate variable values over the study period (1897-2007). Black
lines indicate mean climate variable values across study stands along with
uncertainty levels equal to two times the standard error (grey shading). The
horizontal red line indicates the estimated climate threshold for each variable (thresholds correspond to the following quantiles, 0.95: summer deficit,
lagged summer deficit; 0.98: fall deficit; 0.85: spring deficit, snow). White
filled points indicate threshold exceedances with no growth response. Red
filled points indicate threshold exceedances for which a growth response was
observed in the five-year period centered on the year of exceedance. Orange
filled points indicate a disturbance response. Variable abbreviations are as
follows: fall deficit (FAL-DEF), spring deficit (SPR-DEF), summer deficit
(SUM-DEF), summer deficit in the previous growing season (SUM-DEFLAG), and mean annual snow pack (SNOW). . . . . . . . . . . . . . . . . . . . 25
Figure 2.6: Violin plots of partial residuals (observed log annual growth increment minus
spline-based estimate) for forest tent caterpillar (FTC) host and non-host
individuals located in stands in the 5th percentile for growth in years a)
1950-1954; b) 1991-1993. Stands in the lowest 5th percentile for growth are
considered likely to have been affected by FTC. Note: Black shading indicates
FTC host trees; grey shading indicates non-FTC host trees. . . . . . . . . . . . 26
Figure 2.7: Cumulative distribution of forest growth responses to all climate variables
under the variable climate effects model as a function of years since initiation
(initiation marks establishment a new cohort of trees within a study stand)
for the threshold/persistent and unknown response categories. Note, the
cumulative growth response functions exclude growth responses that occurred
early in the study period when growth data from less than 20 study stands
were available (1987-1929). Grey shading highlights responses 20 to 45 years
following initiation when understory stem density and inter-tree competition
are high. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
Figure 3.1: (A) Location of study stands within western (Alberta: 34 stands) and eastern
(Quebec: 14 stands) study regions relative to mean summer (June-August) climate moisture index values (1970-2000, Fick and Hijmans, 2017). (B) Mean
summer 3-month Standardized Precipitation-Evapotranspiration Index values
for western and eastern study stands (1902-2015, BeguerĂ­a and Serrano, 2017).

40

Figure 3.2: Diagram of hierarchical model structure including tree- and stand-level submodels and the estimation of antecedent water deficit and antecedent insect
defoliation variables. Rectangles delineate observations, while ovals delineate estimated quantities. Note t0 indicates the beginning of the study period
and L indicates the maximum number of years before present for which
ecological memory to past climatic water deficit is estimated. . . . . . . . . . . 45

xi

Figure 3.3: Ecological memory functions for climatic water deficit (A) and insect defoliation (B) based on annual radial tree growth in eastern and western study
regions. Lag indicates the number of years before current year growth. . . . . . 48
Figure 3.4: Estimated antecedent climatic water deficit (A) and insect defoliation (B)
values for simulated water deficits of 10 mm and insect defoliation events
lasting one, two, and five years. Dashed lines highlight the year for which
maximum antecedent values are estimated. . . . . . . . . . . . . . . . . . . . . 49
Figure 3.5: Estimated effects of antecedent insect defoliation, antecedent climatic water
deficit, and their interaction on tree growth in eastern and western study
regions by species (defoliator host vs. non-host) and size (large- vs. smalldiameter) categories. Points represent posterior median antecedent variable
values based on study data. Relative response surfaces correspond to mean
tree growth under antecedent conditions relative to regional mean tree growth
over the study period (East: 1968-1998, West: 1968-2010) after controlling
for tree size. Response surfaces were generated by imposing a dense grid
over the range of modeled antecedent variable values. . . . . . . . . . . . . . . 51
Figure 4.1: Overview of state space model to approximate forest dynamics driven by
FVS-LSlite and constrained by tree-ring and forest inventory data. . . . . . . . 66
Figure 4.2: Graphical depiction of state space model including (A) the estimation of
individual tree diameter at breast height, and (B) the estimation of tree density.
Arrows indicate dependence structure. . . . . . . . . . . . . . . . . . . . . . . 68
Figure 4.3: Posterior summary of stand-scale variables including trees per acre (A), basal
area per acre (B), and quadratic mean diameter (C) relative to forest inventory
observations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
Figure 4.4: Posterior summary of annual basal area increment and basal area mortality
relative to forest inventory estimates equal to the average annual basal area
growth or mortality between inventory observations. . . . . . . . . . . . . . . . 73
Figure 4.5: Posterior summary of diameter at breast height state space estimates relative
to tree-ring and forest inventory diameter observations for a selection of five
red pine trees. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74

xii

Figure 5.1: Illustration of theoretical charcoal deposition to a lake if charcoal particles
arising from regional fires were distinguishable from particles arising from
local fires (in practice, charcoal particles from different sources are indistinguishable, i.e., we observe only y(τ j,i ) = yb (τ j,i ) + y f (τ j,i )). The figure also
depicts the formation of a sediment charcoal record from deposited charcoal
and the structure of the charcoal data (i.e., counts of charcoal over discrete,
non-overlapping time intervals. Note, the figure does not depict charcoal
arising from secondary sources such as surface water runoff or sediment mixing. 85
Figure 5.2: Example creation of new temporal support τ ∗ from three individual lake
records. Note the temporal misalignment across lakes and τ ∗ . . . . . . . . . . . 90
Figure 5.3: Univariate model results for simulated sediment charcoal record generated
by CharSim. Upper panel indicates observed charcoal counts along with the
posterior mean charcoal count (blue line). Second panel illustrates posterior
mean foreground (orange line) and background (black line) intensities. Third
panel plots posterior mean probability of fire estimates for each observed time
interval (black line) along with the upper and lower bounds of the 95 percent
credible interval (gray shading) and the optimal threshold (red line). The
points in the third panel correspond to true fire events occurring during the
sample interval with black dots delineating correctly identified fires and red
open dots delineating missed fires. The arrow highlights the single falsely
identified fire during the simulated study period. Lower panel indicates
observed charcoal counts with the color and shape indicating whether the
count was correctly identified as a true fire (black points), no fire (gray
points), missed true fire (red open circles), or falsely identified fire (red points).

98

Figure 5.4: Location of study lakes within the Yukon Flats region of Alaska relative to
areas burned in wildland fire since 1940 Common Era. . . . . . . . . . . . . . . 99
Figure 5.5: Univariate model results for Chopper Lake (a) and Screaming Lynx Lake
(b). Upper panel indicates observed charcoal counts along with the posterior
mean charcoal count (blue line). Middle panel illustrates posterior mean
foreground (orange line) and background (black line) intensities. Lower
panel plots posterior mean probability of fire estimates for each observed
time interval (black line) along with the upper and lower bounds of the 95
percent credible interval (gray shading) and the optimal threshold (red line). . . 102
Figure 5.6: Multivariate model results for Chopper Lake (a) and Screaming Lynx Lake
(b). Upper panel indicates observed charcoal counts along with the posterior
mean charcoal count (blue line). Middle panel illustrates posterior mean
foreground (orange line) and background (black line) intensities. Lower
panel plots posterior mean probability of fire estimates for each observed
time interval (black line) along with the upper and lower bounds of the 95
percent credible interval (gray shading) and the optimal threshold (red line). . . 102
xiii

Figure 5.7: Background charcoal deposition intensity for each lake in the Yukon Flats
data set over the entire study period (10,680 to -59 YBP, relative to 1950
CE) based on the multivariate point process model. The bold line in the
upper panel is a regional loess smooth function reflecting mean changes in
background charcoal deposition across all lakes (fit using a span of 0.15). The
lower panel indicates the number of lake records used to estimate background
charcoal deposition over time. . . . . . . . . . . . . . . . . . . . . . . . . . . . 104
Figure A.1: Graphical depiction of climate coefficient evolution and variable climate effects (VCE) model dependencies with arrows indicating direct dependence. . . . 118
Figure A.2: Results of application of the Bayesian Lasso to full set of climate variables
excluding seasonal aggregations of actual evapotranspiration (AET). Results
are similar when seasonal aggregations of AET are included and seasonal
aggregations of potential evapotranspiration (PET) are excluded. Points represent the posterior median coefficient value for each climate variable with
95 percent credible intervals (CIs) indicated with solid lines. A dashed line
corresponding to a coefficient value of zero is provided for reference. Points
and lines are color coded according to whether 95 percent CIs overlap zero:
green - the lower bound of the 95 percent CI is greater than zero, red - the
upper bound of the 95 percent CI is less than zero, grey - the 95 percent CI
includes zero. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128
Figure B.1: Antecedent climatic water deficit values for the 14 eastern (Quebec) study
stands over the study period (1968-1998) estimated by applying posterior
samples of antecedent weights to modeled climatic water deficit. . . . . . . . . 134
Figure B.2: Antecedent climatic water deficit values for the 34 western (Alberta) study
stands over the study period (1968-2010) estimated by applying posterior
samples of antecedent weights to modeled climatic water deficit. . . . . . . . . 135
Figure B.3: Antecedent defoliation values for the 14 eastern (Quebec) study stands over
the study period (1968-1998) estimated by applying a spherical decay function
using posterior samples of the decay parameter (φ) to observed defoliation events.136
Figure B.4: Antecedent defoliation values for the 34 western (Alberta) study stands over
the study period (1968-2010) estimated by applying a spherical decay function
using posterior samples of the decay parameter (φ) to observed defoliation events.137
Figure C.1: Posterior summary of stand-scale variables including trees per acre (A), basal
area per acre (B), and quadratic mean diameter (C) relative to forest census
observations relative to an unconstrained run of the Lake States Variant of the
USDA Forest Service Forest Vegetation Simulator (FVS) growth and yield
model for a red pine plantation started from bare ground. . . . . . . . . . . . . 139

xiv

Figure D.1: Uncertainty in local fire identification for Chopper Lake based on the univariate point process model. Upper panel presents the coefficient of variation
(CV) for the mean fire return interval as a function of the probability of fire
threshold. The probability of fire threshold corresponding to the minimum
CV is selected as the optimal threshold value. Lower panel presents the
proportion of posterior samples (1,000 samples, post burn-in) which identify
a local fire in year t applying the optimal probability of fire threshold. The
number of fires refers to the posterior median number of local fires identified over the length of the sediment charcoal record (the 95 percent credible
interval is provided in parentheses). . . . . . . . . . . . . . . . . . . . . . . . . 145
Figure D.2: Uncertainty in local fire identification for Screaming Lynx Lake based on
the univariate point process model. Upper panel presents the coefficient of
variation (CV) for the mean fire return interval as a function of the probability
of fire threshold. The probability of fire threshold corresponding to the
minimum CV is selected as the optimal threshold value. Lower panel presents
the proportion of posterior samples (1,000 samples, post burn-in) which
identify a local fire in year t applying the optimal probability of fire threshold.
The number of fires refers to the posterior median number of local fires
identified over the length of the sediment charcoal record (the 95 percent
credible interval is provided in parentheses). . . . . . . . . . . . . . . . . . . . 146
Figure D.3: Univariate model results for lakes in Yukon Flats network. Upper panel
indicates observed charcoal counts along with the posterior mean charcoal
count (blue line). Middle panel illustrates posterior mean foreground (orange
line) and background (black line) intensities. Lower panel plots posterior
mean probability of fire estimates for each observed time interval (black line)
along with the upper and lower bounds of the 95 percent credible interval
(gray shading) and the optimal threshold (red line). . . . . . . . . . . . . . . . 147
Figure D.4: Multivariate model results for lakes in Yukon Flats network. Upper panel
indicates observed charcoal counts along with the posterior mean charcoal
count (blue line). Middle panel illustrates posterior mean foreground (orange
line) and background (black line) intensities. Lower panel plots posterior
mean probability of fire estimates for each observed time interval (black line)
along with the upper and lower bounds of the 95 percent credible interval
(gray shading) and the optimal threshold (red line). . . . . . . . . . . . . . . . 150

xv

CHAPTER 1
INTRODUCTION

1.1

Factors Driving Forest Change

Forests are impacted by factors operating at a variety of spatio-temporal scales. Climate
affects regional species composition at centennial-to-millennial timescales (Davis and Shaw, 2001).
The occurrence of climate extremes, such as drought events, affects forests over shorter time
frames, altering demographic processes and in severe cases contributing to large-scale tree mortality
(Breshears et al., 2005; Allen et al., 2015). Disturbance events impact forest systems at a range of
scales. Wind, for example, blows down trees within a forest stand creating irregularly-sized canopy
gaps that may be filled by advanced regeneration or newly-recruited individuals increasing diversity
across the forested landscape on a decadal-to-centennial timescale (Bormann and Likens, 1979;
Lorimer and White, 2003). Stand-replacing fire, in contrast, kills the majority of trees within a
stand, releases soil nutrients, and resets successional pathways for centuries to come (Turner et al.,
2007). Finally, the relative abundance of regional species within a given forest stand and the stand’s
structure are shaped on annual-to-decadal scales by forest demographic processes such as growth,
mortality, dispersal, and recruitment (Pacala et al., 1996; Hansen et al., 2001; Clark et al., 2010).
The different factors impacting forest systems are not independent. Rather, there are complex
interactions between climate, disturbance, and forest demographic processes, which likely involve
feedbacks, threshold effects, and non-linear forest responses (Allen et al., 2015). As an example,
consider a lodgepole pine (Pinus contorta Dougl.) forest growing in the Rocky Mountain West.
A mountain pine beetle (Dendroctonus ponderosae Hopk.) outbreak within the forest may impact
the xylem structure of individual trees (Safranyik and Carrol, 2006). Xylem-damaged trees are
more likely to suffer cavitation (an embolism within the xylem) during a subsequent drought event
potentially leading to mortality if drought severity surpasses some unknown physiological threshold
(McDowell et al., 2008, 2011). The existence of standing dead lodgepole pine trees increases the
1

risk of catastrophic fire due to the prevalence of ladder fuels (Bigler and Veblen, 2011, but see
Simard et al. (2011) for evidence of reduced probability of crown fire due to bark beetle mortality
in lodgepole pine stands). Worse, if drought conditions continue following the fire, they may
contribute to regeneration failure of lodgepole pine thereby altering the composition and function
of the forest well into the future (Stephens et al., 2013; Harvey et al., 2016).

1.2

Adaptive Forest Management

Forests face an uncertain future as rapidly changing environmental conditions lead to altered
climate and shifts in historic disturbance regimes (Allen et al., 2015). Adaptive forest management is
focused on maintaining healthy, productive, and fully-functional forest ecosystems under changing
environmental conditions (Millar et al., 2007; D’Amato et al., 2011). In the short-term, this includes
managing forests to generate characteristics (e.g., composition, density, structure) which promote
resistance and resilience to changing conditions (Millar et al., 2007). For example, thinning
treatments have been shown to reduce forest growth responses to drought (resistance) and/or speed
return to pre-drought growth rates (resilience) in a variety of forest ecosystem types (Laurent et al.,
2003; McDowell et al., 2006; D’Amato et al., 2013; Sohn et al., 2016; Bottero et al., 2017). Over
the long-term, the goal of adaptive management is to facilitate the formation of robust, complex
adaptive systems capable of maintaining ecosystem function under uncertain and variable future
conditions (Puettmann, 2011).
The development of both short- and long-term adaptive management approaches requires mechanistic understanding of the joint effects of the different factors impacting forest demographic processes and their interactions (Purves and Pacala, 2008; Vanderwel and Purves, 2014; Clark et al.,
2014). Such understanding enables the prediction of forest change and allows for the development
and testing of adaptive management strategies. Factors affecting demographic processes can be
partitioned into exogenous factors such as climate and disturbance, and endogenous forest factors
including species composition, density, canopy structure, and tree size/age distribution (Bormann
and Likens, 1979). Changes in endogenous factors over time are collectively referred to as forest

2

dynamics and are particularly important in adaptive management contexts because they can be
modified through management and may buffer forest responses to climate and disturbance (Oliver
and Larson, 1996). Thus, if we can identify forest conditions promoting robust ecosystems, we
may be able to develop management strategies to maintain critical forest processes in the face of
environmental change (Millar et al., 2007; D’Amato et al., 2011; Puettmann, 2011).

1.3

Ecological Process & Bayesian Hierarchical Models

As discussed above, we seek mechanistic understanding of joint forest demographic responses
to a variety of factors operating at different spatio-temporal scales. This is a common challenge in
ecology where we would like to make inference about ecological processes and the factors that drive
them based on noisy, often indirect observations. Bayesian hierarchical methodologies are ideally
suited to meet this challenge as they allow disparate observations at varying spatio-temporal scales
to be combined within an integrated statistical framework in order to learn about an ecological
process of interest. This is an area of much research and has been developed and described in a
wide range of statistical and ecological literature (Berliner, 1996; Wikle et al., 2001; Wikle, 2003;
Calder et al., 2003; Clark, 2005; Hooten et al., 2007; Ogle and Barber, 2008; Cressie et al., 2009;
Cressie and Wikle, 2011; Dietze et al., 2013; Hobbs and Hooten, 2015). Berliner (1996) offers
a particularly elegant summary of the power of Bayesian hierarchical methods in the context of
time series analysis. Specifically, the mechanics of a Bayesian hierarchical model in which the
inferential goal is advancing understanding of an underlying process can be written as
(i)

(ii)

[process,parameters|data] ∝ [data|process,parameters]×
[process|parameters] × [parameters],
(iii)

(1.1)

(iv)

where in the current context, process refers to an ecological process of interest, parameters are a
collection of unknown functional parameters describing the process and its variability, and data
comprises noisy, potentially indirect observations of the process (Berliner, 1996). In Equation
1.1, distribution (i) represents the posterior distribution providing estimates of the process and its
functional parameters conditional on a set of noisy observations. Distribution (ii) corresponds to
3

the likelihood or data model conditional on the latent process of interest and unknown functional
parameters. Distribution (iii) describes the latent process of interest and importantly, variability in
the process given its imperfect representation and uncertainty in unknown functional parameters.
Distribution (iv) is the prior distribution (or set of prior distributions) representing prior knowledge
of all unknown functional and variance parameters.
Central to Bayesian hierarchical models such as the one presented in Equation 1.1 are functions
describing the ecological process of interest (i.e., process models). Process models represent
current scientific understanding of the drivers or mechanisms underlying an ecological process—
they are the mathematical expression of a scientific hypothesis (Wikle et al., 2001). Frequently in
ecology, these models describe changes in ecological processes over time and space. Statistical
data models are integrated alongside process models to account for the structure and variability
of different observations (i.e., likelihood functions as defined by distribution (ii) in Equation 1.1).
The combination of process and data models allows for inference about a latent ecological process
of interest based on noisy observations. In particular, we are able to partition process error (i.e.,
error in scientific understanding of a process) from observation error (i.e., noise due to imperfect
sampling).
The partitioning of process and observation error is extremely powerful for the advancement
of scientific understanding. For example, we can use process-based, Bayesian hierarchical models
to identify sources of uncertainty indicating whether estimates of an ecological process are limited
by poor scientific understanding (a mis-specified process model) or lack of sufficient data. Such
frameworks also allow for competing scientific hypotheses to be tested—does one model significantly reduce process error relative to an alternative model? Often overlooked in such frameworks
is the use of prior distributions to further integrate scientific understanding (Hobbs and Hooten,
2015). Specifically, we might use a semi-informative prior to constrain a well-studied parameter
to a scientifically-plausible range. The structure of process-based, Bayesian hierarchical models
(as represented in Equation 1.1) are the foundation for a variety of advanced statistical methods
including Bayesian data assimilation, model-data fusion, state space models, and a wide array of

4

spatio-temporal models (Cressie and Wikle, 2011).

1.4

Statistical Models for Forest Ecology and Management

The goal of the research presented here is to advance understanding of the factors impacting
forest demographic processes over time and space in order to inform forest management to maintain
these processes in the face of changing environmental conditions. Bayesian hierarchical models
following the conceptual framework presented in Equation 1.1 are developed and applied to a range
of forest data types. The models synthesize noisy observations of forest demographic processes
such as growth, mortality, and regeneration with mathematical approximations of these processes.
Observations are imperfect and subject to measurement error. Models of demographic processes
may be empirical and serve as the first step toward building a mechanistic model, or they may be
mechanistic and represent current scientific understanding of the demographic process of interest.
Whether they are empirical or mechanistic, the demographic process models are also imperfect
and subject to potential mis-specification and process error. Understanding of the different factors
driving changes in forest demographic processes resulting from the application of the Bayesian
hierarchical models described herein can be used to identify forest characteristics which promote
resistance and resilience to changing environmental conditions. These characteristics serve as
targets for forest management, couching ecological results within an adaptive management context.

1.5

Overview of Chapters

The following dissertation includes four chapters. The first three chapters focus on understanding interactions between climate, non-stand replacing disturbance (e.g., insect defoliation, wind),
and forest dynamics. The chapters are related with the methods evolving and advancing from one
chapter to the next. The fourth chapter develops a model to reconstruct regional fire regimes. It is
distinct from the first three chapters in its focus on paleo-fire reconstruction, although the modeling
approach has many similarities to the models applied in the first three chapters. A brief summary
of each chapter is provided below.

5

Chapter 2 - Time varying effects of climate on tree growth: A state space model was developed,
which allows the effects of seasonal water deficits on inter-annual radial tree growth to vary over
time and in relation to climate extremes, insect defoliation, and forest dynamics. The model was
applied to tree ring-width increment data from a collection of stands of varying composition,
structure, and development stage in northeastern Minnesota. Potential interactions between tree
growth responses to extreme water deficits, insect defoliation events, and endogenous forest factors
(e.g., density) were considered, but not explicitly tested.
Chapter 3 - Ecological memory of boreal tree growth to drought and insect defoliation and
their interaction: A hierarchical model was developed to quantify the ecological memory of tree
growth to past climatic water deficit and insect defoliation, derive antecedent variables reflecting
the persistent and cumulative effects of these stressors on tree growth, and test for their interactive
effects. The model builds on the approach applied in Chapter 2, but importantly allow for the
interactive effects of water deficit and insect defoliation on tree growth to be tested. The model
was applied to extensive tree growth, weather, and defoliation survey data from western and eastern
regions of the Canadian boreal forest characterized by contrasting tree compositions, climates, and
insect defoliators.
Chapter 4 - Assimilation of tree rings and forest census data to model past forest dynamics: An
approach to quantify interactions between climate extremes, disturbance, and forest dynamics was
formalized. Specifically, a state space model was developed to reconstruct past forest dynamics and
their interaction with climate extremes and non-stand replacing disturbance events. The state space
framework is driven by a forest growth and yield model and assimilates individual tree ring-width
series and forest inventory data to constrain estimates of forest dynamics and advance understanding
of the impacts of climate extremes and disturbance events on forest growth and mortality. The
model was applied to a collection of red pine plantations in northern Minnesota defined by different
thinning regimes.
Chapter 5 - Model-based approach to wildland fire reconstruction: A point process model
to reconstruct wildland fire regimes using sediment charcoal records was developed. The model

6

explicitly estimates the probability of local fire occurrence associated with sediment charcoal
counts and allows for charcoal records from multiple regional lakes to be pooled in order to reduce
uncertainty in regional fire history. The point process model was applied to sediment charcoal
records from a network of 13 boreal lakes located in the Yukon Flats region of interior Alaska to
reconstruct its fire regime over the past 10,000 years.
A synthesis of the key findings from the four chapters and a discussion of unanswered questions
and future research directions is provided in the final chapter (Chapter 6).

7

CHAPTER 2
VARIABLE EFFECTS OF CLIMATE ON FOREST GROWTH IN RELATION TO
CLIMATE EXTREMES, DISTURBANCE, AND FOREST DYNAMICS

2.1

Abstract

Changes in the frequency, duration, and severity of climate extremes are forecast to occur under
global climate change. The impacts of climate extremes on forest productivity and health remain
difficult to predict due to potential interactions with disturbance events and forest dynamics—
changes in forest stand composition, density, size and age structure over time. Such interactions
may lead to non-linear forest growth responses to climate involving thresholds and lag effects.
Understanding how forest dynamics influence growth responses to climate is particularly important
given stand structure and composition can be modified through management to increase forest
resistance and resilience to climate change. To inform such adaptive management, we develop
a hierarchical Bayesian state space model in which climate effects on tree growth are allowed to
vary over time and in relation to past climate extremes, disturbance events, and forest dynamics.
The model is an important step toward integrating disturbance and forest dynamics into predictions
of forest growth responses to climate extremes. We apply the model to a dendrochronology data
set from forest stands of varying composition, structure, and development stage in northeastern
Minnesota that have experienced extreme climate years and forest tent caterpillar defoliation events.
Mean forest growth was most sensitive to water balance variables representing climatic water deficit.
Forest growth responses to water deficit were partitioned into responses driven by climatic threshold
exceedances and interactions with insect defoliation. Forest growth was both resistant and resilient
to climate extremes with the majority of forest growth responses occurring after multiple climatic
threshold exceedances across seasons and years. Interactions between climate and disturbance
were observed in a subset of years with insect defoliation increasing forest growth sensitivity to
water availability. Forest growth was particularly sensitive to climate extremes during periods of

8

high stem density following major regeneration events when average inter-tree competition was
high. Results suggest the resistance and resilience of forest growth to climate extremes can be
increased through management steps such as thinning to reduce competition during early stages of
stand development and small-group selection harvests to maintain forest structures characteristic
of older, mature stands.

2.2

Introduction

Understanding the effects of climate on forest productivity is integral to predicting the response
of forest ecosystems to global climate change. Changing climatic conditions have important
implications for sustainable forest management, a fundamental goal of which is to maintain healthy
and productive forests in perpetuity. Projected consequences of climate change, in addition to
global warming trends, include changes in the frequency, severity, and duration of extreme climate
or weather events (IPCC, 2013). These extreme events have the potential to profoundly alter the
productivity and health of forest ecosystems (Allen et al., 2010, 2015).
The increased frequency of droughts combined with warmer temperatures in the southwestern
US, for example, is projected to lead to growth declines in dominant coniferous species and in severe
cases large-scale forest mortality (Breshears et al., 2005; Williams et al., 2010). Other analyses
predict similar growth declines and mortality in mixed oak and pine forests of the southeastern
US and elsewhere around the globe (Klos et al., 2009; Berdanier and Clark, 2016). The impact
of droughts on forest health and productivity is compounded with potential interactions between
drought and other abiotic and biotic disturbances such as wildfire and forest damaging insects (Dale
et al., 2001). In particular, the occurrence of droughts may facilitate increased insect populations,
or insect damage may exacerbate the effects of drought leading to more severe forest mortality
(McDowell et al., 2008; Anderegg et al., 2015). Forest responses to interactions between climate
extremes and disturbance are expected to be complex and exhibit non-linear behavior involving
lags and threshold effects (Betancourt et al., 2004; Williams et al., 2010; Macalady and Bugmann,
2014).

9

Forest sensitivity to climate extremes has been shown to vary depending on endogenous forest stand characteristics such as stem density, age, species composition, and developmental or
successional stage (Laurent et al., 2003; Klos et al., 2009; D’Amato et al., 2013). Changes in
these characteristics are driven by forest dynamics (e.g., Oliver and Larson, 1996), altering forest
responses to climate extremes over time and space. Further, stand characteristics can be modified through forest management to increase forest resistance and resilience to changing climatic
conditions in the short term and facilitate forest adaptation to climate change in the long term
(Millar et al., 2007; Puettmann, 2011). Effective forest management in the face of global climate
change requires understanding and predicting changes in the complex interactions between climate
extremes, disturbance, and forest dynamics (Dale et al., 2001). New analytical approaches capable
of dealing with non-linear forest responses to climate that change over time and space are needed
to inform adaptive forest management (Betancourt et al., 2004).
Tree rings are valuable data for understanding the effects of climate and disturbance on forest
productivity. In particular, tree rings can be used to infer relationships between inter-annual
climate variability and tree growth and to identify past extreme climate/weather and disturbance
events (Cook, 1987; Cook and Kairiukstis, 1990; Babst et al., 2014). In the current study, we are
interested in the combined effects of past climate extremes, disturbance, and forest dynamics as
observed in tree rings. Our goal was to develop a statistical model that allows for inference regarding
if and how disturbance and stand characteristics driven by forest dynamics modify forest growth
responses to climate extremes. Following from our goal, we developed a hierarchical Bayesian state
space model that allows the effects of climate variables on radial tree growth to vary over time. We
hypothesized tree growth responses to climate would change over time in relation to the occurrence
of climate extremes, disturbance, and variation in stand characteristics. The model was motivated
by the need to identify forest conditions that promote resistance and resilience to climate extremes,
which can be used to inform forest management to minimize the impact of future extreme events
on forest productivity.
The hypothesis that forest growth responses to climate vary over time is not new in forest

10

ecology. Indeed, a number of dendrochronology methods exist to estimate time-varying climate
response functions, most notably, applications of the Kalman filter to tree-ring analysis (Visser,
1986; Visser and Molenaar, 1988; Van Deusen, 1989), and moving correlation analysis (Biondi,
1997, 2000; Carrer and Urbinati, 2006). The Kalman filter approach, in particular, has been applied
to identify changes in the effects of climate on tree growth attributable to air pollution (Innes and
Cook, 1989), and to interactions with forest dynamics (Van Deusen, 1987). The model developed
herein nests the Kalman filter within a Bayesian hierarchical state space framework with several
important model properties. First, the hierarchical model structure allows for changes in climate
effects over time to be explicitly modeled as a function of past disturbance and forest dynamics
(assuming suitable data are available). Secondly, the model can accommodate non-linear functions
to model changing climate effects through time as well as non-normal error structures. Third, the
hierarchical Bayesian approach provides explicit error quantification from posterior distributions
as well as tractable error propagation across model components. We applied the state space
model to a dendrochronology data set from northeastern Minnesota to demonstrate its potential
to identify changes in forest growth responses to climate driven by interactions between climate
extremes, disturbance, and forest dynamics. The data set includes radial growth measurements from
individual trees located in 35 forest stands of varying age and size structure, species composition,
and development stage.

2.3

Modeling Approach

We apply two models to estimate tree growth as a function of climate. The first model moves
the analytical steps involved in a dendrochronology response function analysis—the objective
of which is to estimate the average effects of a set of climate variables on annual tree growth
over a fixed period—into an integrated hierarchical Bayesian model allowing for explicit error
propagation across steps (Schofield et al., 2016). The second model extends the first, allowing
climate effects on tree growth to vary annually over the study period using a hierarchical Bayesian
state space approach. Throughout, we refer to the first model as the fixed climate effects model,

11

and the second model as the variable climate effects model. Sources of uncertainty under the fixed
climate effects model include steps to detrend/standardize individual tree growth records, generate
composite growth records or chronologies, and random variability in individual and composite
growth records. The variable climate effects model shares the same sources of uncertainty as the
fixed climate effects model with additional uncertainty associated with the evolution of climate
effects through time. An overview of the fixed and variable climate effects models are included in
the following two subsections. Figure 2.1 provides a schematic representation of the two models.
The statistical details and discussion of Bayesian inference applied to estimate model parameters
are provided in Appendix A.

2.3.1

Fixed Climate Effects

We begin by decomposing individual tree growth measured in terms of radial growth increment
(hereafter growth increment) into component sources of variability to explain the different submodels included in the fixed climate effects (FCE) model. Specifically,
Individual Tree Growth = Long-term Trend + Inter-annual Variability + Random Error
where the long-term trend captures low-frequency changes in growth due to tree size and age, the
inter-annual variability captures high-frequency changes in growth potentially driven by climate,
and the random error term captures residual variability.
Long-Term Trend: A smoothing spline is commonly applied in dendrochronology to model
the long-term, low-frequency trend in individual growth attributable to tree size and age (Cook and
Peters, 1981). Smoothing splines are a highly flexible method to model natural phenomena using a
set of polynomial basis functions (Wood, 2006). We apply a penalized spline regression submodel
in the FCE model to capture long-term size and age effects, although a variety of smoothing
approaches can be used here (see Appendix A for a complete discussion).
Inter-Annual Variability: High-frequency, inter-annual variability in tree growth records is
frequently attributed to the effects of climate and forest dynamics (Cook, 1987). We model mean
12

inter-annual growth variability with an additive annual stand effect that represents the average
deviation of each tree within a stand from the tree’s long-term growth trend within a given year
(i.e., the stand effect is equivalent to the mean inter-annual variability across all trees in a stand
after controlling for tree size and age). The additive stand effect is modeled as a function of
observed climate. We model the mean stand-level variation rather than tree-level variation because
individual trees within a stand are likely to exhibit differential responses to climate variability and
forest dynamics.
Random Error: The residual error is modeled using an autoregressive process to explicitly
account for temporal autocorrelation in annual tree growth increments. Specifically, we apply a
first order autoregressive (AR1) model for the residual error. We note the residual error can also
be modeled as a function of the mean growth increment to account for heteroscedasticity. We did
not choose to do so here.
Log Transformation: We model growth increments on the log scale to ensure positive growth
estimates consistent with Clark et al. (2007). Log transforming growth increments results in
multiplicative errors equivalent to modeling non-transformed growth increments using a negative
exponential model (Schofield et al., 2016). Log transforming growth increments also reduces the
heteroscedasticity frequently observed in tree-ring records.
Combining the submodels, we model annual tree growth increments in a hierarchical Bayesian
model as follows. Let i index individual trees (i = 1, . . . , n), j index stands ( j = 1, . . . , k), and
t index years (t = 1, . . . , T) where n, k, and T are the total number of trees, stands, and years,
respectively. Let j(i) indicate the stand j in which the ith tree is located (e.g., j(i) = 3 indicates the
ith tree is located in stand 3). Finally, define y to be the observed growth increment, such that yi,t
is the observed growth increment for tree i in year t. Individual tree growth is modeled as,
log(yi,t )
| {z }

log-transformed
growth increment

= x0i,t βi +
|{z}
long-term
trend

Îą j(i),t
|{z}

inter-annual
variability

+ i,t
|{z}

(2.1)

random
error

where xi,t includes tree age covariates, βi is a set of tree-specific regression coefficients, ι j(i),t is
the additive effect of being located in stand j during year t, and i,t is the residual error modeled
13



2 /(1−φ2 ) ,
as an AR1 process. We assume the error term is normally distributed, i,t ∟ N 0, σpe
2 is the pure error variance and φ is the temporal autocorrelation coefficient. Stand
where σpe

effects reflecting mean inter-annual growth variability across all trees in a stand are modeled using
observed climate as,
ι j,t = f 0j,t θ + v j,t

(2.2)

where f j,t includes observed, standardized climate covariates, θ is a set of stand-level regression
coefficients, and v j,t is a random error term assumed to be independent with respect to time both


within and across stands, v j,t ∼ N 0, τ 2 .

2.3.2

Variable Climate Effects

The variable climate effects (VCE) model extends the FCE model to allow climate regression
coefficients (θ in Equation 2.2) to vary over time. The climate regression coefficients are treated as
state variables in the VCE model and evolve over time such that a unique set of climate coefficients
is estimated for each year in the study period (θ 1, θ 2, . . . , θ T ). Annual climate coefficient estimates
are updated using the Kalman filter and are informed by coefficient values for the previous and
subsequent years (t − 1, t + 1) and annual stand effect estimates (Figure A.1). Annual climate
coefficients are estimated using stand effects for a five-year period centered on the current year
(t − 2 through t + 2). The use of a five-year moving window allows for partial temporal pooling
of tree growth data similar to a moving correlation analysis (Biondi, 1997), and provides increased
sample size to estimate annual climate effects.
The tree-level model (Equation 2.1) is unchanged in the VCE model. The stand-level model
(Equation 2.2) is updated to integrate time-varying climate coefficients.
ι j,t = f 0j,t θ t + v j,t

(2.3)

The evolution of climate coefficients over time is modeled using a random walk.
θ t = θ t−1 + wt
14

(2.4)

Parameters

Variable Climate
Effects

θ, τ 2

θ 0:T ,
τ 2 , η, Σ(η)

βi

Stand
Effects

Sub-Models

Model
Components

Fixed Climate
Effects

0
Îąj,t = fj,t
θ + vj,t

Long-Term
Trend
x0i,t β i

Data

Stand
Effects
0
Îąj,t = fj,t
θ t + vj,t
θ t = θ t−1 + wt

Inter-Annual
Variability
Îąj,t

2
φ, σpe

Auto-Regressive
Process
i,t = φi,t−1 + ei,t

Residual
Error
i,t

Individual Tree
Growth
yi,t

Figure 2.1: Schematic of fixed climate effects and variable climate effects tree growth models.
Note: Fixed and variable climate effects sub-models are applied separately; 0:T subscript indicates
all time points from 0 to T.

We assume both the stand effect error and random walk error terms follow normal distributions,
v j,t ∼ N (0, τ 2 ), and wt ∼ N (0, Σ θ ). Equations 2.3 and 2.4 define a state space or dynamic linear
model framework with Equation 2.3 serving as the observation equation and Equation 2.4 the
process or state equation (West and Harrison, 1997). Details on the state space modeling approach
and the numerical methods used to estimate model parameters including the application of the
Kalman filter are provided in Appendix A.
Criteria for evaluating variable climate effects: Time-varying climate coefficients can be
difficult to interpret. Understanding the factors that contribute to a significant forest growth
response to a climate variable in a given year is especially challenging, particularly when multiple
climate variables are considered. We applied the following criteria to evaluate time-varying climate
15

coefficients and facilitate interpretation of VCE model results. A climate variable was considered
to have no effect on mean annual growth at the stand scale in a given year if the 95 percent credible
interval for the climate variable coefficient included zero. A climate variable was considered to
have a weak effect in a given year if the 95 percent credible interval for its coefficient did not
include zero, but the climate regression model (i.e., f 0j,t θ t in Equation 2.3) explained less than 25
percent of the variability in the five years of annual stand effects centered on the current year (i.e.,
annual r 2 < 0.25). Mean annual growth at the stand scale was considered sensitive to a climate
variable in a given year if the 95 percent credible interval for its coefficient did not include zero,
and the climate regression model explained at least 25 percent of the variability in the five years of
annual stand effects centered on the current year (in the remainder we use forest growth response
to describe variable by year combinations for which forest growth sensitivity was observed).
We partitioned years during which there was evidence of growth sensitivity for one or more
climate variables into four categories. First, we defined climatic thresholds as the upper quantiles of
observed mean annual climate variable values across stands (quantile values used varied depending
on the climate variable to ensure identification of at least one year with a growth response; Figure
2.5). Growth sensitivity to a climate variable in a given year was attributed to a threshold exceedance
if the threshold for a variable was exceeded within the five-year moving window of annual stand
effects used to estimate climate coefficients. Growth sensitivity to a climate variable was attributed
to a persistent exceedance effect if there were continued years of growth sensitivity following the
five-year period centered on a climatic threshold exceedance. Growth sensitivity was attributed
to interactions with disturbance if years of growth sensitivity coincided with forest tent caterpillar
outbreak years for known host species. Finally, growth sensitivity to a climate variable in a given
year was attributed to unknown sources if it did not meet any of the previous criteria.

16

2.4
2.4.1

Data
Tree Growth Data

We applied the FCE and VCE models to previously-published tree growth data collected in 2010
from 35 forest stands in and around Superior National Forest in northeastern Minnesota (Figure 2.2;
Foster et al., 2014, 2016). The current analysis differs from previous applications of the northeastern
Minnesota tree growth data set in that the VCE model estimates time-varying climate coefficients
based on stand-level growth to understand the effects of past forest disturbance and dynamics on
growth responses to climate extremes. The study region has a continental climate defined by cold
winters (mean January temperature -15 °C) and short summers (mean July temperature 19 °C).
Mean annual precipitation is 600 to 800 mm with much of the total precipitation falling as snow.
Stands were selected to represent the predominant forest communities in the broader geographical
area based on National Forest Inventory and Analysis data from 2004 to 2008 and included
a mixture of species compositions, age structures, and development stages. The study region
spans the temperate-boreal forest ecotone and sampled forest types reflected this biogeographic
setting ranging from common temperate forest types, such as Acer saccharum-dominated northern
hardwood forests, to boreal forests dominated by Pinus banksiana, Populus tremuloides, and Picea
spp. Three replicate 400-m2 circular plots were established within each stand. Increment cores
were collected at breast height (1.3 m) from all live trees with a diameter at breast height (DBH)
larger than 10 cm. Increment cores were measured using a Velmex measuring stage and crossdated
according to standard dendrochronological techniques (Holmes, 1983; Yamaguchi, 1991). The
most-recent year in which growth data was available for all study trees was 2007.
The DBH and species of sample trees were recorded and tree locations were mapped (relative
to plot center coordinates). Tree age was estimated by defining the pith as the recruitment year
and counting the number of growth rings from recruitment to present. Data suitable for climate
modeling exist for 2,291 trees representing 15 unique species located across 105 forest plots (Table
2.1). Sampled trees are assumed to have established as a new cohort following a large disturbance

17

49°0'0"N

¸

93°0'0"W

92°0'0"W

Longitude

91°0'0"W

90°0'0"W
Canada

Superior Nat'nl Forest

Minnesota

Wisconsin

Michigan

!

Latitude

!
!
!
! ! !

!

48°0'0"N

!

!
!

!
!

0

25

!
!

!!
!
!

!

47°0'0"N

!

!

!

!
!

!
!
!
!
!
!

!

!

!

!

Lake Superior

50 Miles

Figure 2.2: Location of 105 forest plots (35 stands, 3 plots per stand) in relation to Superior National
Forest in northeastern Minnesota, USA.
event (e.g., timber harvest or fire). The year of new cohort establishment (or initiation) was set
equal to the 25th percentile of the tree recruitment year distribution for each stand to account for
the presence of trees from older cohorts (additional details provided in Foster et al., 2014). The
start of the study period was set to 1897 to be consistent with the earliest available climate data,
although the growth records for a subset of trees date back before 1897 (study period = 1897 to
2007). As is the case with nearly all dendrochronology data sets, the number of trees and stands
observed each year increases with time reflecting trees that established between the start and end
of the study period (Figure 2.4).

18

Table 2.1: Summary of species sampled.
Species
Code
Abies balsamea
ABBA
Acer rubrum
ACRU
Acer saccharum
ACSA
Betula papyrifera
BEPA
Fraxinux nigra
FRNI
Larix laricina
LALA
Picea glauca
PIGL
Picea mariana
PIMA
Pinus banksiana
PIBA
Pinus resinosa
PIRE
Pinus strobus
PIST
Populus grandidentata POGR
Populus tremuloides
POTR
Quercus rubra
QURU
Thuja occidentalis
THOC
Species

Number of
Trees
365
87
175
273
132
10
96
400
383
33
56
23
93
118
47

Number of
Plots
59
30
16
50
9
6
30
36
23
9
15
4
25
11
12

First Year
Observed
1897
1899
1899
1897
1897
1903
1929
1897
1919
1903
1898
1927
1925
1919
1897

Mean % Relative
Basal Area
9 (0, 62)
3 (0, 23)
9 (0, 97)
13 (0, 70)
7 (0, 97)
1 (0, 16)
5 (0, 50)
13 (0, 98)
13 (0, 93)
3 (0, 47)
5 (0, 77)
2 (0, 33)
6 (0, 66)
7 (0, 91)
6 (0, 68)

Notes: Number of plots is the number of plots in which each species is found out of a
maximum of 105 plots. First year observed is the first year in the study period a growth
record exists for a tree of the corresponding species (several trees have records that date
back prior to 1897). Mean percent relative basal area is calculated across all study stands
based on tree diameters in 2007; the minimum and maximum percent relative basal areas
across study stands are provided in parentheses.
2.4.2

Forest Tent Caterpillar

The forest tent caterpillar (Malacosoma disstria) is an important native defoliating insect in eastern
North America. There have been several forest tent caterpillar (FTC) outbreaks in the study
region between 1987 and 2007. Most notably, FTC outbreaks resulted in significant defoliation
of susceptible trees during the following periods: 1951-1959; 1964-1972; 1989-1995; 2000-2006
(Reinikainen et al., 2012). FTC defoliation in 2001 was particularly severe with greater than 7.5
million acres of susceptible hardwood forests in the state suffering defoliation (Albers et al., 2014).
Study species that are known FTC hosts include: Acer saccharum, Betula papyrifera, Populus
grandidentata, P. tremuloides, and Quercus rubra.

19

Table 2.2: Summary of seasonal aggregations of climate variables. Bullets indicate that a seasonal
aggregation is calculated for a given variable.
Variable
Mean Tmin
Mean Tmean
Mean Tmax
Total AET
Total PET
Total DEF
Mean SNOW

Fall
(Sep - Nov)t−1
•
•
•
•
•
•

Winter
(Dec - Feb)t
•
•
•

•

Spring
(Mar - May)t
•
•
•
•
•
•

Summer
(Jun - Aug)t
•
•
•
•
•
•

Summer Lag
(Jun - Aug)t−1
•
•
•
•
•
•

Notes: Tmin, Tmean, Tmax indicate minimum, mean, and maximum temperature, AET and
PET indicate actual and potential evapotranspiration, DEF indicates climatic water deficit (PET
- AET), SNOW indicates snow pack. Subscripts, t = year of growth, t − 1 = year preceding
growth.
2.4.3

Climate Data

Mean monthly temperature and precipitation estimates were obtained for the study period at a 4-km
resolution from PRISM (PRISM Climate Group, 2013). Climate data were assigned to individual
stands by intersecting plot centroids with the PRISM grid and averaging climate observations
across the three stand plots if they fell within different grid cells (occurs for 3 out of 35 stands).
A number of studies have shown that temperature and precipitation are poorly correlated with
plant distribution and growth in comparison to water balance metrics that translate raw climate
observations into variables with direct physiological relevance to plant function (Stephenson, 1998;
Dyer, 2004; Lutz et al., 2010). We derived monthly values of potential evapotranspiration (PET),
actual evapotranspiration (AET), climatic water deficit (DEF: PET - AET), and mean snow pack
using a modified Thornthwaite-type water balance model (Lutz et al., 2010). We calculated
seasonal aggregations for each variable where relevant as detailed in Table 2.2 for a total of 28
climate variables.
Bayesian Variable Selection: We applied the Bayesian Lasso to select a reduced set of climate
variables with the greatest effect on annual tree growth (Park and Casella, 2008). While not a
formal model-based variable selection technique, the Bayesian Lasso shrinks the coefficient values

20

for unimportant variables to zero in a regression model (Hooten and Hobbs, 2015). We chose to
apply the Bayesian Lasso given its ability to accommodate collinear variables since many of our
climate variables were correlated. Additional details on the Bayesian Lasso and its implementation
are provided in Appendix A. Applying the Bayesian Lasso, the 28 climate variables were pared
down to a final set of five climate variables: fall deficit (FAL-DEF), spring deficit (SPR-DEF),
summer deficit (SUM-DEF), summer deficit in the previous growing season (SUM-DEF-LAG),
and mean annual snow pack (SNOW). All model results are restricted to this final set of climate
variables. Intuitively, we expected water deficit variables to have a negative effect on tree growth.
Snow pack, reflecting spring soil water recharge given the parametrization of the water balance
model, was expected to have a positive effect on tree growth, though reduced growing season length
due to prolonged snow cover could negatively affect growth.

2.5
2.5.1

Results
Fixed Climate Effects (FCE) Model

Mean annual stand-level tree growth was sensitive to all five water balance variables in the model as
indicated by the 95 percent credible intervals not overlapping with zero (Figure 2.3). Specifically,
the four variables representing seasonal climatic water deficit (FAL-DEF, SPR-DEF, SUM-DEF,
SUM-DEF-LAG), where larger values indicate greater water deficit, were negatively related to
mean annual growth at the stand level, while snow pack was positively related to mean annual
growth. Climatic deficit in the fall, summer, and summer before the year of growth were most
related to mean annual growth at the stand level (credible interval bounds are farthest from zero).
The climate coefficient estimates in the FCE model represent the average effects of the five climate
variables over the study period.
Posterior variance estimates for the FCE model are provided in Table A.1. Notably, the firstorder autocorrelation coefficient was roughly 0.37 indicating that tree-ring records were moderately
autocorrelated even after detrending. The individual-tree variance was approximately 0.29, while
the inter-annual variance (capturing stand-level variability) was approximately 0.05, both on the log
21

FAL−DEF

SPR−DEF

SUM−DEF

SUM−DEF−LAG

SNOW

−0.025

0.000

0.025

Std. Coefficient Value

Figure 2.3: Standardized coefficient values for fixed climate effects model. Note: Points represent
posterior median coefficient estimates; black lines indicate 95 percent credible interval bounds; a
dashed grey line at a coefficient value of zero is provided for reference. Variable abbreviations are
as follows: fall deficit (FAL-DEF), spring deficit (SPR-DEF), summer deficit (SUM-DEF), summer
deficit in the previous growing season (SUM-DEF-LAG), and mean annual snow pack (SNOW).
scale. The individual tree growth variance was roughly six times the stand-level growth variance
on the log scale.

2.5.2

Variable Climate Effects (VCE) Model

Estimates of annual effects for each of the five water balance variables were obtained applying the
VCE model. The evolution of each variable over the study period is presented in Figure 2.4. While
the FCE model demonstrated that mean annual growth at the stand level was sensitive to all five
water balance variables, there is evidence under the VCE model that the sensitivity of mean annual
growth to each variable changed in strength and, in some cases, direction over the study period
following the sensitivity criteria described in Section 2.3.2 (Figure 2.4). Table 2.3 summarizes
the results shown in Figure 2.4 partitioning sensitive years for each climate variable into different
response categories. The most common source of growth sensitivity to any climate variable was a

22

threshold exceedance, either during the exceedance year (∟41 percent of all sensitive growth years)
or in the years following an exceedance, i.e., a persistent response (∟13 percent of all sensitive
growth years; Figures 2.4 and 2.5). Growth sensitivity to one or more climate variables coincided
with forest tent caterpillar defoliation of host species in 20 percent of all sensitive growth years.
There is evidence that trees driving large stand-level growth decreases (indicative of defoliation)
during periods of regional forest tent caterpillar defoliation and growth sensitivity to climate were
forest tent caterpillar hosts (Figure 2.6). Finally, growth sensitivity to one or more climate variables
was due to unknown sources (i.e., could not be attributed to a climatic threshold exceedance or
forest tent caterpillar defoliation event) for roughly 26 percent of all sensitive growth years.
Table 2.3: Summary of tree growth sensitivity to climate variables (reference for Figure 2.4).
Variable
SPR-DEF
SUM-DEF
SUM-DEF-LAG

Threshold
Exceedance
1934-1937
1950-1954
1908-1912
1934-1938
1963
1933
1937-1939
1963

Persistent Response
Disturbance
Threshold Exceedance

Other

NA

NA

NA

1913
1939-1943

1954
1991-1993

1902
1907

1940-1943

1953-1954
1991-1993

1946-1947
1950, 1975

FAL-DEF

1975-1979

NA

NA

SNOW

1975-1977

NA

Percent of
Growth Responses

1950-1954
1992-1993

1901-1902
1908-1913
1940-1941
1947
1909, 1933
1947

41.25

12.5

20

26.25

Notes: Variable abbreviations are as follows: fall deficit (FAL-DEF), spring deficit (SPRDEF), summer deficit (SUM-DEF), summer deficit in the previous growing season (SUMDEF-LAG), and mean annual snow pack (SNOW).
Posterior variance estimates for the VCE model are provided in Table A.1. The variance
estimates were consistent with the FCE model, except inter-annual stand-level variance, which
was slightly smaller under the VCE model (0.042 vs 0.05). As in the FCE model, individual-tree
23

Number of Stands

35.0

SPR−DEF

0.1
0.0
−0.1
−0.2

30.0
25.0
20.0
15.0
10.0
5.0

Zero Response
Weak Response
Threshold Response
Disturbance Response
Persistent Response
Unknown Response

SUM−DEF−LAG

SUM−DEF

0.1
0.0
−0.1
−0.2
−0.3

0.0
−0.1
−0.2

0.2

SNOW

FAL−DEF

0.0

0.1

−0.5

0.1
0.0

−1.0
−0.1
1897 1907 1917 1927 1937 1947 1957 1967 1977 1987 1997 2007

1897 1907 1917 1927 1937 1947 1957 1967 1977 1987 1997 2007

Year

Figure 2.4: Evolution of climate coefficient values for each climate variable over the study period
(1897-2007). Solid black line and points indicate posterior mean coefficient values. Dashed lines
delineate 95 percent credible intervals. Points are colored to indicate different response categories.
Zero Response: credible interval includes zero; Weak Response: credible interval does not contain
zero, but annual r 2 < 0.25; Threshold Response: strong response to climate (credible interval
does not contain zero and annual r 2 > 0.25) within two years of threshold exceedance; Persistent
Response: strong response to climate in years immediately following a threshold exceedance;
Disturbance Response: strong response to climate in years where forest tent caterpillar is present
in study region; Unknown Response: strong response to climate not attributable to threshold
exceedance, persistent response, or disturbance. The upper right panel indicates the number of
study stands from which tree growth data exists in relation to the study period. The grey shading
in each panel indicates the period during which tree growth data from fewer than 20 study stand
are available (1897-1929). Variable abbreviations are as follows: fall deficit (FAL-DEF), spring
deficit (SPR-DEF), summer deficit (SUM-DEF), summer deficit in the previous growing season
(SUM-DEF-LAG), and mean annual snow pack (SNOW).

24

SPR DEF [mm]
SUM DEF [mm]
SUM DEF LAG [mm]
FAL DEF [mm]

25

Exceedance − No Response

Exceedance − Response

Disturbance Response

20
15
10
5
0

40
30
20
10
0

40
30
20
10
0
30
20
10
0

SNOW [cm]

20
15
10
5
1897

1907

1917

1927

1937

1947

1957

1967

1977

1987

1997

2007

Year

Figure 2.5: Observed climate variable values over the study period (1897-2007). Black lines
indicate mean climate variable values across study stands along with uncertainty levels equal to
two times the standard error (grey shading). The horizontal red line indicates the estimated climate
threshold for each variable (thresholds correspond to the following quantiles, 0.95: summer deficit,
lagged summer deficit; 0.98: fall deficit; 0.85: spring deficit, snow). White filled points indicate
threshold exceedances with no growth response. Red filled points indicate threshold exceedances
for which a growth response was observed in the five-year period centered on the year of exceedance.
Orange filled points indicate a disturbance response. Variable abbreviations are as follows: fall
deficit (FAL-DEF), spring deficit (SPR-DEF), summer deficit (SUM-DEF), summer deficit in the
previous growing season (SUM-DEF-LAG), and mean annual snow pack (SNOW).

25

a)
Host

Non−Host

b)
Host

Non−Host

−5.0

−2.5

0.0

Partial Residuals

Figure 2.6: Violin plots of partial residuals (observed log annual growth increment minus splinebased estimate) for forest tent caterpillar (FTC) host and non-host individuals located in stands
in the 5th percentile for growth in years a) 1950-1954; b) 1991-1993. Stands in the lowest 5th
percentile for growth are considered likely to have been affected by FTC. Note: Black shading
indicates FTC host trees; grey shading indicates non-FTC host trees.
growth variance was roughly six times the stand-level growth variance on the log scale. The large
tree-level variance relative to stand-level variance in both the FCE and VCE models is consistent
with previous analyses (Foster et al., 2016).

2.6

Discussion

Climate change and associated extreme drought events are expected to fundamentally alter the
structure and functioning of forest ecosystems across wide portions of the globe (Clark et al., 2016).
The localized impacts of these events on forest processes, such as productivity, are likely to vary as
a function of tree- and stand-level characteristics including species, size, age, and density leading to
differential effects across a landscape and over time. Most approaches to modeling climate effects

26

on forest growth have focused on “average species” responses limiting our understanding of how
differences in forest conditions may affect the severity of climate impacts. This study presents
a modeling framework that underscores the importance of forest dynamics in predicting forest
growth responses to climate extremes and disturbance, and highlights the potential for management
regimes focused on manipulating stand structure and density to increase resistance and resilience to
future climate change. The current analysis focuses on the interactive effects of climate extremes,
disturbance, and forest dynamics on growth in mesic forests in northeastern Minnesota. We note
the modeling approach developed herein may prove even more useful for understanding forest
growth responses to drought and disturbance in drier ecosystems such as in the southwestern US.
We begin this section with a discussion of the FCE and VCE model results. We then discuss the
broader implications of the VCE model results to advance understanding of forest growth responses
to climate extremes and potential forest management applications.

2.6.1

Fixed Climate Effects

We defined a set of 28 climate variables that may affect tree growth in northeastern Minnesota
indicative of temperature, precipitation, evaporative water demand, and climatic water deficit
(Table 2.2). The Bayesian Lasso provides an objective method to identify the subset of variables
to which tree growth is most sensitive by shrinking the coefficient values of unimportant climate
variables to zero (Hooten and Hobbs, 2015, see Figure A.2). The value of such a tool in analyses
of the effects of climate on ecological processes is great and has not been previously employed in
tree-ring analyses. We found water balance variables (climatic water deficit and snow pack) had
the largest impact on inter-annual tree growth in northeastern Minnesota. Results indicate that tree
growth was sensitive to all five water balance variables selected by the Bayesian Lasso over the
study period with climatic water deficit exhibiting negative growth effects regardless of season,
and snow pack (a measure of spring soil water recharge) exhibiting positive growth effects (Figure
2.3). Water availability is important in the study region where summers can be dry and soils are
generally shallow and formed largely from glacial till with poor water retention. The climatic water
27

deficit and snow pack variables reflect the interaction between temperature, precipitation, and soil
water holding capacity. The FCE model results indicate that tree growth in the study region was
more sensitive to these interactions than to raw temperature and precipitation values (Figure A.2)
underscoring the importance of translating climate data into physiologically-relevant variables in
studies of tree growth as noted in previous studies (Stephenson, 1998).

2.6.2

Variable Climate Effects

We developed the VCE model to better understand the role of past disturbance and forest dynamics
in shaping forest growth responses to climate extremes. Application of the VCE model to tree
growth data from northeastern Minnesota indicates tree growth was sensitive to water balance
variables in punctuated intervals of one to several years. We partitioned periods of tree growth
sensitivity into four categories to identify climate extremes and elucidate potential interactions
between climate extremes and past forest tent caterpillar defoliation (as defined in Section 2.3.2).
Climatic Threshold Exceedance: The thresholds set for each climate variable in the model
were used to identify extreme climate values (large climatic water deficits) which might lead to a
forest growth response (Figure 2.5). In particular, we sought to identify a threshold for each climate
variable above which a forest growth response is always observed. Forest growth responses to
climate variable threshold exceedances occurred in several ways. There were strong responses to
singular exceedances, e.g., the observed growth response to summer deficit in 1910, a pronounced
drought year in the region (Clark, 1989), and the response to large fall deficit in 1977 (Figures 2.4
and 2.5). There were also responses to an exceedance that closely followed several exceedances
in a short period, e.g., the observed response to summer deficit exceedances in 1933 and 1936
following an exceedance in 1930 (Figure 2.5, Table 2.3). Finally, there were several instances of a
response to an exceedance of one climate variable when exceedances of multiple climate variables
occurred coincidentally, e.g., the negative response to fall deficit in 1977 and the positive response
to a large snow pack in 1975 coincided with large spring deficits in the same year (Figure 2.5).
Figure 2.5 demonstrates that we cannot establish a threshold for each climate variable, above
28

which, there is a high probability of a growth response. Rather, there were a number of threshold
exceedances for each variable (except fall deficit) with no forest growth response (Figure 2.5). This
may be due to the study period (111 years) being too short to observe sufficient variability in water
balance variables to establish meaningful thresholds. It also suggests the stands in the study area
are relatively resistant to the effects of isolated climate extremes. As highlighted above, responses
to threshold exceedances coincided with multiple exceedances in a short period or exceedances of
multiple variables.
The persistent response category provides a measure of forest growth resiliency to the exceedance of a climate variable threshold. A persistent growth response indicates study stands
were sensitive to a climate variable threshold exceedance for several years following the year of
exceedance. There were only three persistent growth responses observed during the study period
despite ten periods where study stands exhibited sensitivity to the exceedance of a climate variable threshold (Figure 2.4, Table 2.3). Specifically, there were persistent growth responses to two
large summer deficits (1910, 1936) and to lagged summer deficit in 1937. The low proportion of
forest growth responses to climate variable threshold exceedances that led to a persistent response
indicates that study stands are relatively resilient to climate extremes.
Forest Disturbance: Contemporary eco-physiological studies note forest growth responses to
large climatic water deficits, as experienced during a drought, may be modified by biotic disturbance
such as forest tent caterpillar defoliation (McDowell et al., 2008; Anderegg et al., 2015). In the
northeastern Minnesota study stands, 20 percent of forest growth responses coincided with forest
tent caterpillar defoliation within the study region (Table 2.3). Only one period of forest tent
caterpillar defoliation coincided with a period of forest growth sensitivity to a climatic threshold
exceedance: response to spring deficit from 1951 to 1954 (Figure 2.4, Table 2.3). Instead, periods of
forest growth sensitivity to one or more climate variables that coincided with forest tent caterpillar
defoliation of host species in the study region occurred when climate variable values were below
set thresholds suggesting defoliation within study stands increased growth sensitivity to water
availability (Figure 2.5).

29

As noted in Anderegg et al. (2015) there can be complex interactive effects of water stress
and insect defoliation on forest productivity. There were two instances in which a forest growth
response to a climate variable occurred in the opposite direction than expected (i.e., positive effect
of climatic water deficit, negative effect of mean snow pack) during periods of forest tent caterpillar
defoliation of host species: a positive response to lagged summer deficit from 1953 to 1954, and
a negative response to snow pack from 1992 to 1993. The positive response to lagged summer
deficit in 1953 and 1954 may have been caused by the presence of drought-weakened trees which
were subsequently killed or further weakened by caterpillar defoliation creating improved growing
conditions for study trees by reducing competition levels. The snow pack variable is indicative of
spring soil water recharge; large snow pack values, however, may cause shortened growing seasons
by delaying leaf flush. The negative response to snow pack in 1992 and 1993, therefore, may have
been due to delayed leaf flush reducing carbon assimilation prior to forest tent caterpillar defoliation
leading to poor growth years. Mean annual snow pack was particularly high in 1992 exceeding
the set threshold value (Figure 2.5). We categorize the growth response to snow pack in 1992 as
a disturbance response, rather than a threshold response given the response coincided with forest
tent caterpillar defoliation of host species in the study region and was in the opposite direction than
expected.
Forest Dynamics: Stand characteristics driven by forest dynamics including age and size structure, stem density, and species composition are likely to impact forest growth responses to climate
extremes. Stem density, in particular, provides a measure of inter-tree competition for light, water,
nutrients, and growing space. A number of studies have demonstrated the resistance and resilience
of forest growth to drought are sensitive to inter-tree competition levels as measured through standlevel basal area (a measure of stand density). Specifically, there is evidence that reducing basal area
via thinning increases forest growth resistance and resilience to drought (Aussenac and Granier,
1988; Laurent et al., 2003; Klos et al., 2009; Martínez-Vilalta et al., 2012; D’Amato et al., 2013;
Sohn et al., 2016; Bottero et al., 2017; but see Floyd et al., 2009 for counter example). The benefit
of thinning, however, may last only a few years and, in some cases, can cause stands to become

30

more sensitive to drought as they mature given an increased presence of large trees with high water
demand due to large leaf area to sapwood area ratios (McDowell et al., 2006; D’Amato et al., 2013).
One goal of the current analysis was to understand the role of stand characteristics, such as
density, in shaping forest growth responses to climate extremes. Unfortunately, past stand density
cannot be inferred from tree rings alone as no growth records exist from previously deceased
trees (Foster et al., 2014). In the current study, we apply knowledge of the fundamental processes
occurring following initiation of a new cohort of trees within a stand based on general models
of stand development (Oliver and Larson, 1996) to analyze the role of past dynamics on forest
growth responses to climate extremes. Specifically, we considered the cumulative growth response
across all climate variables for the threshold/persistent and unknown response categories as a
function of mean years since initiation excluding the initial study period (1897 to 1929) when
growth data from less than 20 study stands were available (Figure 2.7). Initiation here defines
a large recruitment event leading to the formation of a new cohort of trees within a study stand.
There is evidence of a large increase in the number of threshold/persistent and unknown growth
responses 20 to 45 years following initiation, on average (Figure 2.7). Study stands 20 to 45
years after initiation are likely to have high stem densities corresponding to high levels of inter-tree
competition and density-dependent mortality. The cumulative threshold/persistent growth response
function suggests study stands are more sensitive to climate extremes (as represented by climate
variable threshold exceedances) during periods of high stem density following large regeneration
events when individual trees experience higher levels of competition, on average, than in older,
mature stands. Over 85 percent of all unknown forest growth responses (i.e., responses that did not
coincide with a climate variable threshold exceedance or forest tent caterpillar defoliation event)
occurred within 33 years of initiation, on average (excluding the first 33 years of the study period),
providing evidence that study stands are more sensitive to water availability when they are young
(Figure 2.7).
Management Applications: The results of the VCE model combined with previous studies of
forest resistance and resilience to drought following thinning suggest that forest managers may be

31

Unknown Response

0.75
0.50
0.25
0.00

Cumulative Growth Response

1.00

Threshold/Persistant Response

10

20

30

40

50

60

70

80

Mean Years Since Initiation

Figure 2.7: Cumulative distribution of forest growth responses to all climate variables under the
variable climate effects model as a function of years since initiation (initiation marks establishment
a new cohort of trees within a study stand) for the threshold/persistent and unknown response
categories. Note, the cumulative growth response functions exclude growth responses that occurred
early in the study period when growth data from less than 20 study stands were available (19871929). Grey shading highlights responses 20 to 45 years following initiation when understory stem
density and inter-tree competition are high.
able to reduce forest sensitivity to climatic water deficit by thinning stands during periods of peak
density and inter-tree competition (i.e., during the stem exclusion phase of development in even-aged
stands). This period corresponds to the stage at which thinning treatments are traditionally applied
to increase resource levels for residual trees and mimic density-dependent mortality. Thinning
from below (removing only trees in intermediate or suppressed canopy positions) may limit the
formation of large canopy crowns with high leaf area to sapwood area ratios reducing sensitivity
to climatic water deficit as stands mature; further, thinning from below may minimize levels of
evaporative demand at the forest floor due to the high levels of canopy cover it maintains relative
to other thinning approaches. In uneven-aged stands where intermediate thinning treatments may
not be applicable, forest managers may be able to increase forest growth resistance and resilience
to water deficit by minimizing forest gap sizes through individual tree or small group selection

32

harvests that limit the amount of forest in the stem exclusion phase of development at a given
time and increase the range of tree sizes and spatial diversity present in a given stand (Churchill
et al., 2013). Moreover, the pronounced influence of forest tent caterpillar outbreaks and their
interaction with climatic water deficit on productivity underscores the importance of maintaining
mixed-species stands with a diversity of host and non-host species to minimize the impact of insect
defoliation on future productivity.
The modeling approach developed here relies on the dynamic nature of forests to advance
understanding of the effects of past disturbance and forest dynamics on forest growth responses
to climate extremes. The model results demonstrate the importance of considering the effects of
disturbance and forest dynamics on forest growth responses to climate extremes if the goal is to
maintain forest productivity under changing climatic conditions. In its current form, the state space
framework applies a random walk to model changes in climate variable coefficients over time,
rather than modeling changes as a function of past disturbance and forest dynamics. Future work
will apply the framework to data sets that provide more specific information on disturbance history
and past forest dynamics allowing climate coefficients to be modeled explicitly as a function of
disturbance and stand characteristics including density and stand age. In particular, the framework
will be applied to data sets with sufficient sample sizes of individuals from forest tent caterpillar
host versus non-host species to explicitly model the effects of insect defoliation on forest growth
responses to climate. Despite the data limitations in the current analysis, the hierarchical state
space framework we have developed provides a novel method to analyze dynamic forest responses
to climate extremes driven by forest disturbance and dynamics.

33

CHAPTER 3
BOREAL TREE GROWTH EXHIBITS DECADAL-SCALE ECOSYSTEM MEMORY TO
DROUGHT AND INSECT DEFOLIATION, BUT NO NEGATIVE RESPONSE TO
THEIR INTERACTION

3.1

Abstract

• Interactions between drought and insect defoliation may dramatically alter forest function under
novel climate and disturbance regimes, but remain poorly understood. We empirically tested
two important hypotheses regarding tree responses to drought and insect defoliation: 1) trees
exhibit delayed, persistent, and cumulative growth responses to these disturbances; 2) physiological
feedbacks in tree responses to these disturbances exacerbate their impacts on tree growth. These
hypotheses remain largely untested at a landscape scale, yet are critical to predicting forest function
under novel future conditions given the connection between tree growth and demographic processes
such as mortality and regeneration.
• We developed a Bayesian hierarchical model to quantify the ecological memory of tree growth to
past water deficits and insect defoliation events, derive antecedent variables reflecting the persistent
and cumulative effects of these stressors on current tree growth, and test for their interactive effects.
The model was applied to extensive tree growth, weather, and defoliation survey data from western
and eastern regions of the Canadian boreal forest characterized by contrasting tree compositions,
climates, and insect defoliators.
• Results revealed persistent tree growth responses to past water and defoliation stress lasting
5-7 and 8-10 years, respectively, depending on study region. Regional differences in ecological
memory highlight the role of climate and insect defoliator dynamics in shaping forest responses
to drought and defoliation. Boreal tree growth was negatively related to antecedent water deficit
and insect defoliation (host trees only) variables. However, there was no evidence of negative
interactions between antecedent variables in either region. Rather, a positive interactive effect was
34

observed for non-host trees likely due to reduced water stress following defoliation.
• Synthesis. Tree growth responses to water and insect defoliation stress were found to last up to a
decade in contrasting regions of the boreal forest. Interactions between water and defoliation stress,
however, were not found to exacerbate their impacts on host tree growth even after accounting
for persistent and cumulative effects. This result, consistent with earlier experiments, suggests
negative feedbacks in host tree responses to drought and insect attack may be weaker than predicted
for defoliator-dominated boreal forest systems.

3.2

Introduction

The boreal forest is one of the largest terrestrial ecosystems (∟1135 Mha) and accounts for
roughly 32 percent (272 Pg C) of all forest-stored carbon on earth (Pan et al., 2011). It is an
important natural resource, contributing to national economies in North America, Europe, and
Asia (Brandt, 2009). The boreal forest is a disturbance-driven ecosystem, defined by the interplay
of droughts, insect damage, pathogens/disease, and stand-replacing fire (Fleming et al., 2000;
Girardin et al., 2013). Increased temperatures and aridity forecast for much of the boreal forest
under global climate change have the potential to fundamentally alter historical boreal disturbance
regimes (Price et al., 2013; Allen et al., 2015). In particular, severe drought events characterized
by high temperatures are already occurring in boreal regions of western North America (Hogg
et al., 2002, 2008; Peng et al., 2011; Michaelian et al., 2011). Further, the severity, extent, and
duration of defoliating insect outbreaks have increased in recent decades (Blais, 1983; Roland,
1993; Pureswaran et al., 2015). Changing climatic conditions are likely an important contributor to
changes in defoliator population dynamics (Cooke et al., 2007; Price et al., 2013). Both drought and
insect defoliation can severely impact boreal forest function causing growth reductions and largescale forest mortality with lasting impacts on the global carbon cycle, regional timber supplies, and
an array of other ecosystem services (Kurz et al., 2008; Hicke et al., 2012; Pothier et al., 2012;
Price et al., 2013).

35

A number of studies have demonstrated spatio-temporal synchrony between drought and defoliating insect damage to the boreal forest (Hogg et al., 2008; Flower et al., 2014; DeGrandprĂŠ et al.,
2017). Several physiological feedbacks in tree responses to drought and insect defoliation have been
proposed to explain the synchrony of such events. First, insect defoliation may reduce tree growth
and xylem formation making trees more susceptible to future drought (Anderegg and Callaway,
2012; Jacquet et al., 2014). Second, limited available carbon and water in drought-stressed trees
may reduce metabolic defenses against defoliating insects (McDowell et al., 2008, 2011). Finally,
reductions in carbon uptake due to the combination of drought and insect defoliation may deplete
non-structural carbohydrates and increase the risk of carbon starvation (Hogg et al., 2008).
Despite the observed synchrony of drought and defoliation events and proposed negative physiological feedbacks, the nature of interactive effects of drought and insect defoliation on tree growth
and mortality remain unclear (Kolb et al., 2016). For example, water stress may increase the
concentration of secondary metabolites in leaves making them less palatable to defoliating insects
(Mattson and Haack, 1987). Further, moderate defoliation may decrease the impact of drought
events by reducing tree density and leaf area (Jacquet et al., 2014). Previous studies indicate both
positive and negative impacts of drought on tree resistance to insect defoliation and vice versa
depending on the severity of the drought or defoliation event (Jactel et al., 2012).
Uncertainty in physiological mechanisms underlying feedbacks between drought and insect
defoliation stress highlight the importance of empirical studies to test for the interactive effects of
drought and insect defoliation on tree growth and mortality (Anderegg et al., 2015). Few such
studies exist. A notable exception is Jacquet et al. (2014), which found no significant interaction
between drought and defoliation stress on experimentally-manipulated maritime pine plantations.
Further, among a number of studies assessing aspen mortality in western North America (sudden
aspen decline [SAD]) due to severe drought and repeated forest tent caterpillar defoliation (Hogg
et al., 2002, 2008; Worrall et al., 2010; Michaelian et al., 2011; Anderegg and Callaway, 2012;
Anderegg et al., 2012, 2013), none explicitly tested for interactions between disturbance agents.
Detecting interactive effects of drought and insect defoliation on tree growth and mortality is

36

complicated given physiological responses to drought and insect defoliation may be delayed or
persist following a disturbance event (Pothier et al., 2005, 2012; Hogg et al., 2008; Worrall et al.,
2010; Anderegg et al., 2013). Further, physiological responses to past disturbance events may
weaken a tree’s ability to respond to future disturbances (Anderegg et al., 2013); that is, there is a
cumulative effect of drought and insect defoliation. Potential delayed, persistent, and cumulative
responses make it difficult to model the interactive effects of drought and insect defoliation using
observations of forest conditions from a single time point. Classical approaches to overcome this
issue involve applying lagged variables (e.g., water availability in the year previous to growth) within
a regression model to estimate growth or mortality. Regression models using lagged variables,
however, can be difficult to interpret, require the researcher to select the number of lags to include,
and often suffer from multicollinearity issues (Heaton and Gelfand, 2012). An alternative approach,
based on recent advancements in statistical ecology, is to quantify antecedent drought and insect
defoliation variables using ecological memory functions reflecting the cumulative effects, relative
importance, and strength of both disturbance types on forest growth and mortality over a specified
time period (Ogle et al., 2015). The use of ecological memory functions allows for novel insights
into the nature of forest responses to past drought and insect defoliation and improves our ability to
detect interactions due to physiological feedbacks.
In the current study, we developed a Bayesian hierarchical model to quantify the ecosystem
memory of trees to water deficit and insect defoliation allowing for the detection of delayed,
persistent, and cumulative growth responses to these stressors. Further, the model tested for
interactive effects of antecedent water and insect defoliation stress on tree growth where derived
antecedent variables reflect potential delayed, persistent, and cumulative tree responses. The model
was applied to extensive tree-ring width, interpolated climate, and aerial defoliation datasets from
the western and eastern regions of the Canadian boreal forest characterized by contrasting climates,
species compositions, and primary defoliating insect species. The hierarchical model allowed us
to combine multi-scale data in an integrated statistical framework with propagation of uncertainty
among scales. The use of two contrasting boreal study systems highlights the generalizability

37

of the model and allows for a robust interpretation of interactions between drought and insect
defoliation. Our analysis focused on tree growth responses alone; however, there is evidence of
strong connections between tree growth and vital forest demographic rates including tree mortality
(Wyckoff and Clark, 2000, 2002; van Mantgem et al., 2003; Das et al., 2007, 2016; Berdanier and
Clark, 2016; Buechling et al., 2017).
Based on current physiological understanding, we hypothesized that in both boreal study regions: 1) tree growth would exhibit persistent responses to water deficit with deficits in the year
previous to the current growing season having the greatest effect on growth; 2) tree growth would
exhibit persistent responses to insect defoliation for multiple years following a defoliation event; 3)
antecedent water and insect defoliation stress would have a negative interactive effect on host tree
growth. In addition to testing each of these hypotheses, we further tested for regional differences in
the ecosystem memory of trees to water deficit and insect defoliation and their antecedent effects
on tree growth. Finally, we tested whether the effects of antecedent water and defoliation stress
vary depending on tree size (large vs. small diameter) and species (defoliator host vs. non-host)
categories within each region.

3.3

Methods

We begin this section by providing information on the forest composition, climate, and defoliating insect population dynamics in the two study regions. We then detail the tree-growth, climate,
and defoliation data applied in the analysis. Finally, we define the Bayesian hierarchical model used
to estimate ecological memory and test for interactive effects of water deficit and insect defoliation
on boreal forest growth.

3.3.1
3.3.1.1

Study Regions
West

The western study region is located in Alberta, Canada and consists of 34 mixed-wood boreal
stands in the western, interior plains extending from 52.0 to 59.0°N and 111.0 to 119.5°W (Figure
38

3.1A). Study stands are predominantly composed of trembling aspen (Populus tremuloides Michx.)
and white spruce (Picea glauca Moench) and established post stand-replacing fire. The higher
growth rate of aspen relative to spruce leads to a bimodal diameter distribution early after stand
establishment until spruce moves into the canopy as a co-dominant species. The region is defined
by a continental climate (cold winters, hot and dry summers) with chronic summer water deficits.
Annual mean monthly temperature and total annual precipitation (1970-2000) ranged from -0.8
to 2.8 °C and 370 to 615 mm (25 to 35 percent falling as snow) across the region (Fick and
Hijmans, 2017). Mean summer (June-August) Standardized Precipitation-Evapotranspiration Index
(SPEI) values, a measure of drought severity relative to average conditions with more negative
values indicative of severe droughts, indicate a period of relative wetness beginning in 1970 and
overlapping the defined study period (1968-2010) with dry years in 1982, 1992, and 2004 (Figure
3.1B). Forest soils in the western study region consist of orthic gray luvisols and brunisols with silty
to clay-loam texture derived from glacial till and glaciolacustrine deposits (Huang et al., 2013).
The forest tent caterpillar (Malascosoma disstria Hub., hereafter “FTC”) is the primary defoliating
insect within western study stands attacking the regionally-abundant trembling aspen (Brandt et al.,
2013). During its larval stage, the FTC, a univoltine lepidopter, defoliates aspen shortly after leaf
out for a 5-6 week period until pupations in mid-to-late June (Parry et al., 1998). Outbreaks can
last for one to several years (Price et al., 2013). Trees are able to produce a second flush of leaves
after the FTC has pupated (Cooke et al., 2007). Outbreaks of the FTC have been linked to regional
growth declines in aspen as well as large-scale mortality when outbreaks co-occur with droughts
or extreme temperatures (Hogg et al., 2005).

3.3.1.2

East

The eastern study region is located in Quebec, Canada and comprises 14 coniferous boreal stands
within the Forest Ecosystem Research and Monitoring Network, extending from 47.0 to 50.1°N
and 66.1 to 75.0°W (Figure 3.1A). Stands are dominated by black spruce (Picea mariana Mill.) or
balsam fir (Abies balsamea L.) accompanied by white birch (Betula papyrifera Marsh.), trembling
39

A

N

B

Mean

Loess Smooth Function

West (Alberta)

Min-Max Range

Study Period

East (Quebec)

Figure 3.1: (A) Location of study stands within western (Alberta: 34 stands) and eastern (Quebec:
14 stands) study regions relative to mean summer (June-August) climate moisture index values
(1970-2000, Fick and Hijmans, 2017). (B) Mean summer 3-month Standardized PrecipitationEvapotranspiration Index values for western and eastern study stands (1902-2015, BeguerĂ­a and
Serrano, 2017).
aspen, white spruce, and jack pine (Pinus banksiana Lamb.). Annual mean monthly temperature
(1970-2000) ranged from -0.3 to 3.2 °C while total annual precipitation ranged from 860 to 1600
mm (Fick and Hijmans, 2017). Mean summer SPEI values for eastern study stands indicate a
period of relative wetness beginning in the 1960s and continuing through the study period (19681998) with dry years in 1968 and 1991 (Figure 3.1B). Soils in the eastern study region consist
mainly of orthic ferro-humic podzols and gleysols derived from glacial till as well as glaciofluvial
and glaciolacustrine deposits with soil depths ranging from 30 to 60 cm (Ouimet et al., 2001).
Defoliation in the eastern study region is caused by outbreaks of the eastern spruce budworm
(Choristoneura fumiferana Clem., hereafter “SBW”), which occur every circa 35 years and can
persist 4-12 years or more in a given stand (Gray, 2008). Balsam fir is the primary SBW host in the
East with white and black spruce also susceptible to attack during severe defoliation events (Nealis
and RÊgnière, 2004). The larval feeding of the SBW, also a univoltine lepidopter, is approximately
6 weeks from early-to-mid May until mid-to-late June depending on the region (RÊgnière et al.,
40

2012). Larva feed primarily on current year foliage and require multiple years to kill a host tree
(Gray, 2008). SBW outbreaks, given their extent and severity, are one of the most damaging natural
disturbance types in the Canadian boreal forest (Fleming et al., 2000).

3.3.2

Data

A total of 34 forest stands were sampled in the western study region using variable-length belt
transects (Figure 3.1A). Approximately ten white spruce trees from a range of size classes and
fifteen dominant or co-dominant aspen trees were sampled within each transect. Either two radial
increment cores were collected on opposite sides of a sampled tree or entire cross-sections were
obtained through harvest. All radial growth sampling was done at breast height (1.3 m). Cores
and cross-sections were dried, mounted, and sanded. Ring-widths were measured using a Velmex
measuring state (increment cores) or WinDendro (cross-sections) to the closest 0.01 mm. A total
of 919 tree growth series were used in the analysis: 471 white spruce and 448 trembling aspen.
Fixed area plots (0.25 ha) were established in each of the 14 eastern study stands in 1986-1998
(Figure 3.1A). Between 1996-1998, radial increment cores were collected at breast height from 25
to 50 healthy, co-dominant trees along the border of each plot. Two increment cores were taken from
opposite sides of each tree. Cores were dried, mounted in wooden blocks, and sanded. Ring widths
of all cores were measured to the nearest 0.01 mm under 40X magnification. A total of 625 tree
growth series were used in the analysis, mostly black spruce (342) and balsam fir (208), as well as
jack pine (35), white spruce (30) and white birch (10). Tree-ring cross-dating in both study regions
was done visually and verified statistically with COFECHA (Grissino-Mayer, 2001). Replicate ring
widths resulting from the collection of multiple increment cores per tree were averaged to produce
a single radial growth record for each sampled tree in the East and West.
We separated trees into large- and small-diameter classes using an 11-cm diameter at breast
height (DBH) threshold in both study regions for use in the tree growth model (see Section 3.3.3).
The 11-cm DBH threshold was determined based on exploratory analysis using detailed plot data
to determine a DBH above which a tree has high probability of being in a co-dominant or dominant
41

canopy position. Each tree’s diameter class was updated annually conditional on its tree-ring
reconstructed DBH. The use of an 11-cm diameter threshold led to relatively equal sample sizes of
annual growth increments for large- and small-diameter trees in both study regions.
Monthly weather data (1950-2012) were generated for each study stand using the BioSIM
interpolation model (RÊgnière, 1996) based on a network of 365 weather stations in the East and
456 in the West. We applied a water balance model adapted from Lutz et al. (2010) to estimate
monthly potential evapotranspiration (PET), actual evapotranspiration (AET), and climatic water
deficit (PET-AET). PET was estimated using a modified Hargreaves equation based on monthly
averages of daily minimum/maximum temperature and incoming solar radiation, total monthly
precipitation, and latitude (BeguerĂ­a et al., 2014). Monthly climatic water deficit estimates were
converted to annual values by summing monthly estimates for growing season months (June through
August). Estimates of soil water holding capacity were necessary to derive AET. Mean soil water
holding capacity was estimated for eastern study stands by converting granulometric content of
unique soil horizons applying equations from Saxton and Rawls (2006) with weights corresponding
to horizon thickness. Granulometric content was estimated based on laboratory analysis of four
1-m2 soil pits dug in each study stand. Less detailed soil data were available for western study
stands. Ordinal estimates of soil water holding capacity for western stands were obtained from the
Canadian Soil Information Service (Soil Landscapes of Canada Working Group, 2010) and were
converted to numeric estimates by taking the median value of each ordinal category.
The provinces of Quebec and Alberta began conducting aerial defoliation surveys in 1967
and 1939, respectively, to estimate the extent and severity of damage caused by major forest
pests (Ministère des Forêts, de la Faune et des Parcs, 2016; C. Whitehouse, Alberta Ministry of
Agriculture and Forestry, personal communication, March 21, 2017). Derived maps indicate annual
insect defoliation according to ordinal severity classes. We derived an annual binary defoliation
variable for each study stand based on defoliation maps. A value of 1 was assigned to the defoliation
variable if a mapped defoliation event of moderate (36-70% loss of foliage in the upper half of the
crown of most trees) to high (71-100% loss of foliage over the entire crown length of most trees)

42

severity intersected a study stand in a given year with 0 assigned otherwise. The start of the study
period in the eastern and western regions was set to 1968 to be consistent with the availability of
defoliation data in both study regions.

3.3.3

Model Specification

Our Bayesian hierarchical model consists of tree- and stand-level submodels (Figure 3.2). The
tree-level submodel estimates annual radial growth increment on the log scale as a function of a
time-varying stand effect reflecting stand-level growing conditions, and a tree’s diameter in the
previous year. The stand-level submodel estimates the time-varying stand effect as a function of
antecedent climatic water deficit and insect defoliation derived using ecological memory functions.
Based on studies demonstrating forest growth responses to drought and insect defoliation are
modified by stand structure and composition (Bergeron et al., 1995; Nealis and RÊgnière, 2004;
McDowell et al., 2006; Jactel and Brockerhoff, 2007; D’Amato et al., 2013; Gleason et al., 2017),
we allow the effects of antecedent water deficit to vary by DBH class (large vs. small) and the
effects of antecedent insect defoliation to vary by tree species categories (defoliator host versus
non-host). Primary inferential interest is in the stand-level submodel including the estimation of
ecological memory functions and regression coefficients for derived antecedent variables. The treelevel submodel is included to control for tree size and its effect on annual radial growth increment.
Combining the different model components results in the following regression equation for the
radial growth increment (y) of the ith tree in year t,
log(yi j hs (t)) = ι j hs (t) + xi j hs (t)β + ei j hs (t),
where j indexes the study stand, h indexes whether the tree is in the host or non-host species
category for the dominant regional defoliating insect, and s indexes the DBH class of the tree in
year t (large: > 11 cm; or small: < 11 cm). A time-varying stand effect (Îą j hs (t)) is uniquely
estimated for each tree species (host/non-host) and size (large/small DBH) category combination,
xi j hs (t) is a tree’s diameter in the previous year, β is an estimated regression coefficient, and ei j hs (t)
43

is a tree-level error term. A single linear term for tree diameter in the previous year was sufficient
to control for tree size in the eastern region. Tree growth in the western region, however, differed
by species (aspen versus spruce) and showed evidence of following a sigmoid growth function
with respect to diameter. As such, the tree-level submodel in the West included species-specific
linear and quadratic terms for tree diameter in the year previous to growth (see Web Supplement).
Consistent with a multiple linear regression model, we assume the tree-level errors are normally
iid

distributed, ei j hs (t) ∼ N(0, σy2 ), where σy2 is a tree-level variance and “iid” stands for independent
and identically distributed. The time-varying stand effect is estimated according to a separate
regression equation, the stand-level submodel, given by,
(h)
(s)
(hs)
α j hs (t) = γ0 + z̃ j (t)γ1 + f˜j (t)γ2 + zff j (t)γ3 + u j hs (t),

where z̃ j (t) is the estimated antecedent defoliation value, f˜j (t) is the estimated antecedent climatic
water deficit value, zff j (t) is an antecedent defoliation-water deficit interaction term ( zff j (t) =
z̃ j (t) × f˜j (t)), the γ’s are estimated regression coefficients, and u j hs (t) is a stand-level error term.
iid

We assume the stand-level errors are also normally distributed, u j hs (t) ∟ N(0, σι2 ), where σι2 is an
inter-annual, stand-level variance.
Ecological memory is quantified by estimating an antecedent weight function where antecedent
weights reflect the relative importance of past environmental conditions on current ecosystem
function (Ogle et al., 2015). We build on the Bayesian framework presented in Ogle et al. (2015) to
estimate ecological memory using splines and spatio-temporal covariance functions rather than a
Dirichlet prior assigned to antecedent weights. The modified approach allows for efficient (i.e., few
model parameters) estimation of complex memory functions. Our approach to estimate ecological
memory is different for climatic water deficit and insect defoliation reflecting differences in forest
growth responses to these disturbance types.
We used penalized regression splines to estimate antecedent weights for past climatic water
deficit. Regression splines provide a flexible method to model complex, non-linear ecological
processes based on linear combinations of piecewise polynomials (Wood and Augustin, 2002).
Regression splines are penalized through the use of a smoothing parameter to control against over
44

Tree-Level Submodel
Log Radial
Growth Increment
(t)

Stand Effect
(t)

=

+

Tree Diameter
(t-1)

+

Tree-Level
Error
(t)

Stand-Level Submodel
Host

Large
Diameter

Non-Host

Antecedent
Defoliation
(t)

Small
Diameter

Antecedent
Water Stress
(t)

+

Large-Host

+

Weight

Defoliation
Time Since Defoliation

Defoliation Events
(t0) to (t)

0

Lag

Antecedent
Interaction
(t)

=

Large-Non-Host

Small-Host

+

Stand-Level
Error
(t)

Small-Non-Host

Antecedent
Defoliation
(t)

L

⇥

Climatic Water
Deficit
(t-L) to (t)

Antecedent
Water Stress
(t)

Figure 3.2: Diagram of hierarchical model structure including tree- and stand-level submodels and
the estimation of antecedent water deficit and antecedent insect defoliation variables. Rectangles
delineate observations, while ovals delineate estimated quantities. Note t0 indicates the beginning of
the study period and L indicates the maximum number of years before present for which ecological
memory to past climatic water deficit is estimated.
fitting. Penalized regression splines have a range of applications in ecology (Wood and Augustin,
2002; Hefley et al., 2017). Use of penalized regression splines in the current analysis allows the
antecedent weight function for water deficit to be entirely determined based on the data (i.e., no
pre-specified functional form). Splines were constructed as a function of time before present up to
a maximum range of years L. We set L = 10 years, but other values can be applied depending on
region and knowledge of forest ecosystem responses to water stress. The antecedent weight value
(w) for a given time lag (`) is estimated as,
nÍ

w` =

o

p
exp i=1 hi (`)Ρi
nÍ
o,
ÍL
p
exp
h
(`)Ρ
i
`=0
i=1 i

45

where hi (`) is the value of the ith spline function evaluated at lag ` and Ρi is its corresponding
regression coefficient (i = 1, 2, . . . , p). The exponential term (equivalent to modeling weights on
the log scale) constrains weights to be positive, while the sum over ` in the denominator imposes
ÍL
a sum-to-one constraint among the weights, `=0
w` = 1. The antecedent water deficit value
(reflecting cumulative water stress) for stand j in year t is then estimated as the weighted sum of
past water deficit observations,
f˜j (t) =

L
Õ

w` f j (t − `),

`=0

where f j (t − `) is the observed climatic water deficit in stand j, ` years before year t. Additional
details on the penalized regression spline approach are provided in Appendix B.
We applied a temporal decay function to model ecological memory to past insect defoliation
given peak defoliation stress levels are expected to occur when a defoliating insect population is
present in a stand. The decay function begins in the year of a defoliation event and decays to zero as
a function of years since defoliation. We applied a spherical decay function in the current analysis
because it achieves a value of zero within a finite interval of time (see Appendix B). The function is
conditional on an estimated temporal range parameter (φ), which reflects the length of ecosystem
memory to past defoliation events. Defoliating insect outbreaks can last more than a single year and,
in some cases, an additional defoliation event can occur before host trees have fully recovered from
a previous defoliation event. To account for defoliator dynamics, we estimate antecedent insect
defoliation as the sum of the decay function values for all defoliation events preceding or coincident
with the current year to capture their cumulative effects. Specifically, antecedent defoliation is
given by,
D j (t)

z̃ j (t) =

Õ

g (t − t(di ); φ) ,

i=1

where D j (t) is the total number of defoliation events in stand j from the beginning of the study
period through year t, g(x; φ) is the spherical decay function, t(d) indicates the year of a defoliation
event, and i indexes individual defoliation events. Reference plots for antecedent climatic water
deficit and insect defoliation variables are provided in Figures B.1 through B.4 in Appendix B.

46

The Bayesian hierarchical model is completed by specifying prior distributions for unknown
parameters. These include tree- and stand-level regression coefficients and variance parameters,
the basis function coefficients for the antecedent water stress weights, and the decay parameter for
the antecedent defoliation memory function. We used Markov chain Monte Carlo (MCMC) to
sample from the joint posterior distribution for all model parameters (Robert and Casella, 2004).
Details on the prior distributions used and the MCMC sampler are provided in Appendix B.

3.4
3.4.1

Results
Ecological Memory

Antecedent weights for climatic water deficit ranged from near zero to 0.4 (Figure 3.3A). Climatic
water deficit in the year previous to growth had the largest antecedent weight in each region (` = 1;
posterior mean: East = 0.23 [0.17,0.29]; West = 0.32 [0.25,0.39]). Weights larger than 0.05 were
also estimated for deficits in the year of growth and 3-5 years previous to growth depending on
the region. Trees within eastern study stands exhibited a more prolonged response to past climatic
water deficit than trees in western study stands. Non-zero antecedent weights were estimated for
climatic water deficit up to 7 years prior to growth in the East compared to 5 years in the West.
Further, the estimated memory function for the West exhibited a faster response to climatic water
deficit than in the East with larger antecedent weights for lags 0-2 and smaller weights for lags
3-10. There was only weak evidence, however, of differences in ecological memory to climatic
water deficit between the two study regions (95 percent credible intervals overlap in Figure 3.3A).
Tree growth in both study regions is estimated to exhibit persistent responses to past climatic water
deficit with multi-year water deficits leading to more persistent responses (Figure 3.4A). Further,
multi-year water deficits are estimated to result in greater antecedent effects reflecting a larger
cumulative tree growth response.
Similar to climatic water deficit, there was weak evidence trees within eastern study stands
exhibited more prolonged memory to past insect defoliation events than trees in western stands
(Figure 3.3B). Based on the posterior distribution of the range parameter (φ), trees in eastern study
47

Antecedent Weight

0.4

A
●

0.3

Posterior Median

●

95% Cred. Int.

●
●

●
0.2

Region

●

●

●

0.1

● East

●

● West

●
●
0.0
0

1

2

3

●

●
●

●
●

●

●

●

5

6

7

8

9

10

4

Lag
1.00

●

Antecedent Weight

B

Region
●
●

0.75

● East
●
●

0.50

● West
●
●

●
●

●
●

0.25

●

●
●

●
●

0.00
0

1

2

3

4

5

6

7

●

●
●

●
●

8

9

10 11 12 13 14 15

●
●

●

●

●

●

Years Since Defoliation Observed

Figure 3.3: Ecological memory functions for climatic water deficit (A) and insect defoliation (B)
based on annual radial tree growth in eastern and western study regions. Lag indicates the number
of years before current year growth.
stands exhibited ecosystem memory to SBW defoliation events occurring up to 11.2 (95% credible
interval: 7.8-14.7) years before present, while trees in western study stands exhibited memory
to past FTC defoliation events up to 8.3 (95% credible interval: 6.9-9.6) years before present.
The ecosystem memory to past defoliation is estimated to lead to persistent and cumulative tree
growth responses to insect defoliation events in both study regions (Figure 3.4B). The length of
persistent responses and magnitude of antecedent effects are estimated to increase the longer an
insect defoliation outbreak lasts.

3.4.2

Antecedent Effects

Mean host tree growth was negatively related to antecedent insect defoliation in both study regions
(Îł1 in Table 3.1). Specifically, the average annual growth of host trees, after controlling for tree
size and assuming zero antecedent water deficit, was 0.63 times smaller per one-unit increase in

48

A

B

EAST

Posterior Median

Disturbance
Endpoints

WEST

Antecedent Insect Defoliation

Antecedent Climatic Water Deficit

95% Cred. Int.

Figure 3.4: Estimated antecedent climatic water deficit (A) and insect defoliation (B) values for
simulated water deficits of 10 mm and insect defoliation events lasting one, two, and five years.
Dashed lines highlight the year for which maximum antecedent values are estimated.
antecedent FTC defoliation in the West, and 0.88 times smaller per one-unit increase in antecedent
SBW defoliation in the East (e.g., 1 + [eγ1 − e0 ]). There was no evidence mean annual growth of
non-host trees was related to antecedent defoliation after controlling for tree size in either region
(small posterior median coefficient values and 95 percent credible intervals overlap with zero in
Table 3.1).
Mean annual growth of large- and small-diameter trees was negatively related to antecedent
water deficit in both study regions with large-diameter trees exhibiting greater sensitivity to water
deficit than small-diameter trees (Table 3.1). There was stronger evidence of increased sensitivity
of large-diameter trees to antecedent water deficit in the West where the magnitude of the posterior
median coefficient for large-diameter trees (-0.012) was four times as large as the posterior median
coefficient for small-diameter trees (-0.003), and the 95 percent credible intervals for the size classes
did not overlap. Specifically, after controlling for tree size effects and assuming no antecedent
insect defoliation, mean annual tree growth was 0.84 times smaller for large-diameter trees and
0.90 times smaller for small-diameter trees in the East experiencing a simulated antecedent climatic
water deficit of 5 mm, equivalent to the regional average over the study period. Under the same
conditions in the West, mean annual tree growth was 0.79 times smaller for large-diameter trees
49

Table 3.1: Posterior summary of stand-level, antecedent variable coefficients. Posterior median
coefficient values are provided for antecedent insect defoliation, antecedent climatic water deficit,
and their interaction for each study region by species (defoliator host vs. non-host) and size (largevs. small-diameter) categories. 95 percent credible intervals are given in parentheses. Coefficients
for which credible intervals do not include zero are bolded.
Parameter
Îł1
Îł2

Îł3

Category

East
West
Antecedent Defoliation
Host
-0.128 (-0.187, -0.077) -0.467 (-0.602, -0.335)
Non-Host
0.026 (-0.055, 0.106)
0.083 (-0.03, 0.201)
Antecedent Water Deficit
Large
-0.034 (-0.043, -0.026) -0.012 (-0.014, -0.01)
Small
-0.02 (-0.029, -0.011) -0.003 (-0.004, -0.001)
Antecedent Defoliation × Antecedent Water Deficit
Host|Large
0.005 (-0.003, 0.012)
0.004 (-0.001, 0.008)
Host|Small
0.007 (-0.001, 0.015)
0.001 (-0.003, 0.005)
Non-Host|Large 0.032 (0.012, 0.067)
-0.001 (-0.005, 0.003)
Non-Host|Small 0.005 (-0.014, 0.029)
0.002 (-0.002, 0.005)

and 0.95 times smaller for small-diameter trees experiencing a simulated climatic water deficit of
20 mm, again equal to the regional average over the study period.
In general, mean annual tree growth in stands from both study regions was not strongly related
to the interaction between antecedent insect defoliation and climatic water deficit (Îł3 in Table 3.1).
The one exception was for large-diameter, non-host trees in the East where there was evidence
mean annual growth was positively related to the interaction between the antecedent variables.
Specifically, mean annual growth among large-diameter, non-host trees experiencing 5 mm of
antecedent water deficit was estimated to be 1.17 times greater for an antecedent defoliation value
of one versus zero after controlling for tree size. Figure 3.5 presents the sum of estimated effects
for antecedent insect defoliation, antecedent climatic water deficit, and their interaction relative to
regional mean annual tree growth over the study period after controlling for tree size. Based on
Figure 3.5, the positive interaction between antecedent variables for large-diameter, non-host trees
in the East is estimated to have the greatest impact at moderate antecedent water deficits and high
antecedent defoliation levels. We observed no evidence of such an interaction in the West where
large-diameter, non-host trees were equally sensitive to antecedent water deficit regardless of the

50

Large-Diameter | Host

Large-Diameter | Non-Host

Small-Diameter | Host

Small-Diameter | Non-Host

EAST
WEST

Antecedent Insect Defoliation (unitless)

Relative
Response

Antecedent Climatic Water Deficit (mm)

Figure 3.5: Estimated effects of antecedent insect defoliation, antecedent climatic water deficit,
and their interaction on tree growth in eastern and western study regions by species (defoliator
host vs. non-host) and size (large- vs. small-diameter) categories. Points represent posterior
median antecedent variable values based on study data. Relative response surfaces correspond to
mean tree growth under antecedent conditions relative to regional mean tree growth over the study
period (East: 1968-1998, West: 1968-2010) after controlling for tree size. Response surfaces were
generated by imposing a dense grid over the range of modeled antecedent variable values.
amount of antecedent defoliation. Additional model coefficients of reduced inferential interest,
including the effects of diameter on tree-level growth and tree- and stand-level variances for both
regions are summarized in Table B.1.

3.5
3.5.1

Discussion
Ecological Memory

The estimated ecological memory functions provide evidence of persistent tree growth responses to
climatic water deficit and insect defoliation in both regions of the Canadian boreal forest supporting
hypotheses 1 and 2, respectively (Figure 3.3). The derived antecedent water deficit and defoliation
variables reflect the cumulative effects of these forest stressors over time (Figure 3.4). The estimated

51

ecological memory functions for climatic water deficit indicate memory to water conditions up to
5 (West) to 7 (East) years in the past (Figure 3.3A). Water deficits in the growing season prior to
the year of growth had the highest antecedent weight in both study regions providing evidence in
support of hypothesis 1. The form of the memory functions indicates strong tree growth response
to water deficits in the current growing season and preceding 1-4 growing seasons followed by a
decline in responsiveness to water deficits more than 4 years in the past (Figure 3.3A). These results
are consistent with previous studies indicating water availability in past years has greater impact
on tree growth and mortality than current growing season conditions (Michaelian et al., 2011;
D’Amato et al., 2013; D’Orangeville et al., 2013). Novel here is the quantification of the relative
importance of water deficits up to 10 years in the past for current tree growth. Estimating the
antecedent effects of multi-year water deficits illustrates the persistent and cumulative effects they
have on boreal tree growth. Specifically, a five-year climatic water deficit of 10 mm is estimated to
have a measurable effect on boreal tree growth for up to 8 (West) to 10 (East) years with maximum
antecedent effects in the fifth consecutive year of water deficit (Figure 3.4A).
There was only weak evidence of regional differences in ecosystem memory to past climatic
water deficit (posterior median weights differ, but 95 percent credible intervals overlap in Figure
3.3A). The form of the estimated ecological memory functions, however, suggests trees within
western study stands are more responsive to water deficits in the current and previous growing
season and less responsive to water deficits 3-5 years in the past than trees in eastern study stands
(Figure 3.3A). The specific drivers underlying regional differences in ecosystem memory to climatic
water deficit are uncertain and warrant future analysis. It is likely regional species composition,
climate, and soils contribute to the different memory functions estimated for the eastern and western
regions. In particular, chronic summer water deficits in the West may lead to drought-adapted stands
in terms of stand density and structure, while favorable growing conditions punctuated by periodic
water deficits in the East may cause structural overshoot contributing to longer-term sensitivity to
past deficits (Jump et al., 2017).
The estimated ecological memory functions for insect defoliation indicate memory to defoliation

52

events 8 (West) to 11 (East) years in the past consistent with hypothesis 2 (Figure 3.3B). Such
memory results in persistent effects of insect defoliation on boreal tree growth, which lengthen the
longer a defoliation event lasts (Figure 3.4B). The magnitude of antecedent effects increases for
multi-year defoliation events (Figure 3.4B), or if the time between defoliation events is less than
the ecosystem memory range (Figures B.3-B.4). These increases in antecedent defoliation values
reflect the accumulation of insect defoliation damage over the course of repeated attacks leading
to reduced tree growth and potential mortality (McDowell et al., 2008; Anderegg and Callaway,
2012). Additionally, the increase in cumulative effects over consecutive years of defoliation leads
to maximum antecedent variable values one or more years after the initiation of a defoliation event
(Figure 3.4B) consistent with previous studies (Pothier et al., 2005, 2012).
While the ecological memory functions for insect defoliation were similar for the two study
regions (as evidenced by overlapping 95 percent credible intervals in Figure 3.3B), there was
evidence the memory of trees in eastern stands to SBW defoliation events was longer than the
memory of trees in western stands to FTC defoliation events (as indicated by the larger range
parameter, φ, in the East). Additional analysis is needed to explore potential mechanisms underlying
regional differences in ecosystem memory to insect defoliation. Still, there are several reasons why
ecosystem memory to insect defoliation may be longer for spruce-fir trees in the East. SBW
outbreaks last an average of 5-15 years (Fleming et al., 2000), and in severe cases can impact forest
stands for up to 20 years (Price et al., 2013). Comparatively, FTC outbreaks last an average of 2-4
years (Price et al., 2013). SBW preferentially attacks new spruce-fir needles immediately following
bud break (Gray, 2008). The needles of balsam fir trees can last up to 4-6 years (Fleming and
Piene, 1992). Thus, an SBW defoliation event occurring in a single year may reduce leaf area in
spruce-fir forests for several years to come with corresponding reductions in growth. While the
FTC is also an early-season defoliator attacking leaves soon after bud break, aspen is able to reflush
leaves several times during a growing season albeit at a significant carbohydrate cost (Anderegg
and Callaway, 2012). While the loss of leaves and the energy requirement to reflush new leaves
reduce annual growth, aspen stands are able to regenerate their leaf area more quickly following a

53

defoliation event than spruce-fir stands.

3.5.2

Antecedent Effects

Antecedent climatic water deficit was negatively related to mean annual growth of both large- and
small-diameter trees in eastern and western study stands (Table 3.1). Further, antecedent insect
defoliation was negatively related to mean annual growth of host trees in both study regions, but
did not have an observable effect on the mean growth on non-host trees. Contrary to hypothesis
3, however, there was no evidence of negative interactive effects of antecedent water and insect
defoliation stress on host tree growth (Table 3.1). These results run counter to current physiological
theory (McDowell et al., 2011; Anderegg et al., 2015), yet are consistent with previous studies
assessing tree responses to interactions between drought and defoliation stress (Jactel et al., 2012;
Jacquet et al., 2014; Kolb et al., 2016). The lack of negative interactive effects observed among
host trees even after accounting for persistent and cumulative tree growth responses to water and
defoliation stress, combined with similar results from previous studies, implies negative interactions
between drought and defoliating insects have minimal effects on tree growth. Instead, insect
defoliation, though damaging to host trees, may offset the impacts of water stress leading to
negative additive effects of water deficit and insect defoliation, but no interactive effects (Table
3.1). The positive interactive effect we observed for non-host trees underscores this potential offset.
Specifically, we found evidence of increased growth of large-diameter, non-host trees in eastern
stands at moderate antecedent water stress when there was high antecedent defoliation (Figure
3.5). This result supports the hypothesis that insect defoliation has a positive effect on non-host
tree growth during droughts by lessening evaporative water demands through leaf area reduction
(Jacquet et al., 2014), and has been demonstrated in previous studies in the eastern Canadian boreal
forest (Duchesne and Ouimet, 2008).
Our Bayesian hierarchical model can be used in a range of forest ecosystems characterized
by different species compositions, climates, and insect disturbance agents (as the current analysis
demonstrates). While offsets in water stress may reduce negative interactions between drought and
54

defoliating insects, such interactions are more likely to be detected among bark beetles given their
impact on tree xylem structure (McDowell et al., 2011; Anderegg et al., 2015). Importantly, our
model can be used to test for persistent and cumulative tree responses to bark beetle attack and the
interactive effects of drought and bark beetles on tree growth and mortality.
Potential non-negative interactions between drought and defoliating insects on tree growth and
mortality warrant further investigation. Several factors, in particular, may have contributed to the
lack of observed negative interactive effects among host trees in the current study.
i) There were no severe droughts in either boreal region during the study period. Figure 3.1B
indicates study periods in the East and West coincide with periods of relative wetness over
the past century. The climate was drier, on average, in both regions prior to circa 1960. More
severe water stress may trigger stronger tree growth responses leading to interactions with
antecedent defoliation not observed in the current study. Indeed, ecological theory indicates
there are likely drought thresholds beyond which tree responses occur, but below which, trees
are able to maintain basic physiological function (Allen et al., 2015).
ii) The current analysis utilizes tree-ring records collected from live trees alone. Trees that
suffered drought- or insect defoliation-induced mortality in the recent past may exhibit different ecosystem memory to these disturbances and stronger responses to their interaction than
surviving trees. The extensive research demonstrating prolonged radial growth suppression
in trees that die during or following a drought relative to surviving trees suggests there is the
potential for stronger interactions among trees suffering mortality (Wyckoff and Clark, 2002;
van Mantgem et al., 2003; Das et al., 2007, 2016; Berdanier and Clark, 2016). Future work
will focus on comparing live versus recently-dead tree responses to antecedent water deficit
and insect defoliation.
iii) The lack of a mechanistic model for interactions between drought and insect defoliation
may also contribute to no observable negative interactive effects among host trees. Once
developed (see Anderegg et al., 2015), such models can be integrated into the Bayesian
55

hierarchical model developed here taking the place of the empirical, linear interaction term
between antecedent variables and may reveal negative interactions not identified in the current
analysis.
If non-negative interactions are confirmed, it would alter forecasts of the impacts of changing
climate and insect disturbance regimes on the boreal forest where the primary biotic disturbance
is insect defoliation, with corresponding impacts to the global carbon cycle, potentially modifying
management approaches within this globally-important ecosystem (Price et al., 2013; Kurz et al.,
2013).

3.6

Conclusions

We found evidence of decadal-scale memory to climatic water deficit and insect defoliation in
the growth of trees from western and eastern regions of the Canadian boreal forest characterized
by different species compositions, climates, and primary defoliating insect populations. Ecological
memory functions for water deficit and insect defoliation indicated slightly different ecosystem
memory between the two study regions with trees from eastern stands exhibiting longer-term
memory to both stressors than trees from western stands. Counter to current physiological theory
and initial hypotheses, we found no evidence of negative interactive effects between antecedent
water and insect defoliation stress even after accounting for persistent and cumulative tree growth
responses to these stressors. This counter-intuitive result, consistent with previous studies, suggests
negative interactions between droughts and insect outbreaks may have minimal effects on tree
growth in defoliator-dominated systems such as the boreal forest due to offsets in water stress
caused by defoliation. Future analysis will apply the Bayesian hierarchical model developed here to
additional data sets including the occurrence of severe droughts, live and recently-dead tree growth
records, and non-defoliating insect damage to further explore interactions between drought and
insect outbreaks. The potential lack of interactive effects between drought and insect defoliation
on boreal tree growth and mortality has important implications for our understanding of future
impacts to the boreal forest under changing climate and insect disturbance regimes and associated
56

future impacts to the global carbon cycle.

57

CHAPTER 4
ASSIMILATION OF TREE-RING AND FOREST INVENTORY DATA TO MODEL
INTERACTIONS BETWEEN CLIMATE AND FOREST DYNAMICS

4.1

Abstract

Rapid changes to long-term temperature and precipitation trends already occurring globally
may contribute to altered climate and disturbance regimes. Managing forest ecosystems in the face
of changing climate and disturbance regimes requires we understand interactions between climate
extremes, disturbance, and forest dynamics—changes in forest composition, density, size, age, and
spatial structure over time. Analyses focused on modeling these interactions face a data challenge.
Tree rings provide sufficiently long records to model climate and disturbance effects on growth,
but provide little information on past forest dynamics. Forest inventory data provide detailed
measurements of forest growth, mortality, composition, and structure over time, but are rarely of
sufficient length or resolution to understand forest responses to climate. We nested a forest growth
and yield model (the Forest Vegetation Simulator [FVS]) within a Bayesian state space framework
to assimilate tree-ring and forest inventory data in order to reconstruct past forest dynamics,
advance understanding of their interactions with climate extremes and disturbance, and identify
forest characteristics promoting resistance and resilience to changing conditions. The state space
model was applied to data from a long-term red pine thinning experiment in northern Minnesota.
The assimilation of tree-ring and forest inventory data constrained estimates of forest growth and
mortality throughout the study period with forest inventory observations reducing uncertainty in
forest density (trees/acre, basal area/acre) 13 to 44 percent, and tree-ring data reducing uncertainty
in individual tree diameter growth by 77 percent (based on 95 percent credible intervals). The
state space model resulted in estimates of forest growth and mortality more closely aligned with
observed values and regional projections of forest dynamics than unconstrained runs of the FVS
model. Future work will integrate climate and disturbance as drivers of forest dynamics within FVS,

58

test for interactions between drought and forest growth and mortality, and quantify the resistance
and resilience of experimental red pine stands to drought stress.

4.2

Introduction

Forests face an uncertain future. Projected shifts in long-term temperature and precipitation
trends have the potential to increase the frequency, severity, and duration of climate extremes such
as droughts, and disturbance events including insect outbreaks, fire, and disease (IPCC, 2013;
Clark et al., 2016). Changes in historical climate and disturbance regimes may have profound
impacts on forest demographic processes including growth, mortality, dispersal, and recruitment
(Hansen et al., 2001; Breshears et al., 2005; van Mantgem and Stephenson, 2007; Price et al.,
2013; Allen et al., 2015; Zhang et al., 2015). Forecasting potential forest responses to changing
environmental conditions requires understanding interactions among the various factors impacting
forest demographic processes (Purves and Pacala, 2008; Dietze and Moorcroft, 2011; Clark et al.,
2014; Vanderwel and Purves, 2014). Factors affecting these processes can be partitioned into
exogenous factors such as climate and disturbance, and endogenous factors including species
composition, density, canopy structure, and tree size/age distribution (Bormann and Likens, 1979).
Endogenous factors (hereafter referred to as “forest characteristics”) are particularly important in
applied forest ecology contexts because they can be modified through management and may buffer
forest responses to climate extremes and disturbance. Thus, if we can identify forest characteristics
promoting ecosystem resistance (little-to-no change in demographic processes due to disturbance)
and resilience (rapid return of demographic processes to pre-disturbance rates), we may be able to
develop adaptive management strategies to maintain vital forest demographic processes in the face
of environmental change (Millar et al., 2007; D’Amato et al., 2011).
Forest dynamics refer to changes in forest stand composition and structure (i.e., the forest
characteristics described above) over time, including stand behavior during and after the occurrence
of disturbance events (Oliver and Larson, 1996). Records of past forest dynamics, when combined
with records of historic meteorological conditions and disturbance chronologies, provide valuable

59

information to advance understanding of interactions between climate, disturbance, and forest
demographic processes. Further, we may use these records to begin to identify forest characteristics
promoting resistance and resilience of demographic processes to climate extremes and disturbance
events. For example, a number of studies have demonstrated reductions in stem density (i.e., forest
thinning treatments) increase the resistance and resilience of forest growth to drought in the short
term (Laurent et al., 2003; Sohn et al., 2016; Bottero et al., 2017), but may lead to the formation of
large canopy trees making forests more sensitive to drought over the long term (McDowell et al.,
2006; D’Amato et al., 2013). As described below, modeling the interactive effects of climate and
disturbance on forest dynamics is hindered by several data challenges.
Common data sources to study interactions between climate, disturbance, and forest dynamics
include forest inventory data and radial ring-width increment records (Clark et al., 2007; Evans et al.,
2017). Forest inventory data (i.e., periodic measurements of fixed sample plots) provide detailed
estimates of forest growth, mortality, regeneration, and composition at a stand scale over time
(Avery and Burkhart, 2002). Inventory records, however, are rarely of sufficient length to infer forest
responses to long-term climatic conditions. Further, plot re-measurements occur periodically, rather
than annually, providing temporally-coarse observations of demographic processes. For example,
annual inventory records from the USDA Forest Service, Forest Inventory and Analysis Program
exist only since 1999 (less than 20-years of data available) with individual plots re-measured on
a 5- or 10-year cycle depending on the region (United States Department of Agriculture, Forest
Service, 2007).
Radial ring-width increment cores (i.e., tree rings) collected for all live trees within an inventory
plot (e.g., all live trees with a diameter at breast height larger than some threshold) provide annual
observations of radial tree growth. Tree-ring records (dating decades to centuries) are sufficiently
long to model climatic effects on forest growth (Babst et al., 2014); however, inferring stand-scale
growth (e.g., basal area increment/acre) from tree rings alone is difficult as no growth records exist
for trees which died before the collection date (Clark et al., 2001; Foster et al., 2014). This is
commonly referred to as the fading record in paleoecology contexts (Swetnam et al., 1999). As a

60

result of the fading record, tree-ring derived estimates of past stand-scale growth are downwardly
biased with the magnitude of the bias increasing as a function of years in the past (Dawson et al.,
Unpublished Data). Further, tree rings alone do not provide estimates of stand-scale mortality or
density (e.g., trees/acre, basal area/acre) critical to reconstructing forest dynamics (but see Lorimer
and Frelich (1989) for a method to reconstruct canopy gap formation and growth release events
from tree rings).
Here we develop a Bayesian state space framework to synthesize tree-ring and repeat forest
inventory data in order to reconstruct past forest dynamics and advance understanding of their
interactions with climate and disturbance. The utility of Bayesian state space frameworks to
advance understanding of complex ecological processes is well documented in recent scientific
literature (Calder et al., 2003; Wikle, 2003; Clark et al., 2007; Dietze et al., 2013). In particular,
Clark et al. (2007) applied such a framework to assimilate tree-ring and diameter-tape measurements
to infer past tree growth. The Bayesian state space framework developed here applies a similar
approach to Clark et al. (2007), but utilizes a forest growth and yield model as a driver of forest
dynamics to jointly estimate forest growth and mortality and incorporate observations of tree
density over time. The framework allows for explicit incorporation of climate and disturbance by
building interaction terms for these variables into growth, mortality, and regeneration functions
included in the growth and yield model. The interaction terms need not be linear and may include
threshold effects consistent with recent ecological hypotheses (Allen et al., 2015). Additional
benefits of the Bayesian state space approach include explicit uncertainty quantification and the
partitioning of process error (uncertainty in growth and yield model projections) from observation
error (uncertainty in tree-ring and forest inventory measurements), which can be used to test
hypotheses regarding interactions between climate, disturbance, and forest dynamics, and refine
growth and yield models (as described in Section 4.5).
Forest growth and yield models vary in their complexity, structure (empirical versus mechanistic), and the level at which forest processes are represented (tree versus stand scale). At a minimum,
growth and yield models predict forest growth and mortality over time as a function of species

61

composition, competition level, tree size and age (or their distributions), and site quality (Weiskittel
et al., 2011). Forest growth and yield models are the primary quantitative tool used by foresters to
forecast future conditions and develop management plans (Avery and Burkhart, 2002). Importantly,
outputs of forest growth and yield models—state variables in the current analysis—are in terms
forest managers understand and in contexts they can use to inform management decisions. As such,
the use of a forest growth and yield model in the current analysis helps bridge the gap between
fundamental ecological research into the effects of climate and disturbance on forest demographic
processes and adaptive forest management.
We apply the Bayesian state space framework to tree-ring and forest inventory data from a
collection of red pine (Pinus resinosa Aiton) plantation stands in northern Minnesota, USA to
achieve the following research objectives: (i) test for the effects of drought on forest dynamics;
(ii) quantify the resistance and resilience of forest growth and mortality to drought; (iii) identify
forest characteristics promoting resistance and resilience to drought (target characteristics); and,
(iv) demonstrate the capability of a Bayesian state space model to assimilate tree-ring and forest
inventory data to reconstruct past forest dynamics with explicit uncertainty.

4.3

Data and Methods

We begin this section by providing background on the structure of the growth and yield model
used to approximate forest dynamics (Section 4.3.1). We then define a Bayesian state space
framework to assimilate tree-ring and forest inventory data in order to reconstruct past forest
dynamics and advance understanding of their interaction with climate extremes (Section 4.3.2).
Finally, we detail the red pine plantation data used to demonstrate the application of the state space
model in the current analysis (Section 4.3.3).

4.3.1

Forest Vegetation Simulator

The USDA Forest Service, Forest Vegetation Simulator (FVS) is an individual-tree, distanceindependent forest growth and yield model (Dixon, 2002). That is, FVS models growth, mortality,

62

and regeneration within a forest stand (the population unit) at the individual-tree scale, but does
not keep track of the relative position of individual trees within a stand. Tracking individual
tree locations is computationally intensive, but allows for the use of spatially-explicit competition
indices, which account for the size and species of trees in a local neighborhood centered on the
target tree (Canham et al., 2004). In lieu of such indices, FVS calculates a competition modifier as
a function of a tree’s basal area relative to the total basal area of all larger trees (Dixon, 2002). This
approach is commonly applied in distance-independent growth and yield models and is thought to
provide a more accurate estimate of competition than stand density (Weiskittel et al., 2011).
A number of FVS variants have been developed for different regions of North America (Dixon,
2002). While the variants have the same general structure (i.e., individual tree, distance independent), the specific functions used to approximate forest growth and yield vary slightly among
variants. Growth and yield functions are parameterized based on detailed tree- and stand-level data
collected within the regional domain of each variant. Currently, 22 variants of FVS exist and are
available for use (see https://www.fs.fed.us/fvs/documents/guides.shtml). We use the Lake States
variant of FVS (FVS-LS) in the current analysis. The FVS-LS variant is applicable to common
tree species located within a geographic range spanning Michigan, Wisconsin, Minnesota, and
eastern portions of North and South Dakota (Dixon and Keyser, 2008). The model was originally
developed in 1993 and updated in 2006.
FVS-LS projects diameter at breast height (DBH) growth as a function of tree-level variables
including current DBH, crown ratio, and the competition modifier discussed above, and stand-level
variables including site index. Tree height is projected as a function of current tree height, age,
the competition modifier, relative height (current height of a tree relative to the average height of
the 40 largest-diameter trees in the stand), as well as site index. FVS-LS grows large- (DBH ≥
5.0) and small-diameter (DBH < 5.0) trees separately with large tree growth driven by diameter
(DBH projected first, then height) and small tree growth driven by height (height projected first,
then DBH). Both density-independent and density-dependent mortality are represented in FVS-LS.
The stand density index (SDI), a derived value based on the empirical relationship between the

63

number of trees per unit area (e.g., trees/acre) and the quadratic mean diameter (QMD = DBH
of tree with average basal area) approximating the maximum stand density for a given mean tree
size (Avery and Burkhart, 2002), is used to determine the type of mortality. Density-independent
mortality is applied when stand density is less than a threshold percentage of the SDI (55 percent
by default) and is estimated as a function of tree diameter. Density-dependent mortality is applied
when stand density exceeds the threshold SDI percentage and includes two components. First,
stand-level mortality is estimated depending on stand density relative to the SDI; second, mortality
is assigned to individual trees as a function of their relative height. A tree list is maintained for
each modeled tree in FVS-LS including its DBH and the number of trees it represents per unit area
(the tree’s expansion factor). The mortality component of FVS-LS acts by reducing the number of
trees per unit area a tree represents—it does not remove individual trees from the tree list as dead
(Dixon and Keyser, 2008). Additional details on the functions used to represent growth, mortality,
and regeneration within FVS-LS and their associated parameter values can be found in Dixon and
Keyser (2008).
We created a reduced version of FVS-LS (hereafter, “FVS-LSlite ”) for use in the state space
framework (Section 4.3.2). FVS-LSlite comprises the growth and mortality components of FVSLS for planted red pine (see Section 4.3.3). Regeneration functions are excluded as regeneration
is not modeled in the current analysis (focus on single-cohort red pine plantations). The base
version of FVS applies an additive random effect term to the projected basal area increment for
each modeled tree to account for stochastic variability in tree growth (Dixon, 2002). We exclude
the random effect term in FVS-LSlite as this variability is captured by process and observation
variance components of the state space model (Section 4.3.2). Further, the state space model
requires deterministic growth and mortality functions. FVS-LSlite is programmed in R providing
a callable function for use in the state space model (see FVS-LSlite in supplementary material; R
Core Team, 2016).

64

4.3.2

State Space Model

A state space model is used to assimilate tree-ring and repeat forest inventory data in order to
reconstruct past forest dynamics and advance understanding of their interaction with climate. FVSLSlite serves as the process model within the state space framework, driving forest dynamics
through time (Figure 4.1). State variables, inputs and outputs of FVS-LSlite at each time point,
include individual tree DBH and forest density expressed as the number of trees per unit area for
each year in a defined study period. The state space framework is setup to accommodate inventory
data collected within a fixed area plot (or complete stand census data) and tree-ring data collected for
a subset of trees within the inventory plot (Section 4.3.3). The DBH of all trees sampled during the
first forest inventory (`) are estimated within the state space model. Each of these trees is assigned a
unique record within the FVS-LSlite tree list and tracked over the study period. Let Di,t denote the
true, unobserved DBH of tree i (i = 1, 2, . . . , `) and Nt denote the true, unobserved tree density for
a target forest stand in year t (t = 0, 1, . . . , T, where T is the last year in the study period). We define
(≡) a state vector θ t including all state variable values at time t (θ t ≡ (D1,t , D2,t , . . . , D`,t , Nt )| ,
where | indicates the transpose of a vector or matrix). FVS-LSlite is used to forecast state variable
values at the next time point subject to a process error (wt ) given their value at the current time
point. Specifically,
θ t = f (θ t−1 |Ω) + wt ,

(4.1)

where f (·|Ω) represents FVS-LSlite conditional on a set of growth and yield functional parameters
(Ω). The Markov property is implicit in Equation 4.1; the current value of θ t depends only on
its value at the previous time point. Importantly, stand-scale variables including basal area per
unit area (BA) and QMD in a given year t can be calculated conditional on the state vector (θ t ) as
follows,
BA(t)|θ t
QMD(t)|θ t


π
=
g(Di,t )(ei,t ) given, g(x) =
(2 ¡ 12)2
i=1
s
!
BA(t)
(y)(2 ¡ 12)2
= g −1 Í
given, g −1 (y) =
,
` e
π
i=1 i,t
x2

Õ̀

65



✓t

f (·|⌦)
FVS-LSlite

1

f (✓

t

✓t

1

+ vt

yt

1; ⌦

)+

1

⌦
Demographic Process Parameters

✓t
wt

1

f (✓

t; ⌦

State Vector (✓ t )

✓ t+1
)+

wt

+1

✓ t + vt

i) Individual tree diameter
ii) Trees per acre

yt

Data (yt )

Error Terms

i) Tree-ring records (annual)

wt
vt

ii) Forest inventory records (periodic)

Process
Observation

✓ t+1 + vt+1

yt+1

Figure 4.1: Overview of state space model to approximate forest dynamics driven by FVS-LSlite
and constrained by tree-ring and forest inventory data.
where ei,t is the per-unit-area expansion factor for tree i in year t, BA is assumed to be expressed in
ft2 , and QMD in inches. Further, the basal area increment per unit area (BAI) in year t, reflecting
the diameter accretion of live trees, and the amount of basal area loss due to tree mortality per unit
area (BAM) can also be derived conditional on the state vector,
BAI(t)|θ t =

Õ̀ 

	
g(Di,t ) − g(Di,t−1 ) (ei,t )

i=1

BAM(t)|θ t =

Õ̀


	
g(Di,t ) ei,t − ei,t−1 ,

i=1

where g(x) is as defined above.
Tree-ring derived diameter estimates (DBH at time of coring minus sum of ring widths between
target year and collection year) exist for a subset of modeled trees (nc 6 ` where nc denotes the
(TR)

number of trees for which tree-ring data exist) on an annual basis. Let di,t

denote the tree-ring

derived DBH for tree i in year t. DBH observations based on diameter-tape measurements collected
during forest inventories exist for each modeled tree on a periodic basis unless a tree dies prior to
66

(FC)

the inventory date. Let di,t

denote the taped DBH for tree i in year t. Tree density is derived

from forest inventory data on a periodic basis by multiplying the within-plot tree count observed in
a given inventory by a stand-level expansion factor (ratio of area of interest [acre, hectare] to plot
(FC)

area). Let nt

denote the inventory-derived forest density estimate in year t. Both tree-ring and

diameter-tape measurements are used to constrain individual tree diameter values, while inventory
data alone are used to constrain tree density values. Let yt denote an observation vector including
all tree diameter and tree density observations in year t.

(TR) (TR)
(TR) (FC) (FC)
(FC) (FC)


 (d1,t
, d2,t , . . . , dnc,t , d1,t , d2,t , . . . , d`,t , nt )| t = inventory year
yt ≡

 (d (TR), d (TR), . . . , d (TR) )|
t = non-inventory year
nc,t
 1,t
2,t
Modeled trees have inventory-based DBH observations alone if they are not cored. Further, trees
that suffer mortality during the study period are missing DBH observations between the year of
death and year T. The DBH of these trees continues to be modeled after death with no data to
constrain estimates. Although a tree dies within the fixed area plot, it still represents trees of the
same species and size within the stand. Observations of individual tree DBH and tree density in
year t are equal to the estimated state values in the same year subject to an observation error (vt ).
Specifically,
yt = At θ t + vt ,

(4.2)

where At is an incidence matrix mapping state variable values to corresponding observations.
Equation 4.2 indicates observations at time t (yt ) are independent of observations at other time
points (t = 1, 2, . . . , t − 1, t + 1, . . . , T − 1, T) conditional on the state vector, θ t .
Any state space model is defined by a process and observation equation. Here, these are given by
Equations 4.1 and 4.2, respectively. Missing from the state space model specification are variances
corresponding to the process and observation error terms. We model separate observational errors
for tree-ring and diameter-tape measurements of tree diameter, and inventory measurements of
tree density. Further, we model separate process errors for tree diameter and density. A normal
distribution is applied to model all observation and process errors. This results in the following

67

Figure 4.2: Graphical depiction of state space model including (A) the estimation of individual
tree diameter at breast height, and (B) the estimation of tree density. Arrows indicate dependence
structure.
error distributions,
iid

2)
w(Di,t ) ∼ N(0, τD
iid

w(Nt ) ∼ N(0, τN2 )
(TR) iid

v(di,t

) ∟ N(0, σd2 )
tr

(FC) iid

v(di,t

) ∟ N(0, σd2 )
fc

(FC) iid

v(nt

) ∟ N(0, σn2 ),
fc

where w(·) and v(·) reference individual components of wt and vt , respectively, “iid” indicates
2 , τ 2 , σ 2 , σ 2 , σ 2 are unknown, scalar variindependent and identically distributed, and τD
nfc
N
dtr
dfc
ances. The state space framework is represented graphically in Figure 4.2, including the temporal

dependence structure of state variables, tree-ring, and forest inventory data.
The state space model is fit using a Bayesian hierarchical approach. Unknown parameters
including the latent state vector at the initial time point (θ 0 ) and scalar variances are assigned prior
distributions as defined in Appendix C. A Metropolis-within-Gibbs Markov chain Monte Carlo
(MCMC) algorithm written in R (R Core Team, 2016) is used to sample from the joint posterior

68

distribution of all model parameters (Robert and Casella, 2004). Details on the MCMC sampler
are provided in Appendix C.

4.3.3

Birch Lake Data

The FVS state space model is applied to tree-ring and forest inventory data from a collection of
red pine (Pinus resinosa Aiton) plantation stands part of the long-term Birch Lake research forest
maintained by the USDA Forest Service in northeastern Minnesota (47°42’N, 91°56’W). Red pine
stands were planted in 1912 at a density of ∟1012 trees per acre. A thinning experiment was
established in 1957 to assess variation in forest stand structure, growth, and mortality as a function
of thinning intensity and method. Forest thinning at five intensity levels (residual basal areas: 30,
61, 91, 122, and 148 ft2 ·ac−1 ) were applied to three replicate 2-acre stands at approximately 5-10
year intervals (thinning years: 1957, 1962, 1972, 1982, 1992, 2003). Three thinning methods were
used (one replicate per intensity level): thinning from above, below, and proportional. Additionally,
three replicate control stands were established in which no thinning occurred. Treatments were
randomly assigned to the 18 study stands (Bradford and Palik, 2009; Powers et al., 2010; D’Amato
et al., 2013).
Experimental stands were sampled using a single ∟0.2 acre fixed area plot (per-acre expansion
factor ≈ 5) located randomly within each stand beginning in 1957 and at ∼5 year intervals thereafter
(10 plot inventories between 1957-2009). The DBH and species of each tree with a DBH greater
than 3.5 inches were recorded within each plot during sampling events. Ring-width increment
cores were collected at breast height from all live trees greater than ∟2 in DBH within each sample
plot in 2009 (D’Amato et al., 2013). Increment cores were prepared, analyzed, and cross-dated
using standard dendrochronological procedures.
The Birch Lake data have been used in previous studies to analyze temporal trends in stand
structure and growth (Bradford and Palik, 2009); tree mortality (Powers et al., 2010); and resistance
and resilience to drought (D’Amato et al., 2013) under different thinning regimes. The Birch
Lake experiment is optimal to test the application of the FVS state space model given it is a
69

simplified system consisting of single cohort, single species stands of known planting density
and establishment year. The frequency (∟5 years) and length (1957-2009) of forest inventory
observations combined with censused tree-ring records make for a unique forest data set. Further,
the experimental data provide an opportunity to test the efficacy of different thinning regimes to
generate target forest characteristics in order to maintain growth and mortality at average rates
during and after climate extreme events.
Here we apply the FVS-driven state space model described in Section 4.3.2 to a single control
stand (9C) from the Birch Lake red pine experiment. Given the lack of data to constrain the seedling
growth stage of FVS, the state space model was initialized 26 years after stand establishment when
the mean tree DBH was estimated to have reached 4.0 inches based on FVS diameter-height growth
relationships. Thus, the study period begins in 1938 (t = 1) and continues through the final forest
inventory in 2009 (t = T).

4.4

Results

The number of trees per acre declined from a posterior mean of 651 (95 percent credible
interval: 606, 713) to 460 (451, 469) over the course of the study period (Figure 4.3A). The basal
area per acre increased from a posterior mean of 113 ft2 (104, 127) to 285 ft2 (272, 303) over
the study period (Figure 4.3B). Finally, the QMD of the study stand increased from 5.64 inches
(5.56, 5.73) to 10.67 inches (10.42, 10.94) over the study period (Figure 4.3C). Forest inventory
observations of tree density reduced uncertainty in all stand-scale variables (note, tree-ring data
was available for all years in the study period). Specifically, mean credible interval widths for tree
density, basal area, and QMD were reduced by 16.3 trees per acre, 2.3 ft2 ·ac−1 , and 0.01 inches,
during forest inventory years relative non-inventory years. The relatively large uncertainty in tree
per acre estimates prior to 1950 reflects the lack of forest density observations to constrain state
value estimates.
Annual basal area increment declined slowly over the model period from a posterior mean of
5.5 ft2 ·ac−1 (5.1, 6.0) in 1938 to 1.9 ft2 ·ac−1 (1.5, 2.1) in 2009 (Figure 4.4). Increased stand-scale

70

Trees/acre

A

Basal area/acre (ft2)

B

Quadratic mean diameter (in)

C

Posterior Mean
95% Credible Interval
Forest inventory observations

Year

Figure 4.3: Posterior summary of stand-scale variables including trees per acre (A), basal area per
acre (B), and quadratic mean diameter (C) relative to forest inventory observations.

71

competition and corresponding tree density reductions (i.e., self thinning) likely contributed to
reductions in basal area increment through time. Annual basal area mortality increased smoothly
from a posterior mean of 0.34 ft2 ·ac−1 (0.31, 0.38) in 1938 to 0.76 ft2 ·ac−1 (0.38, 0.80) in 1983
(Figure 4.4). The smooth increase in mortality from 1938 to 1983 was due to density-independent
mortality, which is driven by tree diameter in FVS—larger diameters lead to greater background
mortality rates. Density-dependent mortality began circa 1984 at a stand age of 71 when the
posterior mean basal area mortality increased to 0.83 ft2 ·ac−1 (0.18, 2.21). The greatest basal area
mortality was estimated to occur in 2009 when the posterior mean was equal to 3.24 ft2 ·ac−1 (0.49,
4.13). The increased uncertainty in basal area mortality following the onset of density-dependent
mortality reflects uncertainty in stand density estimates relative to the SDI.
Although the primary inferential goal of the current analysis is at the stand scale, we obtain
estimates of individual-tree diameter growth from the state space model (Figure 4.5). Posterior
mean tree DBH ranged from 5.29 inches (90 percent credible interval: 0.68, 7.58) at the start of the
study period (1938) to 10.17 inches (5.80, 13.34) at the end of the study period (2009). Similar to
stand-scale variables, forest inventory and tree-ring based observations of DBH reduced uncertainty
in diameter estimates. Specifically, the mean 95 percent credible interval width for DBH in years
with no data available was 3.53 inches compared to 0.80 inches during years where only tree-ring
data were available, and 0.75 inches during years where both tree-ring and forest inventory data
were available.
The posterior mean variance for both tree-ring and diameter-tape observations of tree diameter
were significantly larger than the posterior mean process variance for tree diameter (σd2 = 0.094,
σd2
fc

2
= 0.370, τD

tr

= 0.071 in Table 4.1). The posterior mean variance for diameter-tape observations

was roughly 4 times larger than the posterior mean variance for tree-ring reconstructed diameter
values. Further, the posterior mean variance for forest inventory observations of tree density
(σn2 = 51.0) was 5 times larger than the posterior mean process variance for tree density (τN2 =
fc

10.6) reported in Table 4.1.

72

Annual basal area mortality (ft2)

Annual basal area increment (ft2)

Posterior Mean
95% Credible Interval
Forest inventory observations

Approximate onset of
density-dependent mortality

Year

Figure 4.4: Posterior summary of annual basal area increment and basal area mortality relative
to forest inventory estimates equal to the average annual basal area growth or mortality between
inventory observations.

73

Diameter breast height (in)

Tree: 1

Tree: 5

Tree: 12

Tree: 59

Tree: 104

Posterior Mean
95% Credible Interval
Tree-ring diameter observations
Forest inventory diameter observations

Year

Figure 4.5: Posterior summary of diameter at breast height state space estimates relative to tree-ring
and forest inventory diameter observations for a selection of five red pine trees.

Table 4.1: Posterior summary of observation and process variance parameters.
Variance
parameter
σd2
tr

σd2

fc
σn2
fc
2
τD
τN2

Posterior Credible interval
mean
2.5%
97.5%
0.094
0.093
0.095
0.370

0.361

0.380

51.0
0.071
10.6

41.9
0.067
8.43

62.0
0.074
13.3

74

4.5

Discussion

Few methods exist to reconstruct forest dynamics. The state space model developed here
assimilates common forest data types including tree-ring and forest inventory observations in order
to reconstruct past forest dynamics with explicit uncertainty and improve understanding of their
interactions with climate extremes and disturbance. Model results indicate the importance of
assimilating tree-ring and forest inventory observations to constrain estimates of forest dynamics
generated by FVS. Uncertainty in all stand-scale variables of interest (trees per acre, basal area
per acre, QMD) was reduced during years in which a forest inventory occurred as evidenced by
narrower credible intervals. Uncertainty in stand-scale variables was further reduced as a function
of the number of trees for which tree-ring diameter estimates existed in a given year. The mean 95
percent credible interval width was reduced by 4.5 ft2 ·ac−1 for basal area and 0.17 inches for QMD
when using tree-ring measurements from 87 versus 80 trees.
Comparing stand-scale growth and mortality estimates for red pine resulting from the state
space model to a bare ground run of FVS indicates that without data to constrain model estimates,
FVS over predicts tree density up to a stand age of roughly 60 (in 1970), then slightly under
predicts densities at older stand ages (Figure C.1). Further, the unconstrained run of FVS predicts
higher diameter growth rates than observed throughout the study period as reflected in the QMD
(Figure C.1). The combination of higher tree density and QMD results in unconstrained FVS
estimates of basal area per acre equal to the maximum basal area for the stand (340 ft2 ·ac−1 ) by
a stand age of roughly 65 years (in 1975) whereas state space estimates never reach the maximum
basal area per acre during the study period (Figure C.1). Mapping unconstrained FVS growth
and density estimates on a red pine stand density management diagram for the northern Great
Lakes region indicates estimates are along or exceed the maximum size-density line given high
tree density and large QMD (Smith and Woods, 1997, results not shown). State space estimates,
conversely, are well contained within the zone of imminent competition and mortality after the onset
of density-dependent mortality. In general, the state space model provides stand-scale estimates
consistent with forest inventory observations and suggests low diameter growth rates delay the
75

onset of density-dependent mortality until an approximate stand age of 70-75 years (1980-1985 in
Figure 4.4).

4.5.1

Model Applications

The Bayesian state space framework partitions process and observation variances for tree density
and diameter. These variance estimates provide a powerful tool to refine the FVS growth and
yield model (or other growth and yield models). Specifically, we can compare the magnitude of
the process variance relative to the observation variance to determine whether estimates of forest
growth and yield may be improved through increased sampling (if observation variance exceeds
process variance) or by modifying FVS to better approximate forest dynamics (if process variance
exceeds observation variance). In the latter case, the state space model provides a mechanism to
test whether changes to FVS improve the model. That is, a significant reduction in process variance
resulting from a modification to FVS relative to its current structure (e.g., implementing a new
diameter growth equation) would provide evidence the modification better approximates observed
forest demographic processes. In the current analysis, the process variance for all parameters
was significantly less than the observation variance (Table 4.1). This suggests model estimates of
growth and mortality may be improved through increased forest inventory sampling and a larger
number of cored trees. This is consistent with expectations given there are only 10 forest inventory
events over the 72-year study period and only 87 trees cored out of 400-500 trees per acre. We
note, however, the necessity of informed prior distributions to ensure the identifiability of variance
parameters in the current analysis (see Appendix C) may limit conclusions about process variance
relative to observation variance.
In addition to variance partitioning, the hierarchical structure of the state space model allows
for tractable propagation of uncertainty to estimates of forest growth and mortality. Here, we found
strong evidence of reduced uncertainty in estimates of forest demographic processes following
the assimilation of tree-ring and forest inventory observations to constrain the FVS model. The
hierarchical structure also allows for the integration of climate variables as drivers of forest dynamics
76

in a straightforward and testable manner. Comparing differences in the magnitude of process
variance parameters based on alternative models for climatic effects on growth and mortality
within FVS provides a tool to test competing hypotheses regarding interactions between climate
and forest dynamics (see Section 4.5.3).

4.5.2

Management Applications

The FVS-driven state space model can be used to identify forest characteristics promoting resistance
and resilience to climate extremes such as drought events. Specifically, if data from a number of
forest stands with different densities and structures exist, the state space framework can be used to
reconstruct forest dynamics within each stand in relation to the occurrence of past climate extremes.
We can then compare the growth and mortality of each stand during and after the occurrence of
a climate extreme event to average growth and mortality rates over the study period to identify
characteristics (i.e., trees per acre, basal area per acre, QMD) which result in reduced responses
to extreme conditions (resistance) or accelerated recovery to pre-extreme growth and mortality
(resilience). The state space framework can also be used to test different management scenarios to
achieve target forest characteristics and increase resistance and resilience to climate extremes. In
particular, alternative thinning regimes may be modeled using the FVS-driven state space model to
determine which regime, if any, generates the desired forest density and size distribution. Further,
the model provides a mechanism to identify which thinning regime, if any, results in increased
resistance and resilience to climate extremes (in terms of forest growth and mortality) relative to
unthinned stands. Importantly, the output of the state space model developed here including trees
per acre, basal area per acre, QMD, and the testing of management regimes, are in terms and
contexts forest managers understand and can apply in the development of long-term management
plans.

77

4.5.3

Future Steps

Future work on the state space model will include steps to ensure individual tree diameters and
growth rates are positive, both within FVS and data models. Further, a selection of critical growth
and mortality functional parameters will be treated as random variables within FVS (currently
these are fixed). We will assign prior distributions to these parameters constraining their values
to scientifically-plausible ranges and allow forest growth and mortality observations to inform
posterior parameter estimates. Most importantly, we will integrate climate variables into FVS to
model interactions between climate extremes and forest dynamics. There are several potential
approaches to integrate climate into FVS. One intuitive approach, based on recent studies of the
interactive effects of climate and competition on tree growth (Rollinson et al., 2016; Buechling
et al., 2017), is to adjust potential tree growth and mortality using a multiplicative factor estimated
as a function of climate variables. Specifically, FVS models diameter growth (X) of a tree i in a
given year t as
∗ ) f (competition) + a,
Xi,t = (Xi,t
∗ is the maximum potential diameter growth for tree i in year t, f (competition) is a
where Xi,t

function of stand-scale competition whose range is zero to one ( f (¡) ∈ [0, 1]), and a is an additive
adjustment factor (Dixon and Keyser, 2008). Similarly, FVS models density-independent percent
mortality (M) for a given tree i in year t as
Mi,t = m(1 + exp(ρ0 + ρ1 Di,t ))−1,
where m is the maximum potential percent mortality, and (ρ0 , ρ1 ) are mortality coefficients (Dixon,
2002). A multiplicative factor estimated as a function of climate (similar to the competition modifier
f (climate)) can be added to both equations to scale diameter growth and percent mortality up or
down depending on climatic growing conditions. We will conduct model selection to identify
meaningful climate variables (e.g., temperature, precipitation, climatic water deficit) and test
different formulations for multiplicative factors. Selection criteria will include changes in process

78

variance relative to observation variance and score functions such as the posterior predictive loss
(Hooten and Hobbs, 2015).
Three known drought events occurred at the Birch Lake research forest during the study period
(1948, 1961, 1988) based historic documents and meteorological records (D’Amato et al., 2013).
The drought events present an opportunity to model interactions between climate extremes and forest
dynamics and quantify resistance and resilience of growth and mortality during and after drought
conditions (see Section 4.5.2). We will apply the updated state space framework with climate
integrated as a driver of forest dynamics (see above) to model growth and mortality responses and
quantify ecosystem resistance and resilience to drought within the Birch Lake experimental stands.
Model results will contribute to advanced understanding of interactions between drought and forest
dynamics by providing functions for climatic impacts on forest growth and mortality as well as
estimates of climatic effects on forest function during and after drought.

4.6

Conclusions

The Bayesian state space model presented here provides a framework to reconstruct forest
dynamics and advance understanding of their interaction with climate extremes including identifying key climate drivers to forest growth and mortality over time. Further, it allows alternative
management scenarios to be tested in terms of their efficacy in promoting resistance and resilience
of vital forest demographic processes to climate extremes over the rotation of a forest stand (i.e.,
the lifetime of a stand). Finally, the state space framework provides a mechanism to refine forest
growth and yield models such as FVS leading to improved estimation and understanding of forest
dynamics. The limiting factor in the applicability of the state space framework to other forest
ecosystems (beyond red pine plantations) is the existence of a suitable growth and yield model (or
any other model of forest dynamics). Assuming a suitable model of forest dynamics exists for a
given region and forest ecosystem type, the state space model may be applied as it is presented here
after replacing FVS-LSlite with the chosen model of forest dynamics. As such, the state space
framework provides an array of opportunities to explore interactions between climate extremes

79

and forest dynamics and test the effectiveness of alternative adaptive management scenarios in a
range of forest ecosystems characterized by different species compositions, regional climates, and
disturbance regimes.

80

CHAPTER 5
A MODEL-BASED APPROACH TO WILDLAND FIRE RECONSTRUCTION USING
SEDIMENT CHARCOAL RECORDS

5.1

Abstract

Lake sediment charcoal records are used in paleoecological analyses to reconstruct fire history
including the identification of past wildland fires. One challenge of applying sediment charcoal
records to infer fire history is the separation of charcoal associated with local fire occurrence and
charcoal originating from regional fire activity. Despite a variety of methods to identify local
fires from sediment charcoal records, an integrated statistical framework for fire reconstruction
is lacking. We develop a Bayesian point process model to estimate probability of fire associated
with charcoal counts from individual-lake sediments and estimate mean fire return intervals. A
multivariate extension of the model combines records from multiple lakes to reduce uncertainty
in local fire identification and estimate a regional mean fire return interval. The univariate and
multivariate models are applied to 13 lakes in the Yukon Flats region of Alaska. Both models
resulted in similar mean fire return intervals (100-350 years) with reduced uncertainty under the
multivariate model due to improved estimation of regional charcoal deposition. The point process
model offers an integrated statistical framework for paleo-fire reconstruction and extends existing
methods to infer regional fire history from multiple lake records with uncertainty following directly
from posterior distributions.

5.2

Introduction

Charcoal particles deposited in lake sediments during and following wildland fires serve as
records of local to regional fire history. Sediment charcoal records are used in paleoecological
analyses to identify individual fire events and to estimate fire frequency and regional biomass
burned (i.e., the amount of organic plant matter consumed due to fire) at centennial to millennial

81

time scales (Clark, 1988, 1990; Whitlock and Millspaugh, 1996; Long et al., 1998; Power et al.,
2008). When combined with sediment pollen records, charcoal deposits can be used to infer
relationships between changing climate, vegetation, and fire regimes including fire frequency, size,
and severity (Clark and Royall, 1996; Clark et al., 1996; Long et al., 1998; Carcaillet et al., 2001;
Higuera et al., 2009; Kelly et al., 2013). In particular, combined sediment charcoal records from
multiple lakes have been used to correlate changes in regional biomass burned with shifts in regional
vegetation and/or climate (Power et al., 2008; Higuera et al., 2009; Marlon et al., 2012; Kelly et al.,
2013).
Charcoal deposits in lake sediments arise from several different sources. Large charcoal particles
(> 100 Âľm) have small dispersal distances and exhibit strong correlation with fire occurrence within
roughly 500 to 1000 meters of lakes (Clark, 1988; Whitlock and Millspaugh, 1996; Gavin et al.,
2003; Lynch et al., 2004; Peters and Higuera, 2007). Small charcoal particles (< 50 Âľm) have
larger dispersal distances (typically 1-20 km) and are indicators of regional biomass burned (Clark,
1988; Clark et al., 1996). Sediment charcoal deposits arise from primary sources—direct transport
during a fire—and secondary sources, including surface transport via wind and water of charcoal
deposited within a lake catchment (Whitlock and Millspaugh, 1996; Higuera et al., 2007). Further,
lake sediments mix over time, redistributing charcoal particles vertically and concentrating charcoal
in the lake center (Whitlock and Millspaugh, 1996). The different depositional sources and sediment
mixing increase the variability in sediment charcoal records and make inference regarding the size
and location of individual fire events difficult (Higuera et al., 2007). Despite the noise present in
sediment charcoal records, the use of such data to accurately identify local fire events has been
consistently demonstrated (Clark, 1990; Gavin et al., 2003; Lynch et al., 2004; Higuera et al., 2007).
Charcoal deposition is often expressed in terms of charcoal accumulation rate to account for
different sedimentation rates over time (CHAR; particles · cm−2 · yr−1 ). Analytical approaches to
identify individual, local fire events based on sediment charcoal records decompose CHAR into
background and peak components. The background component captures low-frequency variability
associated with time-varying charcoal production rates (e.g., changes in biomass burned), sec-

82

ondary charcoal deposition, sediment mixing, and charcoal arising from regional sources. The
peak component captures high-frequency variability associated with local fire events as well as
measurement and random error (Clark et al., 1996; Long et al., 1998). Approaches to estimate
background accumulation include low-pass filters applied to Fourier-transformed CHAR (Clark
et al., 1996; Carcaillet et al., 2001) and locally-weighted regression models (Long et al., 1998;
Gavin et al., 2006; Higuera et al., 2009). Charcoal peaks are defined as the residuals resulting
from raw CHAR series minus background CHAR or the ratio between raw and background CHAR.
A threshold is used to distinguish charcoal peaks indicative of local fire events from false peaks
attributable to elevated background deposition (Clark and Royall, 1996). Optimal thresholds, in
terms of correct identification of local fires, are estimated using sensitivity analysis (Clark and Royall, 1996) or upper quantiles of a Gaussian mixture model (Gavin et al., 2006), lacking independent
fire records to identify and validate optimal threshold values (Higuera et al., 2009).
While methods to identify local fire events based on sediment charcoal records have been
well developed over the past 30 years, an integrated statistical framework for fire identification is
still lacking (Higuera et al., 2010). We build upon existing charcoal analysis methods to develop
a hierarchical Bayesian point process model for fire identification and estimation of fire return
intervals (FRIs). The point process model offers a fully model-based approach to charcoal analysis
with several important properties. The model operates on charcoal counts directly, using an offset
term to control for sedimentation rate. We generate an explicit estimate of the probability of fire
for each charcoal count. The hierarchical Bayesian approach makes for tractable error propagation
allowing for a complete treatment of uncertainty sources in sediment charcoal records including
uncertainty associated with sediment age models. The model is easily extended to multivariate data
sets, allowing for pooling of sediment charcoal records among lakes. While methods currently exist
to pool charcoal records (Power et al., 2008), the point process model requires no transformation
or interpolation of charcoal counts, improving interpretability of results and avoiding potential
introduction of non-quantifiable error to charcoal data sets. The modeling approach objectively
identifies parameter values controlling the decomposition of sediment charcoal into background

83

and peak components via regularization. Most importantly, our hierarchical Bayesian point process
model provides an integrated probabilistic framework to identify local fires and estimate FRIs across
multiple lake records with explicit uncertainty quantification.
The remainder of this chapter is organized as follows. In Section 5.3, we develop a point
process model for charcoal deposition using data from a single lake (Section 5.3.1) and a regional
network of lakes (Section 5.3.2). Estimation of local fire probability and mean FRIs are described
in the context of developing the single-lake model and the multiple-lake model. Implementation
of the two models applying Bayesian inference is described in Section 5.3.3. We demonstrate
the application of the single-lake and multiple-lake models to both simulated data and observed
data from a regional network of lakes (Section 5.4). We conclude with a discussion of modeling
properties and results (Section 5.5).

5.3

Bayesian Point Process Model for Charcoal Deposition

We construct a Bayesian point process model that relates charcoal deposition in lake sediments
to local and regional fire occurrence. Central to the model is the separation of charcoal arising
from regional and secondary sources from charcoal attributable to local fires (as depicted in Figure
5.1). In practice, charcoal particles arising from different sources (i.e., regional and secondary
sources versus local fire events) are indistinguishable in sediment charcoal records (apart from the
size distinction noted earlier). We separate background from peak deposition by assuming charcoal
particles are generated by independent processes in time: a smooth background process, exhibiting
low-frequency changes in charcoal deposition rates over time, and a highly variable foreground (or
peak) process, exhibiting high-frequency changes in charcoal deposition rates associated with local
fire events. Total charcoal deposition is proportional to the sum of the background and foreground
processes (Clark and Royall, 1996). The separation of total charcoal deposition into background
and foreground processes provides the necessary analytical mechanism to identify local fire events
from noisy sediment charcoal records. We begin this section by defining a univariate point process
model to identify local fires events. We then extend the univariate model to accommodate sediment

84

local fire|
j,f

Pj ⌘ Pr{local fire}

j,f

=

j

= foreground intensity

j,b

+

regional charcoal|

j,f

j,b

j,b

Random Labeling

Pj

j

(1

Pj )

= background intensity

j

Local Fire Event

Regional Fire Event

j,b

j,f

Time

ion

Time

it
os

La

Sediment Charcoal Record

⌧j,i

{

Dj

⌧j,i = time interval

Observed Count

p
De
ke

yf (⌧j,i )
yb (⌧j,i )

Dj = temporal domain

Foreground Count
Background Count

Time

Figure 5.1: Illustration of theoretical charcoal deposition to a lake if charcoal particles arising
from regional fires were distinguishable from particles arising from local fires (in practice, charcoal
particles from different sources are indistinguishable, i.e., we observe only y(τ j,i ) = yb (τ j,i ) +
y f (τ j,i )). The figure also depicts the formation of a sediment charcoal record from deposited
charcoal and the structure of the charcoal data (i.e., counts of charcoal over discrete, non-overlapping
time intervals. Note, the figure does not depict charcoal arising from secondary sources such as
surface water runoff or sediment mixing.
charcoal records from a regional network of lakes using a multivariate model.

5.3.1

Univariate Model

Charcoal counts are observed over time intervals spanned by the bottom and top ages of a sediment
(b)

(a)

(b)

(a)

core section: τ j,i = t j,i −t j,i , where t j,i and t j,i are the bottom and top ages of sediment core section
i (i = 1, . . . , n j ) from lake j ( j = 1, . . . , k). The τ j,i correspond to non-overlapping time intervals
Ðn j
such that i=1 τ j,i = D j where D j is the temporal domain of lake j. Throughout, we use τ j to denote
the set of observed time intervals, the temporal support, for lake j (i.e., τ j = (τ j,1, τ j,2, . . . , τ j,n j )0).
Let y(τ j,i ) equal the observed charcoal count for τ j,i defining yτ j = (y(τ j,1 ), y(τ j,2 ), . . . , y(τ j,n j ))0.
We model the accumulated charcoal particles over each observed time interval using a Poisson
likelihood conditional on a latent, continuous intensity process mapped to the observed temporal

85

support. Specifically, we let

n

yτ j |λ j ∼

j
Ö

Poisson(µ(τ j,i )),

(5.1)

i=1

where µ(τ j,i ) arises as the temporal aggregation of the continuous intensity process,
∍
µ(τ j,i ) =
Îť j (t)dt,
t ∈ τ j,i,
τ j,i

for i = 1, . . . , n j and j = 1, . . . , k. We assume the latent continuous intensity process Îť j can
be decomposed into additive continuous background and foreground intensity processes: Îť j =
Îť j,b + Îť j, f , where Îť j,b and Îť j, f denote the background and foreground intensities, respectively.
This assumption allows us to model the arrival of charcoal particles from regional/secondary
sources and local fire events as separate, independent Poisson point processes conditional on Îť j,b
and Îť j, f (Figure 5.1).
The background and foreground intensities both vary with changes in environmental conditions,
but at different frequencies. The background intensity varies with climate and regional vegetation
patterns on a centennial-to-millennial scale; the foreground intensity varies with local fuel loads
and weather patterns on an annual-to-decadal scale. We model changes in the background and
foreground intensities at the observed temporal support τ j , assuming intensity values are constant
within each observed time interval. This represents a discretization of the continuous intensity
functions. The length of observed time intervals, however, are extremely short relative to the span
of sediment charcoal records ( |τ j,i |/|D j | ≈ 1.0 × 10−3 where |·| represents the length) such that the
approach mimics an inhomogeneous Poisson point process over D j . Importantly, the environmental
conditions driving changes in the intensities are unlikely to vary substantially in the observed time
intervals.
We apply a basis expansion in time to approximate the temporal dynamics of the background
and foreground intensities as a proxy for changing environmental conditions. Specifically, we
model the background and foreground intensity processes at temporal support τ j on the log scale

86

as,


(b)
(b)
log λ j,b (τ j,i ) = β0, j + x(τ j,i )0 β j


(f)
(f)
log λ j, f (τ j,i ) = β0, j + x(τ j,i )0 β j
(f)

(b)

(b)

where β0, j and β0, j are intercept terms, β j

(f)

and β j

(5.2)

are p-dimensional vectors of regression

coefficients, and x(τ j,i ) is a p-dimensional set of known covariate values corresponding to basis
function values at p knots. We assume the latent intensity processes Îť j,b and Îť j, f are independent
conditional on their respective regression coefficients. We apply natural cubic splines as our basis
functions in (5.2), although alternative spline or predictive process basis functions may also be used.
Following from the previous discussion of the environmental factors driving changes in the two
intensities over time, we assume the background intensity process is smooth, while the foreground
(b)

process is highly variable. We regularize β j

(f)

and β j

to control the relative smoothness and

volatility of the two processes (Hooten and Hobbs, 2015). Specifically, we apply separate scalar
(b)

(f)

penalties to β j and β j

(equivalent to informed prior variances, see Section 5.3.3.1) to constrain

the background intensity to be smooth and allow the foreground intensity to be sufficiently flexible
to capture short periods with high charcoal counts attributable to local fires. Optimal penalties are
identified based on out-of-sample prediction (Section 5.3.3.3). Lacking observations of background
and foreground charcoal counts (as depicted in Figure 5.1), the coefficients in (5.2) would be
unidentifiable given identical covariate values, without the above assumption and corresponding
regularization.
Following from (5.2), we can express the mean background (Âľb ) and foreground (Âľ f ) charcoal
counts at temporal support τ j as
µb (τ j,i ) =

∍

µ f (τ j,i ) =

∍

τ j,i

τ j,i

(b)
(b)
β +x(τ j,i )0 β
j |τ |
λ j,b (τ j,i ) dt = e 0, j
j,i

λ j,b (τ j,i ) dt =

(f)
(f)
β +x(τ j,i )0 β
j |τ
0,
j
e

(5.3)
j,i |.

The mean in (5.1) can then be expressed as µ(τ j,i ) = µb (τ j,i ) + µ f (τ j,i ), defining the total charcoal
count in τ j,i as a function of independent, homogeneous Poisson point processes conditional on
87

λ j,b (τ j,i ) and λ j, f (τ j,i ). The point process model applied is similar to a Cox process in that charcoal
counts are observed over disjoint time intervals (τ j,i ) that span D j conditional on the sum of
variable intensity functions. Unlike a Cox process, however, the latent background and foreground
intensities are non-stochastic. For example, adding normally-distributed random error terms to the
background and foreground intensity functions in (5.2) results in two independent log-Gaussian
Cox processes. The variances associated with such random error terms are unidentifiable in the
current application given we do not observe background and foreground charcoal counts.

5.3.1.1

Probability of Fire

Charcoal influx at time t arising from local fire events is distinguished from regional/secondary
sources according to P j (τ j,i ) ≡ Pr{fire event local to lake j at time t} for t ∈ τ j,i . That is, a charcoal
particle arriving at lake j at time t generated conditional on λ j (τ j,i ) = λ j,b (τ j,i ) + λ j, f (τ j,i ) for
t ∈ τ j,i , is labeled as a background or foreground particle according to independent Bernoulli
trials with probability P j (τ j,i ) (Diggle, 2014). It follows that λ j, f (τ j,i ) = P j (τ j,i )λ j (τ j,i ), and
λ j,b (τ j,i ) = (1 − P j (τ j,i ))λ j (τ j,i ) (Ross, 2010, see Figure 5.1), so that
P j (τ j,i ) =

λ j, f (τ j,i )
λ j, f (τ j,i ) + λ j,b (τ j,i )

.

(5.4)

While the coefficients in (5.2) may not all be statistically identifiable, even with the constraints
placed on the background and foreground intensities, the resulting probability of fire as defined in
(5.4) is identifiable (see Appendix D).

5.3.1.2

Mean Fire Return Interval

A chronology of local fire events is derived after fitting the point process model by establishing
a probability of fire threshold, which indicates a local fire event when exceeded. Specifically, we
transform P j (τ j,i ) into a binary variable Z j (τ j,i ) with one indicating a local fire event according to

88

f

P j (τ j,i ) −→ Z j (τ j,i ) where


 1 if P j (τ j,i ) > ξ,

f P j (τ j,i )|ξ =
,

 0 if P j (τ j,i ) 6 ξ





for a defined probability of fire threshold Ξ. The fire chronology generated depends on the specific
probability of fire threshold applied. We discuss selection of an optimal threshold in Section
5.3.3.4.
The series of estimated fire events resulting from the application of the probability of fire
threshold constitute a temporal Poisson process defined by a rate parameter (α−1 ). A property
of Poisson processes is the interarrival times, in this case, the time intervals between local fire
events are independently and identically distributed as exponential random variables with mean
Îą. The exponential mean represents the mean FRI, while its inverse, the Poisson rate parameter
represents the frequency of local fires. The maximum likelihood estimate (MLE) for Îą is equal to
the observation period divided by the number of fire events observed for a given lake, for example,
α̂ j = Ín
j

|D j |

.

Z (τ )
i=1 j j,i

We apply Bayesian inference to estimate the exponential mean parameter (Îą) as described in
Section 5.3.3.4, rather than calculating the MLE, to allow for estimation of a regional mean FRI (see
Section 5.3.2.1). It is common in fire ecology to apply a Weibull likelihood function to estimate
the mean FRI. Applying a Weibull likelihood function allows for the probability of a fire event to
increase as a function of time elapsed since the last fire event. We did not apply a Weibull likelihood
function in the current analysis because sediment charcoal records were of relatively short length
for a number of study lakes, leading to poor estimation of the two Weibull likelihood parameters.

5.3.2

Multivariate Model

The multivariate model follows directly from the univariate model and allows for joint estimation
of fire probabilities and mean FRIs across multiple lakes. We combine observations from all lakes

89

τ1,1

τ1,n1

τ2,1

Lake 2

τ2,n2

τ3,1

Lake 3

τ3,n3

τ1∗

Lake 1

τn∗∗

Combined

Figure 5.2: Example creation of new temporal support τ ∗ from three individual lake records. Note
the temporal misalignment across lakes and τ ∗ .
yτ ≡

(y0τ , y0τ
1
2

. . . , y0τ )0 and define new temporal support τ ∗
k

the temporal misalignment among individual lake records

0
∗
∗
∗
= τ1 , τ2 , . . . , τn∗ to accommodate
(Figure 5.2). τ ∗ is defined by non

overlapping intervals of equal length |τi∗ | = |τi∗0 |, ∀ i, i0 = 1, . . . , n∗ , that span D∗ , the temporal
Ð n∗ ∗
domain spanned by all lakes combined D∗ = i=1
τi .
The background intensity process for each lake is defined at temporal support τ ∗ allowing
charcoal counts to be pooled across lakes to estimate a regional mean background intensity. Comparable with the univariate model, we model the background intensity process for each lake on the
log scale, but apply a new set of cubic regression splines defined for p∗ knots corresponding to
temporal support τ ∗ . Specifically,


(b)
(b)
log λ j,b (τi∗ ) = β0, j + x(τi∗ )0 β j ,
where x(τi∗ ) is a p∗ -dimensional set of known cubic regression spline covariate values equal to the
ith row of the n∗ × p∗ matrix X.
The background intensity process defined at temporal support τ ∗ is mapped back to the observed
Í
temporal support for each lake τ j by an N × n∗ matrix A where N = kj=1 n j . Specifically, given
A ≡ (A01, . . . , A0k )0 where each A j is an n j × n∗ dimensional matrix,


(b)
(b)
λ j,b (τ j,i ) = a0j,i exp β0, j 1 + Xβ j
where a j,i is a n∗ -dimensional vector equal to the ith row of A j and 1 is a n∗ -dimensional vector of
90

ones. The lth entry of a j,i is equal to
a j,i (τl∗ ) = |τ j,i ∩ τl∗ |
such that,

l = 1, . . . , n∗

∗

n
Õ

a j,i (τl∗ ) = |τ j,i |, i = (1, . . . , n j ), j = (1, . . . , k).

l=1

A joint, regional background intensity would ideally be modeled as a spatio-temporal process across
multiple lakes at temporal support τ ∗ . Given the time and challenges associated with collecting
sediment charcoal records, however, sediment samples from regional lake networks are rarely of
sufficient size (in terms of the number of lakes sampled) to allow estimation of a spatio-temporal
background intensity process. In the current analysis, charcoal counts are pooled across lakes to
(b)

(b)

estimate a regional mean background intensity process by assigning the β0, j and β j exchangeable
normal prior distributions, as described in Section 5.3.3, in lieu of estimating a spatio-temporal
background intensity process. The foreground intensity process is modeled exactly as in the
univariate model. Specifically, the foreground process is modeled at the observed temporal support
(f)

for each lake (τ j ) using lake-specific foreground coefficients (β j ) corresponding to a set of p
knots. Probability of fire estimates are calculated for each lake independently according to (5.4).

5.3.2.1

Regional Mean Fire Return Interval

We seek inference regarding the regional mean FRI in addition to individual-lake mean FRIs
under the multivariate model. We apply a partial pooling approach to estimate the regional mean
FRI across lakes. Specifically, individual-lake mean FRI values (Îą j ) are assigned exchangeable
2 . The
log-normal priors centered on the log of a regional mean FRI (log α∗ ) with variance σfri

partial pooling approach allows charcoal records from each lake to inform the regional average,
but penalizes lakes with large uncertainty in their mean FRI value estimate. The regional mean
2 ) quantifies the deviation of individual-lake mean FRI values from the
FRI variance parameter (σfri

regional average.

91

5.3.3

Bayesian Implementation

The univariate and multivariate models are completed by specifying prior distributions for remaining unknown parameters. These include background and foreground regression coefficients and
mean FRI parameters.

5.3.3.1

Univariate Model

We assigned normal priors to regression coefficients under the univariate model: β0, j ∟ N(0, σ02 ),
(b)

2
β j ∼ N(0, σb,2 j S−1
j ), β0, j ∟ N(0, σ0 ), and β j
(f)

(b)

(f)

2
∼ N(0, σ 2f , j S−1
j ). In our specification, σ0 is fixed

at a large value defining a diffuse normal prior. The σb,2 j and σ 2f , j terms represent scalar penalties (or
regulators) that ensure the background process is smooth while the foreground process is sufficiently
flexible to capture irregular charcoal counts subject to the constraint σb,2 j < σ 2f , j (identification of
optimal penalties is described in Section 5.3.3.3). The p× p matrix S j consists of known coefficients
defined as a function of the selected knot values (Wood, 2006). Note that S j is not full column
rank, rather its rank is p − 2 given that the second derivative of the boundary knots are equal to
(b)

(f)

zero for natural cubic splines. Thus, the priors for β j and β j

are improper, but can be shown to

result in proper posterior distributions. Combining the likelihood from (5.1) with the priors for the
regression parameters, the joint posterior distribution for a single lake under the univariate model,
using notation similar to Gelman et al. (2014), is proportional to:
n

j
Ö
i=1

Pois(y(τ j,i )| µ(τ j,i )) × N(β0, j |σ02 ) × N(β j |σb,2 j , S j ) × N(β0, j |σ02 ) × N(β j |σ 2f , j , S j ). (5.5)

5.3.3.2

(b)

(b)

(f)

(f)

Multivariate Model

We specified identical normal priors for the foreground regression coefficients under the multivariate model as in the univariate model, and exchangeable normal priors for the background


(b)
(b)
(b) (b)
(b) 0
coefficients: β 0 ∼ N(µ0,b 1, τb2 I k ) and β (b) ∼ N(Rµ b, Σ b ) where β 0 ≡ β0,1, β0,2 . . . , β0,k ,


(b)0 (b)0
(b)0 0
β (b) ≡ β 1 , β 2 , . . . , β k
, µ0,b and τb2 are univariate mean and variance parameters, 1

92

is a k-dimensional vector of ones, I k denotes a k-dimensional identity matrix, µ b is a p∗ dimensional mean vector, Σ b is a k p∗ × k p∗ covariance matrix, and R is a k p∗ × p∗ incidence
matrix equal to 1 ⊗ I p∗ . The univariate variance τb2 quantifies inter-lake variation in the intercept
of the background intensity. The covariance matrix ÎŁ b can be decomposed into inter-lake and
within-lake covariance in background coefficients. Defining a k p∗ × k p∗ block diagonal matrix
2 S∗−1, σ 2 S∗−1, . . . , σ 2 S∗−1 ) where S∗ is a p∗ × p∗ matrix of known coefficients
Sb ≡ Diag(σb,1
b,2
b,k

associated with the p∗ knots defined for τ ∗ , we can express the background coefficient covariance
matrix as Σ b = LSb L0 where LL0 is the Cholesky decomposition of H ⊗ I p∗ where H is a kdimensional covariance matrix. The matrix Sb accounts for covariance in background coefficients
within lakes, while H captures covariance among lakes. We apply a spatial covariance function to
construct H, although any valid covariance function can be used. Specifically, H ≡ Hs (θ b ) where
(Hs (θ b )) j, j 0 = c(||s j − s j 0 ||; θ b ) given s j indicates the geographic location of lake j and θ b are
unknown spatial covariance parameters.


In the current analysis, we assigned normal priors to the mean intercept ¾0,b ∟ N 0, σ02 and


mean regression coefficients ¾ b ∟ N 0, ψ 2 I where σ02 is fixed at a large value, while ψ 2 is set to
(b)

an appropriate order of magnitude for the β j based on univariate model results. The univariate
among-lake standard deviation parameter τb was assigned a uniform prior: τb ∼ Unif(aτ, bτ ).
The scalar penalty on the background deposition process for each lake (σb,2 j ) was set equal to the
2 , σ 2 , . . . , σ 2 ). We apoptimal background penalty value from the univariate model σ 2b ≡ (σb,1
b,2
 b,k


plied an exponential spatial covariance function to form H: (Hs (θ b )) j, j 0 = σs2 exp −φ||s j − s j 0 ||

where θ b ≡ (σs2, φ)0. The partial sill (σs2 ) represents spatial variance in background regression
coefficients among lakes and has the potential, if it is large, to generate highly-variable, unconstrained background deposition processes for individual lakes. To avoid the generation of overly
flexible background deposition processes, we fixed the partial sill at one (σs2 = 1). The spatial
decay parameter was treated as a free parameter and estimated applying a diffuse uniform prior:
φ ∼ Unif(aφ, bφ ). Combining the joint likelihood for charcoal counts from all lakes with the priors
for the multivariate regression model parameters, the joint posterior distribution for all lakes under

93

the multivariate model is proportional to:
n

j
k Ö
Ö
j=1 i=1

Pois(y(τ j,i )| µ(τ j,i )) × N(β 0 | µ0,b, τb2 ) × N(β (b) |µ b, σ 2b, S∗, φ)×
(b)

k
Ö
j=1

(f)
N(β0, j |σ02 ) ×

k
Ö

(5.6)

N(β j |σ 2f , j , S j ) × N(µ0,b |σ02 )×
(f)

j=1

N(µ b |ψ 2 ) × Unif(τb |aτ, bτ ) × Unif(φ|aφ, bφ ).
We use a Metropolis-within-Gibbs MCMC algorithm (Robert and Casella, 2004) to sample
from the posterior distributions in (5.5) and (5.6).

5.3.3.3

Identification of Optimal Penalties

Optimal background (σb,2 j ) and foreground (σ 2f , j ) penalties (prior variances) were identified based
on out-of-sample model prediction. Specifically, we conducted a grid search over a range of
penalty values based on initial exploratory modeling with five unique penalty values assigned to the
background and foreground: 25 total combinations. The gridded search was conducted applying
the univariate point process model to a single validation data set for each lake with 25 percent
of observations held out. Prediction of held-out charcoal counts was carried out via composition
(b)

(b)

(f)

(f)

sampling using posterior samples of background and foreground coefficients (β0, j , β j , β0, j , β j )
to generate λ j,b (τ j,i ) and λ j, f (τ j,i ). We sampled yho (τ j,i ) ∼ Pois(µ(τ j,i )) in a one-for-one fashion,
where yho (τ j,i ) indicates a held-out observation. The resulting samples of yho (τ j,i ) represent
the posterior predictive distribution of yho (τ j,i ). Optimal penalty terms were identified as the
background and foreground variances that minimized the posterior predictive loss calculated for
the hold-out data (Gelfand and Ghosh, 1998). The posterior predictive loss rewards accuracy of
predictions with a penalty for large variance in predictions indicative of over parameterization. For
the multivariate model, we applied the optimal background and foreground penalty values identified
for each lake under the univariate model.

94

5.3.3.4

Mean FRI

The mean FRI under the univariate and multivariate models is estimated using FRIs calculated
after applying a probability of fire threshold to sampled posterior probability of fire values at each
iteration of the Gibbs sampler. Specifically, for the `th iteration of the Gibbs sampler, we obtain

 (`)
(fire) (fire)
(fire)
a set of fire event times t j,1 , t j,2 , . . . , t j,m
, where m j is the total number of fire events
j

observed for lake j, by conditioning samples on Z(τ j,i )(`) = 1. A given FRI is equal to the elapsed
time between two consecutive fire events. For example, for lake j and iteration `, the rth FRI is


(`)
(fire)
(fire) (`)
given by FRI j,r = t j,r+1 − t j,r
, for r = (1, 2, . . . , m j − 1).
We assigned a semi-informative, conjugate inverse-gamma prior to the mean FRI parameter
ι j ∟ InvGamma(aι, bι ) under the univariate model centering its density over the possible range
of FRIs for study lakes based on previous analyses (Kelly et al., 2013). Combining the prior with
the exponential likelihood of the FRIs (see Section 5.3.1.2), we obtain an inverse-gamma posterior
distribution for Îą j conditional on the derived set of FRIs. Specifically, for the `th iteration of

0
(`)
(`)
(`)
the Gibbs sampler, FRI j = FRI j,1, FRI j,2, . . . , FRI j,m j −1 and α j |FRI j ∼ InvGamma(aα +
(`)
mj

− 1, bα +

(`)
Ím j −1
r=1

(`)

FRI j,r ).

We seek inference regarding the regional mean FRI in addition to individual-lake mean FRIs
under the multivariate model. We assigned diffuse uniform priors to the regional mean FRI (α∗ )
and the inter-lake standard deviation (σfri ). Combining the priors and the exponential likelihood
of the FRIs, the joint posterior distribution conditional on the model parameters defined in (5.6) is
proportional to:
j −1
k mÖ
Ö
j=1 r=1

Exp(FRI j,r |α j ) ×

k
Ö

2 ) × Unif(α∗ |a , b ) × Unif(σ |a , b ).
N(log α j |log α∗, σfri
∗ ∗
fri σ σ

j=1

We estimated mean FRI values under the univariate and multivariate models using a range of
probability of fire thresholds: Ξ = (0.50, 0.55, . . . , 1.00)0. Ideally, we would select the probability
of fire threshold yielding the greatest accuracy in terms of identifying local fire events (i.e., the
most-accurate mean FRI estimate). A common challenge with the use of sediment charcoal records,
however, is the lack of independent fire history data to conduct model validation. In the absence
95

of independent fire history data, we sought an optimal probability of fire threshold in the sense
of providing a precise mean FRI estimate based on a consistent fire chronology for each posterior
sample of P j (τ j ) ≡ (P(τ j,1 ), P(τ j,2 ), . . . , P(τ j,n j ))0. We estimated the coefficient of variation for
posterior mean FRI samples (σ̃/¾̃ where σ̃ and ¾̃ are the posterior sample standard deviation and
mean, respectively). We selected the probability of fire threshold that minimized the coefficient
of variation as the optimal threshold. The optimal threshold provided the most-consistent mean
FRI estimate based on posterior samples. This approach is similar to previous paleoecological
approaches relying on sensitivity analysis to identify an optimal threshold (Clark et al., 1996).

5.4

Model Application

We apply the point process model to both simulated data and to sediment charcoal records from
a 13-lake network in interior Alaska. The use of simulated data allows for proper model validation,
which is difficult using sediment charcoal records due to the lack of observed fire data to compare
with probability of fire estimates.

5.4.1

Simulation Study

We tested the accuracy of the point process model by applying it to simulated sediment charcoal
records generated using CharSim (Higuera et al., 2007). CharSim is a semi-mechanistic model that
generates fires on a landscape and maps the subsequent deposition of charcoal particles to a target
lake. The amount of charcoal deposited in the target lake is proportional to the size of the fire, its
proximity to the target lake, the atmospheric injection height of charcoal particles during the fire,
and secondary charcoal deposition and sediment mixing within the lake (see Higuera et al., 2007,
for additional details).
We applied the univariate point process model to a simulated sediment charcoal record generated
by CharSim to mimic a fire regime consistent with historic fire regimes in the Alaskan boreal forest
(as described in Table 2 of Higuera et al., 2007). We defined a binary variable for each sample
interval of the simulated record with one indicating a local fire event within 100 m of the target lake

96

within an interval, and 0 indicating no local fires within an interval. Probability of fire estimates
from the point process model were converted to binary fire occurrences as described in Section
5.3.1.2 and compared with true fire occurrence values to determine the percentage of fires correctly
identified.
The simulated record was 4760 years long divided into 238 equal-length sample intervals (each
interval was 20 years). There were 41 fires within 100 m of the target lake leading to a true mean
FRI of 116 years. The point process model correctly identified 38 of the 41 fires occurring over the
simulated study period applying an optimal threshold value of 0.95, corresponding to a 93 percent
fire identification rate (Figure 5.3). The model was accurate in its identification of local fires with
only a single falsely identified fire between 1880 to 1920 years before present (YBP). Despite the
accuracy of the point process model in identifying true local fires, the mean FRI was over estimated:
posterior mean FRI equaled 191 years (95 percent credible interval: 135 to 271 years). The upward
bias of the point process model in estimating the mean FRI was due to the model’s inability to
separate fires occurring in close temporal proximity as unique fire events. For example, there were
40 sample intervals that included a true fire in the simulated record (note, two fires occurred within
a single sample interval), but only 25 unique threshold exceedances (or peaks) were estimated. A
more in-depth discussion of source of bias is provided in Section 5.5. Model results were sensitive
to the definition of a local fire event. Specifically, if a local fire is defined as occurring within 1000
m of the target lake (rather than 100 m), the fire identification rate drops to 75 percent.

5.4.2

Yukon Flats

We applied the univariate and multivariate point process models to previously-published sediment
charcoal records from 13 lakes in the Yukon Flats region of Alaska (Figure 5.4; Kelly et al., 2013).
The Yukon Flats region is dominated by boreal forests and has a fire regime characterized by
stand-replacing fires with return intervals of several decades to centuries.

97

100

Fitted Charcoal Count

Count

75
50
25
0
4000

3000

2000

1000

0

5
Background
Foreground

Intensity

4
3
2
1
0
4000

3000

2000

1000

0

2000

1000

0

2000

1000

0

1.00

P(T)

0.75
0.50

False
False
Peak
Peak

0.25
0.00
4000

3000

4000

3000

100

Count

75
50
25
0

Time (YBP)

Figure 5.3: Univariate model results for simulated sediment charcoal record generated by CharSim.
Upper panel indicates observed charcoal counts along with the posterior mean charcoal count (blue
line). Second panel illustrates posterior mean foreground (orange line) and background (black line)
intensities. Third panel plots posterior mean probability of fire estimates for each observed time
interval (black line) along with the upper and lower bounds of the 95 percent credible interval (gray
shading) and the optimal threshold (red line). The points in the third panel correspond to true fire
events occurring during the sample interval with black dots delineating correctly identified fires
and red open dots delineating missed fires. The arrow highlights the single falsely identified fire
during the simulated study period. Lower panel indicates observed charcoal counts with the color
and shape indicating whether the count was correctly identified as a true fire (black points), no fire
(gray points), missed true fire (red open circles), or falsely identified fire (red points).

98

Yukon Flats Region

66°10'0"N

Jonah

!
(
!
( Screaming Lynx
!
( Granger
!
( Windy
!
( Lucky

66°0'0"N

Epilobium

!
( Reunion
!
( Chopper

!
( Latitude

!
(
!
( Robinson

!
(

Landing

65°50'0"N

¸

West Crazy
!
( Picea

!
(

Area Burned (1940 CE to Present)
0

146°0'0"W

10

20 Km

145°30'0"W

Figure 5.4: Location of study lakes within the Yukon Flats region of Alaska relative to areas burned
in wildland fire since 1940 Common Era.
5.4.2.1

Univariate Model

We applied the univariate model to each of the 13 lakes in the Yukon Flats data set and the
multivariate model to all lakes jointly. Mean FRI values from the univariate model varied from
roughly 134 years for fires local to Chopper Lake to roughly 356 years for fires local to Screaming
Lynx Lake. Table 5.1 provides mean FRI value estimates for each of the 13 lakes. Optimal threshold
values providing the most-precise mean FRI estimate for each lake varied from 0.50 to 0.75 (Table
5.1). Results of the univariate model for Chopper and Screaming Lynx Lakes are provided in Figure
5.5 (similar plots are provided for all lakes in Appendix D along with additional plots for Chopper
and Screaming Lynx Lakes illustrating the identification of an optimal threshold and uncertainty

99

in local fire chronologies; Chopper and Screaming Lynx Lakes are selected to illustrate the lakes
with the shortest and longest mean FRI, respectively).
Geospatial fire perimeter data for the state of Alaska date back to 1940 Common Era (CE;
Alaska Fire Service, 2016). Restricting the fire perimeter data to fires local to lakes within the
Yukon Flats data set, there were five fires since 1940 within 100 m of one or more of the 13 study
lakes (a threshold of 100 m from lake edge was used to distinguish local fire events based on results
of the simulation study) the earliest of which occurred in 1985. Although the local fire record is
insufficient in length to conduct proper model validation (i.e., roughly 70 years of local fire data
for a study period of over 10,000 years is less than 1 percent data coverage), we can compare the
true occurrence of local fires to local fires identified by the point process model to assess model
performance. Two of the five fires local to at least one study lake occurred after the end of the
sediment charcoal records for local lakes: Big Creek Fire in 2009, Discovery Creek Fire in 2013.
The Preacher Creek Fire in 2004 occurred local to Picea and Epilobium Lakes, however, the point
process model did not identify a local fire for either of these lakes in the most-recent 50 years.
There were two unnamed fires, the first in 1985 and the second in 1988, local to several study lakes.
A fire event was identified by the point process model within the past 25 years in four out of seven
lakes local to the 1985 fire and two of five lakes local to the 1988 fire.

5.4.2.2

Multivariate Model

Joint probability of fire estimates generated using the multivariate model varied in magnitude from
the univariate model results. Figure 5.6 presents the multivariate model results for Chopper and
Screaming Lynx Lakes. Comparing the results presented in Figure 5.6 to those from the univariate
model (Figure 5.5), the probability of fire estimates for Chopper Lake are slightly higher under the
multivariate model than the univariate model, while the probability of fire estimates for Screaming
Lynx Lake are roughly consistent between the two models. In general, probability of fire estimates
were higher under the multivariate model than the univariate model. The different magnitudes
of probability of fire estimates under the multivariate model necessitated calculating new optimal
100

Table 5.1: Summary of univariate and multivariate Poisson process model results for each lake
in the Yukon Flats data set. Mean FRI is equal to the posterior mean fire return interval with 95
percent credible intervals in parentheses.
Univariate Model
Lake
Optimal
Mean
Cred. Int.
Threshold
FRI
Width
Chopper
0.60
134 (91,197)
106
Epilobium
0.50
197 (133,292)
159
Granger
0.55
286 (187,436)
249
Jonah
0.50
159 (111,228)
117
Landing
0.70
303 (206,445)
239
Latitude
0.75
158 (106,235)
129
Lucky
0.75
136 (92,200)
108
Picea
0.50
318 (227,448)
221
Reunion
0.70
244 (167,354)
186
Robinson
0.55
142 (94,213)
119
Screaming Lynx
0.60
356 (255,497)
242
West Crazy
0.50
204 (133,315)
182
Windy
0.50
162 (110,239)
128

Multivariate Model
Optimal
Mean
Cred. Int.
Threshold
FRI
Width
0.85
144 (97,216)
119
0.50
201 (136,293)
157
0.50
246 (165,368)
203
0.55
148 (107,207)
100
0.70
291 (201,419)
218
0.75
137 (96,194)
98
0.75
138 (93,206)
114
0.50
312 (223,429)
206
0.70
237 (165,345)
180
0.80
124 (87,180)
93
0.55
357 (255,502)
247
0.70
200 (129,306)
177
0.65
155 (107,226)
119

fire thresholds for each lake (Table 5.1). Optimal threshold values ranged from 0.50 to 0.85 for
the multivariate model consistent with slightly higher probability of fire estimates compared to the
univariate model. Mean FRI values estimated using the multivariate model were consistent with the
mean FRI values estimated for each lake using the univariate model; however, the credible interval
widths for the multivariate model were narrower than under the univariate model. Specifically, 10
out of 13 lakes had narrower credible intervals under the multivariate model than the univariate
model (Table 5.1). The average credible interval width for the mean FRI was 156 years under the
multivariate model versus 168 years under the univariate model. In addition to lake-specific, local
mean FRIs, we applied the joint probability of fire estimates and optimal thresholds for each lake
to estimate a joint regional mean FRI as described in Section 5.3.2.1. We estimated a regional
mean FRI over the study period (10,680 to -59 YBP, relative to 1950 CE) of roughly 187 years (95
percent credible interval: 136 to 261 years) for the Yukon Flats region based on the lake network
data.
The multivariate model provides inference on regional background charcoal deposition. Ap-

101

300
150

Count

Count

Fitted Charcoal Count
200
100

50

0

0
1500

1000

500

0

9000

Intensity

Intensity

Background
Foreground

40
20
0

3000

0

9000

6000

3000

0

9000

6000

3000

0

Background
10

Foreground

5
0

1500

1000

500

0

1.00

1.00

0.75

0.75

P(T)

P(T)

6000

15

80
60

Fitted Charcoal Count

100

0.50

0.50
0.25

0.25
0.00

0.00
1500

1000

500

0

Time (YBP)

Time (YBP)

(a)

(b)

Figure 5.5: Univariate model results for Chopper Lake (a) and Screaming Lynx Lake (b). Upper
panel indicates observed charcoal counts along with the posterior mean charcoal count (blue line).
Middle panel illustrates posterior mean foreground (orange line) and background (black line)
intensities. Lower panel plots posterior mean probability of fire estimates for each observed time
interval (black line) along with the upper and lower bounds of the 95 percent credible interval (gray
shading) and the optimal threshold (red line).

150

Fitted Charcoal Count

200

Count

Count

300

100

Fitted Charcoal Count

100
50

0

0
1500

1000

500

0

9000

6000

3000

0

9000

6000

3000

0

9000

6000

3000

0

15
Background
40

Intensity

Intensity

60
Foreground

20
0

Background
Foreground

5
0

1500

1000

500

0

1.00

1.00

0.75

0.75

Prob(Fire)

Prob(Fire)

10

0.50
0.25
0.00

0.50
0.25
0.00

1500

1000

500

0

Time (YBP)

Time (YBP)

(a)

(b)

Figure 5.6: Multivariate model results for Chopper Lake (a) and Screaming Lynx Lake (b). Upper
panel indicates observed charcoal counts along with the posterior mean charcoal count (blue line).
Middle panel illustrates posterior mean foreground (orange line) and background (black line)
intensities. Lower panel plots posterior mean probability of fire estimates for each observed time
interval (black line) along with the upper and lower bounds of the 95 percent credible interval (gray
shading) and the optimal threshold (red line).
102

plying an exponential spatial covariance function to describe spatial correlation in background
regression coefficients (as described in Section 5.3.3.2) resulted in an estimated effective spatial
range of approximately 16.5 km (95 percent credible interval: 14.4 to 19.1 km) where the effective
spatial range is defined as the distance at which the correlation drops to 0.05. This suggests the
parameters describing local, background charcoal intensity for individual lakes are similar for lakes
within 20 km of each other. Finally, the multivariate model provides an estimate of the background
charcoal deposition in each lake over the entire study period, although the sediment charcoal records
for most lakes are shorter than the full study period. The background charcoal intensities for each
lake in the Yukon Flats network are plotted together in Figure 5.7 along with a regional loess smooth
function. The background charcoal deposition for most lakes exhibited a similar pattern, with a
long-term increase in background charcoal deposition from 6000 YBP to present, a sharp increase
in background deposition roughly 3000 YBP, and a secondary increase 1000 YBP followed by a
decrease in background deposition roughly 500 YBP.

5.5

Discussion

The use of sediment charcoal records to reconstruct past fire regimes is challenging given
charcoal counts rather than past fire occurrences are observed. Further, observed charcoal counts
include charcoal generated during local fires as well as charcoal stemming from regional fire activity
and secondary sources. The goal of our analysis was to construct an integrated statistical framework
for local fire identification and the estimation of mean FRIs based on sediment charcoal records
from individual lakes building on previous approaches to paleo-fire reconstruction. We further
sought to advance existing approaches to reconstruct regional fire history through the development
of a multivariate model, which combines sediment charcoal records from multiple lakes to identify
local fires and jointly estimate mean FRIs at individual-lake and regional scales. Here, we discuss
the key results of the application of the univariate and multivariate point process models to the
Yukon Flats data set and connect the results of the current analysis to previous studies in the same
region.

103

1.5

Intensity

1.0

0.5

Number of Records

0.0

12
9
6
3
10000

8000

6000

4000

2000

0

Time (YBP)
Figure 5.7: Background charcoal deposition intensity for each lake in the Yukon Flats data set over
the entire study period (10,680 to -59 YBP, relative to 1950 CE) based on the multivariate point
process model. The bold line in the upper panel is a regional loess smooth function reflecting mean
changes in background charcoal deposition across all lakes (fit using a span of 0.15). The lower
panel indicates the number of lake records used to estimate background charcoal deposition over
time.

104

Mean FRI estimates from the univariate and multivariate point process models applied to the
Yukon Flats data set ranged from 100 to 350 years (Table 5.1). The lakes with the largest mean
FRI estimates have sediment charcoal records dating back the longest among study lakes: Reunion,
Granger, Landing, Picea, and Screaming Lynx Lakes all have records that date back at least 5000
YBP. This pattern in consistent with previous interpretations of Holocene fire history in Alaskan
boreal forests, which highlight increased fire activity over the last several thousand years beginning
with the local arrival of black spruce between 6000 to 3000 YBP (Lynch et al., 2004; Higuera
et al., 2009; Kelly et al., 2013).We observe a similar pattern of increased fire activity in plots of
the background charcoal deposition intensity for each lake over the study period derived from the
multivariate model (Figure 5.7). The increase in fires across the Yukon Flats region from 6000 to
3000 YBP reflected in the background intensity suggests that lakes with sediment charcoal records
dating back prior to 3000 YBP should have longer mean FRI values than lakes with relatively short
records. A secondary peak in the background charcoal deposition intensity is observable around
500 YBP coincident with the Medieval Climate Anomaly, a period of increased temperatures
and drought frequency (1000-500 YBP), followed by the Little Ice Age, a period of cooler and
wetter climatic conditions (500-80 YBP). Finally, modeled background charcoal intensities for
individual lakes indicate an increase in biomass burned in recent decades, although the increase
is not reflected in the regional loess smoother (Figure 5.7). The modeled background charcoal
intensities are consistent with composite CHAR records (i.e., mean charcoal accumulation rate
among lakes) calculated using the Yukon Flats data (Kelly et al., 2013).
The univariate Bayesian point process model provides a model-based approach to estimate
probability of fire values associated with sediment charcoal records from a single lake and convert
those probabilities into mean FRI estimates. The multivariate model allows for correlation among
lakes in the parameters used to estimate the background charcoal deposition intensity. The background intensity reflects regional charcoal sources and exhibits low-frequency changes over time
associated with factors such as species composition and climate. As such, background charcoal
deposition should be similar among lakes in the same region with a high potential for correlation

105

in background deposition process parameters. As expected, mean FRI estimates for each lake
are similar based on the univariate and multivariate models (Table 5.1). The multivariate model,
however, resulted in mean FRI estimates with reduced uncertainty. Specifically, the 95 percent
credible interval for the mean FRI was narrower in 10 out of 13 lakes in the Yukon Flats network
with a mean credible interval width of 156 years for the multivariate model versus 168 years for
the univariate model. The reduced uncertainty in mean FRI estimates under the multivariate model
provides some evidence that background charcoal deposition is indeed correlated among lakes
located in the same region and that we can reduce uncertainty in estimates of background charcoal deposition by accounting for such correlation. The effective spatial range for the background
(b)

process regression parameters (β j ) was estimated to be 16.5 km and provides some indication of
the distance within which background charcoal deposition is similar among lakes within the Yukon
Flats region. This estimate is consistent with the previous analysis using the Yukon Flats data set,
which found significant correlation between composite CHAR and regional area burned within a
20-km radius (Kelly et al., 2013).
We estimated a regional mean FRI of roughly 187 years (95 percent credible interval: 136
to 261 years) for the Yukon Flats over the study period applying the partial pooling approach
described in Section 5.3.2.1. The partial pooling approach also produces estimates of mean FRI
values for individual lakes similar to the univariate and multivariate results presented in Table 5.1.
However, we do not see the same reduction in uncertainty in individual-lake mean FRI estimates
when conducting partial pooling. Specifically, credible interval widths were narrower in only 6
out of 13 lakes with the remaining intervals comparable to univariate model results. The partial
2 ) and combines uncertainty in mean FRI
pooling approach adds two additional parameters (α∗, σfri

values across lakes. As such, it is not surprising the partial pooling approach does not lead to
the same reductions in uncertainty as generating individual-lake mean FRI estimates based on the
multivariate model. We envision the partial pooling approach being applied only in the setting
where a researcher is interested in estimating a regional mean FRI, otherwise, the individual-lake
approach is preferred.

106

The univariate point process model we developed achieves our goal of developing an integrated
statistical framework for local fire identification and estimation of mean FRIs based on sediment
charcoal records for individual lakes. The Bayesian hierarchical model structure allows for tractable
propagation of additional uncertainty sources in paleo-fire reconstructions. In particular, uncertainty in sediment age models can be integrated by treating the ages of sediment core sections as
unobserved, latent variables in the point process model. The multivariate extension of the point process model provides a novel approach for paleo-fire reconstruction applying multiple lake records
to make inferences at both individual-lake and regional scales. Specifically, the multivariate model
provides estimates of individual-lake mean FRIs, a regional mean FRI, and background charcoal
deposition intensity indicative of regional biomass burned. When applied to the Yukon Flats data
set, pooling of individual-lake records under the multivariate model led to reduced uncertainty in
individual-lake mean FRIs. We expect the multivariate model to provide even greater reductions
in uncertainty in individual-lake mean FRI values, and improved estimates of regional parameters
including the regional mean FRI, when applied to larger regional lake networks (i.e., networks with
20 plus lakes), assuming all lakes share a common regional fire regime.

107

CHAPTER 6
CONCLUSIONS

6.1

Research Synthesis

The preceding chapters demonstrate a range of spatio-temporal modeling methodologies applied to forest growth, mortality, and disturbance data. Chapter 5 develops a stand alone model to
reconstruct regional fire regimes using sediment charcoal records. Chapters 2-4 build upon one
another with a common research goal in mind: to advance understanding of interactions between
climate extremes, disturbance, and forest dynamics. Understanding these interactions is a necessary step toward predicting forest responses to future environmental conditions and is critical to
informing adaptive forest management aimed at maintaining forest function under changing climate
and disturbance regimes.
There was evidence in Chapter 2 that tree growth responses to water deficits varied over time
depending on the severity of water stress, their co-occurrence with insect defoliation events, and
forest density. In particular, there was evidence water deficit thresholds might exist below which
tree growth is unaffected, but above which water deficit triggers strong growth reductions. This is
consistent with contemporary physiological theory (McDowell et al., 2011; Allen et al., 2015), and
is important given tree growth is correlated with a range of forest demographic processes including
mortality and is an indicator of forest health (Buechling et al., 2017). Further, the results of Chapter
2 underscored the dependence of tree growth on past water deficits and their potential interaction
with insect defoliation events.
Building on the results of Chapter 2, we developed a modeling framework in Chapter 3 to
quantify the persistent and cumulative effects of water deficit and insect defoliation on current tree
growth applying the concept of ecosystem memory. The framework allowed for the interactive
effects of cumulative water and defoliation stress on tree growth to be explicitly tested. The
model framework was applied to two contrasting regions of the Canadian boreal forest with results
108

revealing persistent tree growth responses to both water deficit and insect defoliation lasting up
to a decade. There was evidence of a positive interactive effect of cumulative water deficit and
insect defoliation among non-defoliator-host trees suggesting these trees benefit from reduced
competition for limited water following defoliation events. There was little evidence, however, of
negative interactive effects among host trees even after accounting for persistent and cumulative
tree responses. These results are counter to the physiological hypothesis that feedbacks in tree
responses to drought and insect defoliation exacerbate their effects on tree growth and mortality
(McDowell et al., 2011; Anderegg et al., 2015). Combined with the lack of evidence for negative
interactions in previous studies (Jactel et al., 2012; Jacquet et al., 2014; Kolb et al., 2016), our
results imply negative interactions between water and defoliation stress on host tree growth may
not be as strong as expected due to insect defoliation offsetting the impacts of water deficits.
Finally, Chapter 4 builds on the methods and results of Chapters 2 and 3 to develop a state space
framework to model the interactions between climate extremes, disturbance, and forest dynamics.
The state space framework is similar to the framework presented in Chapter 2, but critically replaces
the random walk for climate effects with an empirical model of forest dynamics providing joint
estimates of demographic processes (growth, mortality, regeneration). The framework provides a
tool to test alternative functions describing the effects of climate extremes and disturbance on forest
dynamics. In particular, based on the results of Chapter 2, we can test whether there is evidence
to support the use of growth and mortality functions incorporating a threshold for climatic water
deficit, and further, quantify the value of such a threshold. The state space framework allows
for the identification of forest characteristics promoting ecosystem resistance and resilience to
climate extremes and disturbance events. Further, it provides a tool to test alternative management
scenarios to maintain forest health and productivity under variable future conditions in terms and
contexts forest managers understand and can implement within long-term management plans. The
framework is a valuable tool for adaptive forest management and connects fundamental forest
ecology research with the applied management goals described in Chapter 1.

109

6.2

Future Research

The analyses presented here represent a step toward improved understanding of interactions
between climate extremes, disturbance, and forest dynamics. Many important questions, however,
remain unanswered and there are myriad opportunities to continue exploring such interactions in
future work. Several of these opportunities are described below.
Ecosystem memory: The concept of ecosystem memory is important because it allows ecologists to quantify the strength, relative importance, and cumulative effects of past environmental
conditions on current ecosystem function (Ogle et al., 2015). A wide range of ecological processes likely exhibit persistent responses to past conditions. Ecological memory functions offer a
tool to quantify and understand the dynamics of these persistent responses. While frameworks to
quantify memory functions are just emerging in ecology, there is a precedent for such frameworks
in spatial statistics where methods exist to estimate spatially-averaged covariates reflecting the response of a process to environmental conditions within a local neighborhood (Heaton and Gelfand,
2012). There is a terrific opportunity to bring these methods to bear in ecological contexts. Future
methodological research will focus on expanding the spline-based approach to quantify ecosystem
memory introduced in Chapter 3 to spatio-temporal processes. That is, develop a flexible approach
to quantify ecological memory functions in time and space reflecting ecosystem responses to environmental conditions in the recent past within a local neighborhood. An ecosystem memory R
package will ultimately be developed to facilitate the use of such methods in ecology.
Analysis of dead trees: The tree growth analyses presented here are limited to tree rings
collected from live trees. Numerous studies have demonstrated trees suffering drought-induced
mortality exhibit reduced radial growth rates relative to the growth of surviving conspecific trees of
similar size for a number of years prior to death (Wyckoff and Clark, 2002; van Mantgem et al., 2003;
Das et al., 2007, 2016; Berdanier and Clark, 2016). Tree death, in general, is thought to be a slow
process occurring over a number of years and resulting from interactions among multiple factors
including tree age/size, water and nutrient availability, insect and wind damage, fungal infection,
and disease (Franklin et al., 1987; McDowell et al., 2011). The approach to quantify the ecological
110

memory of tree growth to past water and insect defoliation stress presented in Chapter 3 provides a
unique opportunity to explore the factors contributing to tree death. Future research will apply the
modeling approach to tree-ring records collected from both live and recently-dead trees following
drought and insect defoliation events to test for differences in their ecological memory to water
and insect defoliation stress. Tree-ring data exist for white spruce trees suffering drought-induced
mortality in Alberta and spruce-fir trees suffering spruce budworm defoliation-induced mortality
in Quebec in addition to the data used in Chapter 3. Ultimately, the tree-ring based approach to
quantify the ecological memory of tree growth to past drought and insect disturbance events can
be expanded to regional scales using data from the International Tree Ring Data Bank to advance
understanding of the mechanisms underlying tree mortality.
Regional modeling of forest dynamics using SORTIE-ND: The state space framework presented in Chapter 4 offers a unique approach to combine disparate tree-ring and forest inventory
data to advance understanding of interactions between climate extremes, disturbance, and forest
dynamics, and inform adaptive forest management. The results of Chapter 4 represent a proof of
concept that past forest dynamics can be reconstructed using a growth and yield model constrained
by common forest growth and mortality data. Future work will replace the FVS model with a
process-based model of forest dynamics, SORTIE-ND (Pacala et al., 1996). SORTIE-ND is an
individual-tree, distance-dependent model of forest dynamics and affords several benefits over the
FVS model. In particular, SORTIE-ND represents competition at the individual tree scale using well known competition indices (Canham et al., 2004). Further, SORTIE-ND accommodates
unevenly aged, multi-species forests. Finally, the SORTIE-ND model offers increased flexibility,
in terms of modifying underlying growth and mortality functions, and improved computational
efficiency relative to FVS. The SORTIE-ND-driven state space model will be applied to regional
forest data sets for which historic meteorological records and disturbance chronologies exist. This
will enable the identification of forest characteristics promoting ecosystem resistance and resilience
to climate extremes and disturbance events and can be used to inform regional silvicultural systems.
Regional data for northeastern North America spanning northern New England to southern Quebec,

111

Canada has been procured and will be analyzed once the SORTIE-ND-driven state space model
has been developed and validated.
Understanding the role of population genetics in forest demographic processes: There
is a range of individual-scale variability in tree growth, mortality, and fecundity. Accounting
for individual-scale variability is integral to understanding how forest demographic processes are
impacted by climate extremes and disturbance (Clark et al., 2011). Variability in individual-tree
responses to similar growing conditions is commonly attributed to unobserved micro-site factors,
genetics, and tree health. Future research will modify growth, mortality, and regeneration functions
from the SORTIE-ND model to allow for individual random effects driven by tree genotypes. The
model will be applied to a mixed, white oak (Quercus alba L.) experimental forest for which
extensive tree-ring, forest inventory, and genetic data exist for a collection of forest stands (the
population unit) to assess the role of population genetics in shaping demographic responses to
drought and insect disturbance. Ecological memory functions will be applied to account for
persistent and cumulative tree responses to past water stress and insect damage. Further, we
will account for competition among trees using spatially-explicit competition indices and forest
stem maps. The analysis will identify within-population variability in growth, mortality, and
fecundity during and following drought events and insect outbreaks. Recent analyses have applied
a similar approach to understand tree growth responses to climate (Housset et al., 2018), but
have not integrated the concept of ecosystem memory, competition, or modeled joint demographic
responses.

6.3

Concluding Remarks

Models are critical to understanding changes in forest processes given the long timescales
(decades to millennia) over which forests persist. Predictions of future changes to forest ecosystems
can be used to inform current management decisions to maintain the productivity, health, and
function of forest systems or to move them in new directions better adapted to future conditions. The
research presented here sought to combine spatio-temporal statistical models, fundamental forest

112

ecology, and applied forest management in order to advance understanding of forest ecosystems
and facilitate their sustainable management. Future research will continue toward these same
fundamental goals building on the methods and results of this dissertation.

113

APPENDICES

114

APPENDIX A
VARIABLE EFFECTS OF CLIMATE ON FOREST GROWTH IN RELATION TO
CLIMATE EXTREMES, DISTURBANCE, AND FOREST DYNAMICS

A.1

Detailed Model Specification

In this section we provide detailed model specifications for the fixed climate effects (FCE)
and variable climate effects (VCE) models. Model specification includes details on the prior
distributions assigned to all model parameters. We discuss specification of hyperparameter values
in the Bayesian Inference section that follows.

A.1.1

Fixed Climate Effects (FCE) Model

We model annual tree growth increments under the FCE model as follows. Let i index individual
trees (i = 1, . . . , n), j index stands ( j = 1, . . . , k), and t index years (t = 1, . . . , T) where n, k, and T
are the total number of trees, stands, and years, respectively. Let j(i) indicate the stand j in which
the ith tree is located (e.g. j(i) = 3 indicates the ith tree is located in stand 3). Finally, define y to
be the observed growth increment, such that yi,t is the observed radial growth increment for tree i
in year t. Individual tree growth is modeled as,
log(yi,t ) = x0i,t βi + ι j(i),t + i,t
i,t = φi,t−1 + ei,t


2
ei,t ∟ N 0, σpe
where xi,t is a p-dimensional vector of known covariate values corresponding to a set of natural
cubic spline basis functions defined for tree age at p knots, βi is a p-dimensional vector of treespecific spline regression coefficients, ι j(i),t is the additive effect of being located in stand j during
year t, and i,t is the residual error. The residual error is modeled as an AR1 process where φ
is an unknown correlation coefficient and ei,t is a pure error term. Growth increments of zero
115

(i.e. missing rings) represented an extremely small portion of the data (less than 1 percent of all
observed growth increments). We applied an additive factor of 0.001 to growth increments of zero
prior to applying the log transformation. We assume the pure error term follows a zero-centered
2 . The AR1 error model and normality assumption
normal distribution with pure error variance σpe

for ei,t induces the follow distribution for i,t ,
i,t ∟ N 0,

2
σpe

!

1 − φ2



where the i,t ’s are temporally correlated. Specifically, Cor i,t , i,t−h = φ h where h is the lag,
h = 0, . . . , Ti , given tree i is observed for a total of Ti years. We assume the growth of individual trees
within and across stands is independent based on exploratory analysis assessing spatial dependence
(i,t ⊥ i0,t 0 ∀ i , i0, t = 1, . . . , Ti , and t 0 = 1, . . . , Ti0 ). Stand effects are modeled using observed
climate as,
ι j,t = f 0j,t θ + v j,t


v j,t ∼ N 0, τ 2
where f j,t is an l-dimensional vector of observed, standardized climate covariates, θ is an ldimensional vector of regression coefficients, and v j,t is a random error term assumed to follow
a zero-centered normal distribution with inter-annual variance τ 2 . We assume stand effects are
independent with respect to time both within and across stands.
We apply conjugate normal priors for the individual tree-level spline and stand-level climate
regression coefficients,


βi ∼ Np µ β, σβ2 S−1


θ ∟ Nl ¾ θ, σθ2 I
where S is a p × p penalty matrix defined based on knot locations (see details in Appendix A), I is
an l-dimensional identity matrix, and ¾ β , σβ2 , ¾ θ , and σθ2 are hyperparameters. We apply uniform
priors for the AR1 correlation coefficient (ensuring the autoregressive process is stationary) and

116

the pure error and inter-annual standard deviation parameters (Gelman, 2006),
φ ∼ Unif (−1, 1)
σpe ∟ Unif (a1, b1 )
τ ∼ Unif (a2, b2 ) .
The hierarchical model is fully specified after setting hyperparameter values for prior distributions
(¾ β, σβ2, ¾ θ, σθ2, a1, b1, a2, b2 ). We discuss hyperparameter values in the Bayesian Inference section.

A.1.2

Variable Climate Effects (VCE) Model

The climate regression coefficients are treated as state variables in the VCE model and evolve
over time such that a unique set of climate coefficients is estimated for each year in the study
period (θ 1, θ 2, . . . , θ T ). The time-varying process is initialized at θ 0 with all subsequent θ t
values informed by the set of stand effects in a five-year moving window centered on the current year αt = (α1,t−2, . . . , αkt−2,t−2, . . . , α1,t+2, . . . , αkt+2,t+2 )0 where kt is the number of stands
with growth observations in year t. Note, the five-year window is truncated for the first two
and last two years in the study period where data does not exist before and after the current year, respectively. The stand effects are informed by individual tree growth observations
yt = (y1,t−2, . . . , ynt−2,t−2, . . . , y1,t+2, . . . , ynt+2,t+2 )0 where nt is the number of trees with observed
growth increments in year t. The time-varying climate coefficient model structure is represented
graphically in Figure A.1.
The tree-level model is unchanged in the VCE model. The stand-level model is updated to
integrate time-varying climate coefficients. We also add a model for the evolution of climate effects
over time. We apply a first order random walk to model changes in the climate coefficients as

117

follows,
ι j,t = f 0j,t θ t + v j,t

(A.1)

θ t = θ t−1 + wt


v j,t ∼ N 0, τ 2

(A.2)

wt ∟ Nl (0, Σ θ )
where wt is an l-dimensional process error vector following a zero-centered multivariate normal
distribution with covariance matrix Σ θ . Equations A.1 and A.2 define a state space or dynamic
linear model framework with (A.1) serving as the observation equation and (A.2) the process or
state equation (West and Harrison, 1997).
θ0

...

θ t−1

θt

θ t+1

Îą1

αt−1

Îąt

y1

yt−1

yt

θ1

...

θT

Climate Coefficients

Îąt+1

ÎąT

Stand Effects

yt+1

yT

Tree Growth Data

Figure A.1: Graphical depiction of climate coefficient evolution and variable climate effects (VCE)
model dependencies with arrows indicating direct dependence.
We apply a conjugate normal prior to initialize the climate coefficient process and a scaled
inverse-Wishart prior for the process error covariance matrix (as described in Gelman et al. 2014),
θ 0 ∟ Nl (m0, C0 )
Σ θ = Diag(Ρ)Σ(Ρ)Diag(Ρ)
Σ(Ρ) ∟ Inv-Wish (ν, V)
iid

Ρs ∟ Unif(a3, b3 )

s = 1, . . . , l

where Ρ is an l-dimensional vector of scale parameters each assigned a uniform prior distribution,
Diag(Ρ) indicates an l-dimensional diagonal matrix with values of Ρ on the diagonal, and Σ(Ρ)
is an l-dimensional correlation matrix. All other priors remain the same as in the FCE model.
Again, the model is fully specified after setting hyperparameter values for prior distributions
(m0, C0, ν, V, a3, b3 ).
118

A.2

Bayesian Inference

We use Markov chain Monte Carlo (MCMC) algorithms to estimate model parameters using
Gibbs sampling (Gelfand and Smith, 1990). Specific implementation of our Gibbs sampler is
discussed in the following section. We apply diffuse normal prior distributions for the FCE
regression parameters setting ¾ β and ¾ θ equal to zero, σβ2 equal to 250, and σθ2 equal to 1e7/k.
The scalar variance hyperparameter for βi serves as a tuning parameter controlling the degree of
smoothness of the penalized spline regression model. We selected a value of 250 based on a gridded
search matching as closely as possible penalized spline regression residuals with residuals from a
relatively stiff smoothing spline representative of those applied in conventional dendrochronology
analyses (i.e. 50 percent frequency response; see Penalized Spline Regression section below).
We apply non-informative uniform priors for the pure error (σpe ) and inter-annual (τ) standard
deviations, but constrain values to an appropriate order of magnitude based on exploratory analysis
setting a1 and a2 to 0.0001 and b1 and b2 to 10.
The same hyperparameter values for σβ2, a1, b1, a2 and b2 are used in the VCE model. The
climate coefficient evolution is initialized using a diffuse normal distribution with m0 = 0 and
C0 = Diag(1e3). The temporally-varying climate coefficient estimates are updated using a Kalman
filter nested within a Gibbs sampler. We set ν = l + 1 and V = I (per Gelman et al. 2014). Finally,
we set a3 = 0 and b3 = 1 constraining the ηs ’s to an appropriate order of magnitude based on
exploratory analysis.

A.3

Gibbs Sampler

We apply Gibbs sampling, a Markov chain Monte Carlo (MCMC) technique, to sample the
joint posterior distribution for all model parameters (Gelfand and Smith, 1990). Here we present
conditional distributions for each parameter and the algorithms used to sample from them. We
begin by defining additional notation used throughout the appendix. All other notation is as defined
in article.

119

A.3.1

Notation
tr.yr s =

T
Õ

nt

t=1

std.yr s =

T
Õ

kt

t=1


0
yi = yi,1, . . . , yi,Ti

i = 1, . . . , n


0
y = y01, . . . , y0n

0
Îą = Îą01, . . . , ÎąT0

0
0
0
β = β 1, . . . , β n
0

Θ = θ 01, . . . , θ T0
A = tr.yr s × std.yr s incidence matrix
n X
X = ⊕i=1
i

where Xi = Ti × p design matrix

FFCE = tr.yr s × l design matrix
FVCE =

T F
⊕t=1
t

where Ft =

t+2
Õ

kt × l design matrix

t−2
n ÎŁ
Σ y = ⊕i=1
i
2
σres

=

where Σi = ni × ni AR1 correlation matrix

2
σpe

1 − φ2

Notes: ⊕ indicates a direct sum forming a block diagonal matrix. The matrix A selects the
corresponding scalar stand effect for each yi,t . The elements of X are known knot values defined for
each tree based on age. The elements of FFCE and FVCE are known standardized climate covariate
values. FVCE is truncated for the first two and last two years in the study period where data previous
to and after the current year do not exist. The ta : tb notation indicates all time points from ta to
tb . Finally, we use Ω to indicate all model parameters and Ω(−) to indicate all model parameters
120

minus the parameter currently being updated.

A.3.2

Climate Coefficients

Under the fixed climate effects (FCE) model, climate coefficients (θ) are sampled from a normal
conditional posterior distribution as follows.
θ |y, Ω(−) ∼ N (V v, V)
where
V=
v=

 0
FFCE FFCE
τ2

F0FCE Îą
τ2

+ C−1
0

 −1

+ m0 C−1
0 .

The state space framework of the variable climate effects (VCE) model sets up a dynamic linear
model for the climate coefficients and stand effects at time t as defined by the following observation
and process equations
Obs:

ιt = Ft θ t + vt



vt ∼ N 0, τ 2 I kt

Proc:

θ t = θ t−1 + wt

wt ∟ N (0, Σ θ )

initialized at t = 0 using θ 0 ∟ N (m0, C0 ). Climate coefficient values are sampled from their
joint conditional posterior distribution N (θ 0:T |y, Σ θ ) using the forward filtering backward sampling (FFBS) algorithm (Carter and Kohn, 1994), which applies the Kalman filter and smoothing
equations. The following steps define the FFBS algorithm (adapted from Petris et al., 2009):
(1) update filtering distribution for θ t (t = 1, . . . , T) using θ t |y1:t ∟ N (mt , Ct ) where
(a)

(b)

mt = at + Rt F0t Qt−1 (αt − ft )

at = mt−1

Ct = Rt − Rt F0t Qt−1 Ft Rt

Rt = Ct−1 + Σ θ
121

(c)
ft = Ft at
Qt = Ft Rt F0t + τ 2 I kt

(a) filtering parameters; (b) one-step ahead prediction parameters for θ t ; (c) one-step ahead
prediction parameters for yt .
(2) sample θ T from its posterior conditional distribution θ T |y1:T , Σ θ ∟ N (mT , CT )
(3) sample θ t from its conditional smoothing distribution θ t |θ t+1, y1:T , Σ θ ∟ N (st , St ) for
(t = T − 1, ..., 1, 0) where
−1 (θ
st = mt + Ct Rt+1
t+1 − at+1 )
−1 C .
St = Ct − Ct Rt+1
t

We use a singular value decomposition (SVD) form for the covariance matrices (Rt , Ct , St ) included
in the FFBS algorithm to improve numerical stability and ensure that each matrix is symmetric
positive-definite. Specifically, we apply the SVD version of the Kalman filter described in Wang
et al. (1992) to complete FFBS steps 1-2. The SVD approach is extended to calculate St relevant
to step 3 as follows:


(L∗ )0 U+ 
0

t 
+ +
(1) form Pt = 
 where L∗ L∗ = Σ−1


θ and Ut , Dt are as defined in Wang et al. (1992);
 D+ −1 

 t



(2) solve for the SVD of Pt retaining the singular values Dt◦ and right singular vectors Vt◦ ;


(3) form Utb = Ut+ Vt◦ and Dtb = Dt◦ −1 ;




(4) then, St = Utb Dtb 2 Utb 0 and st = mt + St Σ−1
θ (θ t+1 − at+1 ).

A.3.3

Stand Effects

Stand effects under the FCE and VCE models (Îą) are sampled from conditional normal distributions
as follows.
α|y, Ω(−) ∼ N (V v, V)

122

where under the FCE model,
A0Σ−1
y A

V=

v=

2
σres

I
+
τ2

! −1

A0Σ−1
y (y − Xβ)
2
σres

+

FFCE θ
.
τ2

FFCE and θ are replaced by FVCE and Θ in the formula for v under the VCE model.

A.3.4

Smoothing Spline Regression Coefficients

The individual tree penalized smoothing spline regression parameters (β) are sampled from the
same conditional posterior normal distribution under the FCE and VCE models as follows.
β|y, Ω(−) ∼ N (V v, V)
where
0 −1

−1

S ÂŞ
Š X Σy X
+
V =­
ÂŽ
2
σres
σβ2
ÂŤ
ÂŹ
v=

A.3.5

X0Σ−1
y (y − Aα)
2
σres

+

¾ βS
σβ2

Autoregressive Error Process Parameters

The pure error standard deviation (σpe ) and correlation coefficient (φ) parameters corresponding to
the first order autoregressive process (AR1) used for the individual tree residual errors are updated
using a Metropolis step applying normal proposal densities for each variable centered on the current
parameter value and proposal variance tuned to achieve an acceptance rate of approximately 44


 x−a 	
percent. We propose a logit-transformed g(x) = log b−x for x ∈ (a, b) value for the pure error
standard deviation parameter such that it has support equal to the whole real line. The normal
proposal density for the correlation coefficient is truncated such that is has support -1 to 1. The log
123

target density function used to update the AR1 parameters is as follows,



 1
tr.yr s
2
l σpe, φ|y, Ω(−) = −
log σres
+ log|Σ−1
y |−
2
2
1
(y − Aα − Xβ)0 Σ−1
y (y − Aα − Xβ) +
2
2σres




log σpe − a1 + log b1 − σpe
with the last two terms equal to the Jacobian correction applied to back transform the proposed
value for σpe .

A.3.6

Inter-annual Standard Deviation

The inter-annual standard deviation (τ) is updated using a Metropolis step with a normal proposal
density centered on the current value of τ and variance tuned to achieved an acceptance rate of
approximately 44 percent. We again propose a logit-transformed value for τ such that is has support
equal to the whole real line. The log target density function used to update τ in the FCE model is
as follows.
 


std.yr s
log τ 2 −
l τ|y, Ω(−) = −
2
1
(α − FFCE θ)0 (α − FFCE θ) +
2
2τ
log (τ − a2 ) + log (b2 − τ)
For the VCE model FFCE and θ are replaced with FVCE and Θ.

A.3.7

Random Walk Process Error

The process error variance for the variable climate effects under the VCE model (Σ θ ) is given a
scaled inverse-Wishart prior. As such, its value is sampled by updating the scale parameters (ηs ’s)
and correlation matrix (Σ(Ρ)) separately. The scale parameters are updated using independent
Metropolis steps proposing logit-transformed values from a normal proposal density centered on
the current parameter value. We apply an adaptation step to achieve a target acceptance rate of
124

44 percent (Roberts and Rosenthal, 2009). The log target density function used to update Ρs for
s = 1, . . . , l is as follows,
"
#
1
1
l ηs |y, Ω(−) = −T log (ηs ) −
SS(θ)s,s Σ(η)−1
s,s −
2 Ρs2


l 

1  1 Õ 1
−1
+
Σ(Ρ)
SS(θ)
0
0
s,s
s,s 
2  Ρs 0
Ρ s0


 s ,s




log (ηs − a3 ) + log (b3 − ηs )
where SS(θ) =

ÍT

t=1 (θ t

− θ t−1 ) (θ t − θ t−1 )0 and Zi, j = Z[i, j], the i, jth element of the matrix Z.

The correlation matrix (Σ(Ρ)) is updated using a Gibbs step sampling directly from an inverseWishart conditional posterior distribution as follows,
Σ(η)|y, Ω(−) ∼ Inv-Wish(ν ∗, V∗ )
where
ν∗ = ν + T
V∗ = V + Diag(η)−1 SS(θ)Diag(η)−1
where the inverse-Wishart distribution is defined as in Gelman et al. (2014).

A.4

Penalized Spline Regression

Smoothing splines are a highly flexible method to model natural phenomena using a set of
polynomial basis functions (Wood, 2006). Alternatives to a smoothing spline include a negative
exponential regression model or a generalized additive model (Cook, 1987; Fajardo and Mcintire,
2012). We apply a penalized spline regression model instead of a smoothing spline as it achieves
equivalent smoothing with fewer model parameters (Wood, 2006). Unlike a smoothing spline in
which knots are defined for every observation, in a penalized spline regression model, a set of
knots less than or equal to the number of observations is specified and regression coefficients are
penalized such that the coefficient values for non-informative knots are forced to zero.
125

We fit a penalized spline regression model using natural cubic regression splines. That is, we
estimate the value for βi that minimizes the following target function:
∍ 1
 00  2
2
||yi − Xi βi || + λ
f (x) dx
0

∍1
where Ν is the penalty term. As noted in Wood (2006), 0 [ f 00(x)]2 dx can be rewritten as β0i Sβi
where S is a matrix of known coefficients defined as a function of the selected knot values for tree
i. The target function is then given by
||yi − Xi βi || 2 + λβ0i Sβi .
Note that S is not full column rank, rather its rank is equal to the total number of knots defined
minus two given that the second derivatives of the boundary knots are equal to zero for natural
cubic splines. Assuming yi |Ω ∼ N (•) and placing a normal prior on βi , we seek a hyperprior
covariance matrix for βi such that the kernel of the conditional posterior distribution has the form
of the above target function. This is achieved by assigning βi ∼ N (0, σβ2 S−1 ). Although S is not
invertible, i.e. the prior placed on βi is improper, we obtain a valid posterior distribution.
The value of σβ2 serves as a tuning parameter controlling the level of smoothing. We determined
an optimal value of 250 for σβ2 based on a gridded search matching as closely as possible the residuals
from a classic dendrochronology smoothing spline. Specifically, we fit 50 percent frequency
response smoothing splines to each individual tree growth record using the dplR package for R
(Bunn, 2008). We then calculated the correlation and spectral coherence between the smoothing
spline and penalized spline regression residuals for each tree varying σβ2 across a range of values.
We selected the value of σβ2 that maximized the mean correlation and mean coherence across all
sample trees.

A.5

Bayesian Lasso

The Bayesian Lasso (short for: least absolute shrinkage and selection operator) is applied
within the stand effects sub-model of the FCE model using the Gibbs sampling approach presented
in Park and Casella (2008) to identify the most important climate variables from a set of 28
126

potential climate covariates. The set of climate variables identified was then used in both the FCE
and VCE models. Implementation of the Bayesian Lasso requires alternative priors for the climate
coefficients (θ) and cannot be applied within the VCE model framework since climate coefficients
vary over time through a first order random walk initialized by specifying a prior for θ 0 . We note
variable selection would not be necessary if the Bayesian Lasso could be implemented within the
VCE model. In this case, the Bayesian Lasso would appropriately shrink the coefficient values of
unimportant climate covariates to zero in both models. This would allow us to compare whether
the same set of climate variables are important under the FCE model versus VCE model. Given
the Bayesian Lasso is not applicable to the VCE model, we apply it to the FCE model, select the
most important climate variables, then use this reduced set of variables in both the FCE and VCE
models throughout the analysis.
Although the Bayesian Lasso can accommodate collinear variables, pairs of variables that are
highly correlated may still lead to erroneous coefficient values (Hooten and Hobbs, 2015). The
climatic water deficit is a linear function of potential and actual evapotranspiration (DEF = PET
- AET). This results in a singular design matrix (F) if all three variables (DEF, PET, AET) are
included in the model. Thus, we apply the Bayesian Lasso twice: once using all climate variables
except PET seasonal aggregations; and, once using all climate variables except AET seasonal
aggregations. In both cases, seasonal aggregations of climatic water deficit and mean annual snow
pack are identified as the most important variables (Figure A.2).

A.6

Variance Partitioning

Posterior variance estimates for the FCE and VCE models are provided in Table A.1. Notably,
the first-order autocorrelation coefficient under both models was roughly 0.37 indicating that tree
ring records were moderately autocorrelated even after detrending. The individual-tree variance
was approximately 0.29 under both models, while the inter-annual variance (capturing stand-level
variability) was approximately 0.05 under the FCE model, and 0.04 under the VCE model, all on
the log scale. The individual-tree growth variance was roughly six times the stand-level growth

127

DEF − Fall
DEF − Spring
DEF − Summer
DEF − Summer Lag
PET − Fall
PET − Spring
PET − Summer
PET − Summer Lag

Climate Variable

SNOW
Tmax − Fall
Tmax − Spring
Tmax − Summer
Tmax − Summer Lag
Tmax − Winter
Tmean − Fall
Tmean − Spring
Tmean − Summer
Tmean − Summer Lag
Tmean − Winter
Tmin − Fall
Tmin − Spring
Tmin − Summer
Tmin − Summer Lag
Tmin − Winter
−0.1

0.0

0.1

Posterior Coefficient Value

Figure A.2: Results of application of the Bayesian Lasso to full set of climate variables excluding
seasonal aggregations of actual evapotranspiration (AET). Results are similar when seasonal aggregations of AET are included and seasonal aggregations of potential evapotranspiration (PET)
are excluded. Points represent the posterior median coefficient value for each climate variable with
95 percent credible intervals (CIs) indicated with solid lines. A dashed line corresponding to a
coefficient value of zero is provided for reference. Points and lines are color coded according to
whether 95 percent CIs overlap zero: green - the lower bound of the 95 percent CI is greater than
zero, red - the upper bound of the 95 percent CI is less than zero, grey - the 95 percent CI includes
zero.
variance on the log scale under both the FCE and VCE models. The large tree-level variance relative
to stand-level variance in both the FCE and VCE models is consistent with previous analyses using
the same dataset (Foster et al., 2016).

128

Table A.1: Posterior distribution summary for model variance terms. Posterior mean values are
provided with 95 percent credible intervals in parentheses. FCE = fixed climate effects model, VCE
2
σpe

.
= variable climate effects model. Note, σy2 =
1−φ2
Parameter
FCE
φ
0.369 (0.362, 0.376)

VCE
0.369 (0.362, 0.376)

2
σpe

0.248 (0.246, 0.250)

0.247 (0.245, 0.249)

σy2

0.287 (0.283, 0.290)

0.286 (0.283, 0.290)

τ2

0.050 (0.047, 0.054)

0.042 (0.039, 0.045)

129

APPENDIX B
BOREAL TREE GROWTH EXHIBITS DECADAL-SCALE ECOSYSTEM MEMORY TO
DROUGHT AND INSECT DEFOLIATION, BUT NO NEGATIVE RESPONSE TO
THEIR INTERACTION

B.1

Model Specification

B.1.1

Tree-Level Submodel

A single linear term for tree diameter in the previous year was sufficient to control for tree size
in the eastern region. Tree growth in the western region, however, differed by species (trembling
aspen versus white spruce) and showed evidence of following a sigmoid growth function with
respect to diameter. As such, the tree-level submodel in the West included species-specific linear
and quadratic terms for tree diameter in the year previous to growth (xi j hs (t)β in the tree-level
submodel → xi j hs (t)β1 + xi j hs (t)2 β2 where h distinguishes trembling aspen [host] from white
(h)

(h)

spruce [non-host]).

B.1.2

Antecedent Weights

Climatic Water Deficit
Penalized regression splines were used to estimate antecedent weights for past climatic water deficit
(Section 3.3.3). We used cubic spline functions defined for k = 10 equally-spaced knots between 0
and L = 10 years in the past calculated using the mgcv package in R (Wood, 2006; R Core Team,
2016). Spline basis functions include an identifiability constraint allowing for the estimation of an
intercept (see Wood, 2006). This constraint results in one fewer spline basis function than there are
knots (p = k − 1). Selection of a smoothing parameter to penalize the regression coefficients for
individual splines is described in Section B.2.

130

Insect Defoliation
We applied a spherical decay function to estimate antecedent weights for past insect defoliation
events. Specifically,
       3



 1 − 32 φx + 12 φx

g(x; φ) =


 0


0<x≤φ

,

x>φ

where x is the number of years since a defoliation event and φ is the estimated range parameter. A
spherical decay function was selected given it achieves a value of zero within a finite interval of
time (i.e., x > φ).

B.2

Bayesian Model Implementation

We assigned normal priors for all regression coefficients within the tree and stand submodels.
(H)

Specifically, each of the following regionally-specific tree-level parameters (β [East], β1
(NH)

[West], β2

(L-H)

, Îł3

β1
Îł3

(H)

(L-NH)

(S-H)

, Îł3

(NH)

[West], β2

(S-NH)

, Îł3

(H)

(NH)

[West]) and stand-level parameters (Îł0 , Îł1 , Îł1

[West],

(L)

(S)

, Îł2 , Îł2 ,

) where “H” indicates host trees, “NH” indicates non-host trees,

“L” indicates large-diameter trees, and “S” indicates small-diameter trees, were assigned a N(0, σ02 )
prior. σ02 is fixed at a large value (1 × 105 ) defining a diffuse normal prior. The cubic spline
basis function coefficients used to model antecedent climatic water deficit weights were assigned
a multivariate normal prior η ∼ N k (0, ση2 S−1 ). Here, S is a k × k matrix of known coefficients
defined as a function of the selected knot values (Wood, 2006). The scalar variance (σΡ2 ) acts as
a smoothing parameter controlling the smoothness of the estimated antecedent weight function.
We assigned σΡ2 a value of 0.1 in the current analysis based on a series of experimental model
runs. The decay parameter for the antecedent defoliation spherical decay function was assigned
a uniform prior distribution φ ∼ Unif(aφ, bφ ) where aφ = 1/15 and bφ = 1. Finally, the tree- and
stand-level standard deviation parameters were assigned identical diffuse uniform prior distributions
σy ∼ Unif(aσ, bσ ) and σα ∼ Unif(aσ, bσ ) where aσ = 1 × 10−4 and bσ = 100.
Combining the likelihood for observed radial growth increments and latent time-varying stand
effects with the above priors, the joint posterior distribution for the combined tree and stand
131

submodels is proportional to (using notation similar to Gelman et al., 2014)

hs (t)
T N(t)
2 Ö
2 n jÖ
Ö
ÖÖ

N(log(yi j hs (t))|α j hs (t), β, σy2 )×

t=1 j=1 h=1 s=1 i=1
T N(t)
2 Ö
2
Ö
ÖÖ

N(α j hs (t)|γ, η, φ, σα2 )×

(B.1)

t=1 j=1 h=1 s=1

N(β|σ02 ) × N(γ0 |σ02 ) ×

2
Ö
h=1

(h)
N(γ1 |σ02 ) ×

2
Ö
s=1

(s)
N(γ2 |σ02 ) ×

2 Ö
2
Ö
h=1 s=1

(hs)

N(Îł3

|σ02 )×

N k (η|ση2, S) × Unif(φ|aφ, bφ ) × Unif(σy |aσ, bσ ) × Unif(σα |aσ, bσ ),
where T is the total number of years in the study period, N(t) is the number of study stands
for data exist in year t, n j hs (t) is the number of trees in stand j of species class h (host/non-host)
and diameter class s (large/small) for which radial growth increment observations exist in year
(H)

(NH)

t, and γ ≡ (γ0, γ1 , γ1

(L)

(S)

(L-H)

, Îł2 , Îł2 , Îł3

(L-NH)

, Îł3

(S-H)

, Îł3

(S-NH) 0
).

, Îł3

All other values are as

previously defined. Note, for the western study region the single β is replaced with species-class
(H)

(NH)

specific linear and quadratic coefficients (see β1 , β1

(H)

(NH)

, β2 , β2

above).

We used a Metropolis-within-Gibbs Markov chain Monte Carlo (MCMC) algorithm (Robert
and Casella, 2004) to sample from the posterior distribution in (B.1) written in R (R Core Team,
2016). For each study region, three chains were run for a total of 5,000 iterations following a 5,000
sample burn-in period. Convergence was assessed visually and using Gelman-Rubin statistics.

132

B.3

Supplementary Results

Table B.1: Posterior summary of tree size effect coefficients and tree- and stand-level variances
for the East and West study regions. Posterior median values are reported along with 95 percent
credible intervals in parentheses. Note a single linear term is used to control for tree size in the
East, whereas host/non-host specific linear and quadratic terms are used in the West.
Parameter
Îł0
β
(H)
β1
(NH)
β1
(H)
β2
(NH)
β2

σy2
σι2

Description

East
West
Stand Effect Coefficients
Intercept
-0.561 (-0.637, -0.483) -0.093 (-0.139, -0.048)
Tree Size Effects
All Trees Linear
0.087 (0.083, 0.091)
NA
Host Tree Linear
NA
0.189 (0.185, 0.194)
Non-Host Tree Linear

NA

0.12 (0.115, 0.126)

Host Tree Quadratic

NA

-0.005 (-0.005, -0.005)

Non-Host Tree Quadratic

NA
Variance
0.324 (0.317, 0.331)
0.236 (0.213, 0.262)

-0.003 (-0.003, -0.002)

Tree-Level
Stand-Level

133

0.25 (0.247, 0.254)
0.35 (0.333, 0.367)

Modeled Climatic Water Deficit
Posterior Mean Antecedent Value
95% Credible Interval

Figure B.1: Antecedent climatic water deficit values for the 14 eastern (Quebec) study stands over the study period (1968-1998) estimated
by applying posterior samples of antecedent weights to modeled climatic water deficit.

134

Modeled Climatic Water Deficit
Posterior Mean Antecedent Value
95% Credible Interval

Figure B.2: Antecedent climatic water deficit values for the 34 western (Alberta) study stands over the study period (1968-2010) estimated
by applying posterior samples of antecedent weights to modeled climatic water deficit.

135

Observed Defoliation Event
Posterior Mean Antecedent Value
95% Credible Interval

Figure B.3: Antecedent defoliation values for the 14 eastern (Quebec) study stands over the study period (1968-1998) estimated by
applying a spherical decay function using posterior samples of the decay parameter (φ) to observed defoliation events.

136

Observed Defoliation Event
Posterior Mean Antecedent Value
95% Credible Interval

Figure B.4: Antecedent defoliation values for the 34 western (Alberta) study stands over the study period (1968-2010) estimated by
applying a spherical decay function using posterior samples of the decay parameter (φ) to observed defoliation events.

137

APPENDIX C
ASSIMILATION OF TREE-RING AND FOREST INVENTORY DATA TO MODEL
INTERACTIONS BETWEEN CLIMATE AND FOREST DYNAMICS

138

C.1

Supplementary Results

A
Unconstrained FVS
Posterior Mean

Trees/acre

95% Credible Interval
Forest inventory observations

Basal area/acre (ft2)

B

Quadratic mean diameter (in)

C

Year

Figure C.1: Posterior summary of stand-scale variables including trees per acre (A), basal area
per acre (B), and quadratic mean diameter (C) relative to forest census observations relative to
an unconstrained run of the Lake States Variant of the USDA Forest Service Forest Vegetation
Simulator (FVS) growth and yield model for a red pine plantation started from bare ground.

139

C.2

Bayesian State Space Model Implementation

The assumption of normally-distributed observation errors (vt ) induces a normal likelihood
function for the data (y1:T where “1 : T” indicates all time points between 1 and T). We assigned
a normal prior for the initial state vector θ with mean m0 and covariance matrix C 0 . Consistent
with the initialization of state space frameworks we set the prior mean m0 equal to a rough
estimate of the mean tree diameter at breast height (DBH) and tree density at the start of the
study period (approximated through linear and spline-based interpolation). The prior covariance
matrix has a diagonal structure with large variance values along the diagonal (1E+05) reflecting low
confidence in the prior mean state values. The assumption of normally-distributed process errors
(wt ) induces a normal likelihood for all state vectors following the initial time point (θ 1:T ). All
2 ,τ 2 ) were assigned inverse gamma prior distributions.
scalar variance parameters (σd2 ,σd2 ,σn2 ,τD
N
tr

fc

fc

Process and observation variance parameters for DBH and tree density are difficult to identify
given similar likelihoods for the data and latent state variables combined with a large amount of
missing observations (unconstrained state estimates) without the use of informed prior distributions.
Moment matching was used to select shape and scale hyperparameters for each variance term
based on prior mean values consistent with previous analyses using tree-ring and diameter-tape
measurements and shape values reflecting the amount of weight to place on observations (Clark
et al., 2007). The shape and scale hyperparameter values for observation and process variances are
presented in Table C.1.
Table C.1: Shape (a x ) and scale (b x ) hyperparameter values for inverse gamma prior distributions
for each observation and process variance parameter (x indicates the variable: dtr , dfc , nfc , D, N).
Variance
Parameter
σd2
tr
σd2
fc
σn2
fc
2
τD
τN2

ax

bx

5 ¡ `dtr

0.10(adtr − 1)

5 ¡ `dfc

10 ¡ `fc
`/500
T

140

0.39(adfc − 1)
50(anfc − 1)
5(a D − 1)
10(a N − 1)

Combining the data and latent state vector likelihood functions with the prior distributions for
the initial state vector and variance parameters, the joint posterior distribution for the state space
model is proportional to (using notation similar to Gelman et al., 2014)
T
Ö
t=1

N(yt |At θ t , Vt ) ×

T
Ö

N(θ t | f (θ t−1 |Ω), W)×

t=1

N(θ 0 |m0, C0 ) × InvGamma(σd2 |adtr, bdtr ) × InvGamma(σd2 |adfc, bdfc )×
tr
fc

(C.1)

2 |a , b )×
InvGamma(σn2 |anfc, bnfc ) × InvGamma(τD
D D
fc

InvGamma(τN2 |a N , b N ),
where Vt and W are observation and process block-diagonal covariance matrices, respectively,
defined as

and

 2

σd I` (t)

 tr dtr



,
Vt = 
σd2 I` (t)

fc dfc




2


σ
d

tr 



τ 2 I

`


W=  D
,

2
τN 



where I p is an identity matrix of dimension p and `¡ (t) is the number of diameter estimates based
on tree-rings, diameter-tape measurements, or tree density estimates in year t.
We used a Metropolis-within-Gibbs Markov chain Monte Carlo (MCMC) algorithm (Robert and
Casella, 2004) to sample from the posterior distribution in (C.1) written in R (R Core Team, 2016).
Details on the sampling algorithm are provided in Section C.3. For each red pine experimental
stand, three chains were run for a total of 500,000 iterations following a 50,000 sample burn-in
period. Convergence was assessed visually and using Gelman-Rubin statistics.

C.3

Sampling Algorithm

Parameter values were updated within the MCMC sampler applying a Gibbs step based on
the conditional posterior distribution for a parameter wherever conjugate priors were a used
141

2 , τ 2 ). All other parameter values (θ
(θ T , σd2 , σd2 , σn2 , τD
0:T−1 ) were updated with MetropoN
fc
tr
fc

lis steps applying a normal proposal density for each parameter centered on its current value and
proposal variance tuned to achieve an acceptance rate of approximately 23 percent. We present full
conditional distributions where they exist and target log-likelihood functions (L) for all remaining
parameters below. We use “|·” to denote conditional on the data and all model parameters not
including the parameter being updated.

Update state values:

θ t |¡

where,










L ∝ −(θ 0 − m0 )| C−1
(θ 0 − m0 ) − (θ 1 − f (θ 0 |Ω))| W−1 (θ 1 − f (θ 0 |Ω))

0











t=0

L ∝ (θ t − Ss)| S−1 (θ t − Ss) − (θ t+1 − f (θ t |Ω))| W−1 (θ t+1 − f (θ t |Ω)) 1 6 t < T















θ T |¡ ∟ N(Ss, S)
t =T





|

S = (At Vt−1 At + W−1 )−1
|

s = At Vt−1 yt + W−1 f (θ t−1 |Ω)

Update variance parameters:
All variance parameters have inverse gamma conditional posterior distributions given normal
likelihood functions and inverse gamma priors. Conditional posterior shape (a∗ ) and scale (b∗ )
(x)

parameter values are provided for each variance parameter below. Note Ht

indicates an incidence

matrix selecting observations and state values corresponding to variable x in year t. Further, fc
denotes an `fc -dimensional vector containing forest census years; fc(i) indexes the ith census year.

142

σd2 |¡
tr

`d
a∗ = adtr + tr
2
T
1 Õ (dtr )
(d )
∗
(Ht (yt − At θ t ))| (Ht tr (yt − At θ t ))
b = bdtr +
2
t=1

σd2 |¡
fc

`d
a∗ = adfc + fc
2
fc(`fc )
1 Õ
(d )
(d )
∗
b = bdfc +
(Ht fc (yt − At θ t ))| (Ht fc (yt − At θ t ))
2
t=fc(1)

σn2 |¡
fc

`
a∗ = anfc + fc
2
b∗

fc(`fc )
1 Õ
(n )
(n )
= bnfc +
(Ht fc (yt − At θ t ))| (Ht fc (yt − At θ t ))
2
t=fc(1)

2 |¡
τD

a∗ = a D +

`
2
T

1 Õ (D)
(D)
b∗ = b D +
(Ht (θ t − f (θ t−1 |Ω)))| (Ht (θ t − f (θ t−1 |Ω)))
2
t=1

τN2 |·
T
a∗ = a D +
2
T

b∗

1 Õ (N)
(N)
= bD +
(Ht (θ t − f (θ t−1 |Ω)))| (Ht (θ t − f (θ t−1 |Ω)))
2
t=1

143

APPENDIX D
A MODEL-BASED APPROACH TO WILDLAND FIRE RECONSTRUCTION USING
SEDIMENT CHARCOAL RECORDS

D.1

Probability of Fire Identification

The probability of fire for a given sample interval τ j,i as defined in (5.4) can be expressed as
(f)
(f)
β +x(τ j,i )0 β
j
e 0, j
(f)
(f)
(b)
(b)
β +x(τ j,i )0 β
β +x(τ j,i )0 β
j + e 0, j
j
e 0, j

which can be simplified to
1


− β∗ +x(τ j,i )0 β ∗
j
0, j



1+e
(f)

(b)

(f)

(b)

∗ = β
where β0,
− β0, j and β ∗j = β j − β j , thereby proving the identifiability of the probability
j
0, j

of fire P j (τ j,i ). The simplified probability of fire is equivalent to the mean response of a logistic
regression model fit using a binary variable indicating the occurrence of a local fire within a given
sample interval. Specifically, we can express the odds of a local fire in a given sample interval as
λ j, f (τ j,i )
λ j,b (τ j,i )

=

(f)
(f)
β +x(τ j,i )0 β
j
e 0, j
(b)
(b)
β +x(τ j,i )0 β
j
e 0, j

∗ + x(τ )0 β ∗ . Thus, there is a direct connection between the
such that the log-odds are given by: β0,
j,i
j
j

mean probability of fire function defined using Poisson count data in (5.4) and a logistic regression
model used to estimate the probability of fire using Bernoulli observations of fire occurrence/nonoccurrence.

144

D.2

Supplementary Figures

Figure D.1: Uncertainty in local fire identification for Chopper Lake based on the univariate
point process model. Upper panel presents the coefficient of variation (CV) for the mean fire
return interval as a function of the probability of fire threshold. The probability of fire threshold
corresponding to the minimum CV is selected as the optimal threshold value. Lower panel presents
the proportion of posterior samples (1,000 samples, post burn-in) which identify a local fire in year
t applying the optimal probability of fire threshold. The number of fires refers to the posterior
median number of local fires identified over the length of the sediment charcoal record (the 95
percent credible interval is provided in parentheses).

145

Figure D.2: Uncertainty in local fire identification for Screaming Lynx Lake based on the univariate
point process model. Upper panel presents the coefficient of variation (CV) for the mean fire
return interval as a function of the probability of fire threshold. The probability of fire threshold
corresponding to the minimum CV is selected as the optimal threshold value. Lower panel presents
the proportion of posterior samples (1,000 samples, post burn-in) which identify a local fire in year
t applying the optimal probability of fire threshold. The number of fires refers to the posterior
median number of local fires identified over the length of the sediment charcoal record (the 95
percent credible interval is provided in parentheses).

146

Figure D.3: Univariate model results for lakes in Yukon Flats network. Upper panel indicates observed charcoal counts along with the
posterior mean charcoal count (blue line). Middle panel illustrates posterior mean foreground (orange line) and background (black line)
intensities. Lower panel plots posterior mean probability of fire estimates for each observed time interval (black line) along with the
upper and lower bounds of the 95 percent credible interval (gray shading) and the optimal threshold (red line).

147

Figure D.3 (cont’d)

148

Figure D.3 (cont’d)

149

Figure D.4: Multivariate model results for lakes in Yukon Flats network. Upper panel indicates observed charcoal counts along with
the posterior mean charcoal count (blue line). Middle panel illustrates posterior mean foreground (orange line) and background (black
line) intensities. Lower panel plots posterior mean probability of fire estimates for each observed time interval (black line) along with
the upper and lower bounds of the 95 percent credible interval (gray shading) and the optimal threshold (red line).

150

Figure D.4 (cont’d)

151

Figure D.4 (cont’d)

152

BIBLIOGRAPHY

153

BIBLIOGRAPHY

Alaska Fire Service (2016). Alaska Fire History: Fire Perimeters Updated February 2, 2016.
Albers, J., Albers, M., and Cervenka, V. (2014). Forest tent caterpillar fact sheet. Minnesota
Department of Natural Resources.
Allen, C. D., Breshears, D. D., and McDowell, N. G. (2015). On underestimation of global vulnerability to tree mortality and forest die-off from hotter drought in the Anthropocene. Ecosphere,
6(8):1–55.
Allen, C. D., Macalady, A. K., Chenchouni, H., Bachelet, D., Mcdowell, N., Vennetier, M.,
Kitzberger, T., Rigling, A., Breshears, D. D., Hogg, E. H. T., Gonzalez, P., Fensham, R., Zhang,
Z., Castro, J., Demidova, N., Lim, J.-h., Allard, G., Running, S. W., Semerci, A., and Cobb, N.
(2010). A global overview of drought and heat-induced tree mortality reveals emerging climate
change risks for forests. Forest Ecology and Management, 259:660–684.
Anderegg, W. R. L., Berry, J. A., Smith, D. D., Sperry, J. S., Anderegg, L. D. L., and Field, C. B.
(2012). The roles of hydraulic and carbon stress in a widespread climate-induced forest die-off.
Proceedings of the National Academy of Sciences, 109(1):233–237.
Anderegg, W. R. L. and Callaway, E. S. (2012). Infestation and hydraulic consequences of induced
carbon starvation. Plant Physiology, 159(4):1866–1874.
Anderegg, W. R. L., Hicke, J. A., Fisher, R. A., Allen, C. D., Aukema, J., Bentz, B., Hood, S.,
Lichstein, J. W., Macalady, K., Mcdowell, N., Pan, Y., Raffa, K., Sala, A., Shaw, D., Stephenson,
N. L., Tague, C., and Zeppel, M. (2015). Tree mortality from drought, insects, and their
interactions in a changing climate. New Phytologist, 208:674–683.
Anderegg, W. R. L., PlavcovĂĄ, L., Anderegg, L. D. L., Hacke, U. G., Berry, J. A., and Field, C. B.
(2013). Drought’s legacy: multiyear hydraulic deterioration underlies widespread aspen forest
die-off and portends increased future risk. Global Change Biology, 19(4):1188–1196.
Aussenac, G. and Granier, A. (1988). Effects of thinning on water stress and growth in Douglas-fir.
Canadian Journal of Forest Research, 18(1):100–105.
Avery, T. E. and Burkhart, H. E. (2002). Forest Measurements. Waveland Press, Inc., Long Grove,
IL, 5 edition.
Babst, F., Alexander, M. R., Szejner, P., Bouriaud, O., Klesse, S., Roden, J., Ciais, P., Poulter,
B., Frank, D., Moore, D. J. P., and Trouet, V. (2014). A tree-ring perspective on the terrestrial
carbon cycle. Oecologia, 176(2):307–322.
BeguerĂ­a, S. and Serrano, S. M. V. (2017). Speibase v.2.5. http://dx.doi.org/10.20350/
digitalCSIC/8508.

154

BeguerĂ­a, S., Vicente-Serrano, S. M., Reig, F., and Latorre, B. (2014). Standardized precipitation
evapotranspiration index (SPEI) revisited: parameter fitting, evapotranspiration models, tools,
datasets and drought monitoring. International Journal of Climatology, 34(10):3001–3023.
Berdanier, A. B. and Clark, J. S. (2016). Multiyear drought-induced morbidity preceding tree death
in southeastern US forests. Ecological Applications, 26(1):17–23.
Bergeron, Y., Leduc, A., Joyal, C., and Morin, H. (1995). Balsam fir mortality following the
last spruce budworm outbreak in northwestern Quebec. Canadian Journal of Forest Research,
25(8):1375–1384.
Berliner, L. M. (1996). Hierarchical Bayesian time series models. In Hanson, K. M. and Silver,
R. N., editors, Maximum Entropy and Bayesian Methods: Santa Fe, New Mexico, U.S.A.,
1995 Proceedings of the Fifteenth International Workshop on Maximum Entropy and Bayesian
Methods, pages 15–22. Springer Netherlands, Dordrecht.
Betancourt, J., D, B., and P, M. (2004). Ecological impacts of climate change. In Report from
a NEON Science Workshop. August 24-25, 2004, Tuscon, AZ. American Institute of Biological
Sciences, Washingtion, D.C., USA.
Bigler, C. and Veblen, T. T. (2011). Changes in litter and dead wood loads following tree death
beneath subalpine conifer species in northern Colorado. Canadian Journal of Forest Research,
41(2):331–340.
Biondi, F. (1997). Evolutionary and moving response functions in dendroclimatology. Dendrochronologia, 15:139–150.
Biondi, F. (2000). Are climate-tree growth relationships changing in north-central Idaho, USA?
Arctic, Antarctic, and Alpine Research, 32(2):111–116.
Blais, J. R. (1983). Trends in the frequency, extent, and severity of spruce budworm outbreaks in
eastern Canada. Canadian Journal of Forest Research, 13(4):539–547.
Bormann, F. H. and Likens, G. E. (1979). Catastrophic disturbance and the steady state in northern
hardwood forests: A new look at the role of disturbance in the development of forest ecosystems
suggests important implications for land-use policies. American Scientist, 67(6):660–669.
Bottero, A., D’Amato, A. W., Palik, B. J., Bradford, J. B., Fraver, S., Battaglia, M. A., and Asherin,
L. A. (2017). Density-dependent vulnerability of forest ecosystems to drought. Journal of
Applied Ecology, 54(6):1605–1614.
Bradford, J. B. and Palik, B. J. (2009). A comparison of thinning methods in red pine: consequences
for stand-level growth and tree diameter. Canadian Journal of Forest Research, 39(3):489–496.
Brandt, J. (2009). The extent of the North American boreal zone. Environmental Reviews, 17:101–
161.
Brandt, J., Flannigan, M., Maynard, D., Thompson, I., and Volney, W. (2013). An introduction to
Canada’s boreal zone: ecosystem processes, health, sustainability, and environmental issues 1.
Environmental Reviews, 21(4):207–226.
155

Breshears, D. D., Cobb, N. S., Rich, P. M., Price, K. P., Allen, C. D., Balice, R. G., Romme,
W. H., Kastens, J. H., Floyd, M. L., Belnap, J., Anderson, J. J., Myers, O. B., and Meyer, C. W.
(2005). Regional vegetation die-off in response to global-change-type drought. Proceedings of
the National Academy of Sciences of the United States of America, 102(42).
Buechling, A., Martin, P. H., and Canham, C. D. (2017). Climate and competition effects on tree
growth in Rocky Mountain forests. Journal of Ecology, 105(6):1636–1647.
Bunn, A. G. (2008). A dendrochronology program library in R (dplR). Dendrochronologia,
26(2):115–124.
Calder, C., Lavine, M., MĂźller, P., and Clark, J. S. (2003). Incorporating multiple sources of
stochasticity into dynamic population models. Ecology, 84(6):1395–1402.
Canham, C. D., LePage, P. T., and Coates, K. D. (2004). A neighborhood analysis of canopy
tree competition: effects of shading versus crowding. Canadian Journal of Forest Research,
34(4):778–787.
Carcaillet, C., Bergeron, Y., Richard, P. J., FrĂŠchette, B., Gauthier, S., and Prairie, Y. T. (2001).
Change of fire frequency in the eastern Canadian boreal forests during the Holocene: does
vegetation composition or climate trigger the fire regime? Journal of Ecology, 89(6):930–946.
Carrer, M. and Urbinati, C. (2006). Long-term change in the sensitivity of tree-ring growth to
climate forcing in Larix decidua. New Phytologist, 170(4):861–872.
Carter, C. K. and Kohn, R. (1994). On Gibbs sampling for state space models. Biometrika,
81(3):541–553.
Churchill, D. J., Larson, A. J., Dahlgreen, M. C., Franklin, J. F., Hessburg, P. F., and Lutz, J. A.
(2013). Restoring forest resilience: From reference spatial patterns to silvicultural prescriptions
and monitoring. Forest Ecology and Management, 291:442–457.
Clark, D. A., Brown, S., Kicklighter, D. W., Chambers, J. Q., Thomlinson, J. R., and Ni, J.
(2001). Measuring net primary production in forests: Concepts and field methods. Ecological
Applications, 11(2):356–370.
Clark, J. S. (1988). Particle motion and the theory of charcoal analysis: Source area, transport,
deposition, and sampling. Quaternary Research, 30(1):67–80.
Clark, J. S. (1989). Effects of long-term water balances on fire regime, north-western Minnesota.
Journal of Ecology, 77(4):989–1004.
Clark, J. S. (1990). Fire and climate change during the last 750 yr in northwestern Minnesota.
Ecological Monographs, 60(2):135–159.
Clark, J. S. (2005). Why environmental scientists are becoming Bayesians. Ecology Letters,
8(1):2–14.

156

Clark, J. S., Bell, D., Chu, C., Courbaud, B., Dietze, M., Hersh, M., HilleRisLambers, J., IbĂĄĂąez,
I., LaDeau, S., McMahon, S., Metcalf, J., Mohan, J., Moran, E., Pangle, L., Pearson, S., Salk, C.,
Shen, Z., Valle, D., and Wyckoff, P. (2010). High-dimensional coexistence based on individual
variation: a synthesis of evidence. Ecological Monographs, 80(4):569–608.
Clark, J. S., Bell, D. M., Kwit, M., Stine, A., Vierra, B., and Zhu, K. (2011). Individual-scale
inference to anticipate climate-change vulnerability of biodiversity. Philosophical Transactions
of the Royal Society of London B: Biological Sciences, 367(1586):236–246.
Clark, J. S., Bell, D. M., Kwit, M. C., and Zhu, K. (2014). Competition-interaction landscapes for
the joint response of forests to climate change. Global Change Biology, 20(6):1979–1991.
Clark, J. S., Iverson, L., Woodall, C. W., Allen, C. D., Bell, D. M., Bragg, D. C., D’Amato, A. W.,
Davis, F. W., Hersh, M. H., Ibanez, I., Jackson, S. T., Matthews, S., Pederson, N., Peters, M.,
Schwartz, M. W., Waring, K. M., and Zimmermann, N. E. (2016). The impacts of increasing
drought on forest dynamics, structure, and biodiversity in the United States. Global Change
Biology, 22(7):2329–2352.
Clark, J. S. and Royall, P. D. (1996). Local and regional sediment charcoal evidence for fire regimes
in presettlement north-eastern North America. Journal of Ecology, 84(3):365–382.
Clark, J. S., Royall, P. D., and Chumbley, C. (1996). The role of fire during climate change in an
eastern deciduous forest at Devil’s Bathtub, New York. Ecology, 77(7):2148–2166.
Clark, J. S., Wolosin, M., Dietze, M., IbĂĄĂąez, I., LaDeau, S., Welsh, M., and Kloeppel, B. (2007).
Tree growth inference and prediction from diameter censuses and ring widths. Ecological
Applications, 17(7):1942–1953.
Cook, E. R. (1987). The decomposition of tree-ring series for environmental studies. Tree-Ring
Bulletin, 47:37–59.
Cook, E. R. and Kairiukstis, L. A., editors (1990). Methods of Dendrochronology: Applications in
the Environmental Sciences. Kluwer Academic Publications, Hingham, MA.
Cook, E. R. and Peters, K. (1981). The smoothing spline: a new approach to standardizing forest
interior tree-ring width series for dendroclimatic studies. Tree-ring bulletin, 41:45–53.
Cooke, B. J., Nealis, V. G., and RÊgnière, J. (2007). Insect defoliators as periodic disturbances
in northern forest ecosystems. In Johnson, E. A. and Miyanishi, K., editors, Plant disturbance
ecology: the process and the response, chapter 15, pages 487–525. Elsevier Academic Press,
Burlington, MA.
Cressie, N., Calder, C. A., Clark, J. S., Hoef, J. M. V., and Wikle, C. K. (2009). Accounting
for uncertainty in ecological analysis: the strengths and limitations of hierarchical statistical
modeling. Ecological Applications, 19(3):553–570.
Cressie, N. and Wikle, C. K. (2011). Statistics for Spatio-Temporal Data. John Wiley & Sons,
Hoboken, NJ.

157

Dale, V. H., Joyce, L. A., McNulty, S., Neilson, R. P., Ayres, M. P., Flannigan, M. D., Hanson, P. J.,
Irland, L. C., Lugo, A. E., Peterson, C. J., et al. (2001). Climate change and forest disturbances.
BioScience, 51(9):723–734.
D’Amato, A. W., Bradford, J. B., Fraver, S., and Palik, B. J. (2011). Forest management for
mitigation and adaptation to climate change: Insights from long-term silviculture experiments.
Forest Ecology and Management, 262(5):803 – 816.
D’Amato, A. W., Bradford, J. B., Fraver, S., and Palik, B. J. (2013). Effects of thinning on
drought vulnerability and climate response in north temperate forest ecosystems. Ecological
Applications, 23(8):1735–1742.
Das, A. D., Battles, J. B., Stephenson, N., and van Mantgem, P. (2007). The relationship between
tree growth patterns and likelihood of mortality: a study of two tree species in the Sierra Nevada.
Canadian Journal of Forest Research, 37(3):580–597.
Das, A. J., Stephenson, N. L., and Davis, K. P. (2016). Why do trees die? Characterizing the drivers
of background tree mortality. Ecology, 97(10):2616–2627.
Davis, M. B. and Shaw, R. G. (2001). Range shifts and adaptive responses to Quaternary climate
change. Science, 292(5517):673–679.
DeGrandprĂŠ, L., Kneeshaw, D., Boucher, D., Marchand, M., Pureswaran, D., and Girardin, M.
(2017). Extreme climatic periods preceded spruce budworm-induced tree mortality in eastern
boreal North America. In Review.
Dietze, M. C., Lebauer, D. S., and Kooper, R. (2013). On improving the communication between
models and data. Plant, Cell & Environment, 36(9):1575–1585.
Dietze, M. C. and Moorcroft, P. R. (2011). Tree mortality in the eastern and central United States:
patterns and drivers. Global Change Biology, 17(11):3312–3326.
Diggle, P. J. (2014). Statistical Analysis of Spatial and Spatio-Temporal Point Patterns. CRC Press,
Boca Raton, FL, 3 edition.
Dixon, G. E. (2002). Essential FVS: A User’s Guide to the Forest Vegetation Simulator. USDA
Forest Service, Forest Management Service Center. Revised 2015.
Dixon, G. E. and Keyser, C. E. (2008). Lake States (LS) Variant Overview – Forest Vegetation
Simulator. USDA Forest Service, Forest Management Service Center.
D’Orangeville, L., Côté, B., Houle, D., and Morin, H. (2013). The effects of throughfall exclusion
on xylogenesis of balsam fir. Tree Physiology, 33(5):516–526.
Duchesne, L. and Ouimet, R. (2008). Population dynamics of tree species in southern Quebec,
Canada: 1970-2005. Forest Ecology and Management, 255(7):3001–3012.
Dyer, J. M. (2004). A water budget approach to predicting tree species growth and abundance,
utilizing paleoclimatology sources. Climate Research, 28(1):1–10.

158

Evans, M. E. K., Falk, D. A., Arizpe, A., Swetnam, T. L., Babst, F., and Holsinger, K. E.
(2017). Fusing tree-ring and forest inventory data to infer influences on tree growth. Ecosphere,
8(7):e01889–n/a. e01889.
Fajardo, A. and Mcintire, E. J. B. (2012). Reversal of multicentury tree growth improvements and
loss of synchrony at mountain tree lines point to changes in key drivers. Journal of Ecology,
100(3):782–794.
Fick, S. and Hijmans, R. (2017). Worldclim 2: New 1-km spatial resolution climate surfaces for
global land areas. International Journal of Climatology.
Fleming, R. A. et al. (2000). Climate change and insect disturbance regimes in Canada’s boreal
forests. World Resource Review, 12(3):521–548.
Fleming, R. A. and Piene, H. (1992). Spruce budworm defoliation and growth loss in young balsam
fir: period models of needle survivorship for spaced trees. Forest Science, 38(2):287–304.
Flower, A., Gavin, D., Heyerdahl, E., Parsons, R., and Cohn, G. (2014). Drought-triggered western
spruce budworm outbreaks in the interior Pacific Northwest: a multi-century dendrochronological record. Forest Ecology and Management, 324:16–27.
Floyd, M. L., Clifford, M., Cobb, N. S., Hanna, D., Delph, R., Ford, P., and Turner, D. (2009). Relationship of stand characteristics to drought-induced mortality in three southwestern piñon–juniper
woodlands. Ecological Applications, 19(5):1223–1230.
Foster, J. R., D’Amato, A. W., and Bradford, J. B. (2014). Looking for age-related growth decline in
natural forests: Unexpected biomass patterns from tree rings and simulated mortality. Oecologia,
175(1):363–374.
Foster, J. R., Finley, A. O., D’Amato, A. W., Bradford, J. B., and Banerjee, S. (2016). Predicting
tree biomass growth in the temperate–boreal ecotone: Is tree size, age, competition, or climate
response most important? Global Change Biology, 22(6):2138–2151.
Franklin, J. F., Shugart, H. H., and Harmon, M. E. (1987). Tree death as an ecological process.
BioScience, 37(8):550–556.
Gavin, D. G., Brubaker, L. B., and Lertzman, K. P. (2003). An 1800-year record of the spatial and
temporal distribution of fire from the west coast of Vancouver Island, Canada. Canadian Journal
of Forest Research, 33(4):573–586.
Gavin, D. G., Hu, F. S., Lertzman, K., and Corbett, P. (2006). Weak climatic control of stand-scale
fire history during the late Holocene. Ecology, 87(7):1722–1732.
Gelfand, A. and Smith, A. F. M. (1990). Sampling-based approaches to calculating marginal
densities. Journal Of The American Statistical Association, 85(410):398–409.
Gelfand, A. E. and Ghosh, S. K. (1998). Model choice: A minimum posterior predictive loss
approach. Biometrika, 85(1):1–11.

159

Gelman, A. (2006). Prior distributions for variance parameters in hierarchical models (Comment
on Article by Browne and Draper). Bayesian Analysis, 1(3):515–534.
Gelman, A., Carlin, J. B., Stern, H. S., and Rubin, D. B. (2014). Bayesian Data Analysis. Chapman
& Hall/CRC, Boca Raton, FL, 3 edition.
Girardin, M. P., Ali, A. A., Carcaillet, C., Gauthier, S., HĂŠly, C., Le Goff, H., Terrier, A., and
Bergeron, Y. (2013). Fire in managed forests of eastern Canada: Risks and options. Forest
Ecology and Management, 294:238–249.
Gleason, K. E., Bradford, J. B., Bottero, A., D’Amato, A. W., Fraver, S., Palik, B. J., Battaglia,
M. A., Iverson, L., Kenefic, L., and Kern, C. C. (2017). Competition amplifies drought stress in
forests across broad climatic and compositional gradients. Ecosphere, 8(7).
Gray, D. R. (2008). The relationship between climate and outbreak characteristics of the spruce
budworm in eastern Canada. Climatic Change, 87(3-4):361–383.
Grissino-Mayer, H. (2001). Evaluating crossdating accuracy: a manual and tutorial for the computer
program COFECHA. Tree-Ring Research, 57:205–221.
Hansen, A. J., Neilson, R. P., Dale, V. H., Flather, C. H., Iverson, L. R., Currie, D. J., Shafer, S.,
Cook, R., and Bartlein, P. J. (2001). Global change in forests: Responses of species, communities,
and biomes: Interactions between climate change and land use are projected to cause large shifts
in biodiversity. BioScience, 51(9):765–779.
Harvey, B. J., Donato, D. C., and Turner, M. G. (2016). High and dry: post-fire tree seedling
establishment in subalpine forests decreases with post-fire drought and large stand-replacing
burn patches. Global Ecology and Biogeography, 25(6):655–669.
Heaton, M. J. and Gelfand, A. E. (2012). Kernel averaged predictors for spatio-temporal regression
models. Spatial Statistics, 2(1):15–32.
Hefley, T. J., Broms, K. M., Brost, B. M., Buderman, F. E., Kay, S. L., Scharf, H. R., Tipton,
J. R., Williams, P. J., and Hooten, M. B. (2017). The basis function approach for modeling
autocorrelation in ecological data. Ecology, 98(3):632–646.
Hicke, J. A., Allen, C. D., Desai, A. R., Dietze, M. C., Hall, R. J., Hogg, E. H. T., Kashian, D. M.,
Moore, D., Raffa, K. F., Sturrock, R. N., and Vogelmann, J. (2012). Effects of biotic disturbances
on forest carbon cycling in the United States and Canada. Global Change Biology, 18(1):7–34.
Higuera, P. E., Brubaker, L. B., Anderson, P. M., Hu, F. S., and Brown, T. A. (2009). Vegetation
mediated the impacts of postglacial climate change on fire regimes in the south-central Brooks
Range, Alaska. Ecological Monographs, 79(2):201–219.
Higuera, P. E., Peters, M. E., Brubaker, L. B., and Gavin, D. G. (2007). Understanding the origin
and analysis of sediment-charcoal records with a simulation model. Quaternary Science Reviews,
26(13–14):1790–1809.

160

Higuera, P. E., Whitlock, C., and Gage, J. A. (2010). Linking tree-ring and sediment-charcoal
records to reconstruct fire occurrence and area burned in subalpine forests of Yellowstone National
Park, USA. The Holocene, 19:996–1014.
Hobbs, N. T. and Hooten, M. B. (2015). Bayesian Models: A Statistical Primer for Ecologists.
Princeton University Press, Princeton, NJ.
Hogg, E. H., Brandt, J. P., and Kochtubajda, B. (2002). Growth and dieback of aspen forests in
northwestern Alberta, Canada, in relation to climate and insects. Canadian Journal of Forest
Research, 32(5):823–832.
Hogg, E. H. T., Brandt, J. P., and Michaelian, M. (2008). Impacts of a regional drought on the
productivity, dieback, and biomass of western Canadian aspen forests. Canadian Journal of
Forest Research, 38(6):1373–1384.
Hogg, E. T., Brandt, J. P., and Kochtubajda, B. (2005). Factors affecting interannual variation
in growth of western Canadian aspen forests during 1951-2000. Canadian Journal of Forest
Research, 35(3):610–622.
Holmes (1983). Computer-assisted quality control in tree-ring dating and measurement. Tree-ring
bulletin, 43(1):69.
Hooten, M. B. and Hobbs, N. T. (2015). A guide to Bayesian model selection for ecologists.
Ecological Monographs, 85(1):3–28.
Hooten, M. B., Wikle, C. K., Dorazio, R. M., and Royle, J. A. (2007). Hierarchical spatiotemporal
matrix models for characterizing invasions. Biometrics, 63(2):558–567.
Housset, J. M., Nadeau, S., Isabel, N., Depardieu, C., Duchesne, I., Lenz, P., and Girardin, M. P.
(2018). Tree rings provide a new class of phenotypes for genetic associations that foster insights
into adaptation of conifers to climate change. New Phytologist. doi:10.1111/nph.14968.
Huang, J.-G., Stadt, K. J., Dawson, A., and Comeau, P. G. (2013). Modelling growth-competition
relationships in trembling aspen and white spruce mixed boreal forests of western Canada. PLOS
ONE, 8(10):1–14.
Innes, J. L. and Cook, E. R. (1989). Tree-ring analysis as an aid to evaluating the effects of pollution
on tree growth. Canadian Journal of Forest Research, 19(9):1174–1189.
IPCC (2013). Climate Change 2013: The Physical Science Basis. Contribution of Working Group
I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change. Cambridge
University Press, Cambridge, United Kingdom and New York, NY, USA.
Jacquet, J.-S., Bosc, A., O’Grady, A., and Jactel, H. (2014). Combined effects of defoliation and
water stress on pine growth and non-structural carbohydrates. Tree Physiology, 34(4):367–376.
Jactel, H. and Brockerhoff, E. G. (2007). Tree diversity reduces herbivory by forest insects. Ecology
Letters, 10(9):835–848.

161

Jactel, H., Petit, J., Desprez-Loustau, M.-L., Delzon, S., Piou, D., Battisti, A., and Koricheva, J.
(2012). Drought effects on damage by forest insects and pathogens: a meta-analysis. Global
Change Biology, 18(1):267–276.
Jump, A. S., Ruiz-Benito, P., Greenwood, S., Allen, C. D., Kitzberger, T., Fensham, R., MartĂ­nezVilalta, J., and Lloret, F. (2017). Structural overshoot of tree growth with climate variability and
the global spectrum of drought-induced forest dieback. Global Change Biology, 23(9):3742–
3757.
Kelly, R., Chipman, M. L., Higuera, P. E., Stefanova, I., Brubaker, L. B., and Hu, F. S. (2013).
Recent burning of boreal forests exceeds fire regime limits of the past 10,000 years. Proceedings
of the National Academy of Sciences, 110(32):13055–13060.
Klos, R. J., Wang, G. G., Bauerle, W. L., and Rieck, J. R. (2009). Drought impact on forest
growth and mortality in the southeast USA: an analysis using Forest Health and Monitoring data.
Ecological Applications, 19(3):699–708.
Kolb, T. E., Fettig, C. J., Ayres, M. P., Bentz, B. J., Hicke, J. A., Mathiasen, R., Stewart, J. E., and
Weed, A. S. (2016). Observed and anticipated impacts of drought on forest insects and diseases
in the United States. Forest Ecology and Management, 380:321–334.
Kurz, W. A., Dymond, C. C., Stinson, G., Rampley, G. J., Neilson, E. T., Carroll, A. L., Ebata, T.,
and Safranyik, L. (2008). Mountain pine beetle and forest carbon feedback to climate change.
Nature, 452(7190):987–990.
Kurz, W. A., Shaw, C. H., Boisvenue, C., Stinson, G., Metsaranta, J., Leckie, D., Dyk, A., Smyth,
C., and Neilson, E. T. (2013). Carbon in Canada’s boreal forest — A synthesis. Environmental
Reviews, 21(4):260–292.
Laurent, M., Antoine, N., and Joël, G. (2003). Effects of different thinning intensities on drought response in Norway spruce (Picea abies (L.) Karst.). Forest Ecology and Management, 183(1):47–
60.
Long, C. J., Whitlock, C., Bartlein, P. J., and Millspaugh, S. H. (1998). A 9000-year fire history
from the Oregon Coast Range, based on a high-resolution charcoal study. Canadian Journal of
Forest Research, 28(5):774–787.
Lorimer, C. G. and Frelich, L. E. (1989). A methodology for estimating canopy disturbance
frequency and intensity in dense temperate forests. Canadian Journal of Forest Research,
19(5):651–663.
Lorimer, C. G. and White, A. S. (2003). Scale and frequency of natural disturbances in the
northeastern US: implications for early successional forest habitats and regional age distributions.
Forest Ecology and Management, 185(1):41–64.
Lutz, J. A., Van Wagtendonk, J. W., and Franklin, J. F. (2010). Climatic water deficit, tree species
ranges, and climate change in Yosemite National Park. Journal of Biogeography, 37(5):936–950.

162

Lynch, J. A., Clark, J. S., and Stocks, B. J. (2004). Charcoal production, dispersal, and deposition
from the Fort Providence experimental fire: interpreting fire regimes from charcoal records in
boreal forests. Canadian Journal of Forest Research, 34(8):1642–1656.
Macalady, A. K. and Bugmann, H. (2014). Growth-mortality relationships in piĂąon pine (Pinus
edulis) during severe droughts of the past century: shifting processes in space and time. PloS
one, 9(5):e92770.
Marlon, J. R., Bartlein, P. J., Gavin, D. G., et al. (2012). Long-term perspective on wildfires in the
western USA. Proceedings of the National Academy of Sciences, 109(9):E535–E543.
Martínez-Vilalta, J., López, B. C., Loepfe, L., and Lloret, F. (2012). Stand-and tree-level determinants of the drought response of Scots pine radial growth. Oecologia, 168(3):877–888.
Mattson, W. J. and Haack, R. A. (1987). The role of drought in outbreaks of plant-eating insects.
BioScience, 37(2):110–118.
McDowell, N., Pockman, W. T., Allen, C. D., Breshears, D. D., Cobb, N., Kolb, T., Plaut, J.,
Sperry, J., West, A., Williams, D. G., and Yepez, E. A. (2008). Mechanisms of plant survival
and mortality during drought: why do some plants survive while others succumb to drought?
New Phytologist, 178(4):719–739.
McDowell, N. G., Adams, H. D., Bailey, J. D., Hess, M., and Kolb, T. E. (2006). Homeostatic
maintenance of ponderosa pine gas exchange in response to stand density changes. Ecological
Applications, 16(3):1164–1182.
McDowell, N. G., Beerling, D. J., Breshears, D. D., Fisher, R. A., Raffa, K. F., and Stitt, M. (2011).
The interdependence of mechanisms underlying climate-driven vegetation mortality. Trends in
Ecology and Evolution, 26(10):523–532.
Michaelian, M., Hogg, E. H., Hall, R. J., and Arsenault, E. (2011). Massive mortality of aspen
following severe drought along the southern edge of the Canadian boreal forest. Global Change
Biology, 17(6):2084–2094.
Millar, C. I., Stephenson, N. L., and Stephens, S. L. (2007). Climate change and forests of the
future: managing in the face of uncertainty. Ecological Applications, 17(8):2145–2151.
Ministère des Forêts, de la Faune et des Parcs (2016). Aires infestÊes par la tordeuse des bourgeons
de l’épinette au Québec en 2016 - Version 1.0. Québec, gouvernement du Québec, Direction de
la protection des forĂŞts, 16 p.
Nealis, V. G. and RÊgnière, J. (2004). Insect host relationships influencing disturbance by the spruce
budworm in a boreal mixedwood forest. Canadian Journal of Forest Research, 34(9):1870–1882.
Ogle, K. and Barber, J. J. (2008). Bayesian data-model integration in plant physiological and
ecosystem ecology. In LĂźttge, U., Beyschlag, W., and Murata, J., editors, Progress in Botany 69,
pages 281–311. Springer Berlin Heidelberg, Berlin, Heidelberg.

163

Ogle, K., Barber, J. J., Barron-Gafford, G. A., Bentley, L. P., Young, J. M., Huxman, T. E.,
Loik, M. E., and Tissue, D. T. (2015). Quantifying ecological memory in plant and ecosystem
processes. Ecology Letters, 18(3):221–235.
Oliver, C. D. and Larson, B. C. (1996). Forest Stand Dynamics, Update Edition. John Wiley &
Sons, New York, NY.
Ouimet, R., Duchesne, L., Houle, D., and Arp, P. A. (2001). Critical loads and exceedances of acid
deposition and associated forest growth in the northern hardwood and boreal coniferous forests
in Québec, Canada. Water, Air and Soil Pollution: Focus, 1(1):119–134.
Pacala, S. W., Canham, C. D., Saponara, J., Silander, J. A., Kobe, R. K., and Ribbens, E. (1996). Forest models defined by field measurements: Estimation, error analysis and dynamics. Ecological
Monographs, 66(1):1–43.
Pan, Y., Birdsey, R. A., Fang, J., Houghton, R., Kauppi, P. E., Kurz, W. A., Phillips, O. L.,
Shvidenko, A., Lewis, S. L., Canadell, J. G., et al. (2011). A large and persistent carbon sink in
the world’s forests. Science, 333(6045):988–993.
Park, T. and Casella, G. (2008). The Bayesian Lasso. Journal of the American Statistical Association, 103(482):681–686.
Parry, D., Spence, J. R., and Volney, W. J. A. (1998). Budbreak phenology and natural enemies mediate survival of first-instar forest tent caterpillar (Lepidoptera: Lasiocampidae). Environmental
Entomology, 27(6):1368–1374.
Peng, C., Ma, Z., Lei, X., Zhu, Q., Chen, H., Wang, W., Liu, S., Li, W., Fang, X., and Zhou, X.
(2011). A drought-induced pervasive increase in tree mortality across Canada’s boreal forests.
Nature Climate Change, 1(9):467–471.
Peters, M. E. and Higuera, P. E. (2007). Quantifying the source area of macroscopic charcoal with
a particle dispersal model. Quaternary Research, 67(2):304–310.
Petris, G., Petrone, S., and Campagnoli, P. (2009). Dynamic Linear Models with R. Use R! Springer
New York.
Pothier, D., Elie, J.-G., Auger, I., Mailly, D., and Gaudreault, M. (2012). Spruce budworm-caused
mortality to balsam fir and black spruce in pure and mixed conifer stands. Forest Science,
58(1):24–33.
Pothier, D., Mailly, D., and Tremblay, S. (2005). Predicting balsam fir growth reduction caused
by spruce budworm using large-scale historical records of defoliation. Annals of Forest Science,
62(3):261–267.
Power, M. J., Marlon, J., Ortiz, N., et al. (2008). Changes in fire regimes since the Last Glacial
Maximum: an assessment based on a global synthesis and analysis of charcoal data. Climate
Dynamics, 30(7):887–907.

164

Powers, M. D., Palik, B. J., Bradford, J. B., Fraver, S., and Webster, C. R. (2010). Thinning
method and intensity influence long-term mortality trends in a red pine forest. Forest Ecology
and Management, 260(7):1138–1148.
Price, D. T., Alfaro, R., Brown, K., Flannigan, M., Fleming, R., Hogg, E., Girardin, M., Lakusta,
T., Johnston, M., McKenney, D., et al. (2013). Anticipating the consequences of climate change
for Canada’s boreal forest ecosystems. Environmental Reviews, 21(4):322–365.
PRISM Climate Group (2013). Oregon State University, http://prism.oregonstate.edu, created 15
Aug 2013.
Puettmann, K. J. (2011). Silvicultural challenges and options in the context of global change: “Simple” fixes and opportunities for new management approaches. Journal of Forestry, 109(6):321–
331.
Pureswaran, D. S., De GrandprÊ, L., ParÊ, D., Taylor, A., Barrette, M., Morin, H., RÊgnière, J.,
and Kneeshaw, D. D. (2015). Climate-induced changes in host tree–insect phenology may drive
ecological state-shift in boreal forests. Ecology, 96(6):1480–1491.
Purves, D. and Pacala, S. (2008). Predictive models of forest dynamics. Science, 320(5882):1452–
1453.
R Core Team (2016). R: A Language and Environment for Statistical Computing. R Foundation
for Statistical Computing, Vienna, Austria.
Régnière, J. (1996). Generalized approach to landscape-wide seasonal forecasting with temperaturedriven simulation models. Environmental Entomology, 25(5):869–881.
RÊgnière, J., St-Amant, R., and Duval, P. (2012). Predicting insect distributions under climate
change from physiological responses: spruce budworm as an example. Biological Invasions,
14(8):1571–1586.
Reinikainen, M., D’Amato, A. W., and Fraver, S. (2012). Repeated insect outbreaks promote multicohort aspen mixedwood forests in northern Minnesota, USA. Forest Ecology and Management,
266:148–159.
Robert, C. and Casella, G. (2004). Monte Carlo Statistical Methods. Springer Texts in Statistics.
Springer-Verlag New York.
Roberts, G. O. and Rosenthal, J. S. (2009). Examples of adaptive MCMC. Journal of Computational
and Graphical Statistics, 18(2):349–367.
Roland, J. (1993). Large-scale forest fragmentation increases the duration of tent caterpillar
outbreak. Oecologia, 93(1):25–30.
Rollinson, C. R., Kaye, M. W., and Canham, C. D. (2016). Interspecific variation in growth
responses to climate and competition of five eastern tree species. Ecology, 97(4):1003–1011.
Ross, S. M. (2010). Introduction to Probability Models. Academic Press, Burlington, MA, 10
edition.
165

Safranyik, L. and Carrol, A. L. (2006). The biology and epidemiology of the mountain pine beetle
in lodgepole pine forests. In Safranyik, L. and Wilson, W. R., editors, The mountain pine beetle:
a synthesis of biology, management, and impacts on lodgepole pine, chapter 1, pages 3–66.
Natural Resources Canada, Canadian Forest Service, Pacific Forestry Centre, Victoria, British
Columbia.
Saxton, K. E. and Rawls, W. J. (2006). Soil water characteristic estimates by texture and organic
matter for hydrologic solutions. Soil Science Society of America Journal, 70(5):1569–1578.
Schofield, M. R., Barker, R. J., Gelman, A., Cook, E. R., and Briffa, K. R. (2016). A model-based
approach to climate reconstruction using tree-ring data. Journal of the American Statistical
Association, 111(513):93–106.
Simard, M., Romme, W. H., Griffin, J. M., and Turner, M. G. (2011). Do mountain pine beetle
outbreaks change the probability of active crown fire in lodgepole pine forests? Ecological
Monographs, 81(1):3–24.
Smith, D. J. and Woods, M. E. (1997). Red pine and white pine density management diagrams
for Ontario. Ontario Ministry of Natural Resources, Southcentral Sciences Section, Technical
Report No. 48.
Sohn, J. A., Hartig, F., Kohler, M., Huss, J., and Bauhus, J. (2016). Heavy and frequent thinning
promotes drought adaptation in Pinus sylvestris forests. Ecological Applications, 26(7):2190–
2205.
Soil Landscapes of Canada Working Group (2010). Soil Landscapes of Canada version 3.2.
Agriculture and Agri-Food Canada.
Stephens, S. L., Agee, J. K., FulĂŠ, P. Z., North, M. P., Romme, W. H., Swetnam, T. W., and Turner,
M. G. (2013). Managing forests and fire in changing climates. Science, 342(6154):41–42.
Stephenson, N. L. (1998). Actual evapotranspiration and deficit: biologically meaningful correlates
of vegetation distribution across spatial scales. Journal of Biogeography, 25:855–870.
Swetnam, T. W., Allen, C. D., and Betancourt, J. L. (1999). Applied historical ecology: Using the
past to manage for the future. Ecological Applications, 9(4):1189–1206.
Turner, M. G., Smithwick, E. A. H., Metzger, K. L., Tinker, D. B., and Romme, W. H. (2007). Inorganic nitrogen availability after severe stand-replacing fire in the Greater Yellowstone ecosystem.
Proceedings of the National Academy of Sciences, 104(12):4782–4789.
United States Department of Agriculture, Forest Service (2007). Forest inventory and analysis
strategic plan: A history of success, a dynamic future, fs-865.
Van Deusen, P. C. (1987). Testing for stand dynamics effects on red spruce growth trends. Canadian
Journal of Forest Research, 17(12):1487–1495.
Van Deusen, P. C. (1989). A model-based approach to tree ring analysis. Biometrics, 45(3):763–779.

166

van Mantgem, P. J. and Stephenson, N. L. (2007). Apparent climatically induced increase of tree
mortality rates in a temperate forest. Ecology Letters, 10(10):909–916.
van Mantgem, P. J., Stephenson, N. L., Mutch, L. S., Johnson, V. G., Esperanza, A. M., and Parsons,
D. J. (2003). Growth rate predicts mortality of Abies concolor in both burned and unburned
stands. Canadian Journal of Forest Research, 33(6):1029–1038.
Vanderwel, M. C. and Purves, D. W. (2014). How do disturbances and environmental heterogeneity
affect the pace of forest distribution shifts under climate change? Ecography, 37(1):10–20.
Visser, H. (1986). Analysis of tree ring data using the Kalman filter technique. IAWA Bulletin n.s.,
7(4):289–297.
Visser, H. and Molenaar, J. (1988). Kalman filter analysis in dendroclimatology. Biometrics,
44:929–940.
Wang, L., Libert, G., and Manneback, P. (1992). Kalman filter algorithm based on singular value
decomposition. Proceedings of the 31st Conference on Decision and Control, pages 1224–1229.
Weiskittel, A. R., Hann, D. W., Kershaw, J. A., and Vanclay, J. K. (2011). Forest Growth and Yield
Modeling. John Wiley & Sons, Ltd., Hoboken, NJ.
West, M. and Harrison, J. (1997). Bayesian Forecasting and Dynamic Models. Springer Series in
Statistics. Springer New York.
Whitlock, C. and Millspaugh, S. H. (1996). Testing the assumptions of fire-history studies:
an examination of modern charcoal accumulation in Yellowstone National Park, USA. The
Holocene, 6(1):7–15.
Wikle, C. K. (2003). Hierarchical Bayesian models for predicting the spread of ecological processes.
Ecology, 84(6):1382–1394.
Wikle, C. K., Milliff, R. F., Nychka, D., and Berliner, L. M. (2001). Spatiotemporal hierarchical
Bayesian modeling tropical ocean surface winds. Journal of the American Statistical Association,
96(454):382–397.
Williams, A. P., Allen, C. D., Millar, C. I., Swetnam, T. W., Michaelsen, J., Still, C. J., and Leavitt,
S. W. (2010). Forest responses to increasing aridity and warmth in the southwestern United
States. Proceedings of the National Academy of Sciences, 107(50):21289–21294.
Wood, S. (2006). Generalized Additive Models: An Introduction with R. Chapman & Hall/CRC
Texts in Statistical Science. Taylor & Francis.
Wood, S. N. and Augustin, N. H. (2002). GAMs with integrated model selection using penalized regression splines and applications to environmental modelling. Ecological Modelling,
157(2):157–177.
Worrall, J. J., Marchetti, S. B., Egeland, L., Mask, R. A., Eager, T., and Howell, B. (2010). Effects
and etiology of sudden aspen decline in southwestern Colorado, USA. Forest Ecology and
Management, 260(5):638–648.
167

Wyckoff, P. H. and Clark, J. S. (2000). Predicting tree mortality from diameter growth: a comparison of maximum likelihood and Bayesian approaches. Canadian Journal of Forest Research,
30(1):156–167.
Wyckoff, P. H. and Clark, J. S. (2002). The relationship between growth and mortality for seven
co-occurring tree species in the southern Appalachian Mountains. Journal of Ecology, 90(4):604–
615.
Yamaguchi, D. K. (1991). A simple method for cross-dating increment cores from living trees.
Canadian Journal of Forest Research, 21(3):414–416.
Zhang, J., Huang, S., and He, F. (2015). Half-century evidence from western Canada shows forest
dynamics are primarily driven by competition followed by climate. Proceedings of the National
Academy of Sciences, 112(13):4009–4014.

168