Twenty years of drought‐mediated change in snag populations in mixed‐conifer and ponderosa pine forests in Northern Arizona

Abstract

Background

Snags (standing dead trees) are important biological legacies in forest systems, providing numerous resources as well as a record of recent tree mortality. From 1997 to 2017, we monitored snag populations in drought-influenced mixed-conifer and ponderosa pine (Pinus ponderosa) forests in northern Arizona.

Results

Snag density increased significantly in both forest types. This increase was driven largely by a pulse in snag recruitment that occurred between 2002 and 2007, following an extreme drought year in 2002, with snag recruitment returning to pre-pulse levels in subsequent time periods. Some later years during the study also were warmer and/or drier than average, but these years were not as extreme as 2002 and did not trigger the same level of snag recruitment. Snag recruitment was not equal across tree species and size classes, resulting in significant changes in species composition and size-class distributions of snag populations in both forest types. Because trees were far more abundant than snags in these forests, the effect of this mortality pulse on tree populations was far smaller than its effect on snag populations. Snag loss rates increased over time during the study, even though many snags were newly recruited. This may reflect the increasing prevalence of white fir snags and/or snags in the smaller size classes, which generally decay faster than snags of other species or larger snags. Thus, although total numbers of snags increased, many of the newly recruited snags may not persist long enough to be valuable as nesting substrates for native wildlife.

Conclusions

Increases in snag abundance appeared to be due to a short-term tree mortality “event” rather than a longer-term pattern of elevated tree mortality. This mortality event followed a dry and extremely warm year (2002) embedded within a longer-term megadrought. Climate models suggest that years like 2002 may occur with increasing frequency in the southwestern U.S. Such years may result in additional mortality pulses, which in turn may strongly affect trajectories in abundance, structure, and composition of snag populations. Relative effects on tree populations likely will be smaller, but, over time, also could be significant.

Background

Understanding the effects of climate change on structure and composition of forest ecosystems presents a major challenge for researchers and managers in the twenty-first century (Millar et al. 2007). Warming climates already have profoundly affected forests throughout the world (Adams et al. 2009; Allen et al. 2010), and climate models predict further warming and drying in many areas (e.g., Seager et al. 2007; Stocker et al. 2013; Garfin et al. 2014) suggesting that changing climates will continue to affect these systems. Forest systems can be affected directly through climate-mediated mortality (Breshears et al. 2005; Mueller et al. 2005; Adams et al. 2009), or indirectly through altered disturbance regimes that result in increased mortality (e.g., McKenzie et al. 2004; Bentz et al. 2010; Wan et al. 2019) or reduced regeneration (Davis et al. 2019). Because drought-mediated mortality may differentially affect tree species and sizes, such mortality may drive significant changes in forest structure and composition (Mueller et al. 2005; Ganey and Vojta 2011; Kane et al. 2014).

Snags, or standing dead trees, are important components of forest systems that are directly affected by drought-mediated trends in tree mortality. These snags serve as biological legacies in forest systems (Thomas et al. 1979; McComb and Lindenmayer 1999; Woldendorp and Keenan 2005), providing habitat for native wildlife (e.g., Bull et al. 1997; Rabe et al. 1998), serving as an important source of coarse woody debris (Harmon et al. 1986; Laudenslayer et al. 2002; Woldendorp and Keenan 2005), and aiding in nutrient cycling and other ecosystem functions. Snags also provide an index of recent tree mortality (Ganey and Vojta 2011; Wu et al. 2017), with recruitment of new snags reflecting temporal changes in tree mortality and relative mortality among species and size classes, potentially providing insight into changes in forest structure and composition.

In the southwestern United States, a pronounced increase in snag recruitment (i.e., tree mortality) followed an extreme climatic year (2002) embedded within a longer-term mega-drought (Breshears et al. 2005; Kopeke et al. 2010; Williams et al. 2020). Several studies reported on short-term trends in tree mortality (which results in snag recruitment) in important forest types following this extreme year (Negron et al. 2009 [through 2004]; Ganey and Vojta 2011 [through 2007]; Kane et al. 2014 [through 2008]), but data documenting longer-term trends in these forests are lacking, and few studies focused specifically on snag populations (but see Ganey and Vojta 2005, 2012, 2014).

Snag populations are governed by the balance between gains and losses in snags. Gains occur when new snags are created by natural tree senescence processes or disturbances such as insect outbreaks, diseases, fire, or droughts. Losses occur when snags are lost to timber or fuelwood harvest, fire, or natural decomposition, processes which are influenced by factors such as snag species and size (Ganey et al. 2015). It is therefore desirable not only to document populations of existing snags but also to understand how those populations change over time and the factors responsible for those changes.

Since 1997, we have sampled snag and tree populations periodically in mixed-conifer and ponderosa pine forests in northern Arizona. Previous papers from this study documented changes in snag populations over 5-, 10-, and 15-year increments within this 20-year period (Ganey and Vojta 2005, 2012, 2014) as well as patterns in: tree mortality from 2002 to 2007 (Ganey and Vojta 2011) and snag longevity from 1997 to 2015 (Ganey et al. 2015). Here, we expanded on this earlier work to summarize drought-mediated changes in snag and tree abundance, and composition of snag and tree populations within these forest types, over the 20-year period from 1997 to 2017, with an emphasis on temporal changes in the snag population, including rates of snag recruitment and loss and the factors driving these processes. Our specific objectives included: (1) Evaluating temporal trends in snag and tree density and snag recruitment and loss rates, (2) Summarizing climate during the study period and evaluating potential relationships between climate patterns and snag recruitment, and (3) Evaluating temporal trends in composition and structure of overall snag populations, as well as in snags that were recruited or lost during the intervals between sampling occasions. These data thus provide information on trends in snag populations during the study period as well as on the factors driving those trends. This information should aid forest managers in understanding the potential effects of future climate patterns on these snag populations, and, to a lesser extent, on the tree populations from which they derive.

Methods

Study area

We established plots for sampling snags and trees within a study area of approximately 73,000 ha within the Coconino and Kaibab National Forests, north-central Arizona. Study plots were randomly located in mixed-conifer (n = 53 plots) and ponderosa pine forests (n = 60 plots) within this area (see Ganey 1999) for details on plot selection). White fir (Abies concolor Lindl. ex Hildebr.), Douglas-fir (Pseudotsuga menziesii [Mirb.] Franco), and ponderosa pine dominated overstories in mixed-conifer forests, accounting for approximately 90 % of total trees in this forest type (as sampled in 2004; Ganey and Vojta 2011). Other relatively common species included Gambel oak (Quercus gambelii Nutt.), quaking aspen (Populus tremuloides Michaux), and limber pine (P. flexilis James), in that order of frequency. Ponderosa pine comprised > 90 % of trees in ponderosa pine forest in 2004 (Ganey and Vojta 2011), with Gambel oak accounting for approximately 8 % of total trees by frequency. Alligator juniper (Juniperus deppeana Steud), Douglas-fir, quaking aspen, limber pine, pinyon pine (P. edulis), and other species of juniper were present in some stands, typically in relatively small numbers.

Study plots included a wide range of topographic conditions and soil types, ranged in elevation from the transition zone between pinyon − juniper woodland and ponderosa pine at lower elevations to the ecotone between mixed-conifer and Engelmann spruce (Picea engelmanni Parry ex Engelm.) − corkbark fir (Abies lasiocarpa var. arizonica [Merriam] Lemmon) forests at higher elevations (Brown et al. 1980), and included both commercial forest lands and administratively reserved lands such as wilderness and other roadless areas. Plots also included areas subject to forest thinning and both prescribed burns and wildfires. Consequently, these plots sampled a wide range of conditions with respect to forest structure and composition, and included areas subjected to recent disturbance events.

Sampling snag and tree populations

We sampled snags in 1-ha plots at five-year intervals beginning in 1997 (i.e., snags were sampled in 1997, 2002, 2007, 2012, and 2017). Within plot boundaries, we sampled all snags ≥ 2 m in height and ≥ 20 cm in diameter at breast height (dbh). The focus on snags ≥ 20 cm in dbh reflected the original study objectives related to wildlife habitat (Ganey 1999) and the assumption that smaller snags were less important to wildlife such as cavity-nesting birds (Scott 1978; Cunningham et al. 1980; Conway and Martin 1993) and/or roosting bats (Rabe et al. 1998; Bernardos et al. 2004; Solvesky and Chambers 2009). Given this minimum diameter, our inference is limited to populations of snags ≥ 20 cm dbh.

We marked all snags with numbered metal tags, allowing us to distinguish pre-existing snags from new snags when re-sampling plots (with some exceptions, see below). We recorded species and dbh (nearest cm) for all snags. We sampled plots from May through August, and did not necessarily sample individual plots on the same date or even in the same month across years. Thus, the elapsed time between consecutive samples for an individual plot could range from slightly < 5 years to slightly > 5 years. We ignored this variability in analysis, and assumed that all intervals between sampling occasions represented a 5-year period.

We sampled live trees ≥ 20 cm in dbh, which were far more abundant than snags, in a 0.09-ha subplot located within each snag plot in 2004 and 2014. We adhered to the minimum 20-cm diameter for consistency with snag sampling. We recorded tree species and dbh (nearest cm) for all trees. We did not mark individual trees, and consequently were not able to determine fates of individual trees. Instead, we focused on overall changes in tree populations, which were driven by the interactions among ingrowth of small (< 20 cm dbh in 2004) trees into our sampled population, growth of trees within that sampled population, and tree mortality. As with snags, our inference is limited to trees ≥ 20 cm dbh.

Analysis of snag and tree populations

We included all sampled plots in our analyses of snag and tree populations, including plots subject to recent disturbances such as wild or prescribed fire, for two reasons. First, our objective was to summarize trends in snag and tree populations on the overall landscape, including areas subject to these disturbances. Second, although the effect of fire could be large on individual plots, the overall effect was relatively small over most 5-year intervals between sampling occasions, because relatively few plots were impacted by moderate-to high severity fire during those intervals (Table 1). Consequently, estimates for all plots were similar to estimates including only plots that did not experience moderate-to high-severity fire. We also included plots subject to forest thinning, but these were even fewer in number than burned plots, and thinning activities therefore had minimal immediate impacts on overall snag or tree populations (although they may impact future patterns of tree mortality and/or wildfire).

Table 1 Number (and percent) of plots in northern Arizona mixed-conifer and ponderosa pine forest that experienced moderate-to high-severity fire by 5-year interval between snag sampling occasions. n = 53 and 60 plots sampled in mixed-conifer and ponderosa pine forest, respectively. Plots were classified as having experienced moderate- to high-severity fire if a visual inspection indicated that considerable tree mortality had occurred within the plot due to fire between sampling occasions

We treated years in which we sampled snag (or tree) populations as sampling occasions and the intervals between those sampling occasions as “sampling periods”. Thus, for snags we had five sampling occasions (hereafter “years”) and four 5-year sampling periods, whereas for trees we had two years and one 10-year sampling period. Years provided estimates of density and composition for snag (or tree) populations at a point in time, and sampling periods provided estimates of change in those parameters between those points in time.

As noted earlier, our plots covered a very wide range of forest structural conditions. Consequently, snag densities varied widely among plots and were so markedly skewed (Ganey and Vojta 2012) that the utility of standard estimators of central tendency such as the mean or median was limited. Therefore, we used Huber’s M-estimator, a generalized maximum–likelihood estimator that provides robust estimates in distributions containing outliers (Huber and Ronchetti 2009), to estimate central tendency in snag density (by year) and change in snag density (by sampling period). We estimated this parameter and associated 95 % bias-corrected confidence intervals using 1,000 bootstrap iterations (Efron and Tibshirani 1993) in IBM SPSS Statistics v 23 (IBM SPSS Statistics, IBM Corp., Armonk, NY, 2015). All subsequent references to mean values for these and other parameters refer to Huber’s M–estimator. For consistency, we used the same methods for tree populations, although those populations were far less skewed than snag populations.

The skewed distributions for snag densities also limited the utility of standard hypothesis tests. Consequently, we used the bootstrapped confidence intervals discussed above to assess significance of observed differences in snag density. For comparisons between years or sampling periods, we assumed that confidence intervals that did not overlap between pairs of years or sampling periods indicated that those years or sampling periods differed significantly from each other. For change during sampling periods, we assumed that confidence intervals that did not overlap zero indicated that snag (or tree) density changed significantly during that period.

Changes in snag and tree density

We estimated snag and tree density within each plot for each year, then summarized snag and tree density across plots within forest type for all years. Because snags were uniquely marked, we also were able to estimate numbers of new snags recruited and existing snags lost during each sampling period. Numbers of new snags may be slightly overestimated for some plots and sampling periods, because the metal tags used to mark snags sometimes melted in plots that experienced moderate- to high-severity fire, making it difficult to distinguish new snags recruited post-fire from snags present before the fire. We suspect that this bias was small, however, for two reasons. First, these disturbances affected relatively few plots during most sampling periods (Table 1). Second, fires hot enough to melt the metal tags also burned most existing snags, meaning that most apparent “new” snags in these areas likely were new. Detection rates for standing snags in unburned plots, estimated using mark-resight methodology, were very high, ranging from 0.983 to 0.994 across major snag species represented (Ganey et al. 2015: Table 3).

For each sampling period, we estimated numbers of snags gained and lost by plot, then summarized these parameters across plots within forest type. Because numbers of snags available at the start of a sampling period varied, we also estimated standardized snag loss as:

Percentage of snag loss (%) = (snags lost from time t to t + 5 years)/(snag density at time t) × 100.

Composition of snag and tree populations

We compared species composition and diameter-class distributions of snag and tree populations among years within forest type using chi-square tests (Conover 1980). We pooled snags or trees across plots within forest type for these comparisons.

We included six major species groups in comparisons of species composition in mixed-conifer forest (white fir, Douglas-fir, quaking aspen, ponderosa pine, Gambel oak, and a sixth group [Other] representing all other species). We included only three species groups (ponderosa pine, Gambel oak, and Other) in tests in ponderosa pine forest. Other species were present in such small amounts in this forest type that including those species as separate categories resulted in multiple cells with expected values < 5, potentially biasing test results (Conover 1980). In contrast, no cells had expected values < 5 after collapsing categories. We recognized five diameter classes in tests involving diameter-class distributions: 20–29, 30–39, 40–49, 50–59, and ≥ 60 cm dbh. We estimated percent change in individual snag species or diameter classes during the 20-year study period as:

Percent change (%) = (2017 value – 1997 value)/(1997 value) × 100.

Climate data

Because climate can strongly affect tree mortality and thus snag creation, we obtained and summarized data on annual precipitation (AP) and cooling degree days (CDD) for the National Weather Service (NWS) weather station at Pulliam airport in Flagstaff, AZ (http://w2.weather.gov/climate/xmacis.php?wfo=fgz; downloaded 14 Dec 2017). We assumed that this station, which was centrally located within the study area at an elevation of 2136 m, provided a valid index to annual variation in broad-scale climate patterns across the study area. We restricted our analysis to the period from 1950 to 2016 because data on CDD were not available for most years prior to 1950. We assumed that this 67-year period, which included the well documented mid-20th century drought prominent throughout the southwest (Hereford 2007), provided a valid index to broad-scale climate patterns across the study area in recent times.

CDD were calculated by NWS using a base temperature of 65°F. Thus, CDD was calculated only for days when the mean temperature was greater than 65°F as:

CDD = mean daily temperature (°F) − 65°F.

We used the data from NWS directly rather than converting temperatures to °C, because the index is useful primarily in a relative context and the actual units were not important. In contrast, we converted AP from inches to cm for consistency with other climate literature. We summarized data for AP and CDD in terms of the difference between annual means for those parameters and the 1950–2016 mean for the indicated parameter, which provides an index showing the magnitude and direction of annual deviation from longer-term normal values.

These measures provided general indices of relative wetness (AP) and warmness (CDD; http://www.weather.gov/key/climate_heat_cool). To provide a third index combining relative wetness and warmness, we used data for Palmer’s Drought Severity Index (PDSI; Palmer 1965; https://www7.ncdc.noaa.gov/CDO/CDODivisionalSelect.jsp#; downloaded 18 March 2019) for the region containing our study area (climate division 2; https://www.esrl.noaa.gov/psd/data/usclimdivs/data/map.html#Arizona). This index uses temperature data and a physical water balance model to estimate potential evapotranspiration and relative drought severity (Palmer 1965, see also https://climatedataguide.ucar.edu/climate-data/palmer-drought-severity-index-pdsi). Index values range from − 10 to 10, with lower and higher values indicating drier and wetter conditions, respectively (Palmer 1965).

Results

Changes in snag and tree density

In all years, mean snag densities in mixed-conifer forest were approximately five times greater than mean densities in ponderosa pine forest (Fig. 1a). Within forest type, mean snag density was similar in 1997 and 2002, increased significantly between 2002 and 2007, and declined slightly after 2007. In mixed-conifer forest, mean snag density remained significantly elevated from 2007 to 2017 relative to earlier years. In contrast, the confidence interval for mean density in 2017 in ponderosa pine forest overlapped with all earlier years, indicating convergence toward pre-2007 levels of snag density from 2012 to 2017. Relative to snag density in 1997, peak mean density in 2007 was 86 % and 79 % greater in mixed-conifer and ponderosa pine forest, respectively, and mean snag density in 2017 was still 83 % and 51 % greater than 1997 density.

Fig. 1

Trends in snag and tree populations in northern Arizona mixed-conifer and ponderosa pine forest, 1997–2017. a Snag density by year; snags were sampled every five years. b Density of newly-recruited snags by 5-year interval between sampling occasions. Each year shown indicates the end of the 5-year interval represented. c Density of snags lost during each 5-year interval between sampling occasions. d Percentage of existing snags lost during each 5-year interval between sampling occasions. e Tree density by year; tree populations were sampled only in 2004 and 2014. Shown in all panels are Huber’s M–estimator (solid line), a generalized maximum–likelihood estimator that yields robust estimates of central tendency in distributions containing outliers (Huber and Ronchetti 2009), and lower and upper bounds (dashed lines) of bias-corrected confidence intervals around M

Mean recruitment of new snags spiked markedly between 2002 and 2007 in both forest types but was similar among the other sampling periods (Fig. 1b). Mean snag recruitment from 2002 to 2007 was more than double recruitment in any other sampling period in mixed-conifer forest, and differed significantly from mean recruitment during all other sampling periods in this forest type. Mean snag recruitment from 2002 to 2007 was almost double any other sampling period in ponderosa pine forest, but due to wide confidence intervals differed significantly only from mean recruitment from 1997 to 2002.

Snag loss rates increased with every consecutive sampling period in both forest types, and this pattern was observed both for absolute numbers of snags lost and for standardized loss rates (Fig. 1c and d). In mixed-conifer forest, mean loss rate was significantly greater from 2007 to 2012 and 2012 to 2017 than from 1997 to 2002. Mean loss rate also was considerably greater in later sampling periods in ponderosa pine forest, but confidence intervals around those loss rates were wide and overlapped for all periods.

Mean tree density declined by 8.3 % and 0.3 % from 2004 to 2014 in mixed-conifer and ponderosa pine forest, respectively, but confidence intervals around mean estimates were wide and overlapped between years in both forest types (Fig. 1e). Mean tree densities in mixed-conifer forest in 2004 and 2014 were approximately ten and five times greater than snag densities in 2002 and 2012, respectively (the closest years in which snags were sampled; compare Fig. 1a and e). In ponderosa pine forest, mean tree density was approximately 20 times greater than snag density in these same among-year comparisons. Consequently, although the observed level of tree mortality was sufficient to significantly alter snag density across time, it did not significantly alter tree density. Note, however, that our initial sample of tree populations in 2004 likely occurred after the major mortality pulse during the study, which appeared to peak from 2002 to 2003 (Negron et al. 2009; Kane et al. 2014).

Snag recruitment and climate

As discussed above, recruitment of new snags was relatively low and similar for three of the four sampling periods, but increased significantly during the period from 2002 to 2007 in both forest types (Fig. 1b). This spike in snag recruitment followed an extremely warm and dry year (2002). In the summer of 2002, PDSI reached its lowest value for our study period (Fig. 2a), as well as its most extreme value in the past 100 years (Koepke et al. 2010). The year 2002 was not the driest year during our study (Fig. 2b), but was by far the hottest year during the study (Fig. 2c), suggesting that the extremely low PDSI in 2002 was driven more by extremely high temperature than by low moisture. PDSI also exceeded the threshold for extreme drought in the summers of 2007 and 2015 (Fig. 2a), but did not reach the extreme level reached in 2002, did not remain low as long as during 2002–2003, and 2007 and 2015 overall were not extremely dry or warm relative to 2002 (Fig. 2). Snag recruitment did not spike following these later years.

Fig. 2

Trends in climate in northern Arizona, 1997–2017. a Palmer’s Drought Severity Index (PDSI; Palmer 1965; Alley 1984; Heim 2002) for the region containing our study area by year. This index ranges from − 10 to 10, with higher and lower values indicating wetter and drier conditions, respectively. Horizontal reference lines indicate thresholds for severe (–3.0; dashed line) and extreme (–4.0; solid line) drought, respectively (Palmer 1965). b Difference between annual precipitation (AP) during the study and the 1950 through 2016 mean at the National Weather Service (NWS) station located at Pulliam airport, Flagstaff, AZ, elevation 2136 m. Values above and below the horizontal reference line indicated years that were wetter or drier than average, respectively. c Difference between annual Cooling Degree Days (CDD) during the study and the 1950 through 2016 mean at the National Weather Service (NWS) station located at Pulliam airport. Values above and below the horizontal reference line indicated years that were warmer or cooler than average, respectively. All data sources are described in text

Composition and structure of snag and tree populations

Species composition of snag populations differed significantly among years in both forest types (Fig. 3a; mixed-conifer: χ² = 905.2, df = 20, P < 0.001, n = 14,857; ponderosa pine: χ² = 24.5, df = 8, P < 0.001, n = 3,287). In mixed-conifer forest, the primary temporal trend was a large increase in white fir snags, with smaller increases in all other species. In ponderosa pine forest, proportional representation decreased for ponderosa pine and increased for Gambel oak snags.

Fig. 3

Trends in species composition of snag populations in northern Arizona mixed-conifer (n = 53 plots) and ponderosa pine (n = 60 plots) forest. a Number of snags by species and year. b Number of newly-recruited snags by species and 5-year period ending in the indicated year. c Number of snags lost within the 5-year period ending in the indicated year, by species. Species acronyms: ABCO = white fir, PIPO = ponderosa pine, POTR = quaking aspen, PSME = Douglas-fir, and QUGA = Gambel oak. Not shown is a composite group comprised of all other species. Snags were pooled across plots within forest type

Much of the change in species composition was driven by patterns in recruitment of new snags. In terms of magnitude, recruitment of new snags was dominated by white fir in mixed-conifer forest and by ponderosa pine in ponderosa pine forest (Fig. 3b). Timing of recruitment of new snags varied among species in mixed-conifer forest, peaking from 2002 to 2007 for white fir and quaking aspen, from 2007 to 2012 for ponderosa pine, and increasing steadily throughout the study for Douglas-fir. In ponderosa pine forest, recruitment of both ponderosa pine and Gambel oak snags peaked between 2002 and 2007. Species composition of snags lost generally was dominated by the same species as newly-recruited snags, but in lower numbers. Increases in snag loss lagged increases in snag recruitment by 5–10 years (Fig. 3c).

Diameter-class distributions of snags also differed significantly among years in both forest types (Fig. 4a, mixed-conifer: χ² = 67.2, df = 16, P < 0.001, n = 14,857; ponderosa pine: χ² = 37.3, df = 16, P = 0.002, n = 3,476). In both forest types, numbers of snags increased across all size classes, with by far the greatest increase occurring in the smallest size classes. These changes largely were driven by recruitment patterns for new snags, with large numbers of snags recruited in the smaller size classes (Fig. 4b). As with species composition, snags lost were dominated by the same size classes as snags recruited, but in lower numbers and with peaks in snag loss lagging peaks in snag recruitment by 5–10 years (Fig. 4c).

Fig. 4

Trends in diameter-class distribution of snag populations in northern Arizona mixed-conifer (n = 53 plots) and ponderosa pine (n = 60 plots) forest. a Number of snags by diameter class and year; b Number of newly-recruited snags by diameter class and 5-year period ending in the indicated year; c Number of snags lost within the 5-year period ending in the indicated year, by diameter class. Snags were pooled across plots within forest type

Species composition of tree populations differed significantly between 2004 and 2014 only in ponderosa pine forest (χ² = 7.105, df = 2. P = 0.028; mixed–conifer forest χ² = 6.640, df = 5, P = 0.359). Ponderosa pine increased from 87.4 to 90.4 % of total trees in ponderosa pine forest over this period, whereas proportions of all other species decreased. Size-class distributions for trees did not differ significantly between 2004 and 2014 in either forest type (mixed-conifer forest χ² = 8.773, df = 4, P = 0.066; ponderosa pine χ² = 5.241, df = 4. P = 0.262), and were numerically dominated by trees in the smaller size classes in both years (results not shown here).

Discussion

Snag numbers increased considerably over our 20-year study period in mixed-conifer forest, with a smaller increase in ponderosa pine forest (Fig. 1a). Increased snag density was driven largely by a pulse of snag recruitment (tree mortality) that occurred between 2002 and 2007, with snag recruitment returning to pre-pulse levels in subsequent sampling periods (Fig. 1b). Thus, increases in snag recruitment appeared more indicative of a short-term tree mortality “event” than of a longer-term pattern of elevated tree mortality (with the possible exception of Douglas-fir, Fig. 3a).

This mortality pulse appeared to be drought-mediated. Our evaluation of the influence of climate on tree mortality was correlative and cannot demonstrate cause and effect relationships. Nevertheless, two aspects of this mortality pulse appear most consistent with drought effects. The first relates to the timing of the mortality pulse, which immediately followed a single extremely warm and dry year in 2002 (Fig. 2; see also Negron et al. 2009; Kane et al. 2014) embedded within a longer-term megadrought (Williams et al. 2020). The second factor suggesting that mortality was drought-mediated was the spatial synchrony of mortality across plots within both forest types plots. This spatial synchrony in mortality appeared more consistent with widespread drought effects than with effects of more spatially limited disturbance events such as fire or insect outbreaks. Those other disturbance agents certainly contributed to mortality in some areas, but their effects were not as widespread, and likely were exacerbated by drought.

The observed mortality pulse appeared to be triggered by the extreme climate year in 2002. PDSI reached its lowest level in over 100 years in 2002 (Kopeke et al. 2010), and a tree-ring reconstruction indicated that moisture stress during 2002 reached one of the highest levels in the past 1400 years (Salzer and Kipfmueller 2005; see also Williams et al. 2020). The extremely low PDSI in 2002 appeared to be driven more by extremely hot conditions than by extremely dry conditions; 2002 was by far the hottest year during the study but was not the driest year (Fig. 2). Thus, conditions in 2002 appeared representative of what Breshears et al. (2005) defined as “global climate change type drought”, characterized by dry conditions coupled with extreme heat. Climate change models suggest that such years will occur more frequently in the southwestern U. S. in the future (Seager et al. 2007; Garfin et al. 2014).

Many later years during the study were either warmer or drier than average, or both, and annual PDSI values suggested some level of drought stress in 16 of 20 years during the study (Fig. 2a). Those later years and periods did not trigger the same levels of tree mortality seen following 2002, however. There are at least four possible explanations for this apparent discrepancy. First, the extreme conditions experienced during 2002, especially with regard to temperature (Adams et al. 2009) may have pushed many trees beyond a physiological threshold for survival (McDowell et al. 2008, 2016; Koepke et al. 2010) that was not reached again in subsequent years. Second, the lack of increased tree mortality following later years that were almost as extreme climatically as 2002 could indicate that the most vulnerable trees died in the post-2002 mortality event, leaving tree populations dominated by trees that were more resilient to drought. Third, remnant trees may not have been more drought-resistant, but drought stress may have been alleviated in subsequent years by previous mortality that reduced tree density and competition for available precipitation (Peet and Christiansen 1987). Fourth, bark beetle populations may have been higher following the extreme climate year in 2002 than in subsequent years, swamping the defenses of many trees (e.g., Negron et al. 2009; Kane et al. 2014).

At present we cannot distinguish between these hypotheses, which are not-mutually exclusive. They clearly have very different implications for future trends in forest structure and composition, however. For example, if hypothesis 1 is correct, we can expect to see similar levels of tree mortality whenever climatic extremes like 2002 occur. In contrast, both hypotheses 2 and 3 suggest that future mortality may be lower even in similarly extreme climate years. Hypothesis 4 is particularly intriguing in the context of understanding future trends, as it may suggest feedback loops that regulate drought-mediated mortality. Bark beetle populations are known to increase in response to drought (e.g., Fettig et al. 2007; Huang et al. 2020), exacerbating mortality in drought-stressed trees (Huang et al. 2020). This increase in bark beetles, along with the associated increase in tree mortality, provides increased food and nest resources for bark insectivores, many of which nest in cavities in snags, resulting in short-term (2–5 years) increases in populations of woodpeckers and other bark insectivores (Koplin 1969; Edworthy et al. 2011; Saab et al. 2019). These populations of bark insectivores then aid in regulating bark beetle numbers (Koplin 1972; Fayt et al. 2005). This raises the possibility that severe drought years following 2002 may have seen lower mortality partly because of these feedback loops operating following earlier spikes in snag recruitment.

Regardless of the explanation for the observed mortality event, that event was large enough to significantly alter species composition and structure of snag populations during the study. All species of snags showed increased recruitment during the study, but the magnitude, timing, and duration of that increase varied among species. For example, white fir and quaking aspen, the species least resistant to drought stress (Ganey and Vojta 2011; after Niinemets and Vallardes 2006), responded strongly and quickly in mixed-conifer forest, with relatively large increases in snag recruitment peaking in the sampling period from 2002 to 2007 and returning to pre-2002 levels in subsequent periods (Fig. 3b). In contrast, the more drought-resistant Douglas-fir did not respond as strongly from 2002 to 2007, but instead snag recruitment increased over subsequent sampling periods. Ponderosa pine also responded fairly strongly in mixed-conifer forest from 2002 to 2007, but snag recruitment for this species did not peak until the period from 2007 to 2012 and remained elevated from 2012 to 2017.

These differences in magnitude and timing of species-specific snag recruitment resulted in significant changes in species composition of snag populations. In mixed-conifer forest, snag populations in 1997 and 2002 were co-dominated by white fir, ponderosa pine, and Gambel oak in approximately equal numbers (Fig. 3a). By 2007, those populations were dominated by white fir snags, and white fir continued to dominate these populations through 2017, although it declined in abundance in later years. In contrast, relative abundance of Douglas-fir increased across all years, peaking in 2017. This pattern may indicate that, unlike some other species, Douglas-fir was responding more to cumulative drought stress than to a one-time mortality event.

In the less species diverse ponderosa pine forest type, snag populations were dominated by ponderosa pine in all years, with Gambel oak the only other species present in considerable numbers. Recruitment of ponderosa pine snags peaked earlier in the warmer and drier ponderosa pine forest (2002–2007; Fig. 3b) than in mixed-conifer forest. That increase largely was offset by increased loss rates of ponderosa pine snags following 2007, however (Fig. 3c).

Diameter distributions of snag populations also changed significantly during the study in both forest types (Fig. 4a). Numbers of snags recruited increased for all size classes from 2002 to 2007, but the magnitude of increase was far greater for the smallest size classes than for the largest size classes. This imbalance was partially but not entirely offset by higher loss rates for snags in the smallest size classes (Fig. 4c). The net result was increased relative abundance of snags in the smallest size classes and decreased relative abundance of snags in the largest size classes.

Note that, although these changes in species composition and diameter distributions were driven primarily by snag recruitment, they were mediated by patterns in snag longevity. Longevity was related primarily to snag species, diameter, and height (Ganey et al. 2015). Standing rates for white fir and ponderosa pine snags were lower than for other major species (Ganey et al. 2015: Table 2), which tended to work opposite the large increases in recruitment for those two species (Fig. 3b). And standing rates for smaller diameter snags were lower than those for larger diameter snags across species (Ganey at al. 2015: Table 5). Again, this suggests that patterns in snag loss worked opposite to patterns in snag recruitment, which was concentrated in the smaller size classes (Fig. 4b). Snag loss rates also may have been influenced by bark beetle activity, especially later in the study. Snags resulting from bark beetle infestations appear to fall relatively quickly (Chambers and Mast 2014; Rhoades et al. 2020), and this also may have contributed to the increase in snag loss rates observed over time in this study (see Fig. 1c and d).

The effect of these changes in composition and structure of snag populations on snag use by native wildlife is unknown. The increase in snag abundance should provide increased foraging opportunities for native birds (Bull et al. 1997), at least in the short term. Effects on nesting birds and/or roosting bats are more unclear. Snag populations in both forest types became increasingly dominated by smaller snags during the study, and populations in mixed-conifer forest became increasingly dominated by white fir snags. These snags are less likely to be used by cavity nesting birds or roosting bats than are larger snags or snags of other species (Rabe et al. 1998; Ganey and Vojta 2004; Solvesky and Chambers 2009). They also generally decay and fall more rapidly than do larger snags or other species of snags (Ganey et al. 2015). The increased abundance of such snags may have contributed to the observed increase in snag loss rates across time in both forest types during the study (Fig. 1c), despite the large number of “new” snags in these populations. Thus, although snag numbers increased during the study, and that increase likely provided a temporary increase in foraging substrates and food resources for birds, many of the recruited snags may represent an ephemeral resource that is not useful to cavity-nesting birds (Ganey and Vojta 2004). Note that these rapidly-decaying and falling snags contribute to increases in log populations, however, and this may benefit native wildlife that feed on, hide in, or den in or under down logs (Bull et al. 1997; Ganey and Vojta 2017).

In contrast to snag populations, tree mortality did not appear to be adequate to significantly alter tree populations during the study. Our ability to evaluate changes in tree populations during the study period was hampered by the more limited sampling of those populations, however, and especially by the timing of sampling. We assessed change in snag populations over a 20-year period, but assessed change in tree populations only over a 10-year period embedded within that longer period. More importantly, we initially sampled tree populations in 2004, after the peak of the tree mortality pulse triggered by extreme drought conditions in 2002 (2002–2003; Negron et al. 2009; Ganey and Vojta 2011; Kane et al. 2014). Thus, tree populations likely showed greater change during the 20-year study period than we were able to demonstrate. Even given these limitations, however, the much greater tree densities relative to snag densities clearly ensure that, for a given mortality level, relative effects will be much lower for tree populations than for snag populations.

Conclusions

Warming climates have profoundly affected forests throughout the world (Adams et al. 2009; Allen et al. 2010), and climate models predict further warming and drying in many areas, including the southwestern United States (e.g., Seager et al. 2007; Stocker et al. 2013; Garfin et al. 2014). Our results demonstrate the dynamic nature of snag populations and document the drivers associated with rapid changes in those populations over time. That is, a drought event in 2002 caused spatially widespread tree mortality. This mortality significantly increased recruitment of new snags, but many of the newly-recruited snags were small in diameter and/or white fir snags. These snags decay and fall faster than larger snags and other species of snags (Ganey et al. 2015), and this contributed to an increasing snag loss rate over time. This rapid loss of snags did not entirely balance the increased recruitment of new snags following the 2002 drought, resulting in higher snag densities after 2007 compared to prior periods, particularly in mixed conifer forests. Significant questions remain, however, including: (1) how long do these changes in abundance, structure, and composition of snag populations persist? (2) how common will similarly extreme climate years be in the future? and (3) will such years produce mortality pulses similar to the one observed in this study, or might future patterns of tree mortality be modified by past mortality events that alter forest structure and composition, and competition for resources? Answering these questions will require long-term studies, but is critical to understanding how changing climate may affect trends in structure and composition of future snag populations and the tree populations from which they derive. The U.S. Forest Service Forest Inventory and Analysis program (Smith 2002), which repeatedly samples plots through time across the landscape, may provide a means to accomplish the longer-term monitoring needed to answer these questions. In the meantime, the mortality observed following the extreme climatic year in 2002, coupled with climate models predicting that such years will become more common in the future (Seager et al. 2007; Garfin et al. 2014; Williams et al. 2020), supports calls for management to increase resilience and drought-resistance in these forest types (Millar et al. 2007; Stephens et al. 2013).