Whiskeytown National Recreation Area is located at the southeastern edge of the Klamath bioregion (Skinner et al., 2006) in Shasta and Trinity Counties, just west of the city of Redding in northern California. WHIS is over 17,000 hectares in size, and is characterized by steep topography and high biodiversity, with a wide variety of forest types including oak woodlands, knobcone pine woodlands, mixed conifer and yellow pine forests, with true fir (red and white) at high elevations (Fry and Stephens, 2006). Elevation in WHIS ranges from 250 to 1,890 meters above sea level. The park also includes large areas of shrubland and sparsely-treed woodlands, especially on the south-facing slopes north of Whiskeytown Lake.

Common conifer species include ponderosa pine (Pinus ponderosa), Jeffrey pine (Pinus jeffreyi), incense cedar (Calocedrus decurrens (Torr.) Floren.), Douglas fir, sugar pine (Pinus lambertiana Dougl.), and white fir (Abies concolor (Gord. & Glend.) Lindl.). Knobcone and gray pine (Pinus sabiniana) are common at low elevations, and some red fir (Abies magnifica var. shastensis) is found in the high elevations of the park. A variety of hardwood species are also found throughout the park, including California black oak (Quercus kelloggii Newb.), canyon live oak (Quercus chrysolepis Liebm.), tanoak (Notholithocarpus densiflorus), golden chinquapin (Chrysolepis chrysophylla), and dogwood (Cornus nuttallii Audubon). See Smith et al. (2021) for a detailed description of vegetation found in the park.

Fires were historically common in the area but decreased in frequency after 1850, with an almost complete absence of fire in the latter half of the 20th century following the widespread adoption of fire suppression (Fry and Stephens, 2006). However, fire has not been completely absent in the last few decades. Starting in the mid-1990s, the National Park Service introduced fuels reduction and other restoration projects, often involving the use of prescribed fire. In addition to these treatments, wildfire has also occurred within the park boundary. Most notably, the Shasta-Trinity Unit Lightning Complex Fire burned through the park in 2008, which burned 4% of the park area at high severity (MTBS, 2018).

The Carr Fire started on July 23rd, 2018, and actively burned for 38 days before containment on August 30th the same year (https://www.fire.ca.gov/incidents/2018/7/23/carr-fire). The fire burned nearly 93,000 hectares, including nearly the entire area of WHIS, destroyed over 1,000 homes in nearby communities, resulted in 8 human deaths, and had an estimated damage cost of >$1.6 billion (USD). The fire occurred under abnormally hot, dry, and windy conditions, and made significant runs when terrain and wind aligned, burning a large portion of WHIS in a single 24-hour period.

Our study focused on the yellow pine and mixed conifer forests of WHIS, which covered over 5,200 ha within the park prior to the Carr Fire ( Figure 1 ). According to data produced by MTBS, over 65% of the conifer forests (3,450 ha) burned at high severity ( Figure 1 , MTBS, 2020). High severity patches were extensive, especially on the steep slopes of the Brandy, Boulder, and Mill Creek drainages.

We leveraged field data from two separate sources for this analysis. We used data from a preexisting network of 0.1 ha fire monitoring (FMH) plots in WHIS, which are designed to monitor changes in fuels and forest structure at random locations following fire (USDI National Park Service, 2003). We selected all plots where seedling density and overstory trees (>15 cm diameter at breast height or 1.37 m, hereafter DBH) were sampled after the Carr Fire, for a total of 23 plots. The FMH plot data captured species, DBH, and status for all overstory trees in the plot based on a 15 cm DBH threshold. The same information was collected for pole size trees (stems ≤ 15 cm DBH) in 0.025 ha subplots. The FMH plots also captured seedling tallies (stems ≤ 2.5 cm DBH) by size class, species, and status in a 50 m² subplot (0.005 ha). The FMH plots also captured vegetation cover using two 50 m point intercept transects, where species and height of vegetation was recorded at 30 cm intervals for a total of 166 points per transect. We averaged the shrub cover estimates from both transects to get a plot level estimate of shrub cover.

In addition to the FMH plots, we established a total of 108 0.1 ha plots in the summers of 2020 and 2021 ( Figure 1 ). We sampled plots over most of the elevation range of the park, with plot elevations ranging from 330 to 1,870 m. For most of these plots (n = 76) we used GIS to randomly select sampling locations within conifer forests according to the following criteria: no slopes > 50%, > 50 m and < 500 m from a road, 100 m from non-vegetated areas (e.g., Whiskeytown Lake), and at least 50 m from one another. The sampling locations were stratified across burn severity and predicted regeneration. Regeneration was predicted using the models developed in (Stewart et al., 2021) using the poscrptR R package (Wright et al., 2020). We used vegetation alliance polygons created for the National Park Service (Fox et al., 2006) to identify conifer forests by selecting polygons with Douglas fir, mixed conifer, ponderosa pine, or red fir forest alliance types. We established 12 of the plots specifically in random locations within unoccupied aerial survey (UAS) sampling areas to assist with future vegetation mapping projects using the same criteria described above, with the exception that plot locations were stratified by burn severity and whether the forest alliance was conifer (described above) or oak, including black oak, blue oak, canyon live oak, interior live oak, Oregon white oak, and tanoak forest and woodland alliances. We established an additional 20 plots to match earlier randomly sampled plot surveys of old-growth conifer forests (Leonzo and Keyes, 2010). Three of the original plots from Leonzo and Keyes (2010) could not be reached due to unsafe conditions, so we sampled three replacement plots using the same random sampling criteria outlined in Leonzo and Keyes (2010).

At each plot, we recorded DBH, species identity, and live and dead status for all standing trees > 15 cm DBH over a 0.1 ha circular plot. We measured sapling (stems 0.1 to 15.0 cm DBH) stem diameter, species location and status in a 0.011 ha subplot at the plot center. We also collected tallies of live seedlings (stems < 1.37 m in height) in the subplot, recording species and height class (0.1–5, 5–10, 11–25, 26–50, 51–75, 76–100, and 100–136 cm). On nine occasions the seedling subplot was reduced in size to 0.001 ha or 0.005 ha because the large numbers of seedlings present made larger subplot sizes impractical. We also visually estimated the cover of shrubs, forbs, and grasses within the plot area. Cover estimates were binned into the following classes: 0–1%, 2–5%, 6–25%, 26–50%, 51–75%, 76–95%, and 96–100%.

We normalized the data for each plot type (FMH and 0.1 ha field plots) by putting basal area and seedlings on a common per hectare scale and binning shrub cover into 25% class bins. Seedlings were not always identifiable beyond the genus level, and we could not determine the species of the lidar mapped trees from aerial imagery. Therefore, we report results for our analyses on seedling density for all conifers combined. We initially performed the analysis at the genus level, with broadly similar results. Plot sizes were similar between plot types, so we expected that any difference in seedling detection probability between plot types due to sampled area is likely to be small.

We used whitebox tools (Wu, 2021) to estimate topographic wetness index from a 1 m-resolution lidar digital elevation model (DEM). First, we breached depressions to remove sinks. We then calculated D-infinity flow direction and flow accumulation for each cell, which we used to calculate topographic wetness index:

We calculated heat load index following methods outlined in McCune and Keon (2002) using the spatialEco package (Evans, 2021). We extracted the average heat load index, topographic wetness index, slope, and elevation from the 1 m DEM for the area of each seedling subplot.

We used the brms package (Bürkner, 2017) in R (R Core Team, 2021) to fit Bayesian generalized linear models estimating the influence of estimated seed availability, shrub cover, topographic wetness index, elevation, and heat load index on postfire conifer regeneration. There were several plots without seedlings, so we used negative binomial hurdle models to estimate both the influence of covariates on the probability of encountering no conifer seedlings as well as the effects on conifer seedling abundance, conditional on their being present (Steel et al., 2013). The model takes the following form:

Where N B is the negative binomial distribution, β 0 is the intercept, and β x i are the regression coefficients for the i t h predictor. The count portion of the model is a zero truncated negative binomial with a log link and shape parameter ϕ :

We used the same predictors in the hurdle and count portions of the model.

We fit separate models with different seed availability variables: three models with seed shadows derived using each of the three dispersal function parameterizations, and an additional model that used the distance to the nearest surviving tree. We standardized all continuous predictors by subtracting the mean and dividing by two standard deviations, which puts continuous and binary predictors on a common scale (Gelman, 2008). We used weakly informative priors in all models, where β 0 is the intercept, β x are the regression coefficients, and ϕ is the shape parameter of the negative binomial distribution:

We ran all models for 2,000 iterations and ensured that all R ^ values did not exceed 1.01. We checked the model fit using posterior predictive checks and leave-one-out cross-validation (LOO, Vehtari et al., 2017). We compared models using the expected log pointwise predictive density ( e l p d ^ ), which estimates the predictive accuracy of the model for each data point held out during LOO (Vehtari et al., 2017). Models can be compared by differencing e l p d ^ estimates ( Δ e l p d ^ ), and the standard error of Δ e l p d ^ can characterize the uncertainty in the model comparison. Generally, models with Δ e l p d ^ of less than four have similar predictive performance (Sivula et al., 2020). Here, we report the on the parameters for the model with the highest e l p d ^ . We calculated mean absolute error (MAE) for each model as an additional measure of model performance. Along with conditional effects for the hurdle and count portions of the model, we report the probability of direction, which describes the probability that the parameter is the same sign as the median of the posterior distribution (Makowski et al., 2019).

Several plots were located in areas that burned multiple times, including 26 in the 2008 Shasta-Trinity Unit Lightning Complex fire. Reburned areas may have lower conifer regeneration than similar areas that have not experienced recent burns if surviving seed trees are killed (Tepley et al., 2017; McCord et al., 2020). We assumed that multiple burns would largely effect any regeneration we observed by removing potential parent trees. We modeled seed availability from the mapped surviving trees, so we elected not to include whether a plot was reburned in any of the models as this was unlikely to add additional information. The data can be found at Wright et al. (2023).

LOO suggested that the models using seed availability extracted from seed shadows modestly outperformed the model using simple distance to live trees ( Supplementary Table 1 , Figure 6 ), with e l p d ^ greater than four at least one standard error from zero. The medium distance model had the highest e l p d ^ , though there was little evidence for major differences between the models using seed availability extracted from seed shadows with e l p d ^ less than four in all cases. The short distance dispersal model had the lowest MAE, but MAE was similar across the seed shadow models. As with LOO, MAE suggested that the distance to live tree model had the worst performance, with Δ MAE of 186 then the next best model. As with LOO, MAE for the models using dispersal kernels were very similar, although the model using the short-distance dispersal had lower MAE than the medium-distance model.

The medium distance model selected by LOO suggested that increased seed availability was associated with increased conifer seedling presence and density ( Figures 7 , 8 ), with probability of direction approaching 1 for both the hurdle and count components of the model. The effects of elevation were more uncertain. The best performing model indicated that increasing elevation resulted in fewer seedling observations (probability of direction 0.96), though there was little strength of evidence for a similar effect on seedling abundance (probability of direction 0.58). There was also little evidence for a consistent effect of topographic wetness index on either seedling presence or abundance, with probability of direction 0.78 and 0.56, respectively). There was some evidence that conifer seedlings were less likely to be found in areas with a high heat load index (probability of direction 0.95). The effect of heat load index on seedling density was more uncertain, but also suggested a positive relationship (probability of direction 0.87).

The effects of shrub cover were mixed. The model suggested that conifer seedlings were more likely to be observed (hurdle portion of the model) with shrub cover in the 26–50% range than in 0–25% (probability of direction 0.91), though the model also suggested that seedling abundance (density portion of the model) was lower in the 26–50% shrub cover class than in lower shrub densities (probability of direction 0.94). Evidence was much weaker for effects in the higher shrub cover classes (probability of direction below 85% in all cases), though it should be noted that the highest two shrub cover classes had very few observations (five and six plots, respectively).

We were unable to directly parameterize the dispersal kernels we used to create the seed shadows, so we chose to limit our analysis to isotropic kernels, even though anisotropic kernels generally outperform simple isotropic kernels where they have been applied (Savage et al., 2011). Because we selected the parameterizations a priori, we were necessarily limited in the possible number of parameterizations that could reasonably be assessed, as well as the number of potential variables that may influence the seed shadows. For example, we elected to ignore the influence of terrain and wind. The effect of wind on seed dispersal is well documented (Greene and Johnson, 1996; Sánchez et al., 2011; van Putten et al., 2012), though the results for the influence of terrain have been more mixed and are likely scale- and species-dependent (Donato et al., 2009; Katul and Poggi, 2012; Peeler and Smithwick, 2020). There is little doubt the rugged, steep terrain in WHIS had some influence on seed dispersal, not the least of which is the effect terrain would have on local wind patterns. We also ignored the effect neighboring trees and other obstacles might have had on seed dispersal, which likely contributed to model error because seeds can disperse further when the parent tree is in or near an open area than when it is surrounded by neighbors (Greene and Johnson, 1996). Our models also attempted to associate established seedlings with seed shadows that were modeled assuming a single dispersal agent (i.e., wind), ignoring any potential additional dispersal agents that may have changed the distribution of seeds, such as the movement of cones or seeds by animals (Rogers et al., 2019).

Additionally, seed availability does not necessarily directly translate to seed establishment and survival. Millerón et al. (2013) found large disparities between seed dispersal kernels and kernels from established saplings, suggesting that dispersal alone is insufficient to capture patterns in recruitment and survival. Although we attempted to include variables that might distinguish between site quality and thus the probability of seedling establishment, we simply could not precisely determine the effects of site quality (i.e., soil moisture or temperature, see Wooten et al., 2022) with our data. For example, the presence of downed logs and shrubs can provide protection and increase soil moisture for regenerating seedlings (Tappeiner and Helms, 1971; Landesmann and Morales, 2018; Marcolin et al., 2019). Furthermore, many conifer seeds are consumed by small mammals and birds before germination (Gashwiler, 1970; Zwolak et al., 2010), which would change the distribution of available seed if the predation pressure was not spatially uniform (Janzen, 1970; Connell, 1971). Beyond predation, the factors influencing seed germination and seedling survival are complex, and include climate, competition, and abiotic factors (Irvine et al., 2009; Tepley et al., 2017; Stewart et al., 2021). We expect that the influence of climatic variables on postfire recruitment will become more apparent in coming years, especially since temperature has been unusually high in the postfire years (data not shown).

We tried to account for variation in fecundity with tree size when calculating estimates of seed availability, and to consider the amount of time available for seeds to disperse by including the number of seasons between the fire and when the plot was sampled. However, there is temporal variation in seed production in most conifer species (Clark et al., 1999b). Masting, the synchronization of periodic seed production between plants (Kelly, 1994), may also have influenced conifer regeneration after the Carr Fire. Masting has been demonstrated in many conifer species, including many of the species included in this study (Wright et al., 2021). Whether or not we sampled following a mast year undoubtedly influenced the total number of seedlings available, and thus the error between real and estimated seed availability. Finally, there was likely substantial model error arising from our inability to identify the species or even survival status of trees using aerial imagery or lidar, which may have obscured the advantages of modeling seed dispersal using individual tree locations.