A new approach for modeling delayed ﬁ re-induced tree mortality

. Global change is expanding the ecological niche of mixed-severity ﬁ re regimes into ecosystems that have not usually been associated with wild ﬁ res, such as temperate forests and rainforests. In contrast to stand-replacing ﬁ res, mixed-severity ﬁ res may result in delayed tree mortality driven by secondary factors such as post- ﬁ re environmental conditions. Because these effects vary as a function of time post- ﬁ re, their study using commonly applied logistic regression models is challenging. Here, we propose overcom-ing this challenge through the application of time-explicit survival models such as the Kaplan-Meier (KM-) estimator and the Cox proportional-hazards (PH-) model. We use data on tree mortality after mixed-sever-ity ﬁ res in beech forests to (1) illustrate temporal trends in the survival probabilities and the mortality hazard of beech, (2) estimate annual survival probabilities for different burn severities, and (3) consider driving factors with possible time-dependent effects. Based on our results, we argue that the combination of KM-estimator and Cox-PH models have the potential of substantially improve the analysis of delayed post-disturbance tree mortality by answering when and why tree mortality occurs. The results provide more speci ﬁ c information for implementing post- ﬁ re management measures.


INTRODUCTION
Climate change will modify survival probabilities of trees due to both changes in average climatic conditions and alterations in disturbance regimes (Allen et al. 2010, Seidl et al. 2017. The past decades illustrated that ongoing changes in climate and land-use may result in increasing burns across all forested biomes (van Lierop et al. 2015), including an expansion of mixedseverity fire regimes into ecosystems where fire is currently rare or absent (Adel et al. 2013, Adámek et al. 2015, Ascoli et al. 2015. In order to develop appropriate silvicultural rehabilitations and conservation measures in forest ecosystems where mixed-severity fires occur or will act as novel disturbance forced by climate change, understanding post-fire mortality processes and related factors is of paramount importance (Scott et al. 2002, Hood et al. 2018. Mixed-severity fires initiate different tree mortality trajectories according to the local burn intensity (Bond andKeeley 2005, Pausas andRibeiro 2017), resulting in spatially heterogeneous stand structures that influence forest recovery and resilience as well as future disturbance dynamics (Stephens et al. 2018). To this purpose, different models describing tree mortality probabilities and trajectories have been developed (for a review see Woolley et al. 2012, Hood et al. 2018), among which logistic regression models are the most commonly applied method.
Logistic regression models always refer to a precise event time point and return a dichotomized (dead/alive) response variable. Predictors are thus unified over a target time interval, potentially ignoring meaningful variation in the mortality process (Singer and Willett 1991) and ignoring possible changes in covariate values over time (Fornwalt et al. 2018). Thus, logistic regression models are well suited for predicting immediate or only slightly delayed tree mortality, which commonly occurs in fire-prone regions and in association with high-severity fires (Hood et al. 2010, Thies and Westlind 2012, Valor et al. 2017, Greyson et al. 2017, Roccaforte et al. 2018, Furniss et al. 2019). However, their dichotomized response variable is unsuited for predicting delayed tree mortality over decades.
Consequently, alternative approaches are needed to account for potential changes in the post-fire effects of secondary mortality factors over time and to improve our understanding of tree mortality associated with mixed-severity fires. The family of time-explicit survival models represents an alternative to logistic regression models by answering both when and why tree mortality occurs. Survival models analyze the time to event occurrence by considering both the event indicator (e.g., death of a tree) and the related timing from baseline (e.g., time since fire). In contrast to logistic regression models, the event is not dichotomized as dead or alive, rather as failure and censored (Fig. 1). Failure occurs when a fire-injured tree dies within the observation period. Censoring arises when the individual has not experienced the event (i.e., death) at the end of the follow-up sequences (time intervals between observations) or at the end of the observation period (right-censoring; see Fig. 1). Trees experiencing death at different time points are thus not merged over a given time interval, and changes in covariate values as time passes can be considered.
Survival analyses rely on various methods spanning from the non-parametric (e.g., the Kaplan-Meier estimator; Kaplan and Meier 1958) over the semi-parametric (e.g., Cox proportional hazards model;Cox 1995) to parametric models (e.g., Accelerated Failure Time Models). Originally developed for medical studies, survival models are becoming increasingly popular in forest science (Staupendahl and Zucchini 2010, Neuner et al. 2014, Brandl et al. 2020 1. Schematic representation of censoring and event happening in survival models. All trees enter the study at the time of fire (baseline) and observed until field assessment (years since fire). At an event occurring at time t observed during field assessment, all trees living equal or longer are integrated in the risk set for estimation. (Fox 2000), but have rarely been applied to describe a fire-induced delayed tree mortality and the related driving factors (Smith and Granger 2017). Since we know that European beech (Fagus sylvatica L.) displays delayed post-fire mortality over decades (up to 20 yr) depending on the burn severity (Maringer et al. 2016), we used this species to explore the suitability of survival models in predicting annual mortality considering secondary factors. Our specific questions are: 1. How does delayed post-fire tree mortality vary over time as a function of burn severity and environmental, climatic, and tree-related characteristics? 2. What are the main factors (predictors) influencing the delayed mortality process and how do their effects vary over time?
To tackle these questions, we use a two-step approach: We first apply the Kaplan-Meier estimator (KM-estimator) to assess the overall tree survival probabilities as a function of single potential mortality-influencing parameters (predictors). We then implement semi-parametric Cox proportional hazards models (Cox-PH model) to estimate the baseline hazards to die as well as the multiplicative impact of predictors on the post-fire tree survival probabilities. Since post-fire beech mortality differ with burn severity (Conedera et al. 2007, Maringer et al. 2016), we implemented three Cox proportional hazards models for different burn severities.

The study case
We sampled 27 beech forests (Fig. 2, Appendix S1: Table S1) across the European Alps, which had experienced a single surface fire of mixed severity in the last 20 yr. Criteria for site selection, data collection, variable assessment in the field, climate variables and data preparation followed the protocol by Maringer et al. (2016) and are described in detail in the supplementary material (Appendix S1). Generally, in the southern Alps wildfires are frequent and develop as surface fires, mostly occurring during the winter months when litter accumulates, grass vegetation is cured and the dry and warm wind (North foehn) drops the relative humidity below 20% (Valese et al. 2014). Fires usually start in the mixed deciduous forests (dominated by oak or chestnut) at lower elevation (below 900 m a.s.l.) and spread into the upper beech belt (900-1700 m a.s.l.). When winter drought conditions are combined with strong winds, extended forest fires may occur in beech stands (Pezzatti et al. 2009, Valese et al. 2014. In contrast, fire frequency is low in the northern Alps and burnt areas rarely exceed 1 ha (Conedera et al. 2018).

The fire ecology of beech
Since fires have historically rarely burnt in beech forests (e.g., Pezzatti et al. 2009), the species has no fire-adaptive traits. Beech does not develop heat-isolating thick barks to protect the vital tissue from lethal temperatures during a fire. Furthermore, it rapidly loses its resprouting capacity with age (Wagner et al. 2010, Packham et al. 2012). Indeed, beech is able to resprout after fire, but the resulting shoots tend to rapidly dieback and do not commonly result in a successful regeneration (van Gils et al. 2010, Maringer et al. 2012, Espelta et al. 2012. Post-fire regeneration in beech forests mostly relies on seed dispersal from surviving seed trees within and around the burn margins (Ascoli et al. 2015, Maringer et al. 2020).

Statistical approach
The family of survival analysis combine three main approaches: the non-parametric estimators, semi-parametric, and parametric models. In the present study, we used a two-step analysis flow, running first the Kaplan-Meier estimator (KMestimator), a non-parametric estimator, and in a second step the Cox proportional hazards model (Cox-PH model) as a semi-parametric model. We implemented the KM-estimator as a preliminary analysis exploring survival times with single variables and looking for possible time variation and significant differences between groups (see Table 1) in low-, moderate-, and high-severity burns, respectively (for the definition of burn severity see Appendix S2). Variables showing a significant (P < 0.05) effects were prioritize in the subsequent applied Cox-PH models (Hosmer et al. 2008). The multiplicative effect of predictors was then calculated using the semi-parametric Cox proportional hazards model (see Appendix S3: Fig. S1). In a last step, the KM-estimator was used again to validate the Cox-PH model (Brandl et al. 2020).

Kaplan-Meier estimator
Survival data are generally modeled as survival probability (S(t)) and mortality hazard (h (t)). The survival probability is the probability that an individual survives from the time of origin (e.g., the date of fire) to a time point (t) in the future (e.g., field assessment). The KM-estimator assumes no mathematical forms of the survival distribution. It multiplies together survival curves for intervals. Hence, it becomes a step function that estimates the probability (S 0 (t)) of  Beers et al. (1966). ‡ Recalculated to the year of fire based on the growth curves provided by Z' Graggen (1992), in case of dead lying trees we used the average diameter.
v www.esajournals.org 4 May 2021 v Volume 12(5) v Article e03458 not experiencing the event at time t according to following survival function: where n i is the number of trees at risk at time t i and d i is the number of trees that died during the period of reference. The KM-estimator thus describes the evolution of the survival probability as function of the time (e.g., years post-fire), which makes it useful for assessing changes in survival probabilities for different groups or treatments.
Since the KM-estimator can only test categorical variables, we divided continuous predictors into ranges below and above their median (Hosmer et al. 2008). Significant differences between two groups were determined by the non-parametric logrank test (Peto et al. 1977).

The Cox proportional hazards model
The Cox-PH model is a semi-parametric model that allows the quantification of predictors on the rate of event incidence (e.g., death) at a particular point in time (e.g., years post-fire). This rate is commonly referred to as the hazard rate (h i (t)that is the hazard rate for unit i at time t).
The Cox-PH model is expressed by the hazard function or force of mortality and can be interpreted as the risk that an event occurs. In our case, it calculates the probability of individual beeches to die after fire at a particular year postfire according to the following equation: hðtÞ ¼ h 0 ðtÞ þ e ðβ 1 x 1 þβ 2 x 2 þ⋯þβ n x n Þ with t representing the survival time, h 0 (t) is the baseline hazard corresponding to the value of the hazard if all the x 1,. . .,n are equal to zero (the quantity exp(0) = 1). The Cox-PH model provides a non-parametrical estimate of the baseline hazard function by assuming that the survival times do not follow any particular distribution (e.g., Weibull distribution). The regression coefficients, β 1,. . .,n , return the effect size of the covariates x 1,. . .,n on the probability of tree mortality. Cox-PH model regression coefficients are loghazard ratios. The exponential coefficients denote the relative change in the hazard of the occurrence of the event of interest (in our case fire-induced mortality) that is associated with a one-unit change of a particular predictor or the change of hazards between groups (e.g., when using a categorical variables). A hazard ratio greater (less) to one indicates that the related covariate is associated with an increasing (decreasing) hazard of death.
Data exploration for each sub-dataset followed the guidelines of Zuur et al. (2010), using Pearson's correlation coefficient and the variance inflation factor (VIF) to test collinearity among continuous variables and the chi-squared tests for the categorical ones.
Cox-PH models were fitted separately for low-, medium-, and high-severity burns. All three Cox-PH models were individual tree-based using the living status (failure/censored) together with post-fire year as response variable. After z-score transformation, single continuous variables were implemented in the Cox-PH models as linear and non-linear terms in order to test for non-linear effects (Keele 2010).
Based on the variable selection procedures as proposed by Glomb (2007), each Cox-PH model was first fitted for single explanatory variables separately. In a next step, we progressively added significant variables into the models until we obtained models with the lowest Akaike information criterion (AIC; Venables and Ripley 2010). Finally, the non-significant variables in the first step were added back in order to confirm or reject the lack of statistical significance. During this process, we additionally tested interactions among variables. The model fit has been assessed for all steps by comparing the AIC (Venables and Ripley 2010) of the nested models and their maximized log-likelihoods.
All statistical procedures were conducted using the statistic software R, version 3.3.3 (R Development Core Team 2014). For survival analysis, we used the survival package (Therneau 2019) and simPH package (Gandrud 2017).
The overall goodness-of-fit of the models were checked with the proportional hazards assumption (PHA) and residual analysis. Since survival analysis contains censored data, there is a different approach for calculating the residuals with respect to logistic regression analysis (Mills 2011). In particular, residuals should refer to the following four different parts of the Cox-PH model. v www.esajournals.org 1. The Cox-Snell residuals, which helps to assess the overall models fit and consists of a residual plot that follows a unit exponential distribution with a hazard ratio of 1 (Cox and Snell 1968). 2. The Schoenfeld residuals that test the fundamental Cox-PH models assumption of constancy of the hazard ratio over time, also known as the Cox proportional hazards assumption (PHA). In our specific case, the best model fit did not meet the PHA when referring to single post-fire years. Therefore, we organized the datasets into time intervals using the simPH package (Gandrud 2017). The underlying assumption when splitting the dataset into time intervals is that the hazard is constant within the time intervals, but can vary across them. Variables violating the PHA were considered as time-dependent and included with a time interaction (f(t)). The hazard rate for unit i with one-time interaction is then estimated based on following model equation: hðtÞ ¼ h 0 ðtÞ þ e ðβ 1 x 1 þβ 2 f ðtÞx 2 þ⋯þβ n x n Þ : 3 The score residuals (Klein and Moeschberger 2010) allow analysis of individual observations that have a large influence on the model. Therefore, score residuals are covariate specific for each observation and each covariate. A high absolute score residual means that the observation has a strong influence on the regression coefficient for the concerned covariate. 4 The Martingale residuals, which are used for evaluating the functional form of the model and consist of the representation of the residuals plotted against each model covariate.

Survival probabilities across burn severities
Comparing the observed survival probabilities by using both the KM-estimator and the logrank test confirmed that the survival probabilities differ significantly among burn severities at the 0.05-level (Fig. 3). The KM-estimator shows that in low-severity burns the survival probability is still 0.9 (SE AE 0.01) seven years post-fire and decreases slowly until it reaches 0.5 (SE AE 0.01) 16 yr post-fire. In moderate-severity burns, the survival probability is lower in the first 15 yr but the mean survival time is also reached at 16 yr post-fire (Fig. 3). In contrast, the survival probability rapidly decreases in high-severity burns, reaching 0.5 (SE AE 0.01) after 11 yr post-fire. During the following seven years (11-18 yr postfire), the survival probability steadily decreases and tends to zero after 18 yr post-fire.

Kaplan-Meier curves for single predictors
The KM-curves for single predictors show post-fire fungi infestation as a significant predictor for beech survival probabilities, indicating a higher mortality risk after fungi infestation (Fig. 4). Further, diameter at breast height (DBH) has a constant significant influence over time, revealing that large-sized trees have a higher probability to survive than small-sized ones regardless of the burn severity class (Appendix S4: Fig. S1). In case of moderate-burn severity, multi-stem beeches display a significant higher survival probability than single stem ones (Appendix S4: Fig. S2).
In addition to tree characteristics, post-fire climate variables also have a significant influence on the survival probabilities of beeches, when tested as single predictors. The logrank test shows that beeches have a significantly higher survival probability in warmer and wetter regions than in cooler and drier ones. This is true for moderate-and high-severity burns, but not for low-severity burns (Fig. 5, Appendix S4: Fig. S3). The influence of the lowest standardized precipitation evapotranspiration index (minSPEI) varies over time and differed significantly for moderate-and high-severity burns. Here, wetter years lead to a lower survival probability within the first decade post-fire, while the effect reversed in the subsequent decade (Appendix S4: Fig. S4).
Site characteristics, like aspect, altitude, and slope, influence the post-fire survival probabilities of beeches when testing the influence as a single predictor. The KM-curves indicate that in low-and high-severity burns, fire-injured beeches growing on south to southwestern exposition have significant higher survival probabilities than those on north to northeastern facing slopes (Appendix S4: Fig. S5). The effect of slope, in contrast, is significant for moderateseverity burn only. Here, tree survival probabilities are higher on steeper slopes (Appendix S4: Fig. S6). The logrank test for altitude indicates significantly lower survival probability with increasing elevation for all burn severity classes. The predictor evolves over the time since fire for all burn severity classes (Appendix S4: Fig. S7).

Concurring factors influencing beech's death
The best Cox-PH models, as indicated by the lowest AIC, include tree, site, and climate parameters for all burn severity classes (Table 2). By v www.esajournals.org holding all variables at their means, the best low-, moderate-, and high-severity models estimated survival probabilities of 0.95, 0.9, and 0.6 at 10 yr post-fire and 0.78, 0.7, and 0.3 at 15 yr post-fire, respectively (Fig. 6).
Tree characteristics such DBH, fungi infestation, and growth habit (mono-vs. polycormic trees) differ in their influence on beech mortality. Regardless of the burn severity, large-sized trees display a higher survival probability than smaller ones. In fact, for each increase in a DBH unit (cm), the hazard to die decreases by 6% (corresponding to a hazard ratio HR = 0.94), 10% (HR = 0.9), and 53% (HR = 0.47) in high-, moderate-, and low-severity models, respectively.
Beech infested by fungi in the post-fire period have a 3.6-times higher risk to die than without any fungal infestation when the burn severity is low to moderate, while according to the model the risk to die is only 84% higher in the highseverity burns (HR = 1.84; Table 2). Beech growth habit is significant for the moderateseverity model only, where it reveals a lower hazard to die for individuals growing as a multiple stem form (HR = 0.9).
Higher annual precipitation lowers the postfire hazard of beech to die in both moderate-and high-severity models, while the variable is not significant for the low-severity burns. Further, higher annual temperatures decrease the hazard to die in low-and moderate-severity burns, v www.esajournals.org whereas wetter springs and summers months (minSPEI) increase the hazard for beech to die in moderate-severity burns only.
Topographical parameters are less important predictors of mortality hazard in all models as revealed by the lower z-values. Aspect plays a significant role in case of low-severity fires, indicating a higher mortality hazard in association with northeastern exposure. Altitude is slightly significant in all severity models but has nearly no effect on changes in the hazard ratio (HR ≈ 1).

DISCUSSION
The survival approach for modeling delayed postfire tree mortality The KM-estimator and the Cox-PH model allowed us to answer questions regarding when and why post-fire delayed mortality occurs in beech forests. The temporal trends were determined by the KM-estimator, whereas the Cox-PH models tested the joint impact of multiple predictors, providing insights on the drivers of the post-fire mortality of beeches.
Similarly to clinical studies, applying survival models to delayed post-fire tree mortality implies that all subjects (trees/patients) have the same initial condition (pre-fire/before treatment) that may change after the application (fire/treatment). The lengths of the survival times are then measured from the initial stage to the event (death) or to the end of the study. However, in contrast to clinical studies, we used a retrospective approach as an alternative to long-term studies (Pickett 1989). Consequently, the time to event was not randomly selected from one target population as in classical medical follow-up studies. Rather, it was the result of the assemblage of wildfire areas that burnt in different years. Hence, all recorded trees were part of the target population, which entered the study at the year of fire (baseline 0, see Fig. 1).
Unfortunately, recent events (≤7 yr post-fire) were underrepresented (N = 34) in our dataset, and in the old burnt sites, trees that rapidly died after the fire may no longer be present due to the fast decay and decomposition rate of beech wood. Both factors may cause an overestimation of the survival probabilities, especially in moderate-and high-severity burns, where the mortality of fire-injured beeches within the first 7 yr after fire is

Notes:
Variables name + linear indicates that the predictor is time-dependent. For abbreviation of the variables see Table 1. Abbreviations are exp(ß), estimated hazard ratio (HR < 1 reduce the hazard to die, HR > 1 increase hazard to die, HR = 1 no changes); z-values as the number of standard errors between ß and 0. Ellipses indicate empty values. ***P < 0.001, **P < 0.01, *P < 0.05, •P < 0.1, n.s. = not significant.
usually higher with respect to low-severity burns (Maringer et al. 2016). Nevertheless, the survival approaches used were confirmed as a valuable method to gain insight on the survival probabilities in eventcaused tree mortality analysis in forest science (Staupendahl and Zucchini 2010, Griess et al. 2012, Neuner et al. 2014, Brandl et al. 2020, even when applied in retrospective studies.

Influence of tree characteristics
We used the KM-estimator to visualize temporal trends and associated violation of the proportional hazard assumption for single predictors (Hosmer et al. 2008) to highlight existing significant differences in the survival probabilities with respect to single parameters such as DBH, fungi infestation and, in case of moderate-severity burns, to growth habit. The results were confirmed by the Cox-PH models, which retained most of such predictors under consideration of their multiplicative effect.
Among variables included in the Cox-PH models, DBH has the strongest impact on tree's survival probabilities (indicated by the z-values), except for high-severity burns, while the relevance of the effect (hazard ratio) decreases faster in low-severity burns than in moderate-and high-severity ones. Low heat intensity during a fire results per definition in minimal (low-severity) effects on trees that mostly survive, while the resulting impact is conversely strong in highseverity burns (Dellasala 2018).
Generally, the relation of mortality as a function of DBH has been reported by several authors for other tree species (McHugh and Kolb 2003, Kobziar et al. 2006, Brando et al. 2012) as well as for beech (Shafiei et al. 2010, Maringer et al. 2016. Small-diameter trees are often burnt around their whole circumference stem, killing all of the vital tissue, while the same fire may only have a minor impact on large-sized trees since most of the vital tissue remains undamaged (Michaletz andJohnson 2006, Lawes et al. 2013). In addition, even if beech does not display marked fire resistance traits (see The fire ecology of beech), larger trees tend to have a slightly thicker bark and deeper root system than smaller individuals (Shekholeslami et al. 2011).
The interaction of individual shoots growing out of a stool (polycormic trees) with the fire front and the related flame and heat transfer into the cambium (Gutesell and Johnson 1996) also influences the survival probability in moderateseverity burns. The residence time of the fire is significantly longer on the leeward side of a stem or of a stool than on the windward side. This increases the heat exposure and lethal damage of the most leeward-sided shoot of a polycormic tree (Gutsell and Johnson 1996), concurrently lowering the impact on the shoots on the windward site. In low-and high-severity burns, the produced low and high heat intensity (Dellasala 2018) and the resulting high and low tree survival probability, respectively, might totally mask any possible effect of the polycormic structure.

Secondary stressors
The duration of heating and the related bark damage may directly affect beech survival by influencing the risk of secondary fungi infestation. Due to its thin bark, beech is known to be susceptible to secondary fungi infestation regardless of the burn severity (Conedera et al. 2007, Maringer et al. 2016. In case of the bark opening, fungi infestation starts within the first couple of years (Conedera et al. 2007), while the compartmentalization processes as a defense reaction last up to three years (Dujesiefken et al. 2005). During this period of defenselessness, wood decaying processes can lead to death.

Influence of climate
In addition to tree characteristics, site-related growing conditions also showed a significant influence on beech survival probabilities after fire. For instance, both the KM-estimator and the Cox-PH models indicated significant higher survival probabilities for beech experiencing moderate and high-severe fires when growing in regions with temperature and precipitation above the mean. Unfortunately, in our study case most of such site-related growing conditions are homogeneous or highly co-varying. For example, climate variables co-vary with geology, so that sites with calcareous bedrock have on average 900 mm less annual precipitation than sites on silicate bedrock. Hence, if beech is stressed during periods of drought on bedrock material with low water storage capacity (Gärtner et al. 2008), post-fire mortality might be also higher than under optimal growing conditions (van Mantgem et al. 2013). Consequently, in our specific case we cannot disentangle climate from other drivers (e.g., geologic and geomorphologic factors), although this reflects the dataset analyzed, rather than the overall suitability of the proposed modeling approach.

CONCLUSION
In our retrospective study, we used the survival analysis approach to model delayed (20 yr post-fire) fire-induced tree mortality by considering a broad combination of driving factors such as tree characteristics, climate, and geomorphological parameters.
With the help of the KM-estimator and the Cox-PH model, we illustrated temporal trends in the survival probabilities and the hazard of beech to die, respectively. In contrast to logistic regressions, the presented survival analyses have the advantage to (1) consider a time line (e.g., years post-fire) together with tree status (e.g., dead) as response variable, (2) estimate the survival probability for each time step, (3) include covariates that may vary over time, and (4) consider censored data. Based on the obtained results in this exploratory retrospective study, we are convinced that both the KM-estimator and Cox-PH models have the potential to substantially improve the modeling performances of delayed tree mortality after fire, thus providing much more specific information for implementing time-explicit restoration measures.

ACKNOWLEDGMENTS
This study was partially supported by the Swiss Federal Office for the Environment (FOEN grand number 00.0137.PZ / L424-1645). We thank Prof. Dippon (University of Stuttgart, Department of Statistic) for providing statistical support. Finally, we acknowledge the helpful comments of two anonymous reviewers.