2.1. Sites and Field Methods
Three fires (Hayman, Jasper, and Egley) in ponderosa pine-dominated forests in the western U.S. were sampled. The Jasper Fire (33,794 ha) burned in the Black Hills of southwest South Dakota in 2000, the Hayman Fire (55,893 ha) in the Front Range of central Colorado in 2002, and the Egley Complex (56,800 ha) in Malheur National Forest of central Oregon in 2007. The three locations have generally similar temperature regimes (30 year summer average of 18 °C for Egley and Jasper, 21 °C for Hayman; winter average −2 °C and −3 °C for Egley and Jasper, respectively, and 0 °C for Hayman) but differ in the amount and timing of precipitation. Precipitation at the Egley Fire primarily falls in winter/spring, whereas the areas around the Jasper and Hayman fires receive the majority of their precipitation in spring/summer. Overall, Egley receives only ~50% and 65% of the precipitation that Jasper and Hayman receive (30 year precipitation averages are 264, 525, and 404 mm, respectively).
Field sites were stratified by three variables based on spatial data layers in raster format: (1) burn severity from the Monitoring Trends in Burn Severity (MTBS; www.mtbs.gov
, accessed on 18 October 2021) one-year post-fire delta Normalized Burn Ratio (dNBR) classified burn severity product (unburned, low, moderate, high), (2) elevation (high or low, based on the elevation spread of each fire), and (3) aspect. Aspect was transformed to TRASP on a 0–1 scale ranging from cool-wet (0 = 30°) to warm-dry (1 = 210°) based on the topographic solar-radiation index transformation [28
]. The total number of sites and year of sampling varied among fires: 19 sites were sampled at Hayman in May/June 2015 (13 years post-fire), 16 sites at Jasper in June 2015 (15 years post-fire), and 41 sites at Egley in May/June 2016 (9 years post-fire).
The MTBS product used to stratify sites is a free product commonly used by both researchers and managers for post-fire ecological studies, impact assessments, and to guide potential restoration. Continuous dNBR values are calculated as pre-fire NBR minus post-fire NBR, where NBR is a normalized ratio calculated from 30 m resolution Landsat near-infrared band (NIR) and short-wave infrared band (SWIR), i.e., (NIR-SWIR)/(NIR + SWIR). Continuous dNBR is then categorized as unburned, low, moderate, or high severity, or increased green based on thresholds determined by the MTBS analyst [29
]. The dNBR index is sensitive to changes in green vegetation, bare soil, and char, where areas with higher char and bare soil will have higher NBR values [30
Within each site, five plots were established 30 m apart in a cross formation, with the azimuth defined by the outer plot established upslope of site center according to the dominant slope. Seedlings (≥15 cm and ≤137 cm in height) and saplings (>137 cm in height with a diameter at breast height (DBH) <10 cm) were tallied by full quarters in 5.6 m radius plots until a minimum of six seedlings had been encountered. Up to six representative seedlings or saplings, if they were determined to represent post-fire regeneration based on estimated age, were measured for the length between terminal bud scars to obtain approximate age and yearly growth increments [32
]. This analysis focused on the most recent seven complete years of growth (2009–2015 for Egley and 2008–2014 for Hayman and Jasper) because field observations and previous work [32
] indicate that bud scars of recent years are more reliably identified.
The distance from each plot to the nearest live mature ponderosa pine assumed to be capable of serving as a seed source (hereafter “seed tree”) was recorded. Fractional cover of green vegetation, non-photosynthetic vegetation, mineral soil, and percent char of soil and non-photosynthetic vegetation was visually estimated in a 1 m2
microplot at the center of each of the five plots. Non-photosynthetic vegetation included woody debris, senesced grass or forbs, tree bark, or leaf and needle litter. Additionally, we measured litter and duff depth, fine woody fuel loads (1, 10, 100 h timelag classes; <1 cm, 1–2.5 cm, and >2.5–7 cm diameter, respectively) were estimated using a photoload guide [33
], and canopy cover of the overstory (trees and shrubs exceeding breast height [1.37 m]) was estimated using a convex spherical densiometer.
At each of the four peripheral plots, standing trees were recorded using a 2 m2
/ha basal area factor prism; for each tree the species, vigor (dead, healthy, unhealthy), and DBH was recorded. At only the center plot, all trees were tallied and measured in a 0.02 ha plot (8 m radius). Percent cover of tall shrubs (>1.37 m height) was ocularly estimated, and large downed woody fuel loads (1000-hr timelag class, >7 cm diameter) were estimated using the photoload guide [34
] within a 0.01 ha subplot (5.6 m radius) at the center plot only.
2.3. Statistical Methods
Of the total sites sampled, we used 36 sites for this analysis (Hayman = 12, Jasper = 11, Egley = 13) by subsetting only sites that had ponderosa pine seedlings present and that were not planted following the fire (to capture natural tree regeneration vs. planted seedlings). To ensure a roughly balanced sample among fires, sites from Egley were chosen at random from the suitable subsample. Each year’s growth was standardized on total seedling height in that year [25
], and a site-level average annual growth was determined by averaging the annual growth of all seedlings at a site. The response variable “seedling height growth” therefore represents the average annual height growth for a seedling at a given site for the past seven years of growth prior to sampling, which was standardized to account for the height of the seedling in a given year of growth. Seedling density was calculated as seedlings per hectare based on the number of seedlings counted divided by the area sampled on a site.
Non-parametric multiplicative regression (NPMR) in HyperNiche v.2 [37
] was used to determine the influence of and potential interactions among predictor variables on both seedling height growth and density. NPMR allows for predictor variables to interact in non-linear, multiplicative ways to influence the response variables [38
]. Each NPMR free search run was done with local mean model and Gaussian weighting, with default medium controls for overfitting, automatic minimum average neighborhood size (number of sites * 0.05), step size of 5, 10% maximum allowable missing estimates, and minimal backtracking search. HyperNiche automatically runs the free search as an iterative process over various combinations of predictor variables, producing thousands of models in the process. To reduce collinearity and duplication among predictor variables, predictor variables with pairwise correlations greater than 0.9 to another predictor variable were dropped in a stepwise approach to retain predictors with the stronger Spearman’s rank correlation to average seedling height growth. This resulted in a final set of 39 variables (Table 1
) that were included in the NPMR free search to examine their influence on seedling growth; 36 predictor variables were included in the seedling density free search, which excluded seedling density variables (total live seedling density, total dead seedling density, and total live sapling density) that were included as predictor variables for the growth models.
The lists of models generated by the NPMR process were sorted by the cross-validated R2 (xR2) to determine what predictor variables appeared in the majority of the strongest 100 models. This process was used to select the best model for seedling height growth and for seedling density, individually, at which point the models were evaluated for tolerances and sensitivity. Tolerances in NPRM are the standard deviations used in the Gaussian smoothing and must be interpreted based on the range of each predictor; higher tolerance indicates that a variable is less important to the model. Sensitivity values range from 1 to 0, where higher sensitivity indicates that a percent change in that predictor will result in a similar percent change to the estimate of the response variable.