Multi-scale assessment of post-fire tree mortality models

Additional keywords: fire mortality modelling; FOFEM; Rim Fire; Sierra Nevada mixed-conifer; Smithsonian ForestGEO; spatial patterns; Yosemite Forest Dynamics Plot.

First-order fire effects models are essential tools used by land managers in fire-prone forest ecosystems. These models are used to estimate the direct effects of fire including tree mortality, soil heating, fuel consumption, and smoke production (Reinhardt and Dickinson 2010). Of particular interest to forest land managers are tree mortality models, which are used to estimate the probability of mortality for individual trees or proportion of mortality for stands (Woolley et al. 2012). These tree mortality models are used in post-fire landscapes to estimate fire severity and assess changes to stand structure, and to develop salvage marking and hazard tree guidelines (Hood et al. 2007). Mortality models are also used by managers in pre-fire applications to predict fire effects, conduct landscape-scale risk assessments, and to develop silvicultural prescriptions and prescribed fire treatments (Reinhardt and Dickinson 2010).

There are two main types of first-order tree mortality models: process-based models that use a mechanistic approach to simulate the processes involved in fire spread, heat-caused injury, and subsequent tree mortality (e.g. Michaletz and Johnson 2006), and empirical logistic regression models that use individual tree-level explanatory variables (e.g. bark thickness, crown base height, crown volume scorched (CVS)) to predict probability of mortality (e.g. Ryan and Reinhardt 1988). While process-based approaches have a strong theoretical basis and contribute to our understanding of exactly how a tree is killed by fire, the complexity of these models has precluded their widespread use among managers (Woolley et al. 2012), and their focus on direct fire damage fails to capture the suite of biological and ecological processes that contribute to delayed mortality 1 to 5 years post-fire (e.g. bark beetles and structural failure, Ryan and Amman 1996). Although empirical logistic regression models lack a mechanistic representation of tree mortality, their simplicity and implicit integration of both immediate and delayed mortality processes (mortality is typically assessed 3 years post-fire) have made them the most practical and widely accessible tool for modelling fire-related tree mortality among both researchers and land managers (e.g. Ryan and Amman 1994; Reinhardt and Crookston 2003; Sieg et al. 2006; Hood et al. 2007; Lutes 2016; Grayson et al. 2017).

Empirical logistic regression models are the foundation for the tree mortality models within widely used fire effects modelling software packages including the First Order Fire Effects Model (FOFEM), BehavePlus, and the Fire and Fuels Extension to the Forest Vegetation Simulator (FFE-FVS). Although over 100 logistic regression models for western North American conifer species have been developed (Woolley et al. 2012), the tree mortality models incorporated within these software packages are based on a single logistic regression model known as the Ryan and Amman (R-A) model. The R-A model was originally developed in the late 1980s by Ryan and Reinhardt (1988), refined by Ryan and Amman (1994), and most recently updated by Hood and Lutes (2017) who parameterised a set of 12 species-specific models that have been incorporated into FOFEM since version 5.7

The R-A model (hereafter FOFEM_RA) is perhaps the single most widely used tree mortality model (Hood et al. 2007; Reinhardt and Dickinson 2010), but it is based on a relatively small sample of trees (n = 2356, Ryan and Reinhardt 1988) and is infrequently validated with independent data (see Hood et al. 2007 and Kane et al. 2017 for validations with many species, and see Ganio et al. 2015 and Ganio and Progar 2017 for validations for Pinus ponderosa and Pseudotsuga menziesii). The recent development of species-specific versions of this model (hereafter FOFEM_SP) improved model performance and increased sample sizes for certain species (average n per species = 1403 trees, Hood and Lutes 2017), but validation of these species-specific models with independent data is still rare (but see Grayson et al. 2017 for a recent validation of FOFEM_SP).

Additionally, the mortality predictions made with these empirical models are inherently reflective of the climatic conditions during which the parameterisation data were collected. As climate influences the susceptibility of trees to fire-related mortality (van Mantgem et al. 2013; Stephens et al. 2018), empirical models developed under past climates may under-predict mortality from fires that burn under hotter and drier conditions.

The FOFEM_RA and FOFEM_SP models were not parameterised with small-diameter stems (<10 cm diameter at breast height, DBH; there are two exceptions – the lower diameter limit was 8 cm for Pseudotsuga menziesii, Ryan and Reinhardt 1988; and 6 cm for ‘yellow pine’, Hood and Lutes 2017), and validation of these models with small stems is rare (but see Engber and Varner 2012 and Kane et al. 2017). Some studies have developed logistic mortality models specifically for small-diameter stems (e.g. Battaglia et al. 2009), but these studies did not conduct a validation of the exact models within FOFEM. Trees <10 cm DBH are more abundant than stems ≥10 cm DBH in many forests, especially if fire has been suppressed for many decades. They influence future fire behaviour and act as ladder fuels, and they are ecologically important as components of the understory (re-sprouters) or advanced regeneration (surviving conifers). Previous studies have found that crown scorch influences mortality differently for small-diameter stems compared with mature trees (Kolb et al. 2007; Engber and Varner 2012), and this interaction may compromise the accuracy with which FOFEM predicts mortality for small-diameter trees.

Previous validations of the FOFEM_RA and FOFEM_SP models have also not assessed how stand-level accuracy may vary across a range of spatial scales. Stand-level accuracy is typically assessed by grouping trees according to their probability of mortality, then calculating the proportion of observed mortality within each group (sensu Hood et al. 2007). Although this approach may be used to estimate stand-level accuracy without requiring a spatially-explicit (i.e. stem-mapped; Lutz et al. 2018a) dataset, it implicitly assumes that fire-related mortality is a spatially homogeneous process. As mechanisms of delayed fire mortality can be spatially auto-correlated (e.g. bark beetle activity is patchy, and trees may survive higher levels of fire damage if they are in a more mesic area), stand-level model accuracy may therefore be modified by the presence or absence of these neighbourhood-level variables. The positive and negative effects associated with spatially correlated mechanisms of delayed mortality may equalise if stand-level accuracy is assessed at a large enough scale, but the scale at which this happens remains unknown. A multi-scale assessment of stand-level accuracy would provide a more robust estimate of model performance, and would enable us to quantify the scale at which local neighbourhoods may modify delayed mortality.

Mortality models are often used in a pre-fire context to assess potential fire effects and plan restoration activity, but a challenge to using the FOFEM models in a pre-fire context is that they require metrics of fire injury as predictor variables. Out of over 100 models that have been developed (Woolley et al. 2012), we found no model that relies exclusively on pre-fire tree attributes to predict mortality. The reason for this is obvious – fire damage is an important determinant of fire-related mortality. However, this limits the utility of these models in pre-fire planning applications because fire-damage attributes must be estimated, and this introduces an additional layer of uncertainty when interpreting the model results and evaluating their predictive accuracy. Metrics of fire damage have been considered necessary to create acceptably accurate mortality models, but these metrics are often tightly correlated with tree attributes that may be measured pre-fire, such as diameter at breast height (Fig. S1, available as Supplementary material to this paper). Given this correlation, a pre-fire model based on diameter, species, and surface fuel loads may predict mortality with an acceptable level of accuracy, and this may be useful to managers seeking to model mortality in a pre-fire context.

A final area for improvement in empirical mortality modelling is to consider the accuracy of mortality models in predicting post-fire tree spatial patterns. Fire is a spatially explicit disturbance process (Meddens et al. 2018) and is an important driver of spatially structured stand dynamics in fire-adapted forests (Larson and Churchill 2012). Spatial pattern is a key element of forest structure (Lutz et al. 2013), and it has an important influence on forest heterogeneity, resilience, and future disturbance dynamics (Stephens et al. 2008). Restoring spatial patterns and heterogeneity characteristic of fire-prone forests has become a central aim of many forest restoration efforts throughout the western United States (e.g. Allen et al. 2002; North et al. 2007, 2009; Churchill et al. 2013), but to our knowledge the accuracy with which logistic mortality models scale up to predict post-fire tree spatial patterns has not been assessed.

Our objective was to advance the science of empirical fire mortality modelling in three ways. First, we conducted a validation of the widely used FOFEM tree mortality models with a sample size of 34 174 stems ranging from 1.0 cm to 200.7 cm DBH, and we assessed both the individual tree- and stand-level accuracy of these models. Second, we developed locally-parameterised logistic regression models based on CVS and DBH to estimate of the maximum accuracy these models may attain, and to compare how different model forms influence model performance across different diameter classes. Third, we developed methods for quantifying model accuracy at a range of spatial scales and assessing the accuracy with which mortality models predict fire-induced change in tree spatial patterns.

This study is unique among existing mortality modelling literature in that the dataset we used is a census rather than a sample. Previous mortality modelling studies generally sample post-fire forests across multiple regions and fire events, selecting a subset of trees to create a reasonably balanced sample across diameter and CVS classes. This approach is well suited to developing mortality models with broad applicability, but it is less optimal for validating those models because the dataset is balanced, but not representative. An optimal validation dataset includes a large number of stems, representing a wide range of diameters, in proportion to their abundance in a particular forest type. A validation with this type of dataset permits mortality models to be tested in a context similar to how they are used by managers – to predict structural and compositional changes to a specific forest following fire. Model accuracy may then be assessed in terms of ‘percentage error’ by diameter class, a metric with very tangible implications for managers using mortality models to estimate fire effects at the stand scale. In this study, we conduct the first validation of FOFEM using a complete census of trees ≥1 cm DBH from a large-scale permanent forest plot.

We conducted this study in the lower montane, mixed-conifer zone of the Sierra Nevada, CA, USA (Fig. 1a). We used data from the Yosemite Forest Dynamics Plot (YFDP; Lutz et al. 2012, 2014b), a 25.6-ha plot affiliated with the Smithsonian ForestGEO network (Anderson-Teixeira et al. 2015; Lutz 2015). The YFDP is located in an old-growth (oldest trees >500 years old) Abies concolor–Pinus lambertiana (white fir–sugar pine) forest between 1774 and 1911 m elevation (Fig. 1b, c) with species composition representative of the Sierra Nevada white fir superassociation (Keeler-Wolf et al. 2012). Within the YFDP, all tree stems ≥1 cm DBH were tagged, identified, mapped in 2009 and 2010 (n = 34 458 live stems; Lutz et al. 2012), and tree status was updated in June 2013, 2 months before the YFDP burned. We considered the five most abundant species within the YFDP: white fir (Abies concolor [Gordon] Lindl. ex Hildebr.; 939 stems ha⁻¹), sugar pine (Pinus lambertiana Douglas; 180 stems ha⁻¹), Pacific dogwood (Cornus nuttallii Audobon; 106 stems ha⁻¹), incense cedar (Calocedrus decurrens [Torr.] Florin; 64 stems ha⁻¹), and California black oak (Quercus kelloggii Newb.; 46 stems ha⁻¹). Although A. concolor was the most abundant species, P. lambertiana had approximately the same pre-fire live basal area (30.6 and 28.8 m² ha⁻¹ respectively) and was the most abundant large-diameter (≥100 cm DBH) stem (n = 343 stems, Table 1).

Fig. 1.  The Yosemite Forest Dynamics Plot (YFDP) is located on the western slope of the Sierra Nevada, CA, USA (a), and was burned in a backfire set in the path of the California Rim Fire in September of 2013 (b). Burn severity classifications based on differenced Normalised Burn Ratio (Miller and Thode 2007). The dimensions of the YFDP are 800 m east–west by 320 m north–south, elevations range from 1774 to 1911 m (c); dots represent 34 458 trees sized by diameter at breast height (DBH; cm) and coloured (shaded in greyscale version) by crown volume scorched (CVS;%). Background: Landsat 8, natural colour image, 10 August 2017.

Table 1.  Number of pre-fire live stems of the five most common species within the Yosemite Forest Dynamics Plot (YFDP)
Immediate mortalities are those stems identified as dead in May 2014, 8 months after the California Rim Fire. Delayed mortalities are stems that survived through May 2014, but died within the following 2 years. Total fire-related mortalities includes all stems that died within 3 years of the fire. Rarefied N is the number of stems retained after the dataset was rarefied by 20% CVS bins. Crown volume scorched (CVS) and diameter at breast height (DBH) (1.37 m) metrics indicate mean values for the stems used in all analyses. Minimum DBH for all species was 1 cm. The range of CVS values was 0 to 100%, and median CVS was 100% for all species. The ‘95th’ column represents the 95th percentile threshold diameter. Bold indicates the total values for all species combined.

The historical (pre-suppression) fire regime in lower mixed-conifer forests of the Sierra Nevada was characterised by frequent (fire-return intervals ranging from 5 to 32 years, Caprio and Swetnam 1995), relatively small (median area = 115 ha, Scholl and Taylor 2010), low- to moderate-severity fires (van Wagtendonk and Fites-Kaufman 2006). In the modern post-suppression era, fires in Sierra Nevada forests have become larger in overall size and have greater proportions of moderate and high severity (van Wagtendonk and Lutz 2007), and this trend is predicted to continue as winter snowpack declines (Lutz et al. 2009). Prior to 1900, the mean fire-return interval in the YFDP was 29.5 years, with the last widespread fire occurring in 1900 (Barth et al. 2015).

The YFDP was burned on 1–2 September 2013 in a management-ignited fire set to control the spread of the Rim Fire, a large wildfire that burned 104 131 ha of Stanislaus National Forest and Yosemite National Park (Fig. 1b; Lydersen et al. 2014). The ignition occurred in the afternoon in a mixed-forest–chaparral vegetation type atop the Crane Flat lookout ~1 km from the plot, and no management action was taken within the YFDP before or after ignition. The fire backed downslope through the western portion of the plot during the night of 1 September, and burned upslope through the remaining eastern portion the following day (Larson et al. 2016). The fire was active within the YFDP for ~30 h, although smouldering continued in some large coarse woody debris through November 2013 (T. Furniss and J. Lutz, pers. obs.). Fire intensity ranged from low- to high-intensity surface fire with some crown torching (based on Fire Behaviour Assessment Team cameras and thermocouples, Lutz et al. 2017a; Fig. 1c). Surface fuel consumption was 95% for litter, 93% for duff, 90% for 1 h fuels, 86% for 10 h fuels, 96% for 100 h fuels, and 61% for ≥1000-h fuels (Larson et al. 2016). Although the fire was management ignited, satellite-derived fire severity (Fig. 1c) was consistent with recent fires in the mixed-conifer zone of the Sierra Nevada (van Wagtendonk 2007; van Wagtendonk and Lutz 2007; Lutz et al. 2009), and the area that was management ignited was found to be indistinguishable from the wildfire-ignited area within the fire footprint (except for that portion that was plume-dominated) using a wide range of remote sensing techniques (Kane et al. 2015).

In May 2014, we revisited every stem in the YFDP and measured CVS, DBH, and live or dead status (hereafter ‘immediate fire mortality’). We also conducted full mortality censuses in the summers of 2015 and 2016 to measure delayed fire-related mortality (i.e. mortality of trees that survived >1 year post-fire). We considered ‘mortality’ as the aboveground death of individual stems, a method that overestimates mortality of some hardwood species that are fire-adapted to re-sprout post-fire. Considering this, we also recorded whether C. nuttallii and Q. kelloggii individuals were sprouting post-fire.

For each species, we extracted the corresponding logistic model forms and regression coefficients used in the FOFEM software (Lutes 2016). Species-specific models based on CVS were available for A. concolor, P. lambertiana, and C. decurrens (FOFEM_SP; Hood and Lutes 2017). Species-specific models were not available for Q. kelloggii and C. nuttallii, so we used the default model (FOFEM_RA; Ryan and Reinhardt 1988). The FOFEM_RA model uses two independent variables, CVS and bark thickness (BT), where BT is calculated according to the function DBH × V_sp, where V_sp is a species-specific coefficient (table 4.76 in Reinhardt and Crookston 2003).

We created locally-parameterised models based on the same independent variables and general model forms used in FOFEM, but with re-parameterised coefficients. For each species, we created 1st, 2nd, and 3rd order polynomial models with CVS as the single independent variable and 3-year post-fire status as the binomial response variable (as with FOFEM_SP), and we chose the model with the minimum Akaike information criterion (AIC) as our final ‘CVS’ model. For Q. kelloggii and C. nuttallii, we created optimised models with both BT and CVS as independent variables, but we eliminated BT from our final models because we attained a better fit using CVS alone. All models we created were generalised linear models with a logit link and the logistic model form:

where P_m is probability of mortality (within 3 years), e is the base of the natural logarithm, β₀ – β_t are regression coefficients, and X₁ – X_t are predictor variables (e.g. CVS, CVS²). We fit each model using the iteratively reweighted least-squares method of maximum likelihood estimation; model coefficients may be found in Table S1 (available as Supplementary material to this paper).

Upon preliminary analysis, we found that the third-order polynomial equations typically used in logistic mortality modelling often resulted in non-monotonic fits to our data (i.e. a tree with 25% CVS had a higher modelled probability of mortality, P_m, than a tree with 50% CVS; Fig. S2, available as Supplementary material to this paper). This may have been caused by the combination of two factors: tree DBH was negatively related to CVS (Fig. S1), and large-diameter trees can be more susceptible to fire damage compared with small trees (Kolb et al. 2007). The higher relative proportion of large-diameter trees at low CVS levels may have inflated P_m (20–30% CVS; Fig. S2), whereas the higher proportion of medium-sized trees may have reduced P_m for intermediate CVS levels (40–60% CVS; Fig. S2). This suggests that the dip in P_m is an artefact related to the interaction between CVS and DBH, and may not reflect the true physiological relationship between CVS and P_m for trees of the same size. Although there are multiple mechanisms through which fire damage may lead to mortality (e.g. reduced photosynthesis, Smith et al. 2016, 2017; increased susceptibility to beetles and pathogens, Parker et al. 2006; Hood and Bentz 2007), these mechanisms are all positively related to fire damage. In other words, a tree with more fire damage will have greater reduction to photosynthetic ability and be more susceptible to biotic mortality agents, and both of these indirect fire effects will increase probability of mortality. We resolved this problem in three different ways: we rarefied our dataset by CVS, we created additional CVS-based models that were constrained to be monotonic (i.e. no dip), and we created CVS-based models with a CVS : DBH interaction term. We rarefied our dataset by 20% CVS bin (10% bins used at high and low limits) by randomly selecting n trees from each CVS bin, where n is the minimum number of trees in any of the bins (n per species reported in Table 1). This rarefaction procedure reduced the severity of the dip, but it did not remove it entirely (Fig. S2). To reduce the dip completely, we developed an alternative CVS-based model form that was constrained to be monotonically increasing (slope was not held constant, but was always positive). This alternative model (CVS_mono) used CVS and CVS² terms to fit the lower range of CVS values, but replaced the CVS³ term with a higher-order polynomial (e.g. CVS⁸). Removing the CVS³ term removed the dip, and the high-order polynomial gave the curve an inflection point and a steep increase in predicted P_m values above the inflection. The power of the polynomial was chosen using an iterative procedure beginning with the lowest value necessary to obtain a monotonic fit (CVS⁴ or CVS⁵) and raising the power by one until minimum AIC was reached. As we increased the power of the polynomial, the inflection point moved to the right and the slope of the curve beyond the inflection point became steeper. We dropped the CVS² term if we obtained a better AIC without it. Although this approach is novel, it was the most tractable way could find to constrain a logistic model to be monotonically increasing. Full R code developed for this model fitting procedure in the Supplementary material.

The final method we used to address the interaction between CVS and DBH was to develop additional CVS-based models which included both DBH and CVS : DBH interaction terms (one interaction term per CVS polynomial level; CVS : DBH, CVS² : DBH, and CVS³ : DBH). In previously developed CVS-based mortality models, DBH has often been omitted because CVS and DBH can be co-linear, and dropping the DBH term may actually improve AIC (this was the case with FOFEM_SP; Hood and Lutes 2017). However, eliminating the DBH term ignores any possible interactions between CVS and DBH (e.g. Kolb et al. 2007; Engber and Varner 2012), and this may bias CVS-based models to be less accurate for large-diameter trees (which are less abundant).

To facilitate comparisons with previous studies, we created one final CVS model based only on trees ≥10 cm DBH (CVS_gt10). Lastly, we generated species-specific mortality models using pre-fire DBH as the single predictor. In total, we evaluated and compared six models for each species: FOFEM (FOFEM_RA for Q. kelloggii and C. nuttallii and FOFEM_SP for all other species), four locally-parameterised CVS-based models (CVS, CVS_mono, CVS + DBH and CVS_gt10) and a locally-parameterised DBH model (DBH). All CVS-based models were parameterised using the rarefied dataset.

We validated all models using 10-fold cross validation (sensu Regelbrugge and Conard 1993; Hood and Lutes 2017) to minimise the bias associated with parameterising and validating models using the same dataset (Kohavi 1995). We performed this procedure by first partitioning our data into 10 random groups. We then used nine groups to parameterise each model and generate predictions for the trees in the 10th group. This process was repeated 10 times, leaving out a different group each time. The resulting model prediction for each tree is therefore based on a model that was parameterised using a different subset of the data. Although models were parameterised using the rarefied dataset, we used the full dataset to validate each model.

We summarised model performance in three ways: individual tree-level accuracy (correctly classifying individual stems), stand-level accuracy (predicting the proportion of mortality per species and diameter class), and accuracy in predicting post-fire tree spatial pattern. A key distinction between these three levels is that individual tree metrics use the model response value (P_m) as a binomial classifier thus requiring a cut-off threshold (e.g. trees with P_m ≥0.5 are predicted to die whereas trees with P_m <0.5 are predicted to survive). Stand-level and spatial pattern accuracy may be summarised by using P_m as continuous probability value (e.g. for trees with P_m = 0.8, 80% will be predicted to die).

We assessed individual tree-level model accuracy by calculating sensitivity (correctly identified mortalities), specificity (correctly identified survival), overall accuracy (percentage correctly categorised), ΔAIC (AIC for each model – the lowest AIC value) and area under the Receiver Operating Characteristics (ROC) curve (AUC). We used a threshold P_m value of 0.5 to classify tress as live or dead.

Although individual tree-level metrics are widely used to evaluate mortality model accuracy (e.g. Grayson et al. 2017), these metrics may not be simply scaled up to evaluate stand-level performance. In other words, a model may have high individual-tree accuracy while simultaneously systematically over- or under-predicting mortality at the population or forest stand level. The strength of stand-level model accuracy assessment is that it permits the evaluation and comparison of mortality models in a way that describes observable changes in forest structure following fire. For example, a model may accurately predict small-diameter mortality while under-predicting large-diameter mortality (due to large-diameter trees having unique vulnerabilities to fire damage, Kolb et al. 2007). This model would have high individual-level accuracy because small-diameter trees are abundant and individual accuracy metrics weight all trees equally, but it would not be reliable for predicting mortality of large-diameter trees. Large-diameter trees represent an ecologically unique element of forest structure (Lutz et al. 2018b), and accurately predicting their fate following fire is important for improved predictions of aboveground biomass and carbon storage (sensu Lutz et al. 2017b).

We used individual-tree model response values to scale up our predictions to the population-level for each species. Although this analysis was based on populations (grouped by species), we refer to this scale as ‘stand level’ to maintain consistency with previous studies (sensu Hood et al. 2007). We assessed model accuracy by grouping trees by species and diameter and calculating percentage error as:

where N_model is the number of modelled mortalities, N_obs is the number of observed mortalities, and N_bin is the number of trees in each diameter class. This approach is similar to the method used by Hood et al. (2007) to assess stand-level accuracy, but we grouped stems by diameter class rather than by P_m (sensu Kane et al. 2017). This allows us to evaluate model performance across a range of diameter classes rather than across the range of P_m values. The number of dead trees predicted were identified for each model by assigning a status of either live or dead based on the continuous value of P_m for that individual stem (e.g. a stem with a P_m of 0.85 has an 85% chance of being identified as a dead tree), and these dead trees were summed for each bin.

We summarised stand-level accuracy across a range of spatial scales by sampling trees within randomly located quadrats within the YFDP ranging in area from 0.04 ha (20 × 20 m) to 2.56 ha (160 × 160 m). We sampled 10 quadrats of each size without replacement and calculated the stand-level error for all trees within each quadrat. There are two potential sources of this stand-level error: (1) small spatial scales include few trees, and there is more stochastic error because there is higher variability between sample means for smaller samples compared with larger samples, and (2) spatially correlated error due to the spatial autocorrelation of agents of delayed mortality. In other words, trees with the same level of fire damage may die at different rates if they are on a ridge v. in a drainage, or if they are in a beetle-kill patch v. a patch with no beetles. As mortality models predict mortality based only on individual tree properties, these neighbourhood-level factors affecting mortality contribute to model error.

To differentiate between these two sources of error, we used a null model designed to capture stochastic error (due to small sample size) but not spatially correlated error. We created this model by pairing each spatially explicit sample with a non-spatial sample of the same number of trees, randomly selected from the entire YFDP. This null model may be interpreted as the amount of error due to small sample sizes, and any excess error is attributable to spatially correlated processes that modify the probability of mortality for all trees within a given area.

We note that spatially correlated error is not the same as patchiness in fire effects: patchy fire behaviour is largely captured by tree-level metrics of fire damage such as CVS (resulting in spatially correlated mortality, but not necessarily model error), whereas spatially correlated error describes model error due primarily to the spatial signature of factors influencing delayed mortality.

We further scaled our assessment of model performance to the community level by evaluating live tree spatial patterns for the entire forest community (all species grouped). We assessed post-fire spatial pattern accuracy by quantifying and comparing pre-fire, post-fire, and modelled post-fire spatial patterns of trees within the YFDP. We used two methods: clump size analysis, a method of local pattern analysis often used to characterise spatial patterns in frequent-fire forest ecosystems (Plotkin et al. 2002; Larson and Churchill 2008); and spatial point pattern analysis, a method that has been used to infer ecological process in a variety of contexts (e.g. Lutz et al. 2014a; Larson et al. 2015; Furniss et al. 2017).

For the clump size analysis, we summarised spatial pattern as the number of clumps of trees as a function of clump sizes (i.e. number of clumps with 2 trees, 3 trees, 4 trees, etc.). We used an inter-tree distance threshold of 6 m to identify distinct clumps (i.e. clumps are groups of trees which are spaced no further apart than 6 m), as this was a found to be an optimal distance by Churchill et al. (2013). We compared each model by generating 99 simulations of mortality (where each tree was assigned a status of live or dead based on P_m) to obtain an estimate of the amount of variability in the modelled spatial pattern for each model.

For the point pattern analysis, we summarised spatial patterns with the pair correlation function (PCF), g(r), a point pattern summary statistic that describes the second-order characteristics of a spatial pattern across a range of scales (Wiegand and Moloney 2004). The PCF is defined as:

where K′(r) is the derivative of Ripley’s K function (Ripley 1976). To summarise observed patterns, we calculated g_pre(r) based on the pattern of trees that were alive pre-fire in 2013, and g_post(r) based on the pattern of trees that survived for at least 3 years post-fire (live as of 2016). We quantified model accuracy in predicting post-fire spatial patterns by generating 99 simulations of mortality for each model (where mortality was re-assigned for each simulation, as with the clump analysis), and g_model(r) was calculated for the set of ‘surviving’ trees for each of these simulations. We generated simulation envelopes (Baddeley et al. 2014) based on the minimum and maximum values from these simulations, and an estimate of the true value of g_model(r) based on the mean value from the simulations. We calculated change in pattern by g_post(r) – g_pre(r). We compared the observed and modelled patterns to the null model of complete spatial randomness (CSR; Wiegand and Moloney 2004) for the static patterns, g_model(r), and the null model of ‘no change’ for the differenced patterns, g_post(r) – g_pre(r).

For the multi-scale accuracy and spatial pattern analyses, we omitted two of the CVS-based models (CVS, CVS_mono, and CVS_gt10 all performed similarly; we retained the best one) to maintain interpretability of the figures. Analyses were conducted in R (ver. 3.4.1, R Foundation for Statistical Computing, Vienna, Austria, https://cran.r-project.org, accessed 7 November 2018) using the spatstat package (ver. 1.52–1; Baddeley et al. 2015; http://www.spatstat.org, accessed 7 November 2018).

Out of 34 174 pre-fire stems (1 cm ≤ DBH < 201 cm), there was immediate mortality of 24 151 stems and delayed mortality of 3818 stems (Table 1). The mean CVS for immediate mortalities was 99.7%, 66% for delayed mortalities, and 23% for surviving trees (Fig. S3, available as Supplementary material to this paper). The mean CVS for trees ≥1 cm DBH was 82%, 19% for trees ≥50 cm DBH, and 12% for trees ≥100 cm DBH. Mortality rates (3-year rates) by diameter class were 82% for trees ≥1 cm DBH, 63% for trees ≥50 cm DBH, and 28% for trees ≥100 cm DBH. Mortality rates for conifers were 84% for A. concolor, 77% for P. lambertiana, and 71% for C. decurrens (Table 1). Mortality rates for the hardwoods were 82% for C. nuttallii and 69% for Q. kelloggii, although 19% of top-killed C. nuttallii and 70% of top-killed Q. kelloggii were re-sprouting post-fire. The mean (and median) DBH was 7.4 cm (5.2 cm) for immediate mortality, 24.8 cm (18.1 cm) for delayed mortality, and 35.1 cm (25.9 cm) for surviving trees.

The locally-parameterised CVS-based models had the best AIC and higher overall accuracy compared with FOFEM and DBH-based models (Table 2). Adding a DBH interaction term (CVS + DBH model) improved AIC for A. concolor, but did not improve AIC for other species. The monotonic CVS model (CVS_mono) improved AIC for all conifers (Table 2). Although the AIC of the FOFEM models was worse than the AIC of the locally-parameterised CVS models, overall accuracy (i.e. total percentage correct) of FOFEM was similar for all three gymnosperms.

Table 2.  Comparison of individual and stand-level model performance for six different species-specific mortality models in predicting fire-related mortality
The six models include a diameter at breast height (DBH) model, three crown volume scorch models (3rd order polynomial, CVS; monotonic polynomial, CVS_mono; and based on trees ≥10 cm, CVS_gt10), one CVS model with a DBH interaction (CVS + DBH), and the models used in the First Order Fire Effects Models (FOFEM) software version 5.7+. Models were parameterised based on rarefied data (n in Table 1) and validated with 10-fold cross validation of the entire dataset. All model equations may be found in Table S1. We quantified model performance for individual trees with mortality percentage correct (Mort.; true positive rate, TPR), survival percentage correct (Surv.; true negative rate, TNR), overall accuracy (Acc.), Akaike information criterion (ΔAIC) and area under ROC curve (AUC). We quantified stand-level model performance using percentage error by diameter class (positive percentage error indicates over prediction of mortality, negative indicates under prediction). The top 5% includes the largest 5% of stems for each species (diameter thresholds for the top 5% may be found in Table 1)

Considering stand-level accuracy, the DBH model was the best model for small- to medium-diameter trees (<20 cm DBH) of most species (Table 2). For A. concolor, the DBH model was the most accurate for stand-level predictions for all diameter classes. For P. lambertiana, the FOFEM_SP model was best for stems <10 cm DBH, the CVS models were best for stems 10 ≤ DBH <100, and the DBH model was most accurate for stems ≥100 cm DBH. Adding the DBH interaction term improved model performance considerably for A. concolor, compared with CVS-based models without the interaction term. The DBH interaction term did not improve performance for most size classes of C. decurrens and P. lambertiana, but the CVS + DBH model did improve accuracy by 2% for large-diameter C. decurrens and P. lambertiana (>100 cm DBH; Table 2). In contrast, the CVS + DBH model decreased model accuracy for the largest stems of both hardwood species. The CVS_mono model was more accurate than the non-monotonic CVS model for stems <20 cm among all species, but the two CVS models were generally equivalent for larger diameter classes (Table 2, Fig. 2). The CVS models based only on stems ≥10 cm (CVS_gt10) were generally equivalent to the base CVS model for stems ≥20 cm DBH, but accuracy was worse for stems <20 cm DBH and all sizes of Q. kelloggii.

Fig. 2.  Probability of mortality (P_m) as a function of crown volume scorched (CVS). Dots represent observed proportion of stems that were killed in each CVS category (10% bins plus 0% and 100% bins), and lines represent species-specific logistic regression models using CVS as the independent variable and binary mortality status as the response. The dots are for visualisation purposes only; the models were parameterised on rarefied, unbinned data. For the CVS + DBH and FOFEM_RA models, shaded areas represent the full range of modelled P_m values (CVS + DBH coloured (light grey in greyscale version), FOFEM_RA in grey (dark grey in greyscale version)), whereas the dotted lines represent P_m while DBH was held constant (using average DBH).

Stand-level accuracy of the FOFEM models was within 15% of observed mortality for A. concolor, C. decurrens and P. lambertiana <40 cm DBH, but mortality was consistently under-predicted (Table 2, Fig. 2). FOFEM under-predicted large-diameter (≥100 cm DBH) P. lambertiana mortality by 17% and over-predicted mortality for Q. kelloggii and C. nuttallii, especially for stems ≥10 cm DBH (18–32%).

Stand-level error was negatively related to spatial scale (Fig. 3a–d). We observed this pattern for all models, as well as for each null model. The decreasing error in the null models indicates that some of the error was due to the small number of trees included in the spatial samples at small scales. As spatial scale increased, sample size grew and error in the null model was reduced. However, we also observed additional error in the spatially explicit samples that was greater than the error captured by the null models, indicating that there was spatially correlated error. The CVS-based models were the most accurate at small scales, predicting mortality within 10% of observed levels at all scales >0.1 ha. Model error for FOFEM and the DBH model fell within 10% at scales >1.4 and >1.2 ha respectively, but the DBH model had much greater error at scales <1 ha.

Fig. 3.  Stand-level accuracy as a function of spatial scale (a–d), and spatial pattern accuracy of mortality models (e–g). Coloured lines and shaded areas represent model predictions for four different mortality models. The y-axis in (a–d) represents the minimum and maximum values of stand-level error from 10 spatially explicit samples of size × within the Yosemite Forest Dynamics Plot (YFDP). The grey envelopes in (a–d) represent null models created by sampling the same number of trees as the spatially explicit sample at scale × contained, but each tree was chosen randomly. The lines in (e–g) represent the pre-fire (grey) and post-fire live spatial pattern (black), and modelled post-fire spatial patterns (colours (shades of grey in greyscale version)). The y-axis in E represents the number of clumps ha⁻¹ as a function of clump size. The y-axis in (f) is the value of the pair-correlation function, g(r), at a range of inter-tree distances, whereas the y-axis in (g) represents change in spatial pattern after the fire. Shaded areas represent a 95% confidence envelope around the model predictions, generated by 99 simulations of mortality. The grey dotted lines represent the expected value of g(r) under the null model of complete spatial randomness (f), or ‘no change’ (g). Shaded regions around the dotted lines represent variation expected under the null model. Values above the line indicate aggregation, and values below indicate hyper-dispersion.

The fire decreased the number of clumps from 3 to 1 ha⁻¹ for clumps of 5–20 trees, and mean clump size decreased from 10 to 7 trees per clump. All models accurately predicted this change (Fig. 3e), although the FOFEM model slightly overestimated the number of clumps at clump sizes ranging from 5 to 26 trees per clump. We were not able to differentiate the CVS- or DBH-based models from the observed pattern using the clump size analysis.

The spatial pattern of live stems within the YFDP became more aggregated after the fire (Fig. 3f). Both the CVS_mono and CVS + DBH models accurately predicted this change in pattern, although the magnitude of the observed increase in aggregation was greater than the models predicted, especially at small scales (0 to 2 m; Fig. 3g). The CVS_mono model had the greatest spatial pattern accuracy: the mean value of g(r) from the 99 simulations was the closest to the observed post-fire pattern, and this was the only model for which the observed value of g_post(r) – g_pre(r) fell within the model’s simulation envelope (at r >2 m; Fig. 3f, g). In contrast to the observed increase in aggregation post-fire, the DBH and FOFEM models both predicted reduced aggregation at small scales post-fire (Fig. 3g). Although FOFEM correctly predicted increased aggregation at spatial scales >1.5 m, the DBH model predicted increased dispersion at all scales 0 to 10 m.

Mortality of individual trees was predicted with a high degree of accuracy by all CVS-based models (CVS, CVS_mono, CVS_gt10, CVS + DBH, and FOFEM; Table 2). Indeed, CVS is a direct measure of physical damage incurred during the fire, and this damage is closely related to tree death (Peterson 1985; Weatherby et al. 1994; Sieg et al. 2006). This result is intuitive and expected, and it supports decades of previous work that have focussed on CVS as a key parameter in empirical mortality models (Ryan and Reinhardt 1988; Woolley et al. 2012; Hood and Lutes 2017; Grayson et al. 2017).

Immediate mortality was composed almost entirely of trees with 100% CVS, whereas delayed mortality was dominated by trees with lower CVS (Fig. S3). Mortality models based on 3-year mortality (such as FOFEM) are not well calibrated to predict immediate mortality, as they will over-predict mortality for trees with CVS <100% (Fig. S3). The FOFEM_SP models more closely reflected the observed patterns in delayed mortality for most species, as to be expected because they were parameterised based on 3-year mortality. We used 3 years as the cut-off for delayed mortality to maintain consistency with previous studies, although recent work has suggested that this timespan may not be sufficient to fully capture all delayed mortality (van Mantgem et al. 2011).

Although the FOFEM_SP models had high individual-tree accuracy, they under-predicted mortality at the stand level for A. concolor, C. decurrens, and P. lambertiana by 2 to 18% (Table 2, Fig. 2). A possible explanation for this is that the Sierra Nevada experienced a multi-year drought from 2012 to 2015 (Belmecheri et al. 2016), and this may have exaggerated mortality (e.g. Stephens et al. 2018). The extent to which fire-related mortality was affected by the drought is impossible to quantify for the YFDP, and a thorough analysis of the interactive effects between bark beetles, drought, and fire on tree mortality in the Sierra Nevada will require a multi-plot analysis that is beyond the scope of this study. Other studies have found that pre-fire climate appears to have a stronger influence on fire-related mortality than does post-fire climate (van Mantgem et al. 2013), suggesting that the timing of the Rim Fire early in the drought may have reduced the potential confound between drought- and fire-induced mortality.

This reveals an important question: will mortality models developed under non-drought conditions be useful in predicting mortality for fires that occur in a future climate characterised by greater drought stress? Although we cannot quantify the degree to which drought influenced mortality, we found that mortality models parameterised for past climates may under-predict mortality when fire events are coupled with multi-year drought. As climate warms and drought becomes more frequent, mortality models may need to be revised to reflect the modified relationship between fire damage and mortality probability under a drier climate. This may be done through re-parameterising models using fires that burned during a drought (as in the present study), or by explicitly including climatic variables as independent variables (more challenging, but more adaptable). Addressing this issue is a high priority for management-oriented modelling research, as the suitability of existing mortality models may decline as the co-occurrence of drought and fire becomes increasingly frequent.

In contrast to the under-prediction of mortality for conifers, the FOFEM_RA model greatly over-predicted stand-level mortality of C. nuttallii and Q. kelloggii by 5 to 32% (Table 2). Even for stems with 0% CVS, the FOFEM_RA model predicted a 72% chance of mortality for Q. kelloggii and 84% chance of mortality for C. nuttallii, but observed mortality was <10% for stems with 0% CVS for both species (Fig. 2). This is consistent with Kane et al. (2017) who found FOFEM over-predicted mortality for three hardwood species (including Q. kelloggii). Our observation that 70% of top-killed Q. kelloggii re-sprouted post-fire indicates that aboveground stem death may not be an appropriate measure of mortality for individuals of re-sprouting species. Re-sprouting individuals have an important effect on post-fire regeneration by rapidly sprouting new stems that grow much faster than conifer seedlings post-fire, and these sprouts maintain spatial pattern of the parent trees. Other studies have observed this discrepancy between actual mortality and top-kill (e.g. Catry et al. 2013), although most mortality models are primarily focussed on top-kill. This result underscores the need to develop more species-specific models for FOFEM, and the need to consider multiple types of mortality for hardwoods.

Surprisingly, the high individual-tree accuracy did not always translate into high stand-level accuracy, and this discrepancy highlights the importance of multiple tests of model performance. This was especially evident with the FOFEM models, which generally performed well at the individual tree level but had high stand-level error for some species and diameter classes. Although this has been considered elsewhere (e.g. Hood et al. 2007; Belote et al. 2015), individual-tree accuracy dominates the literature as the most widely used method for model parameterisation, validation, and comparison (e.g. Hood and Lutes 2017; Grayson et al. 2017). These results support the inclusion of stand-level accuracy metrics in future model development as this permits the explicit consideration of model accuracy for rare but important subpopulations such as large-diameter trees. The disparity between individual and stand-level accuracy should also be considered when using mortality models for stand-level inference and management decisions such as post-fire salvage harvest.

Pre-fire models based on DBH alone were often more accurate than the CVS-based models at the stand-level, most notably in predicting mortality of A. concolor and large-diameter P. lambertiana (Table 2). Although this may be partially explained by the tight correlation between CVS and DBH (Fig. S1), DBH should not perform better than CVS if it is primarily acting as a proxy for CVS. An alternative explanation is that there could be an interaction between CVS and DBH at large diameters that is not captured by using CVS alone (McHugh and Kolb 2003; Kolb et al. 2007). In other words, a small-diameter stem with 80% CVS may not have the same P_m as a large-diameter stem with 80% CVS. The presence of a CVS : DBH interaction also explains why stand-level accuracy was generally higher for the CVS + DBH model compared with CVS alone.

Interestingly, despite the improved stand-level performance of the CVS + DBH, AIC was best for CVS alone (all species except A. concolor). Standard model selection procedures are generally based only on AIC, and would therefore select the CVS model without a DBH term (sensu Hood and Lutes 2017). Although this maximises AIC, our results demonstrate AIC does not necessarily correlate with stand-level accuracy, especially for large-diameter trees. To maximise model accuracy across a broad range of tree diameters, we suggest evaluating models based on stand-level accuracy in addition to AIC, and including a DBH interaction term to account for the known interactive effects between CVS and DBH.

The high accuracy of the DBH model at the stand level suggests that post-fire data may not be necessary to develop an acceptable stand-level mortality model. However, as DBH does not implicitly capture fire behaviour as CVS does, site-specific fuel information (e.g. surface fuel loads) would have to be incorporated to calibrate the model for use in different forests. The resulting model would be similar to a CVS-based model used in a pre-fire context, which relies on a separate model to predict scorch height based on fuel levels. This would enable users to assess mortality as a function of DBH rather than CVS, and this may be useful for managers seeking to minimise (or maximise) mortality for specific diameter classes.

A problem we encountered with the DBH models is that standard logistic model forms did not represent the data well. First-, second-, and third-order polynomial models all predicted low P_m values (0 to 0.1%) at the large end of the diameter range for all conifer species, whereas the data show the true proportion of mortality was closer to 50% (Fig. S4, available as Supplementary material to this paper). We expect the true probability of mortality to be minimised for intermediate diameters and to increase for trees at the upper diameter range (e.g. McHugh and Kolb 2003). There are a few possible explanations for this: large-diameter trees often have cavities and scars from previous fires that can serve as an entry point for flames and embers and can lead to mechanical failure (7 out of 20 immediate mortalities >100 cm DBH were due to mechanical failure at a fire scar); they often have large-duff mounds around the base that can smoulder and cause prolonged soil heating at the root crown (Varner et al. 2009), and they may be more susceptible to beetle attack and drought post-fire (Hood and Bentz 2007; Kolb et al. 2007). Alternative approaches such as using constrained linear regression to predict the proportion of mortality (based on binned data) offer more control over the shape of the model response curve (compared with logistic regression), and this may result in models that more closely reflect the observed shape of the DBH–mortality response curve.

The standard third-order polynomial CVS models predicted individual-tree mortality with a high degree of accuracy, but the CVS_mono models improved individual and stand-level accuracy for all conifers in this study (Table 2). Graphical representations of these models demonstrate differences in the fundamental shape of these two model forms (Fig. 2, S2). The third-order polynomial models had a dip in P_m values in the middle of CVS range whereas the CVS_mono models had a steady increase (e.g. P. lambertiana) or a slight plateau (e.g. A. concolor). Although there may be a biological basis for a plateau in P_m at intermediate CVS levels, a negative relationship is biologically implausible. Fitting monotonic models solved this problem while still capturing the subtleties of the underlying data that would have been lost by using a simple linear formula (e.g. low slope at low CVS, steep slope at upper CVS). Additionally, CVS_mono models maintained high P_m for trees with 100% CVS whereas rarefied CVS models reduced P_m, causing the CVS models to under-predict mortality for trees with 100% CVS. Although the exact monotonic models we generated may not be broadly applicable to other sites as the position of the inflection point (determined by the power of the polynomial and the coefficients) may be highly susceptible to specific sites and fire events, we suggest a monotonic model form as an alternative if the common third-order polynomial model form results in a non-monotonic fit. Monotonic CVS models were more accurate at all scales, and they bear more relevance to the underlying processes determining mortality.

Stand-level accuracy was dependent on spatial scale for all models. We were able to partition this error into two categories: error due to small sample size (this error is represented by the null models; grey regions in Fig. 3a–d) and spatially correlated error (this error is represented by the coloured areas outside of the grey regions in Fig. 3a–d). We were not able to identify the mechanisms driving the spatially correlated error, but there are a few plausible explanations based on the spatially explicit nature of ecological processes that mediate delayed mortality. For example, imagine two stands that were burned identically; FOFEM will predict the same probability of mortality for each tree, and the same proportion of mortality at the stand level. If one of these stands experienced a bark beetle attack immediately post-fire whereas the other stand did not, the FOFEM mortality prediction will be too low for the stand that was attacked by beetles (and too high for the beetle-free stand), resulting in spatially correlated error. One could imagine a similar scenario based on other mechanisms that mediate delayed mortality, such as spatial variability of soil depth or moisture holding capacity, or patchily distributed pathogens. Mortality models do not yet include these mechanisms as input variables, but they implicitly consider them by predicting mortality based on average mortality rates. If a prediction is made for a large-enough area, these spatially correlated factors will be integrated together and the average mortality rates become a reasonable approximation. However, the scale at which this happens must be fairly large as we found spatially correlated error at all spatial scales up to the maximum scale of our analysis (2.5 ha).

Restoring structural heterogeneity and spatial patterns to reflect pre-fire suppression reference conditions is a frequent goal of prescribed fire restoration efforts (North et al. 2009), but the ability of mortality models to accurately predict change in spatial pattern remains unstudied. We found that both point pattern analysis and clump size analysis detected a change in pattern post-fire, but point pattern analysis was more sensitive to differences among modelled post-fire patterns. A likely explanation for this is that our choice of a 6-m threshold for the clump size analysis limits the scale of pattern detection to spatial scales >6 m, and the greatest pattern dynamics we observed using the point pattern analysis were at spatial scales 0 to 2 m. The clump size analysis may be better suited to larger-scale pattern dynamics, which all mortality models predicted reasonably well.

The point pattern analysis showed that the spatial pattern of live trees became more aggregated post-fire, and the CVS_mono model predicted this with the greatest accuracy (Fig. 3f, g). The FOFEM model and the DBH model were inaccurate at small spatial scales, predicting increased dispersion rather than increased aggregation post-fire. FOFEM accuracy improved at scales >1.5 m, whereas the DBH model was inaccurate at all scales 0–10 m.

The DBH model predicted a more dispersed (i.e. regular) post-fire pattern because the model predicted that large-diameter trees would survive and small-diameter trees would die (Fig. S4). The resulting spatial pattern is therefore closely based on the pre-fire spatial pattern of large-diameter trees, which was more regular than the pre-fire pattern of small-diameter trees. The CVS models, in contrast, use CVS, which is an indirect measure of fire intensity and it therefore reflects the spatial auto-correlation of fire behaviour and subsequent mortality. In other words, small scale patterns in fire intensity are driven by small scale heterogeneity in fuels, stand structure, and wind (Thaxton and Platt 2006; Hiers et al. 2009; Loudermilk et al. 2012); and measurements of CVS implicitly captured this patchiness. The mortality we predicted by the CVS models was therefore spatially structured, resulting in a predicted post-fire spatial pattern that closely reflects the observed post-fire pattern.

The FOFEM models had poor spatial pattern accuracy at small scales even though they are based on CVS. This was because the FOFEM_RA model greatly over-predicted mortality of the two hardwood species (15% for C. nuttallii, 19% for Q. kelloggii), and those species tend to be highly aggregated (Fig. S5a, available as Supplementary material to this paper). This resulted in a modelled post-fire spatial pattern that included very few C. nuttallii or Q. kelloggii and was therefore much more dispersed (Fig. S5b). The CVS-based models predicted hardwood survival more accurately, resulting in modelled spatial patterns that more closely reflected the observed post-fire pattern (Fig. 3f, g). This finding highlights the need to develop species-specific mortality models for hardwoods which will improve the spatial pattern accuracy of FOFEM.

The YFDP dataset is unique among the datasets used in previous mortality model validations, and many of the novel contributions of this study may be attributed to the representative nature of this dataset. As the population of each species and diameter class was a known quantity, we were able to quantify percentage error – a metric that easily translates to population-level predictions for a given forest stand. Although the data were numerically dominated by small- and medium-diameter trees, we were able to quantify percentage error across the full range of diameter classes, and this demonstrated the importance of using stand-level metrics in addition to individual-tree accuracy when comparing models. The spatial nature of the dataset allowed us to determine the scales at which there was spatially correlated error in model predictions, and to differentiate this from error due to small sample sizes. We were also able to use mortality models to predict change in spatial patterns, and we demonstrated how this can be used to compare models.

The large sample size in this study allowed us to conduct a detailed validation of FOFEMs for five species, and to evaluate logistic mortality models in a few novel ways. The FOFEM_SP models had high individual tree accuracy, but they systematically under-predicted mortality at the stand level. Some of this error is likely due to the timing of the Rim Fire in the middle of a severe, multi-year drought. We showed that existing mortality models may not be sufficient to accurately predict mortality under future climates, and this may require re-parameterising mortality models to capture the modified relationship between CVS and mortality under drought conditions. The models we developed in this study accomplish this for five species, although we recommend validating these models with independent data before they are applied in a management context. The FOFEM_RA model was inadequate for C. nuttallii and Q. kelloggii, and accuracy may be improved by 15 to 21% (individual accuracy, Table 2) by using the CVS models we created (Table S1). Incorporating these species-specific models into FOFEM will also improve community-level spatial pattern predictions.

We compared five different locally-parameterised models and found that there was no clear winner. Although CVS models had high individual-tree accuracy, DBH models were often more accurate at the stand level. Incorporating a DBH interaction term into CVS-based models improved stand-level performance, but the interaction term also introduced a considerable decline in accuracy for stems at the upper limit of the DBH range (e.g. A. concolor >100 cm DBH, Q. kelloggii >20 cm DBH, C. nuttallii >5 cm DBH). The CVS_mono models were generally the most accurate according to both AIC and stand-level metrics, but the monotonic model fitting procedure we developed created models with inflection points that may be less accurate when validated with independent data. These results emphasise the need to develop a stronger theoretical foundation for model forms; model parameters should be chosen based on known processes and interactions, and model forms should be constrained to biologically plausible relationships.

We found that stand-level accuracy improved with increasing spatial scale, although 2.5 ha was not large enough to eliminate spatially correlated error entirely. The spatial analysis revealed the importance of spatial processes to fire-related mortality, and it demonstrated that CVS captures some of this spatial variability. This study not only provides robust estimates of multi-scale mortality model accuracy for five species – it offers a re-evaluation of the fundamental approach (i.e. individual-level accuracy metrics and third-order polynomial model forms) that have dominated empirical mortality modelling for 30 years, and it reveals that existing models may be inadequate for fires that occur during drought. Logistic mortality models are well suited to individual tree predictions, but evaluating accuracy with multi-scale performance metrics and exploring alternative modelling approaches will enhance our ability to model population and community level fire effects.

The authors declare that there are no conflicts of interest.