UNLOCKING THE FOREST INVENTORY AND ANALYSIS DATABASE:
APPLICATIONS TO NATION-WIDE FOREST HEALTH MONITORING

By

Hunter Stanke

A THESIS

Submitted to

Michigan State University

in partial fulfillment of the requirements

for the degree of

Forestry – Master of Science

2020

ABSTRACT

UNLOCKING THE FOREST INVENTORY AND ANALYSIS DATABASE:
APPLICATIONS TO NATION-WIDE FOREST HEALTH MONITORING

By

Hunter Stanke

Forest Inventory and Analysis (FIA) is a US Department of Agriculture Forest Service pro-

gram that aims to monitor changes in forests across the US. FIA hosts one of the largest

ecological datasets in the world, though its complexity limits access for many potential users.

rFIA is an R package designed to simplify the estimation of forest attributes using data

collected by the FIA Program. Specifically, rFIA improves access to the spatio-temporal es-

timation capacity of the FIA Database via space-time indexed summaries of forest variables

within user-defined population boundaries. The package implements multiple design-based

estimators, and has been validated against official estimates and sampling errors produced by

the FIA Program. The package has been made open-source is freely available for download

from the Comprehensive R Archive Network.

In recent decades, forests of the western US have experienced unprecedented change in cli-

mate and forest disturbance regimes, and widespread shifts in forest composition, structure,

and function are expected in response. However, large-scale, comprehensive assessments of

tree population performance have yet to be conducted in the region. We develop an index of

forest population performance based on repeated censuses of field plots, and apply this index

to assess the status of the most abundant tree species in the western US. Our study provides

empirical evidence to suggest tree species in the western US are exhibiting strong divergence

in population performance, with over half (70%) of species experiencing range-wide popu-

lation decline. We found spatial variation in population performance across the ranges of

all species, indicating range shifts are already underway. Our results further indicate that

species decline can seldom be attributed to a single forest disturbance agent, highlighting

the importance of considering multiple risks factors in broad-scale forest management.

Copyright by
HUNTER STANKE
2020

ACKNOWLEDGEMENTS

I would like to thank my adviser Dr. Andrew Finley for his academic guidance over the past

year. I’d also like to acknowledge my committee members Dr. Aaron Weed and Dr. David

Macfarlane for their invaluable comments and criticisms of my thesis work.

iv

TABLE OF CONTENTS

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

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

1.2.1
1.2.2
1.2.3

CHAPTER 1 RFIA: AN R PACKAGE FOR ESTIMATION OF FOREST AT-
TRIBUTES WITH THE US FOREST INVENTORY AND ANAL-
YSIS DATABASE . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.1
Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Software design . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Sampling and estimation procedures
. . . . . . . . . . . . . . . . . .
Software testing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.3 rFIA package features . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.3.1 Loading FIA data . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.3.2
. . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . .
1.3.3 Estimating forest attributes
1.3.3.1 Grouped Estimates . . . . . . . . . . . . . . . . . . . . . . .
1.3.3.2 Unique populations of interest . . . . . . . . . . . . . . . . .
1.3.3.3 Grouping by user-defined areal units . . . . . . . . . . . . .
1.3.3.4 Alternative estimators . . . . . . . . . . . . . . . . . . . . .
1.3.3.5 Parallelization . . . . . . . . . . . . . . . . . . . . . . . . .
1.4 Case study - Michigan ash decline . . . . . . . . . . . . . . . . . . . . . . . .
1.5 Extensions/limitations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.7 Acknowledgements
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Spatial and temporal subsets

CHAPTER 2 OVER HALF OF TOP TREE SPECIES IN DECLINE IN THE

WESTERN UNITED STATES . . . . . . . . . . . . . . . . . . . . .
2.1
Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2 Materials and Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2.1 Field observations . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2.2 Climate data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2.3 Forest stability index . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2.4
. . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Statistical analysis

2.4.1 Range-wide population performance and regional shifts in forest

composition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.4.2 Forest disturbances as drivers of tree population performance . . . . .
2.4.3 Early indications of species range shifts . . . . . . . . . . . . . . . . .
2.4.4
Intra-specific divergence in population performance among size classes

v

1
1
2
2
3
7
8
8
10
10
11
13
14
14
15
16
19
20
20

21
21
24
24
25
26
29
31
39

39
41
42
43

2.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

44

APPENDIX . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

45

BIBLIOGRAPHY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

55

vi

LIST OF TABLES

Table 1.1: Comparison of estimates produced by rFIA and EVALIDator for select
forest attributes. Estimates from both tools were produced using the
“temporally-indifferent” estimator. . . . . . . . . . . . . . . . . . . . . . .

Table 1.2: List of core functions available within rFIA. . . . . . . . . . . . . . . . . .

7

8

Table 1.3: List of core arguments available in rFIA estimator functions.

. . . . . . .

11

Table 2.1: Scientific and common names of ten most abundant species across our
study area, listed in order of increasing abundance (i.e., bottom being
most abundant; abundance defined by estimated total number of trees).
Sample size (number of re-measured FIA plots where the species occured)
provided for small diameter (2.54-12.6 cm DBH), large diameter (DBH
≥ 12.7 cm DBH), and combined populations (DBH ≥ 2.54 cm). As large
diameter and small diameter populations of a species may exist on the
same plot, the sample size of the combined population is not equal to
the sum of large diameter and small diameter sample sizes.

. . . . . . . .

Table 2.2:

Intra-specific differences in range-wide mean Forest Stability Index (FSI)
between tree size classes. Significance (i.e., p-values) derived from Kruskall-
Wallace test using plot-level FSI values.
. . . . . . . . . . . . . . . . . . .

32

36

vii

LIST OF FIGURES

Figure 1.1: Software availability . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 1.2: FIA ground plot design (left) and an example map of forested condition

classes on a subplot (right).

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

Figure 1.3: Schematic diagram of the features of rFIA: the blue boxes represent
decision points for the user; green boxes represent data or results that
are produced and/or processed by rFIA functions, which are represented
in the orange boxes.

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

3

4

9

Figure 1.4: Example code for producing basic population and plot-level estimates

using rFIA.

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

12

Figure 1.5: Example code to produced population estimates of forest attributes

grouped by discrete categories using rFIA.

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

13

Figure 1.6: Example code to produce population estimates of forest attributes within

a user-defined population of interest using rFIA. . . . . . . . . . . . . . .

14

Figure 1.7: Example code to produce population estimates of forest attributes for
multiple values of λ using the exponential moving average estimator in
rFIA. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 1.8: Changes in live ash TPA (left) and annual mortality rates (right), with
associated sampling errors (bottom; 67% confidence), across counties in
the Lower Peninsula of Michigan following establishment of emerald ash
borer. Missing values are shaded in gray. All plots were produced using
the plotFIA function.

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

Figure 1.9: Changes in ash sawlog volume (top; board feet (BF) per acre) and
associated sampling errors (bottom; 67% confidence) for each design-
based estimator, across the Lower Peninsula of Michigan from 2000 to
2018. All plots were produced using the plotFIA function.

. . . . . . .

15

17

18

viii

Figure 2.1: Example relationship of change in TPA (S, x-axis) and change in BAA
(B, y-axis) for subalpine fir (Abies lasiocarpa) in Colorado (all plots
remeasured from 2015-2018). Small diameter circles represent plot level-
indices of S and B, and the single large diamond represents the region-
wide mean for each dimension. Small squares distributed on the 1:1 line
(gray-dashed line) represent plot-level FSI values corresponding to six
example plots in the upper-right of the graph (connected by by black
segments). All points are colored by their respective FSI value. Points
in the bottom-left of the graph are characterized by losses in both TPA
and BAA, indicating the species is likely declining within these stands.
Conversely, points in the upper-right of the graph are characterized by
increases in TPA and BAA, indicating species is likely expanding in the
stand (i.e., increasing in density).

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

Figure 2.2: Spatial variation in forest disturbance severity and long-term climate
patterns across the study region. All variables standardized to a com-
mon scale for visual interpretation of spatial variation.

. . . . . . . . . .

Figure 2.3: Region-wide prevalence (% of total tree population across the study
region) of the 10 most abundant tree species in the western US in 2018.
Estimates provided for combined (left), large diameter (center), and
small diameter populations (right). Species listed in order of increasing
abundance with respect to the combined population (i.e., bottom being
most abundant).

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

Figure 2.4: Range-wide mean Forest Stability Index (FSI) of the 10 most abundant
tree species in the western US. Estimates provided for combined (left),
large diameter (center), and small diameter populations (right). Popu-
lation decline (red) occurs when the FSI is negative and the associated
confidence interval does not include zero. Conversely, population expan-
sion (blue) occurs when the FSI is positive and the associated confidence
interval does not include zero.

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

27

31

33

34

ix

Figure 2.5: (Top) Spatial variation in population performance of species across their
ranges. Forest Stability Index (FSI) summarized within ecoregion sub-
sections using the simple moving average estimator. Population expan-
sion occurs if the mean FSI significantly exceeds zero within an area (net
increase in abundance), and population decline occurs if the mean FSI
is significantly less than zero (net decrease in abundance). Populations
are determined to be stable within an area if the confidence interval of
estimated mean FSI includes zero. (Bottom) Spatial variation in intra-
specific differences in population performance between tree size classes.
Areas where small diameter populations exhibit significantly higher FSI
than large diameter populations shaded in purple. Conversely, areas
where small diameter populations exhibit significantly lower FSI than
large diameter populations shaded in yellow. Otherwise populations are
determined to be changing at equal rates (green).

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

Figure 2.6: Standardized coefficients of the linear mixed-effect model predicting
plot-level FSI as a function of forest disturbance severity and long-term
climate patterns. All variables were scaled at the species level, hence
the relative importance of each predictor in explaining intra-specific
variation in population performance may be directly inferred from the
magnitude of each coefficient. That is, coefficient values can be directly
compared across tree size classes (i.e., vertically) but not across species
(i.e., horizontally). Blue shaded backgrounds indicate the population
is expanding across the study region (e.g., large diameter grand fir;
top left). Red shaded backgrounds indicate the population is declin-
ing across the study region (e.g., large diameter Engelmann spruce; top
right). Error bars represent the 95% confidence intervals in the estimate
of each coefficient.

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

Figure A.1: Example code to download and process data used in software validation.

35

37

46

Figure A.2: Example code to estimate live tree abundance . . . . . . . . . . . . . . .

47

Figure A.3: Example code to estimate live tree basal area . . . . . . . . . . . . . . .

47

Figure A.4: Example code to estimate live tree merchantable volume. . . . . . . . . .

48

Figure A.5: Example code to estimate live tree saw log volume.

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

49

Figure A.6: Example code to estimate live tree aboveground biomass . . . . . . . . .

49

Figure A.7: Example code to estimate live tree carbon . . . . . . . . . . . . . . . . .

50

Figure A.8: Example code to estimate live tree biomass growth.

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

51

x

Figure A.9: Example code to estimate tree mortality. . . . . . . . . . . . . . . . . . .

51

Figure A.10: Example code to estimate tree removals.

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

52

Figure A.11: Example code to estimate coarse woody debris volume.

. . . . . . . . . .

53

Figure A.12: Example code to estimate coarse woody debris biomass.

. . . . . . . . .

53

Figure A.13: Example code to estimate coarse woody debris carbon.

. . . . . . . . . .

54

Figure A.14: Example code to estimate total forested area.

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

54

xi

CHAPTER 1

RFIA: AN R PACKAGE FOR ESTIMATION OF FOREST ATTRIBUTES

WITH THE US FOREST INVENTORY AND ANALYSIS DATABASE

1.1 Introduction

Forest Inventory and Analysis (FIA) is a United States Department of Agriculture For-

est Service program with the goal of monitoring and projecting changes in forests across

the United States (US) [USDA Forest Service, 2019a]. The program collects, publishes, and

analyzes data describing the extent, condition, volume, growth, and use of trees from all

land ownerships in the nation, with some records dating back to the early 1930s [Tinkham

et al., 2018, Smith, 2002]. In 1999, the FIA program established a systematic grid of per-

manent ground plots that are evenly divided into panels measured in a continuous cycle,

allowing spatially unbiased estimates of forest attributes to be computed on an annual basis

[Gillespie, 1999, Smith, 2002]. The flexibility inherent to the FIA inventory system creates

unprecedented opportunities to assess forest change across space, through time, and within

unique populations of interest [Tinkham et al., 2018, Gray et al., 2012].

The research significance of the FIA program has increased in recent decades, with pri-

mary applications in carbon cycling, forest growth, forest health, and remote sensing [Tin-

kham et al., 2018]. The FIA program hosts one of the largest ecological datasets in the world

by spatial and temporal extent [Tinkham et al., 2018], including records from over 5.8 million

trees measured at least twice and encompassing a range of ecological diversity unmatched

by any other large-scale national forest inventory system [Tomppo et al., 2010]. The breadth

of attributes currently measured by the FIA program, ranging from forest composition and

structure to soil chemistry and invasive species abundance, makes it a vast and powerful

resource for monitoring the status and trends in forest attributes from the scale of individual

trees to subcontinents [Tinkham et al., 2018].

1

The primary limitation for individuals using FIA data in their own analyses is the com-

plex sample design, database structure, and Structured Query Language used by the FIA

program [Tinkham et al., 2018, Kromroy et al., 2008]. FIA data are publicly available in

several formats (e.g., FIA DataMart) and estimation is facilitated through online tools (e.g.,

EVALIDator) [USDA Forest Service, 2019b, Pugh et al., 2018] and a non-public R package,

FIESTA [Frescino et al., 2015]. That said, FIA data are difficult for many non-FIA users to

interpret and understand [Tinkham et al., 2018] and there is potential to increase the use of

these incredibly rich data among industry professionals, academic scientists, and the general

public. To promote use of this valuable public resource and extend the reach of the FIA

program and the publicly available data, there is a need for a flexible, user-friendly tool that

simplifies the process of working with FIA data for experts and novices alike. To this end,

we developed rFIA, an add-on package for R [R Core Team, 2018b].

rFIA implements FIA’s design-based estimation procedures [Bechtold et al., 2005] for

over 50 forest attributes using a simple, yet powerful design. rFIA greatly improves access

to the spatio-temporal estimation capacity inherent to the FIA program by allowing space-

time indexed summaries of forest attributes to be produced within user-defined population

boundaries (e.g., geographic, temporal, biophysical). The package enhances the value of

FIA for temporal change detection and forest health monitoring by implementing five design-

based estimators that offer flexibility in a balance between precision and temporal smoothing.

Our intention in developing rFIA is to provide a versatile, user-friendly software that allows

all R users to unlock the value of the FIA program.

1.2 Methods

1.2.1 Software design

We designed rFIA to be intuitive to use and support common data representations by di-

rectly integrating other popular R packages into our development. We achieve efficient joins,

queries, and data summaries with the dplyr package [Wickham et al., 2015]. Specifically,

2

Figure 1.1: Software availability

Name of software rFIA
Type of software Add-on package for R
First available 2019
Program languages R
License GPL 3
Code Repository https://cran.r-project.org/web/packages/rFIA/index.html
Installation in R install.packages("rFIA")
Developers Hunter Stanke, Andrew O. Finley
Contact Address Department of Forestry,

Michigan State University, East Lansing, MI, USA

we leverage dplyr for joining and filtering FIA tables, facilitating hierarchical grouping of

summary attributes, and implementing non-standard evaluation in rFIA core functions. We

achieve efficient space-time query and summary within user-defined population boundaries

(i.e., spatial polygons) with the sf package [Pebesma, 2018]. Parallel processing is im-

plemented with the parallel package [R Core Team, 2018a]. Parallel implementation is

achieved using a snow type cluster [Tierney et al., 2018] on any Windows OS, and with

multi-core forking [R Core Team, 2018a] on any Unix or Mac OS.

1.2.2 Sampling and estimation procedures

The design-based estimation procedures used by the FIA program and implemented by rFIA

have been widely described in the literature [Bechtold et al., 2005, McRoberts et al., 2005a,

Hoffman et al., 2014], and hence will only be briefly described here. All estimators for pop-

ulation totals, ratios, and associated variances are derived from the Horvitz-Thompson esti-

mator and hence incorporate design information via inverse-probability weighting [Horvitz

and Thompson, 1952]. Bechtold et al. [2005] describes in detail the theoretical basis for

estimators used by the FIA program and implemented by rFIA.

The FIA program conducts forest inventories using a multi-phased sampling procedure

designed to reduce variance through stratification [Bechtold et al., 2005, McRoberts et al.,

2005a]. In the pre-field phase, remotely sensed imagery are used to stratify land area by

3

Figure 1.2: FIA ground plot design (left) and an example map of forested condition classes
on a subplot (right).

determining dominant land use in each pixel within the population of interest. In the core

phase, permanent ground plots are systematically distributed across the US at a rate of

approximately 1 plot per 6000 acres using a hexagonal sampling frame [Bechtold et al.,

2005, McRoberts et al., 2005a]. Each plot is assigned to a single stratum based on the the

pre-field stratification of the plot center. If any portion of a plot is determined to contain a

forest land use from pre-field stratification, a field crew will visit the site and measure core

FIA variables [Bechtold et al., 2005]. In the intensive phase, additional forest and ecosystem

health variables are measured on 5-15% of established core plots (approximately 1 plot per

96,000 acres).

FIA permanent ground plots consist of clusters of four subplots (Figure 1.2), where

tree attributes are measured for all stems 5.0 inches diameter at breast height (DBH) and

larger. Within each subplot, a microplot is established where tree attributes are measured

for saplings (1.0-4.9 inches DBH). Each subplot is surrounded by a macroplot, on which rare

events such as large trees are optionally measured. In addition to tree attributes, data are

collected that describe the area where trees are located on each subplot and macroplot. If

area attributes (e.g., ownership group, forest type, stocking) vary substantially across a plot,

4

the land area within the plot is divided into distinct domains referred to as condition classes

so that tree data can be properly associated with area classifications (Figure 1.2).

Plot-level estimates of forest attributes (e.g., tree biomass, forested area) are obtained

by summing observations or estimates of the attribute, which are within the population of

interest (e.g., live trees, private land), across the plot. Observations are multiplied by 1 if the

attribute is within the population of interest, and 0 otherwise. Thus if all observations on a

plot fall outside the population of interest, observations will sum to 0 for the plot. For ratio

estimates (e.g., tree biomass per unit area), it is possible to specify different populations of

interest for the numerator and denominator attributes.

The FIA program uses an annual panel system to estimate both current inventory and

change. Individual panels are represented by a subset of ground plots that are measured

in the same year and represent complete spatial coverage across the population of interest

[McRoberts and Miles, 2016]. In the eastern US, inventory cycles consist of 5 or 7 annual

panels, with 20 or 15 percent of ground plots measured in each year, respectively. In the

western US, inventory cycles include 10 annual panels and thus 10 percent of ground plots

are measured each year [Bechtold et al., 2005]. Variance reduction can often be achieved by

combining current panels with data from previous panels, although FIA does not prescribe

a core procedure for panel combination because a single approach is unlikely to be suitable

for all estimation objectives or across a wide variety of spatial, temporal, and population

conditions [Bechtold et al., 2005].

In rFIA, users may choose from one of four unique estimators that combine data from

multiple panels or choose to return estimates from individual annual panels (ANNUAL),

thereby forgoing panel combination entirely. The “temporally-indifferent” (TI) estimator is

used by default, essentially pooling all panels within an inventory cycle into a large periodic

inventory. The TI estimator is commonly used by other FIA estimation tools, like EVAL-

IDator [USDA Forest Service, 2019b, Pugh et al., 2018]. Alternatively, individual panels

may be weighted by employing a moving-average estimator, including the simple moving

5

average (SMA), linear moving average (LMA), or exponential moving average (EMA). The

SMA applies equal weight to all annual panels (Eq. 1.1), while the LMA and EMA apply

weights that decay linearly or exponentially as a function of time since measurement, respec-

tively (Eq. 1.2-1.3). In each case of the moving average, panel weights sum to one across an

inventory cycle. Panel weights are computed as follows:

1
N

wp,LM A =

wp,SM A =

p(cid:80)N
i=1 pi
(cid:80)N
λ(1 − λ)1−p
i=1 λ(1 − λ)1−pi

wp,EM A =

(1.1)

(1.2)

(1.3)

where wp is the constant positive weight for annual panel, p, described as the sequential

index of the panel within an inventory cycle (i.e., p = 1, 2, . . . , N ), N is the total number

of annual panels in the inventory cycle, and λ is a decay parameter, ranging between zero

and one, that controls the behavior of the exponential weighting function (Eq. 3). As λ

approaches 1, panel weights become nearly evenly distributed across the inventory cycle,

and estimates produced by the EMA will approach those of the SMA. As λ approaches 0,

weights are skewed towards the most recent panel and estimates of the EMA will approach

those of individual annual panels. The inclusion of λ makes the EMA the most versatile

estimator offered in rFIA.

In general, sample variance is minimized when panel weights are evenly distributed across

the inventory cycle (e.g., SMA), with the trade-off of introducing temporal lag-bias by weight-

ing recent measurements the same as older measurements. Hence, equal weighting schemes

could be undesirable in settings where variable values change over time. Alternatively, ap-

plying lower weights to less recent panels reduces the effective sample size of the inventory,

effectively increasing variance but reducing temporal smoothing by favoring more recent

measurements. We advise users of rFIA to be aware of this trade-off between precision and

temporal smoothing when considering various estimators offered in the package.

6

Table 1.1: Comparison of estimates produced by rFIA and EVALIDator for select forest
attributes. Estimates from both tools were produced using the “temporally-indifferent”
estimator.

Forest Attribute

rFIA

Live tree abundance (trees/acre)
Live tree basal area (f t2/acre)
Live tree merchantable volume (f t3/acre)
Live tree sawlog volume (f t3/acre)
Live tree aboveground biomass (tons/acre)
Live tree aboveground carbon (tons/acre)
Annual net biomass growth (tons/acre/year)*
Annual mortality (trees/acre/year)*
Annual removals (trees/acre/year)*
Coarse woody material volume (f t3/acre)
Coarse woody material biomass (tons/acre)
Coarse woody material carbon (tons/acre)
Total forest area (acres x 10−3)

Estimate
432.63
121.19
2625.99
1648.77
75.99
37.99
1.06
1.47
0.36
299.87
3.08
1.52
1789.61

EVALIDator
SE
4.46
2.13
2.65
3.56
2.38
2.38
6.39
6.93
31.09
23.92
25.25
25.14
2.29

Estimate
432.63
121.19
2625.99
1648.77
75.99
37.99
1.06
1.47
0.36
299.87
3.08
1.52
1789.61

SE
4.46
2.13
2.65
3.56
2.38
2.38
6.39
6.93
31.09
23.92
25.25
25.14
2.29

1.2.3 Software testing

We have conducted extensive validation for all estimated attributes for small and large

areas against EVALIDator [USDA Forest Service, 2019b]. Here we present an abbreviated

version of a validation for the state of Connecticut in the year 2018 using the “temporally-

indifferent” estimator (Table 1.1). Current status estimates were produced for live trees
with DBH ≥ 1.0 inch and annual change estimates were produced for all stems with DBH
≥ 5.0 inches. Down woody material estimates were produced for the 1000 hour fuel class

(coarse woody debris). Total forest area estimates were produced using plots containing a

forested condition. All estimates and associated sampling errors produced by rFIA match

those produced by EVALIDator. A detailed description of the validation described in Table

1.1, including code used to produce estimates using rFIA, can be found in Appendix A.

7

Table 1.2: List of core functions available within rFIA.

Description
Estimate land area in various classes
Estimate volume, biomass, and carbon stocks of standing trees
Spatial and temporal queries for FIA data
Estimate species diversity
Estimate volume, biomass, and carbon stocks of down woody material
Lookup Evaluation IDs (EVALIDs) by year and evaluation types
Download FIA data, load into R, and optionally save to disk
Estimate recruitment, mortality, and harvest rates
Estimate areal coverage of invasive species
Produce static and animated plots of FIA summaries
Load FIA database into R environment from disk
Estimate abundance of seedlings

Function
area
biomass
clipFIA
diversity
dwm
findEVALID
getFIA
growMort
invasive
plotFIA
readFIA
seedling
standStruct Estimate forest structural stage distributions
tpa
vitalRates
writeFIA

Estimate abundance of standing trees
Estimate live tree growth rates
Write in-memory FIA database to disk

1.3 rFIA package features

rFIA is capable of estimating more forest attributes from FIA data than any other publicly

available tool and offers unmatched flexibility in defining unique populations of interest

and producing space-time indexed summaries of forest attributes. Users can install the

released version of rFIA from CRAN (v0.1.0, 28 October 2019), or alternatively install the

development version from Github (https://github.com/hunter-stanke/rFIA) using devtools

[Wickham et al., 2019]. Table 1.2 depicts a list of core functions available in rFIA. A

schematic diagram for using rFIA to produce population estimates of forest attributes can

be found in Figure 1.3.

1.3.1 Loading FIA data

Users can automatically download state subsets of the FIA database [USDA Forest Service,

2019b], load these data into R, and optionally save to disk using getFIA. Simply specify

the state abbreviation code for the state(s) of interest as a character vector in the states

8

Figure 1.3: Schematic diagram of the features of rFIA: the blue boxes represent decision
points for the user; green boxes represent data or results that are produced and/or
processed by rFIA functions, which are represented in the orange boxes.

argument of getFIA (e.g., states = "CT" for Connecticut). Alternatively, getFIA can be

used to download, load, and optionally save select table(s) contained in the FIA database

by specifying the names of the tables of interest in the tables argument (e.g., tables =

c("TREE", "PLOT") for the TREE and PLOT tables). All FIA data are downloaded from the

FIA DataMart [USDA Forest Service, 2019b] and saved as comma-delimited text files in a

local directory.

If FIA data are already saved on disk as comma-delimited text files, readFIA can be used

to load these data into R. readFIA requires that FIA data files maintain original FIA naming

conventions (as downloaded from the FIA DataMart) for referential integrity. All FIA data

files should be stored in a single directory with no sub-directories for proper loading with

readFIA. We recommend using getFIA to download and save new FIA data and readFIA to

reload these data in future R sessions. All data loaded with getFIA or readFIA are stored

9

as a modified list object called a FIA.Database. All common list operations are valid for the
FIA.Database object class, and individual FIA tables can be accessed using the $ operator.

1.3.2 Spatial and temporal subsets

Spatial and/or temporal subsets of a FIA.Database may be implemented with clipFIA

if users are interested in producing estimates for select areal regions and/or time periods

contained within the spatial-temporal extent of the FIA.Database. Such subsets are not

required to use rFIA estimator functions (listed in Table 1.3) to produce estimates of for-

est attributes. However, limiting the spatial and temporal extent of the query region will

conserve memory and decrease processing time.

Users may subset a FIA.Database to the boundaries of a spatial polygon object by

specifying the name of the spatial object in the mask argument of clipFIA. All spatial

polygon classes from the sp and sf packages are supported. To obtain the most recent

subset of a FIA.Database by reporting year, specify mostRecent = TRUE in clipFIA (TRUE

by default). If a FIA.Database contains multiple states with different reporting schedules,

setting mostRecent = TRUE will return the data necessary to produce estimates for the most

recent reporting year in each state. Alternatively, users may specify matchEval = TRUE to

obtain a subset of data associated with reporting years which are common among all states

in the FIA.Database.

1.3.3 Estimating forest attributes

rFIA includes a range of estimator functions designed to produce population and plot-level

estimates directly from FIA.Database objects. All estimator functions share a similar design,

although minor nuances exist between functions due to variation in the forms of estimates

they invoke (Table 1.3). We will demonstrate the functionality of rFIA estimator functions

using biomass, a function designed to estimate volume, biomass, and carbon of standing

trees.

10

Table 1.3: List of core arguments available in rFIA estimator functions.

bySpecies bySizeClass treeDomain areaDomain tidy method

area
biomass
diversity
dwm
growMort
invasive
seedling
standStruct
tpa
vitalRates

*

*

(default)

*

*
*

*
*

*

*
*

*
*

*

*

*
*

*
*
*
*
*
*
*
*
*
*

*
*
*
*
*
*
*
*
*
*

*

*

To produce population estimates for the region contained within the spatial extent of

a FIA.Database, specify the name of the FIA.Database as the db argument of any rFIA

estimator function (Figure 1.4). All rFIA estimator functions return population estimates

and associated sampling errors for each reporting year in the FIA.Database by default. To

return estimates at the plot-level, specify byPlot = TRUE. In this case, estimates will be

returned for each occasion that an individual plot was measured.

1.3.3.1 Grouped Estimates

Often it is useful to produce estimates grouped by discrete categories, such as species, forest

type, ownership group, or diameter class. rFIA estimator functions can produce estimates

grouped by fields contained in the PLOT, COND, and TREE tables of a FIA.Database. To

produce grouped estimates, specify the name(s) of fields representing the grouping variable(s)

as the grpBy argument of any estimator function (Figure 1.5). If more than one grouping

variable is provided to grpBy, grouping will occur hierarchically based on the order variable

names are listed. In addition to grpBy, some rFIA estimator functions include bySpecies and

bySizeClass arguments for convenience. Set either of these arguments as TRUE to produced

estimates grouped by species and/or 2.0 inch DBH classes, respectively. Alternatively, users

may specify grpBy = SPCD to group estiamtes by species and use the makeClasses function

11

Figure 1.4: Example code for producing basic population and plot-level estimates using
rFIA.

## Load state subset for Rhode Island,
## included with rFIA
data("fiaRI")
## Estimates for entire state
biomass(db = fiaRI)

# A tibble: 5 x 19

<dbl>
2396.
2438.
2491.
2500.
2491.

<dbl>
1356.
1385.
1419.
1422.
1419.

YEAR NETVOL_ACRE SAWVOL_ACRE BIO_AG_ACRE BIO_BG_ACRE BIO_ACRE
<dbl>
81.6
82.9
84.6
84.8
84.4

<int>
1 2014
2 2015
3 2016
4 2017
5 2018
# ... with 13 more variables: CARB_AG_ACRE <dbl>, CARB_BG_ACRE <dbl>,
#
#
#
#

CARB_ACRE <dbl>, NETVOL_ACRE_SE <dbl>, SAWVOL_ACRE_SE <dbl>,
BIO_AG_ACRE_SE <dbl>, BIO_BG_ACRE_SE <dbl>, BIO_ACRE_SE <dbl>,
CARB_AG_ACRE_SE <dbl>, CARB_BG_ACRE_SE <dbl>, CARB_ACRE_SE <dbl>,
nPlots_VOL <dbl>, nPlots_AREA <dbl>

<dbl>
68.0
69.1
70.6
70.8
70.4

<dbl>
13.5
13.7
14.0
14.1
14.0

## Plot-level estimates
biomass(db = fiaRI, byPlot = TRUE)

# A tibble: 352 x 11

PLT_CN YEAR NETVOL_ACRE SAWVOL_ACRE BIO_AG_ACRE BIO_BG_ACRE BIO_ACRE
<dbl>

<dbl>

<dbl>

<dbl>

<dbl>

0
0
0

84.0

103.

<dbl> <int>
2009
2009
2009
2009
2009

0
0
0

1 1.45e14
2 1.45e14
3 1.45e14
4 1.45e14
5 1.45e14
# ... with 347 more rows, and 4 more variables: CARB_AG_ACRE <dbl>,
#

CARB_BG_ACRE <dbl>, CARB_ACRE <dbl>, nStems <int>

2335.
2760.

1481.
1690.

70.4
86.3

13.6
16.7

0
0
0

0
0
0

0
0
0

to define their own diameter classes. More information on variable definitions in the FIA

Database can be found in the FIA Database Description and User Guide for Phase 2 [Burrill

et al., 2018].

12

Figure 1.5: Example code to produced population estimates of forest attributes grouped by
discrete categories using rFIA.

## Grouping by Forest Type
biomass(db = fiaRI, grpBy = FORTYPCD)

# A tibble: 120 x 20

<int>
103
104
105
167

<dbl>
4030.
3467.
3876.
2545.

<dbl>
3362.
2311.
3403.
1589.

<dbl>
77.8
74.3
85.2
60.6

YEAR FORTYPCD NETVOL_ACRE SAWVOL_ACRE BIO_AG_ACRE BIO_BG_ACRE BIO_ACRE
<dbl>
94.8
90.4

<int>
1 2014
2 2014
3 2014
4 2014
# ... with 116 more rows, and 13 more variables: CARB_AG_ACRE <dbl>,
#
#
#
#

CARB_BG_ACRE <dbl>, CARB_ACRE <dbl>, NETVOL_ACRE_SE <dbl>,
SAWVOL_ACRE_SE <dbl>, BIO_AG_ACRE_SE <dbl>, BIO_BG_ACRE_SE <dbl>,
BIO_ACRE_SE <dbl>, CARB_AG_ACRE_SE <dbl>, CARB_BG_ACRE_SE <dbl>,
CARB_ACRE_SE <dbl>, nPlots_VOL <dbl>, nPlots_AREA <dbl>

<dbl>
17.0
16.1
18.1
13.2

103.

73.8

1.3.3.2 Unique populations of interest

In many cases, the population of interest is a subset of that represented by the full FIA.Database.

For example, we may be interested in producing estimates for live stems greater than 12.0

inches DBH on state-owned land. Within biomass (and most other estimator functions), we

can use the treeDomain and areaDomain arguments to describe our population of interest in

terms of variables contained in the PLOT, COND, and TREE tables of a FIA.Database (Figure

1.6).

Here treeDomain describes the population of interest for the numerator (tree biomass)

and areaDomain the population of interest for the denominator (state-owned forest land

area). In each case, the argument takes the form of a logical predicate that is defined in terms

of the variables in PLOT, TREE, and/or COND tables of the FIA.Database. Multiple conditions

can be combined within either argument using & or | symbols (and/or, respectively).

13

Figure 1.6: Example code to produce population estimates of forest attributes within a
user-defined population of interest using rFIA.

## Live tree biomass (DBH > 12") on state land
biomass(fiaRI, treeDomain = DIA > 12, areaDomain = OWNCD == 31)

# A tibble: 5 x 19

<dbl>
1368.
1396.
1487.
1430.
1441.

<dbl>
1095.
1094.
1145.
1101.
1111.

YEAR NETVOL_ACRE SAWVOL_ACRE BIO_AG_ACRE BIO_BG_ACRE BIO_ACRE
<dbl>
41.7
42.7
45.6
44.4
44.7

<int>
1 2014
2 2015
3 2016
4 2017
5 2018
# ... with 13 more variables: CARB_AG_ACRE <dbl>,
#
#
#
#
#

CARB_BG_ACRE <dbl>, CARB_ACRE <dbl>, NETVOL_ACRE_SE <dbl>,
SAWVOL_ACRE_SE <dbl>, BIO_AG_ACRE_SE <dbl>,
BIO_BG_ACRE_SE <dbl>, BIO_ACRE_SE <dbl>, CARB_AG_ACRE_SE <dbl>,
CARB_BG_ACRE_SE <dbl>, CARB_ACRE_SE <dbl>, nPlots_TREE <dbl>,
nPlots_AREA <dbl>

<dbl>
34.8
35.6
38.0
37.0
37.3

<dbl>
6.96
7.13
7.61
7.39
7.45

1.3.3.3 Grouping by user-defined areal units

To produce estimates grouped by unique areal units, specify the name of a spatial polygon

object (class sp or sf) defining the areal units of interest as the polys argument of any rFIA

estimator function. All FIA data are automatically re-projected to match the projection

of the input spatial polygon object prior to initiating estimation procedures. All fields

originally contained in the input spatial polygon object will be preserved in the estimates

output by the estimator function. To return estimates as an sf spatial polygon object,

specify returnSpatial = TRUE.

1.3.3.4 Alternative estimators

All rFIA estimator functions use the “temporally-indifferent” estimator by default (method

= "TI") for consistency with other FIA estimation tools, like EVALIDator [USDA Forest

Service, 2019b, Pugh et al., 2018]. As an alternative, users may set the method argument to

14

Figure 1.7: Example code to produce population estimates of forest attributes for multiple
values of λ using the exponential moving average estimator in rFIA.

## Most recent year in Rhode Island
## Multiple lambda values
biomass(clipFIA(fiaRI), method = "EMA" ,

lambda = seq(from=0.1, to=0.9, by=0.1)

# A tibble: 9 x 20

<dbl>
2451.
2455.
2449.
2432.
2405.

<dbl>
1387.
1397.
1400.
1393.
1376.

<dbl>
68.8
69.2
69.2
68.9
68.4

<dbl> <int>
0.1 2018
0.2 2018
0.3 2018
0.4 2018
0.5 2018

lambda YEAR NETVOL_ACRE SAWVOL_ACRE BIO_AG_ACRE BIO_BG_ACRE
<dbl>
13.7
1
13.7
2
13.7
3
13.7
4
5
13.6
# ... with 4 more rows, and 13 more variables: BIO_ACRE <dbl>,
#
#
#
#
#

CARB_AG_ACRE <dbl>, CARB_BG_ACRE <dbl>, CARB_ACRE <dbl>,
NETVOL_ACRE_SE <dbl>, SAWVOL_ACRE_SE <dbl>, BIO_AG_ACRE_SE <dbl>,
BIO_BG_ACRE_SE <dbl>, BIO_ACRE_SE <dbl>, CARB_AG_ACRE_SE <dbl>,
CARB_BG_ACRE_SE <dbl>, CARB_ACRE_SE <dbl>, nPlots_TREE <dbl>,
nPlots_AREA <dbl>

"SMA", "LMA", "EMA", or "ANNUAL", to use the simple moving average, linear moving average,

exponential moving average, or annual estimator, respectively.

If using the exponential

moving average, users may also modify the exponential decay parameter, λ (Eq. 1.1), with

the lambda argument. By default, lambda is set to 0.5, although can be set to any value on

the interval (0, 1). If multiple values are specified as the lambda argument, one unique set of

estimates will be returned for each unique value of lambda (Figure 1.7). For example, lambda

= seq(from=0.1, to=0.9, by=0.1) will produce nine unique sets of estimates, grouped by

lambda.

1.3.3.5 Parallelization

Parallel processing is available in all rFIA estimator functions using the nCores argument.

nCores indicates the number of physical cores to be used and may be set to any positive

15

integer up to the number of physical cores available on a given machine. Serial processing

is implemented by default (nCores = 1). Parallelization may substantially reduce memory

during processing. Thus, users should consider implementing serial processing if computa-

tional resources are limited. If implementing parallel processing, we recommend users set

nCores to one less than the number of physical cores available on their machine to ensure

computational resources are available for other processes (e.g., OS, browsers).

1.4 Case study - Michigan ash decline

We demonstrate the utility of rFIA for forest resource monitoring by assessing the decline

of ash (Fraxinus spp.) populations across the Lower Peninsula of Michigan following the

establishment of emerald ash borer (Agrilus planipennis Fairmaire). Emerald ash borer is an

invasive forest insect that was first discovered in southeastern Michigan in 2002. The insect

spread quickly across the state and is considered one of the most destructive and costly forest

insects to invade the United States [Poland and McCullough, 2006, Aukema et al., 2011].

Using the “temporally indifferent” estimator in rFIA, we estimated annual changes in live

ash trees per acre (TPA) and tree mortality rates (mortality TPA per year) by county during

the interval 2006-2018. To highlight differences among alternative design-based estimators

implemented in rFIA, we use each estimator to assess changes in ash sawlog stocks (board feet

per acre) across the entire Lower Peninsula from 2000 to 2018 and compare their relative

performance (lambda = 0.5 used for EMA estimator). All estimates were produced for

white ash (Fraxinus americana L.), green ash (Fraxinus pennsylvanica Marshall ), and black
ash (Fraxinus nigra Marshall ) trees ≥ 5 inches DBH.

A rapid decline of ash populations across Michigan counties is evident from our analysis.

Live ash TPA decreased by 61% and annual mortality increased by 583% across the study

region from 2006 to 2018 (“temporally-indifferent” estimator). Elevated mortality rates

and population decline are apparent in the southeastern portion of the state soon after the

establishment of emerald ash borer in the region (Figure 1.8, 2006 and 2010).

In years

16

Figure 1.8: Changes in live ash TPA (left) and annual mortality rates (right), with
associated sampling errors (bottom; 67% confidence), across counties in the Lower
Peninsula of Michigan following establishment of emerald ash borer. Missing values are
shaded in gray. All plots were produced using the plotFIA function.

following, ash population decline and elevated mortality became evident across the Lower

Peninsula, likely associated with the rapid expansion of emerald ash borer across the state

(Figure 1.8, 2014 and 2018).

Sampling error is directly related to the number of non-zero observations (i.e., FIA ground

plots) used to compute estimates of each variable, and cannot be determined if less than two

non-zero observations are available. Decreasing precision (i.e., increased sampling errors)

associated with live TPA (Figure 1.8; bottom-left) and sawlog volume (Figure 1.9; bottom)

estimates indicates that the number of live ash trees located on FIA ground plots decreased

over time. Ash mortality appeared to be a relatively rare event prior to the expansion of

emerald ash borer across Michigan, as many counties recorded fewer than two observations of

ash mortality per sampling period (gray shaded counties in bottom right of Figure 1.8). The

increased frequency of ash mortality observations across Michigan counties, and associated

reduction in sampling error, provides support for the trend of increased mortality rates that

17

Figure 1.9: Changes in ash sawlog volume (top; board feet (BF) per acre) and associated
sampling errors (bottom; 67% confidence) for each design-based estimator, across the
Lower Peninsula of Michigan from 2000 to 2018. All plots were produced using the
plotFIA function.

is evident following the establishment of emerald ash borer across the state.

A trade-off between temporal smoothing and precision is evident in comparing the behav-

ior of design-based estimators implemented in rFIA (Figure 1.9). For example, the annual

estimator (ANNUAL) appears to indicate a sharper and earlier decline (lower temporal

smoothing) in the ash sawlog resource than the “temporally-indifferent” (TI) and simple

moving average (SMA) estimators, though the annual estimator consistently produces higher

sampling errors (lower precision) (Figure 1.9). Estimators that give more weight to recent

observations, like the annual (100% of weight on most recent observation), linear moving

average (LMA), and exponential moving average (EMA) tend to exhibit lower temporal lag-

bias but at the cost of higher variance. In contrast, estimators that distribute weight more

evenly across observations with respect to time, like the “temporally-indifferent” (no annual

weighting, treated as periodic inventory) and simple moving average tend to exhibit higher

18

temporal smoothing but lower variance. In practice, we recommend users carefully consider

this trade-off and choose an estimator that achieves both their desired precision standards

and an acceptably low level of temporal lag-bias.

1.5 Extensions/limitations

In the future, we hope to expand rFIA to include model-based and model-assisted tech-

niques, allowing users to leverage various forms of auxiliary data to improve estimation of

forest attributes. Specifically, we hope to improve estimation within small domains (e.g.,

spatial/ temporal extents, rare forest attributes) by implementing a number of indirect and

composite small area estimators. Further, we hope to improve the value of the FIA Database

for detection and analysis of rapid changes in forest health through the implementation of

Kalman filters as time-series estimators.

Program R is an “in-memory” application by design, and rFIA requires all input data to

be held in RAM. Hence, users may experience memory management challenges when com-

putational resources are limited. Memory constraints may inhibit some users from loading

or processing large subsets of FIA data (relative to available computational resources), and

high memory usage may result in decreased computational efficiency. Future development

seeks to leverage dplyr [Wickham et al., 2015] back ends for external database engines such

as SQLite and Apache Spark, thus reducing memory limitations and potentially offering

substantial improvements in computational efficiency for large data.

Exact coordinates of FIA ground plots are not available in the public version of the FIA

Database to protect privacy rights of private land owners and preserve the ecological integrity

of ground plots [Tinkham et al., 2018, McRoberts et al., 2005a]. Plot locations are randomly

displaced up to 1 km from their true locations (i.e., “fuzzed”) and the coordinates of up to

20% of plots on private land in each county are exchanged (i.e., “swapped”) [Bechtold et al.,

2005, Tinkham et al., 2018]. The effects of using “fuzzed and swapped” plot coordinates are

thought to be negligible for design-based estimation across large regions (800,000 acres and

19

larger) [McRoberts et al., 2005b], however uncertainty in true plot locations may produce

results within unknown amounts of error when used in conjunction with other spatially

explicit data layers for model-based estimation [Sabor et al., 2007]. Hence, we encourage rFIA

users to consider this uncertainty when using design-based estimation for small populations,

and when producing plot-level estimates to be used with spatially explicit auxiliary data in

a model-based or model-assisted estimation framework.

1.6 Conclusion

The FIA database is among the most valuable ecological datasets in the world, though its

complexity limits access for many potential users. We developed rFIA to simplify the estima-

tion of forest attributes using FIA data, intending to provide a flexible, yet powerful toolset

which is accessible to all R users. Ultimately, we hope that rFIA will improve the accessibility

and relevance of the FIA database and promote use of FIA data among a larger, more diverse

audience. We encourage users to apply rFIA widely and report any issues and/or desired

extensions on our active issues page (https://github.com/hunter-stanke/rFIA/issues).

1.7 Acknowledgements

This work was supported by National Science foundation grants DMS-1916395, EF-

1241874, EF-1253225, National Aeronautics and Space Administration Carbon Monitoring

System program, and United States National Park Service and Forest Service.

20

CHAPTER 2

OVER HALF OF TOP TREE SPECIES IN DECLINE IN THE WESTERN

UNITED STATES

2.1 Introduction

Persistent shifts in forest composition, structure, and function depend largely on the

demographic response of tree species to changing environmental drivers and disturbance

regimes [Johnstone et al., 2016]. Many forests in the western United States (US) have

experienced recent increases in the extent, severity, and frequency of wildfire [Turner et al.,

2019, Stevens-Rumann et al., 2018], drought [Clark et al., 2016, Cook et al., 2004], and insect-

pest outbreaks [Wong and Daniels, 2017, Bentz et al., 2010], due in part to changing climate

and past forest management (i.e., fire suppression). Likewise, large-scale tree mortality

events [Van Mantgem et al., 2009, Breshears et al., 2005] and recruitment failures [Davis

et al., 2019, Harvey et al., 2016] indicate that widespread forest change is already underway

in the region. In the face of such dramatic change, efforts to quantify the joint demographic

response of tree species to anthropogenic and natural stressors are urgently needed to inform

forest management and improve predictions of shifts in forest ecosystem services across the

western US.

Many previous studies have quantified individual-tree demographic responses to novel

stressors (e.g., global change type drought) in the western US [Van Mantgem et al., 2009,

Ma et al., 2012, Py et al., 2006, Adams and Kolb, 2005, 2004], though few have examined

such responses at higher levels of biological organization [Clark et al., 2016] and most only

consider a single facet of tree demography. It is difficult to use knowledge inferred from trees

to predict stand-level or landscape-level response to stress due to inter-tree interactions (i.e.,

competition, shading) [Clark et al., 2014]. Further, individual facets of tree demography

(e.g., mortality, growth, recruitment) are known to interact to influence forest population

21

dynamics, and hence no single facet can fully explain variation in such dynamics [Franklin

et al., 2002, Connell et al., 1984]. That is, increased tree mortality rates do not imply

population decline if increased mortality is matched by a subsequent increase in productivity

(i.e., growth, recruitment) [Assmann, 2013]. As such, a great deal of uncertainty remains

regarding the current status of tree populations (e.g., expanding, declining) and how multiple

drivers may interact to influence the stability of tree populations across broad spatial domains

in the western US.

Multiple studies have suggested changes in the spatial distribution of tree species may

be underway in the western US [Monleon and Lintz, 2015, Serra-Diaz et al., 2016, Bell

et al., 2014a], driven by shifts in population dynamics (e.g., failed recruitment) at leading

and trailing range margins. However, most efforts to quantify species range shifts rely on

“life stage for time” substitution, wherein the potential for range shifts is measured as the

magnitude of divergence between adult and juvenile tree distributions across geographic

and/or bioclimatic spaces at a single point in time [Lenoir et al., 2009, Zhu et al., 2014,

2012]. Trees of each life stage established at different times and potentially under different

environmental settings, and divergent distributions between tree life stages may arise from

a changing environmental conditions over time. However, realized species niches often differ

across tree life stages and change through the ontogeny of trees [Grubb, 1977, Miriti, 2006,

Bertrand et al., 2011]. Hence it is possible that divergence in distributions across tree

life stages could reflect ontogenetic differences rather than the effects of exogenous drivers

(e.g., climate), and additional efforts that explicitly account for population dynamics in the

estimation of species range shifts are warranted [M´aliˇs et al., 2016].

Tree population performance (i.e., net change in abundance over time; population-level

expansion, decline) emerges from the joint demographic response of individual trees to en-

dogenous (e.g., inter-tree competition) and exogenous drivers (e.g., wildfire). Disparities in

tree population performance across species, life stages, and environmental gradients may

provide early indications of large-scale forest change, including species range shifts [Jackson

22

et al., 2009, Woodall et al., 2013, Fei et al., 2017]. For example, increased tree population

performance at high elevations relative to low elevations may indicate upslope migration

over time (i.e., upslope shift in abundance distribution), and inter-specific divergence in tree

population performance may indicate shifts in forest composition within a region. Further,

intra-specific divergence in population performance across tree life stages may predict future

shifts in forest structure (i.e., shifts in abundance distribution across life stages within a

region). However, comprehensive assessments of tree population performance have yet to be

conducted in the western US, likely because such assessments require detailed, temporally

replicated data spanning large spatial domains, which were unavailable in the region until

recently [Goeking, 2015].

We sought to quantify the population performance of 10 tree species and two size classes

across their ranges in the western US, and identify potential drivers of population instabil-

ity (i.e., expansion or decline) where it exists. We define forest population performance in

terms of shifts in abundance over time, where population stability is achieved by a balance

between mortality, recruitment, and growth. That is, when mortality significantly exceeds

recruitment and growth, abundance shifts downward and population decline is evident, char-

acterized by decreased density and/or spatial extent (i.e., range contraction). In contrast, if

recruitment and growth significantly exceed mortality, abundance shifts upward and popu-

lation expansion occurs, characterized by increased density and/or spatial extent (i.e., range

expansion). We adapt a method presented in Lintz et al. [2016] to account for transient

dynamics in tree demography associated with forest structural development and succession,

developing an index of forest population performance that is independent of forest composi-

tion, stand age, and ecological setting.

We draw upon over 24,000 repeated censuses of US Forest Service Forest Inventory and

Analysis (FIA) plots to answer the following questions: (1) Does range-wide population

performance differ among the 10 most abundant tree species in the western US, and what does

spatial variation in population performance indicate about shifts in species distributions?

23

(2) Do tree size classes exhibit intra-specific variation in population performance, and if so,

which regions experience the greatest differences? (3) How do major forest disturbances and

climate influence the population performance of tree species and size classes, and what do

these relationships suggest about the drivers of species decline where it exists?

2.2 Materials and Methods

2.2.1 Field observations

Since 1999, the FIA program has operated an extensive, nationally-consistent forest inventory

designed to monitor changes in forests across all lands in US [USDA Forest Service, 2019a].

We used FIA data from 10 states in the western US (Washington, Oregon, California, Idaho,

Montanta, Utah, Nevada, Colorodo, Arizona, and New Mexico) to quantify forest population

performance, excluding Wyoming due to a lack of repeated censuses. This region spans a

wide variety of climatic regimes and forest types, ranging from temperate rain forests of the

coastal Pacific Northwest to pinyon-juniper woodlands of the interior southwest [Eyre, 1980].

Although the spatial extent of the FIA plot network represents a large portion of the current

range of all species examined in this study, substantial portions of some species ranges (e.g.,

western hemlock) extend beyond the study region into Canada and/or Mexico and therefore

were not fully addressed here.

The FIA program measures forest attributes on a network of permanent ground plots that

are systematically distributed at a rate of approximately 1 plot per 2500 hectares across the

US [Smith, 2002]. For trees 12.7cm diameter at breast height (DBH) and larger, attributes

(e.g., species, DBH, live/dead, mortality agent) are measured on a cluster of four 168m2
subplots [Bechtold et al., 2005]. Trees 2.54 − 12.7cm DBH are measured on a a microplot
(13.5m2) contained within each subplot, and rare events such as very large trees are measured

on an optional macroplot (1012m2) surrounding each subplot [Bechtold et al., 2005]. In the

western US, one-tenth of ground plots are measured each year, with re-measurements first

occurring in 2011.

24

We characterized forest disturbance severity at each remeasured FIA plot as the annual

fraction of mortality, in terms of tree basal area per unit area (BAA), attributed by FIA

field crews to one of four mortality agents during the re-measurement interval: fire, insects,

disease, and harvest/ landclearing. Specifically, we estimate disturbance severity for each

mortality agent as:

dj =

2
r

BAAmort,j

BAA2,live + BAA1,live

(2.1)

where dj is disturbance severity for mortality agent j, BAAmort,j is the the total BAA killed

by mortality agent j during the re-measurement interval r (including stems that recruited

and subsequently died within the measurement interval), and BAA1,live and BAA2,live is the

total live BAA at initial and final measurements, respectively. Hence, we represent annual

change in BAA due to each mortality agent with respect to the average live BAA between

measurements, yielding a symmetric index of change that accommodates initial measurement

equal of (i.e., no trees present at initial measurement) [T¨ornqvist et al., 1985].

2.2.2 Climate data

We acquired 30 year precipitation and mean temperature normals from the 4-km resolution

PRISM climate dataset [PRISM Climate Group, 2010] and extracted data for each FIA plot

location to generate long-term climate variables. While forests may respond to variation

in climate at relatively fine spatial scales [Millar et al., 2018, Case and Peterson, 2005],

particularly in regions with complex topography, coarse resolution climate data are often

adequate to characterize broad-scale climatic controls on forest processes [Ashcroft, 2010,

Pearson and Dawson, 2003], such as those addressed here. Furthermore, our ability to

account for fine scale climatic microrefugia that may not be represented by broad scale

climate data was limited by precision in FIA plot locations, which are randomly displaced

25

up to 1 km from their true locations to protect privacy rights of private land owners [Tinkham

et al., 2018].

We used the Standardized Precipitation Evapotransporitive Index (SPEI) to characterize

drought severity between repeated measurements on each FIA plot. The SPEI is a multi-

scalar drought index widely used for assessing the impacts of drought on forests and other

terrestrial ecosystems [Wu et al., 2018, Vicente-Serrano et al., 2013, Jiang et al., 2019, Huang

et al., 2015, Greenwood et al., 2017]. The SPEI has advantages over other drought indices

(e.g., Palmer Drought Severity Index) as it allows direct comparison of drought severity

across geographic space and can be used to identify drought on a variety of time-scales

[Vicente-Serrano et al., 2010]. We computed the SPEI using 4-km monthly precipitation

and mean temperature estimates from the PRISM climate dataset [PRISM Climate Group,

2010]. Specifically, we computed the SPEI on a monthly basis with an 18 month time-

lag for the period 1999-2018 using the SPEI R package [Beguer´ıa et al., 2014, Beguer´ıa

and Vicente-Serrano, 2017], estimating potential evapotranspiration using the Thornthwaite

method [Thornthwaite, 1948]. We chose to compute the SPEI with an 18 month time lag as

many forest processes (i.e., mortality, growth) commonly exhibit lagged responses to drought,

ranging from months to several years [Vanoni et al., 2016, Bigler et al., 2007, Klockow et al.,

2018, Clark et al., 2016]. We defined drought severity at each FIA plot as the mean monthly

SPEI experienced at each location during the re-measurement interval.

2.2.3 Forest stability index

We adapt a method developed by Lintz et al. [2016] to produce an index of forest population

performance for each remeasured FIA plot. We first estimate annual rate of change in live

trees per unit area (∆T P A) and live basal area per unit area (∆BAA) as change in each

variable between measurements divided by the re-measurement interval on each plot. We

then scale each rate of change by dividing by the respective sample standard deviation:

26

Figure 2.1: Example relationship of change in TPA (S, x-axis) and change in BAA (B,
y-axis) for subalpine fir (Abies lasiocarpa) in Colorado (all plots remeasured from
2015-2018). Small diameter circles represent plot level-indices of S and B, and the single
large diamond represents the region-wide mean for each dimension. Small squares
distributed on the 1:1 line (gray-dashed line) represent plot-level FSI values corresponding
to six example plots in the upper-right of the graph (connected by by black segments). All
points are colored by their respective FSI value. Points in the bottom-left of the graph are
characterized by losses in both TPA and BAA, indicating the species is likely declining
within these stands. Conversely, points in the upper-right of the graph are characterized by
increases in TPA and BAA, indicating species is likely expanding in the stand (i.e.,
increasing in density).

Si =

∆T P Ai

sd(∆T P A)

Bi =

∆BAAi

sd(∆BAA)

(2.2)

(2.3)

where Bi is scaled annual change in BAA on plot i, Si is scaled annual change in T P A on

plot i, and sd(x) denotes the standard deviation of variable x. We represent change on a

unit area basis to account for differences in sampling areas between tree size classes. We

use the sample standard deviation to scale ∆T P A and ∆BAA so that each rate has equal

variance (σ = 1), and hence one unit change in S is equivalent in magnitude to one unit

27

change in B.

We project plot-level S and B perpendicularly onto the 1:1 line to produce a single metric

of population performance, herein referred to as the Forest Stability Index (FSI) (Figure 2.1;

small squares). Reducing S and B to a single dimension (i.e., the FSI) facilitates summary

across geographic, temporal, and bioclimatic spaces, and simplifies interpretation of the

direction and magnitude of population change. Summarized across a large sample of plots

(i.e., a meta-population), the FSI supports the following premise: if net changes in BAA and

TPA are significantly greater than 0 (Figure 2.1, top-right), the underlying population is

likely expanding in the region (i.e., increasing in density and/or spatial extent). Conversely,

if net changes in BAA and TPA are significantly less than 0 (Figure 2.1, bottom-left),

population decline is likely and further examination to identify potential drivers is warranted.

Finally, if net changes in BAA and TPA are not significantly different from zero (Figure 2.1,

center), change in the underlying population is negligible (i.e., the population is stable).

Standardized changes in TPA (S) and BAA (B) may be directly compared among species

and across site conditions under the premise that mortality, recruitment, and growth rates,

which implicitly define S and B, will be balanced for any stable population (i.e., net change

of zero).

Inherent differences among species life-history traits (e.g., shade tolerance) and

spatial variation in resource availability (e.g., light, moisture, nutrients) maintain a demo-

graphic trade-off between tree growth and mortality rates [Assmann, 2013, Wright et al.,

2004, Grubb, 1977, Pacala et al., 1996]. That is, variation in mortality rates are offset by

variation in growth rates across species and site conditions, with fast-growing species and

high productivity sites tending to experience higher mortality rates than slow-growing species

and low productivity sites.

Independently, S and B are confounded by transient dynamics associated with stand

structural development and succession, however we argue that these dynamics are offsetting

when both indices are considered simultaneously using the FSI. For example, mortality in

young forests may exceed recruitment due to inter-tree competition (negative S), albeit this

28

loss is offset by growth on surviving trees (neutral or positive B) and is therefore sustainable

(Figure 2.1, top-left) [Franklin et al., 2002]. Furthermore, mortality of individual large trees

is likely to exceed recruitment and growth among surviving stems when measured in terms of

basal area in old forest (negative B), where basal area is often distributed disproportionately

among large trees [Parobekov´a et al., 2018]. However in the same scenario, mortality may be

far exceeded by recruitment in terms of trees per acre, as the loss of an individual large tree

may coincide with recruitment of many young trees (positive S) (Figure 2.1, bottom-right).

2.2.4 Statistical analysis

We computed the FSI by species and by tree size class for all remeasured FIA plots in the

western US (N = 24, 229). We included plots on both public and private lands and included
all live stems (DBH ≥ 1 inch). All plot measurements included in our analysis occurred
from 2001-2018, with an average re-measurement interval of 9.78 years (± 0.005 years). For

brevity, we restricted our analysis to consider the 10 most abundant tree species in the

western US and simplified tree size classifications to include two classes: small diameter
(2.54-12.6 cm DBH) and large diameter (DBH ≥ 12.7 cm). We identified the most abundant

tree species across our study area using the rFIA R package [Stanke et al., 2020], defining
abundance in terms of estimated total number of trees (DBH ≥ 2.54 cm) in the year 2018. We

excluded species that exhibit non-tree (i.e., shrub, subshrub) growth habits across portions

of the study region (e.g., Quercus gambelii ).

We applied stratified random estimators described by the FIA program [Bechtold et al.,

2005] to estimate the range-wide mean and variance of the FSI for each population defined

by species and tree size class. Specifically, we used the simple moving average estimator

implemented in the rFIA R package [Stanke et al., 2020] to compute estimates from a series of

eight annual panels (i.e., sets of plots re-measured in the same year) ranging from 2011-2018.

The simple moving average estimator combines information from annual panels with equal

weight (i.e., irrespective of time since re-measurement), thereby allowing us to characterize

29

long-term patterns in population performance. We determine populations to be stable if

the 95% confidence intervals for estimated range-wide FSI included zero. Alternatively,

if confidence intervals of estimated range-wide FSI do not include zero, we determine the

population to be expanding when the FSI is positive and declining when the FSI is negative.

We assessed geographic variation in species population performance, independent of tree

size class, by estimating the FSI for each species within ecoregion subsections [Moore et al.,

2016] using the simple moving average estimator. As a direct measure of changes in abun-

dance, spatial variation in the FSI is indicative of spatial shifts in species distributions during

the re-measurement interval (i.e., range expansion/contraction and/or within-range abun-

dance shifts). That is, the distribution of any population shifts towards regions of high

population performance and away from regions of low population performance during the

temporal frame of sampling. We map estimates of the FSI for each ecoregion subsection

to assess spatial patterns in species population performance across their range and identify

regions where widespread geographic shifts in species distributions may be underway.

We determined if range-wide population performance of individual species diverged sig-

nificantly between tree size classes by applying Kruskall-Wallis tests to plot-level FSI values

grouped by species and size class. We assessed geographic variation in this intra-specific

divergence in population performance by applying Kruskal-Wallis tests to the same plot-

level FSI values within each ecoregion subsection. Hence for each species, we assessed if

mean FSI values differed significantly between tree size classes within each ecoregion sub-

section included in its range. We then map significant differences to assess spatial patterns

of divergence in species population performance between tree size classes.

We sought to determine the relative importance of forest disturbance processes and long-

term climate patterns (Figure 2.2) in shaping the population performance of species and

tree size classes across our study region. We used linear mixed-effects models to predict

plot-level FSI as a function of forest disturbance severity rates (i.e., insects, fire, disease,

harvest/landclearing, and drought) and long-term climate normals (i.e., precipitation and

30

Figure 2.2: Spatial variation in forest disturbance severity and long-term climate patterns
across the study region. All variables standardized to a common scale for visual
interpretation of spatial variation.

mean temperature). We included species and tree size as random effects and allowed slopes

of each predictor to vary across groups. We ensured minimal collinearity between predictor

variables using variance inflation factors. We standardized all predictors at the species level

prior to analysis, and compared coefficients estimated by the model to interpret the relative

importance each predictor in explaining variation in the FSI.

2.3 Results

We sought to quantify the population performance of the 10 most abundant tree species

across their ranges in the Western US, and identify drivers of instability where it exists.

31

Table 2.1: Scientific and common names of ten most abundant species across our study
area, listed in order of increasing abundance (i.e., bottom being most abundant; abundance
defined by estimated total number of trees). Sample size (number of re-measured FIA plots
where the species occured) provided for small diameter (2.54-12.6 cm DBH), large diameter
(DBH ≥ 12.7 cm DBH), and combined populations (DBH ≥ 2.54 cm). As large diameter
and small diameter populations of a species may exist on the same plot, the sample size of
the combined population is not equal to the sum of large diameter and small diameter
sample sizes.

Number of plots

Scientific Name
Tsuga heterophylla
Abies grandis
Juniperus osteosperma

Common name
Western hemlock
Grand fir
Utah juniper
Engelmann spruce Picea Engelmannii
Quaking aspen
Common pinyon
Ponderosa pine
Subalpine fir
Lodgepole pine
Douglas-fir

Populus tremuloides
Pinus edulis
Pinus ponderosa
Abies lasiocarpa
Pinus contorta
Pseudotsuga menziesii

Large diameter Small diameter Combined
3262
2734
3446
3079
1723
3076
7309
3174
4556
12284

3149
2580
3425
2965
1563
2998
7177
3003
4388
12086

1617
1439
891
1287
935
1423
2352
1964
2276
4688

We identified the 10 most abundant tree species by their estimated total number of stems,

including 7 distinct genera and 3 families (Table 2.1). Together these top 10 species accounted

for 68.9% of all trees across the study region (64.4% of total basal area) (Figure 2.3). As

expected, small diameter stems were more abundant than large diameter stems for most

species in terms of tree number. The relative abundance distributions provided in Figure 2.3

are intended to aid in interpretation of assessments of population performance that follow,

providing a reference upon which we can base assessments of population change.

Differences in range-wide mean FSI illustrate substantial variation in population perfor-

mance among species. Irrespective of tree size class, we found 70% of species populations

in decline and the remaining 30% expanding across the study region (Figure 2.4, left). The

highest rates of decline appeared in Engelmann spruce and subalpine fir populations, char-

acteristic species of moist, high-elevation regions of the Rocky Mountains. In contrast, we

observed the highest rates of population expansion in grand fir and western hemlock, which

32

Figure 2.3: Region-wide prevalence (% of total tree population across the study region) of
the 10 most abundant tree species in the western US in 2018. Estimates provided for
combined (left), large diameter (center), and small diameter populations (right). Species
listed in order of increasing abundance with respect to the combined population (i.e.,
bottom being most abundant).

occur predominantly in low-mid elevation forests of the Pacific Northwest (i.e., Washington,

Oregon, Idaho, and western Montana).

Qualitatively, we observed considerable variation in the spatial patterns of population

performance among species (Figure 2.5 (Top)). Populations of Engelmann spruce, Douglas-

fir, and ponderosa pine appear to be generally expanding in the Pacific Northwest and Sierra

Nevadas, while subalpine fir and lodgepole pine are generally declining in the region. Grand

fir and western hemlock appear to be expanding predominantly in the Northern Rockies

and Interior Pacific Northwest, and declining in the coastal Pacific Northwest. Widespread

decline of lodgepole pine is evident in the Central Rocky Mountains, while no clear pattern

emerges for Douglas-fir, ponderosa pine, Engelmann spruce, and subalpine fir in the region

(i.e., both expansion and decline occur frequently). Decline appears to be wide-spread among

Engelmann spruce, subalpine fir, and ponderosa pine in the Sourthern Rocky Mountains.

Quaking aspen populations generally appear to be in decline in the southern portions of its

33

Figure 2.4: Range-wide mean Forest Stability Index (FSI) of the 10 most abundant tree
species in the western US. Estimates provided for combined (left), large diameter (center),
and small diameter populations (right). Population decline (red) occurs when the FSI is
negative and the associated confidence interval does not include zero. Conversely,
population expansion (blue) occurs when the FSI is positive and the associated confidence
interval does not include zero.

range (Southern Rocky Mountains) and expanding in the northern portions (Northern Rocky

Mountains, Interior Pacific Northwest). It appears Utah juniper tends to be declining near

the spatial margins of its range (i.e., decline appears more frequently in outermost areas of

its range) and expanding in the central portion of its range. The opposite pattern appears for

common pinyon, with expansion prevalent near the spatial margins of its range and decline

prevalent in the central portion.

We found evidence of intra-specific divergence in range-wide population performance

between tree size classes for 80% of species (Figure 2.4, Large Diameter vs. Small Diameter).
Specifically, we found significant differences (P ≤ 0.05) in mean FSI between small diameter

and large diameter populations for all species except Utah juniper and Engelmann spruce

34

Figure 2.5: (Top) Spatial variation in population performance of species across their
ranges. Forest Stability Index (FSI) summarized within ecoregion subsections using the
simple moving average estimator. Population expansion occurs if the mean FSI
significantly exceeds zero within an area (net increase in abundance), and population
decline occurs if the mean FSI is significantly less than zero (net decrease in abundance).
Populations are determined to be stable within an area if the confidence interval of
estimated mean FSI includes zero. (Bottom) Spatial variation in intra-specific differences
in population performance between tree size classes. Areas where small diameter
populations exhibit significantly higher FSI than large diameter populations shaded in
purple. Conversely, areas where small diameter populations exhibit significantly lower FSI
than large diameter populations shaded in yellow. Otherwise populations are determined
to be changing at equal rates (green).

35

Table 2.2: Intra-specific differences in range-wide mean Forest Stability Index (FSI)
between tree size classes. Significance (i.e., p-values) derived from Kruskall-Wallace test
using plot-level FSI values.

FSI ± 95% CI

Species
Grand fir
Western hemlock
Quaking aspen
Utah juniper
Douglas-fir
Ponderosa pine
Common pinyon
Lodgepole pine
Subalpine fir
Engelmann spruce

Large diameter Small diameter Difference P Value
0.36 ± 0.004
≤ 0.001
0.05 ± 0.004
≤ 0.001
≤ 0.001
-1.44 ± 0.002
0.04 ± 0.001
0.413
0.25 ± 0.001
≤ 0.001
-0.15 ± 0.001
≤ 0.001
-0.22 ± 0.001
0.014
≤ 0.001
-2.38 ± 0.003
-0.62 ± 0.002
≤ 0.001
-0.85 ± 0.002
0.281

2.36 ± 0.020
1.58 ± 0.024
2.63 ± 0.020
-0.22 ± 0.004
-0.87 ± 0.004
-0.39 ± 0.008
-0.44 ± 0.006
3.61 ± 0.020
-0.39 ± 0.008
-0.44 ± 0.005

-2.00
-1.53
-4.07
0.26
1.12
0.24
0.22
-5.99
-0.23
-0.41

(Table 2.2). The FSI of small diameter populations significantly exceeded that of large

diameter populations of grand fir, western hemlock, quaking aspen, lodgepole pine, and

subalpine fir. In contrast, the FSI of large diameter populations significantly exceeded that

of small diameter populations of Douglas-fir, ponderosa pine, and common pinyon.

We used Kruskall-Wallace tests to identify regions (i.e., ecoregion subsections) with sig-

nificant intra-specific differences in mean population performance between tree size classes.

We found no significant divergence in population performance across the majority of each

species range (2.5 (bottom); Equal rates), however general spatial patterns of divergence do

appear for some species. Specifically, the FSI of large diameter Douglas-fir populations ex-

ceeds that of small diameter populations frequently in the coastal Pacific Northwest. Other

areas where the FSI of large diameter populations exceeds that of small diameter populations

were observed sporadically across the ranges of Douglas-fir, Utah juniper, western hemlock,

grand fir, ponderosa pine, lodgepole pine, and Engelmann spruce. The opposite, where the

FSI of small diameter populations exceeds that of large diameter populations, was observed

sporadically across the ranges of all species except Utah juniper.

Results of our linear mixed model reveal that forest disturbance processes are generally

36

Figure 2.6: Standardized coefficients of the linear mixed-effect model predicting plot-level
FSI as a function of forest disturbance severity and long-term climate patterns. All
variables were scaled at the species level, hence the relative importance of each predictor in
explaining intra-specific variation in population performance may be directly inferred from
the magnitude of each coefficient. That is, coefficient values can be directly compared
across tree size classes (i.e., vertically) but not across species (i.e., horizontally). Blue
shaded backgrounds indicate the population is expanding across the study region (e.g.,
large diameter grand fir; top left). Red shaded backgrounds indicate the population is
declining across the study region (e.g., large diameter Engelmann spruce; top right). Error
bars represent the 95% confidence intervals in the estimate of each coefficient.

much more important than long-term climate patterns in shaping the population perfor-

mance of species and tree size classes (Figure 2.6). We found a negative relationship be-

tween the FSI and all forest disturbance severity variables except drought, which exhibited

both positive and negative relationships with the FSI dependent upon species and tree size

class. That is, increases in severity of all forest disturbances, with the exception of drought,

coincided with decreased FSI values (i.e., higher rates of decline or lower rates of expansion)

for all species and tree size classes. In contrast, we found that long-term climate variables

(i.e., long-term precipitation and mean temperature) vary considerably in their relationship

to the FSI across species and tree size classes.

Among large diameter stems, we found fire severity was among the two most important

predictors of the FSI for all species except western hemlock (Figure 2.6). We found harvest/

37

landclearing to be of high importance (i.e., top two predictors) for large diameter grand

fir, western hemlock, Douglas-fir, Utah juniper, and ponderosa pine, all of which except

ponderosa pine exhibited range-wide expansion. Insect outbreak severity was of high relative

importance for large diameter common pinyon, lodgepole pine, subalpine fir, and Engelmann

spruce, all of which we determined to be undergoing range-wide population decline. We found

disease severity to be the most important predictor of large diameter quaking aspen FSI,

but less important for all other species. We found drought to be among the least important

forest disturbance processes associated with large diameter population performance across

all species.

We found the relative importance of forest disturbance processes in relation to the pop-

ulation performance to be similar for each species across size classes, although some notable

differences did appear (Figure 2.6). We found that insect outbreak severity was substan-

tially more important in predicting the FSI of small diameter populations of all species

except western hemlock, Douglas-fir, and Engelmann spruce, where Engelmann spruce ex-

hibited the opposite relationship (i.e., insect outbreak more important for large-diameter

than small diameter). Similarly, we found that disease severity was substantially more im-

portant in predicting the FSI of small diameter populations of grand fir, quaking aspen,

ponderosa pine, lodgepole pine, subalpine fir, relative to that of large diameter populations

of the same species. We found the importance of harvest/ landclearing to be markedly lower

in small diameter populations of species for which it was identified as an important predic-

tor of FSI in large diameter populations. Finally, we found drought to be a more important

predictor of small diameter populations of all species except ponderosa pine, Utah juniper,

Douglas-fir, and Engelmann spruce, where Douglas-fir and Engelmann spruce exhibited the

opposite relationship.

While generally less important than disturbance processes, we found significant relation-

ships between long-term climate variables and population performance for multiple species

and tree size classes. The FSI of large and small diameter populations of grand fir and

38

western hemlock, as well as populations of small diameter Douglas-fir and lodgepole pine,

exhibited significant negative relationships with long-term mean precipitation, indicating

improved population performance in dry regions. We observed the opposite relationship

in large diameter Douglas-fir and small diameter quaking aspen (i.e., improved population

performance in wet regions). Further, we observed significant positive relationships between

the FSI and large diameter populations of western hemlock, Douglas-fir, and Engelmann

spruce, as well as small diameter populations of quaking aspen and lodgepole pine, indicat-

ing improved population performance in warm regions. Neither large and small diameter

populations of Utah juniper, ponderosa pine, common pinyon, and subalpine fir exhibited a

significant relationship with long-term climate variables.

2.4 Discussion

2.4.1 Range-wide population performance and regional shifts in forest compo-

sition

Our results indicate a majority (7 of 10) of the most abundant tree species in the western US

are undergoing significant range-wide population decline (Figure 2.4). This result coincides

with recent reports of widespread tree mortality across the region [Van Mantgem et al., 2009,

Breshears et al., 2005, Bigler et al., 2007, Lintz et al., 2016], suggesting that heightened tree

mortality is likely driving forest decline in the western US. However, our results indicate that

population performance varies considerably between species, and that some species (grand

fir, western hemlock, and quaking aspen) are experiencing significant range-wide population

expansion. Such dramatic divergence in range-wide population performance between tree

species may provide early indications of broad-scale shifts in forest composition across the

western US, where high performing species increase in abundance relative to lower performing

species.

We observed the highest rates of range-wide population decline in Engelmann spruce,

subalpine fir, and lodgepole pine, all of which tend to occur at high-elevation near the upper

39

forest ecotone across the western US [Peet, 1981]. Combined, these species represent 25.3%

of all trees across the region (Figure 2.3). Previous efforts have indicated that subalpine tree

species are among the most vulnerable to future changes in climate and forest disturbance

regimes [Bell et al., 2014b, Conlisk et al., 2017, Wong and Daniels, 2017, Bigler et al., 2007,

Malone et al., 2018]. Our results indicate this heightened vulnerability may already be

manifesting across the western US, serving as an early warning sign of widespread, pervasive

decline of subalpine forests.

In contrast, we observed the highest rates of range-wide population expansion in grand

fir and western hemlock (Figure 2.4), together representing 7.3% of all trees across the

western US. Both grand fir and western hemlock are extremely shade tolerant, late seral

species occurring predominately on cool, mesic sites in the northern portion of our study

region [Burns, 1990, Peet, 1981]. Population expansion of shade tolerant species is consistent

with compositional shifts expected under accelerated secondary succession [Franklin et al.,

2002]. Across the dry domain of the western US, decades of fire exclusion and selective

harvesting have created crowded, light-competitive stand conditions that favor the expansion

of shade tolerant tree species [Hessburg et al., 2005, Hessburg and Agee, 2003, Naficy et al.,

2010, Gallant et al., 2003]. Hence, we argue the population expansion of grand fir and

western hemlock emerges as an artifact of past forest management (i.e., fire suppression,

selective silvicultural systems), and presents a prominent challenge for contemporary dry

forest restoration.

We found population performance was most stable among characteristically montane

(i.e., Douglas-fir, ponderosa pine, quaking aspen) and woodland species (i.e., Utah juniper,

common pinyon) in the western US, indicated by FSI values closest to zero. Significant range-

wide decline was evident in all such species except quaking aspen, which was determined to be

expanding slightly across its range (Figure 2.4). Our results concur with those presented by

Lintz et al. [2016], who reported higher levels of density-independent mortality in subalpine

species relative to montane and woodland species.

It is important to note that the high

40

frequency of decline among montane and woodland species, albeit marginal, may not be

undesirable if management objectives aim to reduce overall tree density, thereby promoting

stand structure more similar to that of pre-suppression eras.

2.4.2 Forest disturbances as drivers of tree population performance

We found insect outbreaks and wildfire to be the most important drivers of population

decline in subalpine tree species (Engelmann spruce, subalpine fir, and lodgepole pine),

although disease was also important for subalpine fir (Figure 2.6). This result is consistent

with previous findings that indicate elevated mortality in lodgepole pine, Engelmann spruce,

and subalpine fir are linked closely with outbreaks of mountain pine beetle (Dendroctonus

ponderosae) [Page and Jenkins, 2007], spruce beetle (Dendroctonus rufipennis) [DeRose and

Long, 2012], and subalpine fir decline (i.e., a complex disorder caused by multiple mortality

agents) [Reich et al., 2016], respectively. Furthermore, increased frequency and extent of

severe fire [Turner et al., 2019] and post-fire drought [Harvey et al., 2016] have been linked to

recruitment declines in subalpine forests across the region. Future shifts in forest disturbance

regimes [Buma et al., 2013, Turner et al., 2019], paired with expected declines in climate

suitability [Dobrowski et al., 2015, Bell et al., 2014b], suggest the future of subalpine tree

species in the western US may be grim.

Though the effects are marginal, we found that population performance of small diameter

grand fir and western hemlock exhibited positive response to drought severity (Figure 2.6),

supporting a narrative of dry-site encroachment by these species. That is, performance of

small diameter populations increases with increasing drought severity, potentially driven

by increased light and/or moisture availability arising from heightened mortality of large

diameter trees that are more susceptible to drought [Bennett et al., 2015, Stovall et al.,

2019]. We found other forest disturbances to be more important than drought in explaining

variation in population performance of grand fir and western hemlock. However the most

important drivers varied substantially between the species, potentially indicating the species

41

have fundamentally different exposures to forest disturbances.

Among montane and woodland species, we found wildfire and insect outbreaks to be

among the most important drivers of population performance, though disease exhibited

strong control on quaking aspen performance, as expected given widespread dieback of the

species due to sudden aspen decline [Worrall et al., 2010, Marchetti et al., 2011, Rehfeldt

et al., 2009]. As expected for important commercial timber species, we found that harvest and

landclearing was also an important driver for Douglas-fir and ponderosa pine. In general,

the importance of major forest disturbances was similar between subalpine and montane

species, however higher rates of decline among subalpine species may further highlight that

these species are more vulnerable to shifts in forest disturbance patterns.

2.4.3 Early indications of species range shifts

Future shifts in tree species distributions in response to changes in climate and forest distur-

bance regimes are widely expected [Davis and Shaw, 2001, Allen et al., 2010, Iverson et al.,

2008]. Previous efforts to detect such shifts have compared the distribution of tree life stages

at a single point in time (i.e., life stage for time substitution), but we present one of the first

to directly incorporate temporally replicated observations of species abundance across the

western US. We observed variation in population performance across the ranges of all species

(Figure 2.5), indicating that significant geographic shifts in species abundance distributions

have occurred over the last two decades. However, we found that the spatial pattern of

species abundance shifts were far more complex than the higher latitudinal–elevation finger-

print generally expected under climate change. The complexity of species abundance shifts

observed in this study are consistent with findings of previous works in the region [Bell et al.,

2014b, Serra-Diaz et al., 2016, 2015, Dobrowski et al., 2015], which attribute such complex-

ity to interactions among climate patterns, topography, and biotic drivers (i.e., competition,

succession).

While generally less important than forest disturbances, we found significant relationships

42

between population performance and long-term climate variables for some tree species (Fig-

ure 2.6), indicating that significant shifts in species abundance distributions have occurred

across climatic gradients. Most notably, we observed a positive association between tree pop-

ulation performance and long-term mean temperature for all species/ size class pairs where a

significant relationship between the variables existed, indicating a shift in abundance toward

historically warmer regions. This result counters the expectation of upslope migration ex-

pected under globally warming temperatures, and indicates that shifts in species abundances

may be driven instead by transient downslope migration mediated by biotic interactions, as

noted in Lenoir et al. [2009]. Relationships were more nuanced between population perfor-

mance of species/ size class pairs and long-term precipitation patterns. However, grand fir

population performance exhibited a strong, negative association with long-term precipita-

tion (i.e., population performance enhanced on historically dry sites), providing additional

support for the narrative of dry-site encroachment by the species.

2.4.4

Intra-specific divergence in population performance among size classes

Variation in population performance across tree size classes can have important ramifica-

tions for forest structure. That is, population expansion among small diameter stems of a

species may offset population decline among large diameter stems, resulting in similar forest

composition, but substantial shift in abundance towards small diameter stems over time

[Dolanc et al., 2013]. We observe significant divergence in range-wide population perfor-

mance between tree size classes for most species (80%), though the magnitude and direction

of divergence varied both by species (Table 2.2) and across space (Figure 2.5; bottom). The

largest divergence appeared in quaking aspen and lodgepole pine, where we observed severe

decline was evident among large diameter stems and rapid expansion among small diame-

ter stems. Both quaking aspen and lodgepole pine have experienced recently elevated rates

of mortality due to sudden aspen decline [Worrall et al., 2010, Marchetti et al., 2011] and

mountain pine beetle [Page and Jenkins, 2007], respectively, and it is likely that these agents

43

contributed to the decline of large diameter populations of each species. However, population

expansion observed in small diameter stems could be attributed to increased regeneration

rates in response to elevated mortality levels (i.e., increased resource availability) within

each species current distribution, an expansion of each species range beyond its current dis-

tribution (i.e., driven by increased regeneration rates at range margins), or a combination of

both. Further efforts to identify the drivers of such intra-specific divergence in population

performance across tree size classes are warranted.

2.5 Conclusions

In recent decades, forests of the western US have experienced unprecedented change in

climate and forest disturbance regimes, and widespread future shifts in forest composition,

structure, and function are expected in response [Williams et al., 2013, Dobrowski et al.,

2015, Liu et al., 2013]. However, large-scale, comprehensive assessments of tree population

performance have yet to be conducted in the region, and uncertainty remains regarding the

current status of tree populations (e.g., expanding, declining) and how multiple drivers may

interact to influence their stability. Our study provides empirical evidence to suggest the

most abundant tree species in the western US are exhibiting strong divergence in population

performance, with over half (70%) of species experiencing range-wide population decline.

We found spatial variation in population performance across the ranges of all species, indi-

cating range shifts are already underway and are pervasive among top species in the region.

Our results further indicate that species decline can seldom be attributed to a single forest

disturbance agent, highlighting the importance of considering multiple, simultaneous risks

factors when interpreting forest decline.

44

APPENDIX

45

APPENDIX

APPENDIX A

This appendix is intended to serve as further reference to section 2.3, Software Testing. We

have conducted extensive validation of rFIA for all available attributes against EVALIDator.

Here we present the code used to produce an abbreviated version of a validation for the state

of Connecticut in the year 2018 (Table 1, section 2.3). First we download the state subset of

the FIA Database for Connecticut and then subset the database to only include observations

necessary to compute estimates for the year 2018 (Figure ).

Figure A.1: Example code to download and process data used in software validation.

## Load rFIA
library(rFIA)

## Download data
ct <- getFIA('CT')

## 2018 Subset
ct <- clipFIA(ct, mostRecent = FALSE, evalid = findEVALID(ct, year = 2018))

A.0.1 Live tree abundance (trees/acre)

EVALIDator

• Land Basis: Forestland
• Numerator: Tree Number
• Denominator: Area
• Estimate: Number of live trees (at least 1 inch d.b.h./d.r.c.), in trees, on forest land
• Evaluation: 092018N CONNECTICUT 2012;2013;2014;2015;2016;2017;2018
• Row Variable: EVALID

46

• Column Variable: EVALID

rFIA

Figure A.2: Example code to estimate live tree abundance

# Trees per acre
tpa_ct <- tpa(ct)
tpa_ct$TPA

A.0.2 Live tree basal area (sq.ft./acre)

EVALIDator

• Land Basis: Forestland
• Numerator: Tree basal area
• Denominator: Area
• Estimate: Basal area of live trees (at least 1 inch d.b.h./d.r.c.), in square feet, on forest

land

• Evaluation: 092018N CONNECTICUT 2012;2013;2014;2015;2016;2017;2018
• Row Variable: EVALID
• Column Variable: EVALID

rFIA

Figure A.3: Example code to estimate live tree basal area

# Basal area per acre (BAA)
tpa_ct <- tpa(ct)
tpa_ct$BAA

47

A.0.3 Live tree merchantable volume (cu.ft./acre)

EVALIDator

• Land Basis: Forestland
• Numerator: Tree volume
• Denominator: Area
• Estimate: Net merchantable bole volume of live trees (at least 5 inch d.b.h./d.r.c.), in

cubic feet, on forest land

• Evaluation: 092018N CONNECTICUT 2012;2013;2014;2015;2016;2017;2018
• Row Variable: EVALID
• Column Variable: EVALID

rFIA

Figure A.4: Example code to estimate live tree merchantable volume.

# Net volume per acre
bio_ct <- biomass(ct)
bio_ct$NETVOL_ACRE

A.0.4 Live tree sawlog volume (cu.ft./acre)

EVALIDator

• Land Basis: Forestland
• Numerator: Tree volume
• Denominator: Area
• Estimate: Net sawlog volume of live trees, in cubic feet, on forest land
• Evaluation: 092018N CONNECTICUT 2012;2013;2014;2015;2016;2017;2018
• Row Variable: EVALID
• Column Variable: EVALID

48

rFIA

Figure A.5: Example code to estimate live tree saw log volume.

# Sawlog volume per acre
bio_ct <- biomass(ct)
bio_ct$SAWVOL_ACRE

A.0.5 Live tree aboveground biomass (tons/acre)

EVALIDator

• Land Basis: Forestland
• Numerator: Tree dry weight
• Denominator: Area
• Estimate: Aboveground biomass of live trees (at least 1 inch d.b.h./d.r.c.), in short

tons, on forest land

• Evaluation: 092018N CONNECTICUT 2012;2013;2014;2015;2016;2017;2018
• Row Variable: EVALID
• Column Variable: EVALID

rFIA

Figure A.6: Example code to estimate live tree aboveground biomass

# Sawlog volume per acre
bio_ct <- biomass(ct)
bio_ct$BIO_AG_ACRE

A.0.6 Live tree aboveground carbon (tons/acre)

EVALIDator

49

• Land Basis: Forestland
• Numerator: Tree carbon
• Denominator: Area
• Estimate: Aboveground carbon of live trees (at least 1 inch d.b.h./d.r.c.), in short

tons, on forest land

• Evaluation: 092018N CONNECTICUT 2012;2013;2014;2015;2016;2017;2018
• Row Variable: EVALID
• Column Variable: EVALID

rFIA

Figure A.7: Example code to estimate live tree carbon

# Sawlog volume per acre
bio_ct <- biomass(ct)
bio_ct$CARB_AG_ACRE

A.0.7 Annual net biomass growth (tons/acre/year)

EVALIDator

• Land Basis: Forestland
• Numerator: Annual net growth dry weight
• Denominator: Area
• Estimate: Average annual net growth of aboveground biomass of trees (at least 5 inch

d.b.h./d.r.c.), in short tons, on forest land

• Evaluation: 092018Y CONNECTICUT 2012;2013;2014;2015;2016;2017;2018
• Row Variable: EVALID
• Column Variable: EVALID

rFIA

50

Figure A.8: Example code to estimate live tree biomass growth.

# Annual biomass growth per acre
vr_ct <- vitalRates(ct)
vr_ct$BIO_GROW_AC

A.0.8 Annual mortality (trees/acre/year)

EVALIDator

• Land Basis: Forestland
• Numerator: Annual mortality number
• Denominator: Area
• Estimate: Average annual mortality of trees (at least 5 inch d.b.h./d.r.c.), in trees, on

forest land

• Evaluation: 092018Y CONNECTICUT 2012;2013;2014;2015;2016;2017;2018
• Row Variable: EVALID
• Column Variable: EVALID

rFIA

Figure A.9: Example code to estimate tree mortality.

# Annual mortality TPA
gm_ct <- growMort(ct)
gm_ct$MORT_TPA

A.0.9 Annual removals (trees/acre/year)

EVALIDator

• Land Basis: Forestland
• Numerator: Annual removal number

51

• Denominator: Area
• Estimate: Average annual removals of trees (at least 5 inch d.b.h./d.r.c.), in trees, on

forest land

• Evaluation: 092018Y CONNECTICUT 2012;2013;2014;2015;2016;2017;2018
• Row Variable: EVALID
• Column Variable: EVALID

rFIA

Figure A.10: Example code to estimate tree removals.

# Annual removals TPA
gm_ct <- growMort(ct)
gm_ct$REMV_TPA

A.0.10 Coarse woody debris volume (cu.ft./acre)

EVALIDator

• Land Basis: Forestland
• Numerator: Down woody material volume
• Denominator: Area
• Estimate: Total volume of CWD, in cubic feet, on forest land
• Evaluation: 092018N CONNECTICUT 2012;2013;2014;2015;2016;2017;2018
• Row Variable: EVALID
• Column Variable: EVALID

rFIA

A.0.11 Coarse woody debris biomass (tons/acre)

EVALIDator

52

Figure A.11: Example code to estimate coarse woody debris volume.

# CWD volume
dwm_ct <- dwm(ct)
dwm_ct$VOL_ACRE[dwm_ct$FUEL_TYPE == '1000HR']

• Land Basis: Forestland
• Numerator: Down woody material biomass
• Denominator: Area
• Estimate: Weight of CWD, in short tons, on forest land
• Evaluation: 092018N CONNECTICUT 2012;2013;2014;2015;2016;2017;2018
• Row Variable: EVALID
• Column Variable: EVALID

rFIA

Figure A.12: Example code to estimate coarse woody debris biomass.

# CWD biomass
dwm_ct <- dwm(ct)
dwm_ct$BIO_ACRE[dwm_ct$FUEL_TYPE == '1000HR']

A.0.12 Coarse woody debris carbon (tons/acre)

EVALIDator

• Land Basis: Forestland
• Numerator: Down woody material carbon
• Denominator: Area
• Estimate: Carbon of CWD, in short tons, on forest land
• Evaluation: 092018N CONNECTICUT 2012;2013;2014;2015;2016;2017;2018
• Row Variable: EVALID

53

• Column Variable: EVALID

rFIA

Figure A.13: Example code to estimate coarse woody debris carbon.

# CWD carbon
dwm_ct <- dwm(ct)
dwm_ct$CARB_ACRE[dwm_ct$FUEL_TYPE == '1000HR']

A.0.13 Total forest area (acres x 10−3)

EVALIDator

• Land Basis: Forestland
• Numerator: Area
• Denominator: None
• Estimate: Area in forestland, in acres
• Evaluation: 092018N CONNECTICUT 2012;2013;2014;2015;2016;2017;2018
• Row Variable: EVALID
• Column Variable: EVALID

rFIA

Figure A.14: Example code to estimate total forested area.

# Total Forest Area
fa_ct <- area(ct)
fa_ct$AREA_TOTAL

54

BIBLIOGRAPHY

55

BIBLIOGRAPHY

Henry D Adams and Thomas E Kolb. Drought responses of conifers in ecotone forests of

northern arizona: tree ring growth and leaf δ 13 c. Oecologia, 140(2):217–225, 2004.

Henry D Adams and Thomas E Kolb. Tree growth response to drought and temperature in
a mountain landscape in northern arizona, usa. Journal of Biogeography, 32(9):1629–1640,
2005.

Craig D Allen, Alison K Macalady, Haroun Chenchouni, Dominique Bachelet, Nate McDow-
ell, Michel Vennetier, Thomas Kitzberger, Andreas Rigling, David D Breshears, EH Ted
Hogg, et al. A global overview of drought and heat-induced tree mortality reveals emerging
climate change risks for forests. Forest ecology and management, 259(4):660–684, 2010.

Michael B Ashcroft. Identifying refugia from climate change. Journal of biogeography, 37

(8):1407–1413, 2010.

Ernst Assmann. The principles of forest yield study: studies in the organic production,

structure, increment and yield of forest stands. Elsevier, 2013.

Juliann E Aukema, Brian Leung, Kent Kovacs, Corey Chivers, Kerry O Britton, Jeffrey
Englin, Susan J Frankel, Robert G Haight, Thomas P Holmes, Andrew M Liebhold, et al.
Economic impacts of non-native forest insects in the continental united states. PLoS one,
6(9):e24587, 2011.

William A Bechtold, Paul L Patterson, et al. The enhanced forest inventory and analysis
program-national sampling design and estimation procedures. Gen. Tech. Rep. SRS-80.
Asheville, NC: US Department of Agriculture, Forest Service, Southern Research Station.
85 p., 80, 2005. doi: 10.2737/SRS-GTR-80.

Santiago Beguer´ıa, Sergio M Vicente-Serrano, Fergus Reig, and Borja Latorre. Standardized
precipitation evapotranspiration index (spei) revisited: parameter fitting, evapotranspira-
tion models, tools, datasets and drought monitoring. International Journal of Climatology,
34(10):3001–3023, 2014.

Santiago Beguer´ıa and Sergio M. Vicente-Serrano. SPEI: Calculation of the Standard-
ised Precipitation-Evapotranspiration Index, 2017. URL https://CRAN.R-project.org/
package=SPEI. R package version 1.7.

David M Bell, John B Bradford, and William K Lauenroth. Early indicators of change:
divergent climate envelopes between tree life stages imply range shifts in the western u
nited s tates. Global Ecology and Biogeography, 23(2):168–180, 2014a.

David M Bell, John B Bradford, and William K Lauenroth. Mountain landscapes offer
few opportunities for high-elevation tree species migration. Global change biology, 20(5):
1441–1451, 2014b.

56

Amy C Bennett, Nathan G McDowell, Craig D Allen, and Kristina J Anderson-Teixeira.
Larger trees suffer most during drought in forests worldwide. Nature Plants, 1(10):15139,
2015.

Barbara J Bentz, Jacques R´egni`ere, Christopher J Fettig, E Matthew Hansen, Jane L Hayes,
Jeffrey A Hicke, Rick G Kelsey, Jose F Negr´on, and Steven J Seybold. Climate change
and bark beetles of the western united states and canada: direct and indirect effects.
BioScience, 60(8):602–613, 2010.

Romain Bertrand, Jean-Claude G´egout, and Jean-Daniel Bontemps. Niches of temperate
tree species converge towards nutrient-richer conditions over ontogeny. Oikos, 120(10):
1479–1488, 2011.

Christof Bigler, Daniel G Gavin, Charles Gunning, and Thomas T Veblen. Drought induces
lagged tree mortality in a subalpine forest in the rocky mountains. Oikos, 116(12):1983–
1994, 2007.

David D Breshears, Neil S Cobb, Paul M Rich, Kevin P Price, Craig D Allen, Randy G
Balice, William H Romme, Jude H Kastens, M Lisa Floyd, Jayne Belnap, et al. Regional
vegetation die-off in response to global-change-type drought. Proceedings of the National
Academy of Sciences, 102(42):15144–15148, 2005.

Brian Buma, Carissa D Brown, Dan C Donato, Joseph B Fontaine, and Jill F Johnstone. The
impacts of changing disturbance regimes on serotinous plant populations and communities.
BioScience, 63(11):866–876, 2013.

Russell M Burns. Silvics of North America: Conifers. Number 654. US Department of

Agriculture, Forest Service, 1990.

Elizabeth A Burrill, Andrea M Wilson, Jeffery A Turner, Scott A Pugh, James Menlove,
Glenn Christiansen, Barbara L Conkling, and Winnie David. The forest inventory and
analysis database: database description and user guide version 8.0 for phase 2. Gen. Tech.
Rep. RMRS-GTR-245. Fort Collins, CO: US Department of Agriculture, Forest Service,
Rocky Mountain Research Station. 946 p., 245, 2018.

Michael J Case and David L Peterson. Fine-scale variability in growth climate relationships
of douglas-fir, north cascade range, washington. Canadian Journal of Forest Research, 35
(11):2743–2755, 2005.

James S Clark, David M Bell, Matthew C Kwit, and Kai Zhu. Competition-interaction
landscapes for the joint response of forests to climate change. Global Change Biology, 20
(6):1979–1991, 2014.

James S Clark, Louis Iverson, Christopher W Woodall, Craig D Allen, David M Bell, Don C
Bragg, Anthony W D’Amato, Frank W Davis, Michelle H Hersh, Ines Ibanez, et al. The
impacts of increasing drought on forest dynamics, structure, and biodiversity in the united
states. Global change biology, 22(7):2329–2352, 2016.

57

Erin Conlisk, Cristina Castanha, Matthew J Germino, Thomas T Veblen, Jeremy M Smith,
and Lara M Kueppers. Declines in low-elevation subalpine tree populations outpace growth
in high-elevation populations with warming. Journal of Ecology, 105(5):1347–1357, 2017.

Joseph H Connell, JG Tracey, and Leonard J Webb. Compensatory recruitment, growth,
and mortality as factors maintaining rain forest tree diversity. Ecological monographs, 54
(2):141–164, 1984.

Edward R Cook, Connie A Woodhouse, C Mark Eakin, David M Meko, and David W Stahle.
Long-term aridity changes in the western united states. Science, 306(5698):1015–1018,
2004.

Kimberley T Davis, Solomon Z Dobrowski, Philip E Higuera, Zachary A Holden, Thomas T
Veblen, Monica T Rother, Sean A Parks, Anna Sala, and Marco P Maneta. Wildfires
and climate change push low-elevation forests across a critical climate threshold for tree
regeneration. Proceedings of the National Academy of Sciences, 116(13):6193–6198, 2019.

Margaret B Davis and Ruth G Shaw. Range shifts and adaptive responses to quaternary

climate change. Science, 292(5517):673–679, 2001.

R Justin DeRose and James N Long. Factors influencing the spatial and temporal dynamics
of engelmann spruce mortality during a spruce beetle outbreak on the markagunt plateau,
utah. Forest Science, 58(1):1–14, 2012.

Solomon Z Dobrowski, Alan K Swanson, John T Abatzoglou, Zachary A Holden, Hugh D
Safford, Mike K Schwartz, and Daniel G Gavin. Forest structure and species traits mediate
projected recruitment declines in western us tree species. Global Ecology and Biogeography,
24(8):917–927, 2015.

Christopher R Dolanc, James H Thorne, and Hugh D Safford. Widespread shifts in the
demographic structure of subalpine forests in the sierra nevada, california, 1934 to 2007.
Global Ecology and Biogeography, 22(3):264–276, 2013.

Francis H Eyre. Forest cover types. Washington, DC: Society of American Foresters, 1980.

Songlin Fei, Johanna M Desprez, Kevin M Potter, Insu Jo, Jonathan A Knott, and Christo-
pher M Oswalt. Divergence of species responses to climate change. Science Advances, 3
(5):e1603055, 2017.

Jerry F Franklin, Thomas A Spies, Robert Van Pelt, Andrew B Carey, Dale A Thornburgh,
Dean Rae Berg, David B Lindenmayer, Mark E Harmon, William S Keeton, David C
Shaw, et al. Disturbances and structural development of natural forest ecosystems with
silvicultural implications, using douglas-fir forests as an example. Forest Ecology and
Management, 155(1-3):399–423, 2002.

Tracey S Frescino, Paul L Patterson, Gretchen G Moisen, and Elizabeth A Freeman. Fi-
esta—an r estimation tool for fia analysts, 2015. URL https://www.fs.usda.gov/
treesearch/pubs/50184.

58

Alisa L Gallant, Andrew J Hansen, John S Councilman, Duane K Monte, and David W Betz.
Vegetation dynamics under fire exclusion and logging in a rocky mountain watershed,
1856–1996. Ecological Applications, 13(2):385–403, 2003.

Andrew JR Gillespie. Rationale for a national annual forest inventory program. Journal of

Forestry, 97(12):16–20, 1999. doi: 10.1093/jof/97.12.16.

Sara A Goeking. Disentangling forest change from forest inventory change: A case study

from the us interior west. Journal of Forestry, 113(5):475–483, 2015.

Andrew N Gray, Thomas J Brandeis, John D Shaw, William H McWilliams, and Patrick
Miles. Forest inventory and analysis database of the united states of america (fia). In:
Dengler, J.; Oldeland, J.; Jansen, F.; Chytry, M.; Ewald, J., Finckh, M.; Glockler, F.;
Lopez-Gonzalez, G.; Peet, RK; Schaminee, J. HJ, eds. Vegetation databases for the 21st
century. Biodiversity and Ecology. 4: 225-231., pages 225–231, 2012. doi: 10.7809/b-e.
00079.

Sarah Greenwood, Paloma Ruiz-Benito, Jordi Mart´ınez-Vilalta, Francisco Lloret, Thomas
Kitzberger, Craig D Allen, Rod Fensham, Daniel C Laughlin, Jens Kattge, Gerhard
B¨onisch, et al. Tree mortality across biomes is promoted by drought intensity, lower
wood density and higher specific leaf area. Ecology Letters, 20(4):539–553, 2017.

Peter J Grubb. The maintenance of species-richness in plant communities: the importance

of the regeneration niche. Biological reviews, 52(1):107–145, 1977.

Brian J Harvey, Daniel C Donato, and Monica G Turner. 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, 2016.

Paul F Hessburg and James K Agee. An environmental narrative of inland northwest united

states forests, 1800–2000. Forest Ecology and Management, 178(1-2):23–59, 2003.

Paul F Hessburg, James K Agee, and Jerry F Franklin. Dry forests and wildland fires of the
inland northwest usa: contrasting the landscape ecology of the pre-settlement and modern
eras. Forest Ecology and Management, 211(1-2):117–139, 2005.

Christ W. Hoffman, Robert L. Mathiasen, Richard W. Hofstetter, Mary Lou Fairweather,
John D. Shaw, John W. Hanna, and Ned B. Klopfenstein. Survey for Armillaria by Plant
Associations in Northern Arizona. Journal of the Arizona-Nevada Academy of Science, 45
(2):76 – 86, 2014. doi: 10.2181/036.045.0204.

Daniel G Horvitz and Donovan J Thompson. A generalization of sampling without replace-
ment from a finite universe. Journal of the American statistical Association, 47(260):
663–685, 1952.

Kaicheng Huang, Chuixiang Yi, Donghai Wu, Tao Zhou, Xiang Zhao, William J Blanford,
Suhua Wei, Hao Wu, Du Ling, and Zheng Li. Tipping point of a conifer forest ecosystem
under severe drought. Environmental Research Letters, 10(2):024011, 2015.

59

Louis R Iverson, Anantha M Prasad, Stephen N Matthews, and Matthew Peters. Estimating
potential habitat for 134 eastern us tree species under six climate scenarios. Forest ecology
and management, 254(3):390–406, 2008.

Stephen T Jackson, Julio L Betancourt, Robert K Booth, and Stephen T Gray. Ecology
and the ratchet of events: climate variability, niche dimensions, and species distributions.
Proceedings of the National Academy of Sciences, 106(Supplement 2):19685–19692, 2009.

Peng Jiang, Hongyan Liu, Shilong Piao, Philippe Ciais, Xiuchen Wu, Yi Yin, and Hongya
Wang. Enhanced growth after extreme wetness compensates for post-drought carbon loss
in dry forests. Nature communications, 10(1):1–9, 2019.

Jill F Johnstone, Craig D Allen, Jerry F Franklin, Lee E Frelich, Brian J Harvey, Philip E
Higuera, Michelle C Mack, Ross K Meentemeyer, Margaret R Metz, George LW Perry,
et al. Changing disturbance regimes, ecological memory, and forest resilience. Frontiers
in Ecology and the Environment, 14(7):369–378, 2016.

Paul A Klockow, Jason G Vogel, Christopher B Edgar, and Georgianne W Moore. Lagged
mortality among tree species four years after an exceptional drought in east texas. Eco-
sphere, 9(10):e02455, 2018.

Kathryn W Kromroy, Jennifer Juzwik, Paul Castillo, and Mark H Hansen. Using forest ser-
vice forest inventory and analysis data to estimate regional oak decline and oak mortality.
Northern journal of applied forestry, 25(1):17–24, 2008.

Jonathan Lenoir, Jean-Claude G´egout, Jean-Claude Pierrat, Jean-Daniel Bontemps, and
Jean-Fran¸cois Dhˆote. Differences between tree species seedling and adult altitudinal dis-
tribution in mountain forests during the recent warm period (1986–2006). Ecography, 32
(5):765–777, 2009.

Heather E Lintz, Andrew N Gray, Andrew Yost, Richard Sniezko, Chris Woodall, Matt
Reilly, Karen Hutten, and Mark Elliott. Quantifying density-independent mortality of
temperate tree species. Ecological indicators, 66:1–9, 2016.

Yongqiang Liu, Scott L Goodrick, and John A Stanturf. Future us wildfire potential trends
projected using a dynamically downscaled climate change scenario. Forest Ecology and
Management, 294:120–135, 2013.

Zhihai Ma, Changhui Peng, Qiuan Zhu, Huai Chen, Guirui Yu, Weizhong Li, Xiaolu Zhou,
Weifeng Wang, and Wenhua Zhang. Regional drought-induced reduction in the biomass
carbon sink of canada’s boreal forests. Proceedings of the National Academy of Sciences,
109(7):2423–2427, 2012.

Frantiˇsek M´aliˇs, Martin Kopeck`y, Petr Petˇr´ık, Jozef Vladoviˇc, J´an Merganiˇc, and Tom´aˇs
Vida. Life stage, not climate change, explains observed tree range shifts. Global change
biology, 22(5):1904–1914, 2016.

60

Sparkle L Malone, Anna W Schoettle, and Jonathan D Coop. The future of subalpine forests
in the southern rocky mountains: Trajectories for pinus aristata genetic lineages. PloS
one, 13(3), 2018.

Suzanne Bethers Marchetti, James J Worrall, and Thomas Eager. Secondary insects and
diseases contribute to sudden aspen decline in southwestern colorado, usa. Canadian
Journal of Forest Research, 41(12):2315–2325, 2011.

Ronald E. McRoberts and Patrick D. Miles. United States of America, pages 829–842.
Springer International Publishing, Cham, 2016. ISBN 978-3-319-44015-6. doi: 10.1007/
978-3-319-44015-6 45.

Ronald E McRoberts, William A Bechtold, Paul L Patterson, Charles T Scott, and Gregory A
Reams. The enhanced forest inventory and analysis program of the usda forest service:
Historical perspective and announcement of statistical documentation. Journal of Forestry,
103(6):304–308, 2005a.

Ronald E. McRoberts, Geoffrey R. Holden, Mark D. Nelson, Greg C. Liknes, Warren K.
Moser, Andrew J. Lister, Susan L. King, Elizabeth B. LaPoint, John W. Coulston, W. Brad
Smith, and Gregory A. Reams. Estimating and Circumventing the Effects of Perturbing
and Swapping Inventory Plot Locations. Journal of Forestry, 103(6):275–279, 09 2005b.
ISSN 0022-1201. doi: 10.1093/jof/103.6.275. URL https://doi.org/10.1093/jof/103.
6.275.

Constance I Millar, David A Charlet, Robert D Westfall, John C King, Diane L Delany,
Alan L Flint, and Lorraine E Flint. Do low-elevation ravines provide climate refugia for
subalpine limber pine (pinus flexilis) in the great basin, usa? Canadian Journal of Forest
Research, 48(6):663–671, 2018.

Maria N Miriti. Ontogenetic shift from facilitation to competition in a desert shrub. Journal

of Ecology, 94(5):973–979, 2006.

Vicente J Monleon and Heather E Lintz. Evidence of tree species’ range shifts in a complex

landscape. PLoS One, 10(1), 2015.

Georgianne W Moore, Christopher B Edgar, Jason G Vogel, Robert A Washington-Allen,
Rosaleen G March, and Rebekah Zehnder. Tree mortality from an exceptional drought
spanning mesic to semiarid ecoregions. Ecological Applications, 26(2):602–611, 2016.

Cameron Naficy, Anna Sala, Eric G Keeling, Jon Graham, and Thomas H DeLuca. Inter-
active effects of historical logging and fire exclusion on ponderosa pine forest structure in
the northern rockies. Ecological Applications, 20(7):1851–1864, 2010.

Stephen W Pacala, Charles D Canham, John Saponara, John A Silander Jr, Richard K
Kobe, and Eric Ribbens. Forest models defined by field measurements: estimation, error
analysis and dynamics. Ecological monographs, 66(1):1–43, 1996.

61

Wesley Green Page and Michael James Jenkins. Mountain pine beetle-induced changes to
selected lodgepole pine fuel complexes within the intermountain region. Forest Science,
53(4):507–518, 2007.

Zuzana Parobekov´a, J´an Pittner, Stanislav Kucbel, Milan Saniga, Michal Fil´ıpek, Denisa
Sedm´akov´a, Jaroslav Vencurik, and Peter Jaloviar. Structural diversity in a mixed spruce-
fir-beech old-growth forest remnant of the western carpathians. Forests, 9(7):379, 2018.

Richard G Pearson and Terence P Dawson. Predicting the impacts of climate change on
the distribution of species: are bioclimate envelope models useful? Global ecology and
biogeography, 12(5):361–371, 2003.

Edzer Pebesma. Simple Features for R: Standardized Support for Spatial Vector Data. The
R Journal, 10(1):439–446, 2018. doi: 10.32614/RJ-2018-009. URL https://doi.org/10.
32614/RJ-2018-009.

Robert K Peet. Forest vegetation of the colorado front range. Vegetatio, 45(1):3–75, 1981.

Therese M Poland and Deborah G McCullough. Emerald ash borer: invasion of the urban
forest and the threat to north america’s ash resource. Journal of Forestry, 104(3):118–124,
2006.

PRISM Climate Group. Prism climate group, 2010. URL http://www.prism.oregonstate.

edu/.

Scott A. Pugh, Jeffery A. Turner, Elizabeth A. Burrill, and Winnie David. The For-
est Inventory and Analysis Database: Population Estimation User Guide. USDA For-
est Service, Washington DC, USA, 2018. URL https://www.fia.fs.fed.us/library/
database-documentation/current/ver80/FIADB.

Camille Py, John Bauer, Peter J Weisberg, and Franco Biondi. Radial growth responses of

singleleaf pinyon (pinus monophylla) to wildfire. Dendrochronologia, 24(1):39–46, 2006.

R Core Team. R: A Language and Environment for Statistical Computing. R Foundation
for Statistical Computing, Vienna, Austria, 2018a. URL https://www.R-project.org/.

R Core Team. R: A Language and Environment for Statistical Computing. R Foundation
for Statistical Computing, Vienna, Austria, 2018b. URL https://www.R-project.org/.

Gerald E Rehfeldt, Dennis E Ferguson, and Nicholas L Crookston. Aspen, climate, and
sudden decline in western usa. Forest Ecology and Management, 258(11):2353–2364, 2009.

Robin M Reich, John E Lundquist, and Kristina Hughes. Host-environment mismatches
associated with subalpine fir decline in colorado. Journal of forestry research, 27(5):1177–
1189, 2016.

Alexia A. Sabor, Volker C. Radeloff, Ronald E. McRoberts, Murray Clayton, and Susan I.
Stewart. Adding uncertainty to forest inventory plot locations: effects on analyses using
geospatial data. Canadian Journal of Forest Research, 37(11):2313–2325, 2007. doi: 10.
1139/X07-067. URL https://doi.org/10.1139/X07-067.

62

Josep M Serra-Diaz, Robert M Scheller, Alexandra D Syphard, and Janet Franklin. Distur-
bance and climate microrefugia mediate tree range shifts during climate change. Landscape
Ecology, 30(6):1039–1053, 2015.

Josep M Serra-Diaz, Janet Franklin, Whalen W Dillon, Alexandra D Syphard, Frank W
Davis, and Ross K Meentemeyer. California forests show early indications of both range
shifts and local persistence under climate change. Global Ecology and Biogeography, 25(2):
164–175, 2016.

W Brad Smith. Forest inventory and analysis: a national inventory and monitoring program.

Environmental pollution, 116:S233–S242, 2002.

Hunter Stanke, Andrew O Finley, Aaron S Weed, Brian F Walters, and Grant M Domke.
rfia: An r package for estimation of forest attributes with the us forest inventory and
analysis database. Environmental Modelling & Software, page 104664, 2020.

Camille S Stevens-Rumann, Kerry B Kemp, Philip E Higuera, Brian J Harvey, Monica T
Rother, Daniel C Donato, Penelope Morgan, and Thomas T Veblen. Evidence for declining
forest resilience to wildfires under climate change. Ecology letters, 21(2):243–252, 2018.

Atticus EL Stovall, Herman Shugart, and Xi Yang. Tree height explains mortality risk during

an intense drought. Nature communications, 10(1):1–6, 2019.

Charles Warren Thornthwaite. An approach toward a rational classification of climate.

Geographical review, 38(1):55–94, 1948.

Luke Tierney, A. J. Rossini, Na Li, and H. Sevcikova. snow: Simple Network of Workstations,

2018. URL https://CRAN.R-project.org/package=snow. R package version 0.4-3.

Wade T Tinkham, Patrick R Mahoney, Andrew T Hudak, Grant M Domke, Mike J
Falkowski, Chris W Woodall, and Alistair MS Smith. Applications of the united states
forest inventory and analysis dataset: a review and future directions. Canadian Journal
of Forest Research, 48(11):1251–1268, 2018.

Erkki Tomppo, Thomas Gschwantner, Mark Lawrence, Ronald E McRoberts, K Gabler,
K Schadauer, C Vidal, A Lanz, G St˚ahl, and E Cienciala. National forest inventories.
Pathways for Common Reporting. European Science Foundation, pages 541–553, 2010.

Leo T¨ornqvist, Pentti Vartia, and Yrj¨o O Vartia. How should relative changes be measured?

The American Statistician, 39(1):43–46, 1985.

Monica G Turner, Kristin H Braziunas, Winslow D Hansen, and Brian J Harvey. Short-
interval severe fire erodes the resilience of subalpine lodgepole pine forests. Proceedings of
the National Academy of Sciences, 116(23):11319–11328, 2019.

USDA Forest Service. Forest inventory and analysis national program, 2019a. URL https:

//www.fia.fs.fed.us/.

USDA Forest Service. Fia tools and data, 2019b. URL https://www.fia.fs.fed.us/

tools-data/.

63

Phillip J Van Mantgem, Nathan L Stephenson, John C Byrne, Lori D Daniels, Jerry F
Franklin, Peter Z Ful´e, Mark E Harmon, Andrew J Larson, Jeremy M Smith, Alan H
Taylor, et al. Widespread increase of tree mortality rates in the western united states.
Science, 323(5913):521–524, 2009.

Marco Vanoni, Harald Bugmann, Magdalena N¨otzli, and Christof Bigler. Quantifying the
effects of drought on abrupt growth decreases of major tree species in switzerland. Ecology
and evolution, 6(11):3555–3570, 2016.

Sergio M Vicente-Serrano, Santiago Beguer´ıa, Juan I L´opez-Moreno, M Angulo, and A El Ke-
nawy. A new global 0.5 gridded dataset (1901–2006) of a multiscalar drought index: com-
parison with current drought index datasets based on the palmer drought severity index.
Journal of Hydrometeorology, 11(4):1033–1043, 2010.

Sergio M Vicente-Serrano, C´elia Gouveia, Jes´us Julio Camarero, Santiago Beguer´ıa, Ricardo
Trigo, Juan I L´opez-Moreno, C´esar Azor´ın-Molina, Edmond Pasho, Jorge Lorenzo-Lacruz,
Jes´us Revuelto, et al. Response of vegetation to drought time-scales across global land
biomes. Proceedings of the National Academy of Sciences, 110(1):52–57, 2013.

Hadley Wickham, Romain Francois, Lionel Henry, Kirill M¨uller, et al. dplyr: A grammar of

data manipulation. R package version 0.4, 3, 2015.

Hadley Wickham, Jim Hester, and Winston Chang. devtools: Tools to Make Developing
R Packages Easier, 2019. URL https://CRAN.R-project.org/package=devtools. R
package version 2.1.0.

A Park Williams, Craig D Allen, Alison K Macalady, Daniel Griffin, Connie A Woodhouse,
David M Meko, Thomas W Swetnam, Sara A Rauscher, Richard Seager, Henri D Grissino-
Mayer, et al. Temperature as a potent driver of regional forest drought stress and tree
mortality. Nature climate change, 3(3):292–297, 2013.

Carmen M Wong and Lori D Daniels. Novel forest decline triggered by multiple interactions
among climate, an introduced pathogen and bark beetles. Global change biology, 23(5):
1926–1941, 2017.

CW Woodall, K Zhu, JA Westfall, CM Oswalt, AW D’amato, BF Walters, and HE Lintz.
Assessing the stability of tree ranges and influence of disturbance in eastern us forests.
Forest Ecology and Management, 291:172–180, 2013.

James J Worrall, Suzanne B Marchetti, Leanne Egeland, Roy A Mask, Thomas Eager, and
Brian Howell. Effects and etiology of sudden aspen decline in southwestern colorado, usa.
Forest Ecology and Management, 260(5):638–648, 2010.

Ian J Wright, Peter B Reich, Mark Westoby, David D Ackerly, Zdravko Baruch, Frans
Bongers, Jeannine Cavender-Bares, Terry Chapin, Johannes HC Cornelissen, Matthias
Diemer, et al. The worldwide leaf economics spectrum. Nature, 428(6985):821–827, 2004.

64

Xiuchen Wu, Hongyan Liu, Xiaoyan Li, Philippe Ciais, Flurin Babst, Weichao Guo, Cicheng
Zhang, Vincenzo Magliulo, Marian Pavelka, Shaomin Liu, et al. Differentiating drought
legacy effects on vegetation growth over the temperate northern hemisphere. Global change
biology, 24(1):504–516, 2018.

Kai Zhu, Christopher W Woodall, and James S Clark. Failure to migrate: lack of tree range
expansion in response to climate change. Global Change Biology, 18(3):1042–1052, 2012.

Kai Zhu, Christopher W Woodall, Souparno Ghosh, Alan E Gelfand, and James S Clark.
Dual impacts of climate change: forest migration and turnover through life history. Global
change biology, 20(1):251–264, 2014.

65