Attributing a Causal Agent and Assessing the Severity of Non-Stand Replacing Disturbances in a Northern Hardwood Forest using Landsat-Derived Vegetation Indices

Abstract Non-stand-replacing disturbances are major drivers of northern hardwood forest dynamics, but are more challenging to characterize using satellite imagery than stand-replacing events. This study proposes a hurdle approach in which disturbance causal agents are first attributed to permanent sample plots that were either partially harvested, had sustained damage from an ice storm or remained undisturbed during the observation period, reaching an overall accuracy of 82.9%. Ordinary least square regression was then used to develop disturbance-specific models to assess the severity of partial harvests and damage from ice storms, with r-squared values of 0.57 and 0.59, respectively. The disturbance-specific models included a different set of predictors, confirming the importance of attributing a causal agent to a disturbance before assessing its severity. The sequence of models was implemented regionally to produce severity maps for two disturbance events, revealing within-stand variability in the severity that could be useful for the planning of future silvicultural actions. Although the proposed models offer acceptable performance, more research is needed to include additional disturbance agents and develop models that better capture the small variations in the spectral reflectance caused by low-severity disturbances, especially in the case of low-intensity partial harvests.

Les perturbations interm ediaires influencent fortement la dynamique des forêts de feuillues nordiques, mais sont plus difficiles a d etecter et a caract eriser a l'aide de l'imagerie satellitaire que les perturbations de forte intensit e. Cette etude propose une approche par etapes dans laquelle les agents causaux des perturbations sont d'abord attribu es a des placettesechantillon permanentes qui ont et e partiellement r ecolt ees, ont subi des dommages dus au verglas ou sont rest ees intactes pendant la p eriode d'observation avec une pr ecision globale de 82.9%.Des mod eles lin eaires distincts permettant l'estimation de la s ev erit e des coupes partielles et des dommages dus au verglas ont ensuite et e d evelopp es, dont les coefficients de d etermination respectifs sont de 0.57 et 0.59.Les mod eles sp ecifiques aux perturbations comprenaient un ensemble diff erent de pr edicteurs, confirmant l'importance d'attribuer un agent causal a une perturbation avant d'en evaluer la s ev erit e. Les mod eles ont et e appliqu es r egionalement pour estimer la s ev erit e de deux ev enements ponctuels, r ev elant une variabilit e dans la s ev erit e a une echelle inf erieure a celle du peuplement forestier.Cette information pourrait permettre une meilleure planification des actions sylvicoles futures.Bien que les mod eles propos es offrent une performance satisfaisante, des recherches suppl ementaires sont n ecessaires pour inclure davantage de types de perturbations interm ediaires et pour d evelopper des mod eles qui capturent mieux les variations subtiles de la s ev erit e des perturbations interm ediaires, en particulier dans le cas des coupes partielles de faible intensit e.

Introduction
Northern hardwood forests are among the most complex and biologically diverse ecosystems in Canada, providing some of its most valued wood products.Their dynamics are driven by disturbances caused by biotic and abiotic agents such as wind, freezing rain, pathogens and herbivory insects, resulting in canopy gaps of varying extent and a characteristic unevenaged structure.The return interval of such non-standreplacing disturbances in these forests is much shorter than that of severe, stand-replacing disturbances, which can reach thousands of years (Payette et al. 1990;Seymour et al. 2002).Although the impact of an individual disturbance is limited, the composite effect of these low-to moderate severity events at the ecosystem scale is a key driver of forest composition, structure and biomass dynamics (Payette et al. 1990;Woods 2000;Woods and Kern 2022).Northern hardwood forests are also subject to anthropogenic disturbances such as harvesting, which become the main driver of forest compositional changes in the most intensively managed regions (Danneyrolles et al. 2019).
Silvicultural systems in these forests generally rely on partial harvesting (Figure 1(a)), which involves removing between 15% and 50% of the basal area (BA, m 2 Áha À1 ) of a stand i.e., the sum of the cross-sectional areas of tree stems at 1.3 m in height over a unit area of land (Merrill 1935).The use of BA in silvicultural prescriptions has the dual advantage of avoiding the imprecision of volume estimates due to the complex architecture of broadleaved trees and being linked to the understory light regime within a stand.Manipulation of BA through partial harvesting helps maintain the uneven-aged structure of forest stands over time by replicating the natural disturbance regime of northern hardwood forests characterized by the dominance of non-stand replacing disturbances (Leak et al. 2014;Seymour et al. 2002).Such silvicultural systems, therefore, favor the regeneration of species with intermediate to high tolerance to shade (Grayson et al. 2012;Jenkins and Chambers 1989).
The projected increase in the frequency and severity of extreme climate events (Peng et al. 2011;Seidl et al. 2017) will inevitably affect the way northern hardwood forests are shaped by disturbances.The extent to which these changes will alter the ecosystem services they provide and their contribution to timber supply remains largely unknown (Gardiner and Moore 2014;Thom and Seidl 2016;Turner 2010).To better forecast these impacts and identify potential remedial actions, it is imperative to get a more accurate picture of the frequency, severity and extent of these more subtle disturbance events.This is usually not achievable from traditional forest inventory approaches, which rely on a limited number of ground plots that are not measured every year and that are not necessarily distributed representatively across the ecosystem (Bowman et al. 2013;Gillis et al. 2005).
Remote sensing from satellite imagery has become an indispensable tool to detect and characterize forest disturbances.Among available products, imagery from the Landsat program is now widely used for the monitoring of forest dynamics, providing continuous data since 1972 with a spatial resolution offering the possibility to study forest change both at global and regional scales, on short and long timescales (Banskota et al. 2014;Coppin and Bauer 1996;Wulder et al. 2022).The validity of conclusions drawn from such studies relies on the accuracy of the methods used to measure change.Detecting and assessing the severity of non-stand-replacing disturbances is much more challenging than it is for severe, stand-replacing disturbances, mostly because the changes they cause are of lower magnitude (Ahmed et al. 2017;Coops et al. 2020;Hermosilla et al. 2015a).
In broadleaved forests, Landsat has proven useful for studying non-stand replacing disturbances such as ice storms and partial harvests.Ice storms are severe events of freezing rain, causing the accumulation of heavy loads of ice on tree branches (Proulx and Greene 2001;Rhoads et al. 2002).This disturbance is common in eastern North America, where climate and physiography favor the development of the required combination of air and precipitation temperatures (Bennett 1959).In the region of Canada dominated by the northern hardwood forest ecosystem, ice accumulation from freezing rain occurs on an annual basis, influencing stand dynamics by causing variable levels of damage to tree crowns (Figure 1(a)), and causing mortality in moderately rare events where ice accumulation exceeds critical levels (Proulx and Greene 2001).In Eastern Ontario, Olthof et al. (2004) used measurements of damage to the canopy taken in 104 ground plots and concluded that changes in the Landsat-derived normalized difference vegetation index (NDVI) could not discriminate between many levels of damage, but could be useful to classify it in a limited number of categories.In the same region and for the same ice storm event, King et al. (2005) used a combination of pre-and post-disturbance values of the Landsat blue, red and shortwave infrared bands as well as biophysical variables to model canopy loss as a continuous variable with moderate success.In the White Mountains in the north-eastern United States, Burnett (2002) classified damage after an ice storm measured in 288 sites in three categories using NDVI and a ratio of shortwave and near-infrared Landsat bands, with accuracy from 78% to 82%.In the Adirondack forests of New York, Millward and Kraft (2004) found a significant and linear relationship between the change in NDVI values and in-situ measurements of damage to the canopy taken on a limited number of plots ( 16) after an ice storm.In the Appalachian mountains, Stueve et al. (2007) also found a statistically significant, but the nonlinear relationship between NVDI changes and the number of downed trees after an ice storm.
Partial harvesting has also been studied using Landsat data.Franklin et al. (2000) concluded that Landsat indices derived from the Tasseled Cap transformation, combining the information contained in six Landsat bands through orthogonal transformation (Crist and Cicone 1984), were moderately accurate for the detection of partial harvests in New Brunswick, while NDVI was easily confused in hardwoods because of the lush understory.In a subsequent study, Franklin et al. (2001) used the Tasseled Cap Wetness (TCW) to classify harvested stands into three categories of intensity.In Maine, Wilson and Sader (2002), Sader et al. (2003) and Jin and Sader (2005) used the normalized difference moisture index (NDMI), based on one of the Landsat shortwave infrared bands in addition to the near-infrared band, to detect and classify clearcut and partial harvests.They concluded that this index performed better than NDVI, and similarly to TCW in detecting harvests.When comparing different harvest practices in British Columbia, Jarron et al. (2017) found that the magnitude of change in the normalized burn ratio (NBR) consecutively to harvesting was representative of the intensity of the studied harvest practices.The NBR was also used by Tortini et al. (2019) to map partial harvests and clearcut in Michigan.Although the detection of partial harvests has received considerable attention, measuring the intensity of the cut has received considerably less attention.Yet, a more precise estimation of the change in harvested stands would be useful for forest managers because the actual basal area removal influences the post-harvest dynamics of the stand (Leak et al. 2014).
One of the main obstacles in studies dedicated to non-stand replacing disturbances is the availability of accurate data on the actual severity of the disturbances.Calibration and validation data are therefore frequently derived from indirect sources such as aerial surveys, high-resolution photos or other remote sensing products (e.g., Tortini et al. 2019).While such sources provide useful estimates of the extent of the damage over a large spatial coverage, they often comprise estimation errors that can limit the possibility to conduct analyses at a resolution that is useful for forest management (Johnson and Ross 2008).Field inventories provide the most accurate validation data but are uncommon because of the prohibitive collection costs and the need for measurements both before and after the disturbance.Yet, when available, such field-validated data allows a clear confirmation of the agent responsible for the disturbance and a precise quantification of the impacts of the disturbance on both the canopy and the sub-canopy layers (Rodman et al. 2021).
Studies on the detection and characterization of non-stand replacing disturbances using satellite imagery have resulted in approaches that rely on different Landsat bands and vegetation indices and have revealed distinct relationships between the change in spectral reflectance and the severity of the disturbances depending on the causal agent.An accurate assessment of the impacts of a disturbance, therefore, requires prior knowledge of its causal agent.Moreover, initial forest conditions are likely to influence the relationship between changes in spectral reflectance and the severity of the event (Harvey et al. 2019).Few studies have so far focused on classifying disturbances according to their causal agent, and those that did were either interested exclusively in standreplacing disturbances (e.g., Coops et al. 2020), were classifying non-stand replacing disturbances versus stand-replacing disturbances (e.g., Hermosilla et al. 2015a), or were conducted in boreal forests in which the initial conditions tend to be more uniform (e.g., Ahmed et al. 2017).
The objective of this study was to identify Landsatderived vegetation indices allowing the correct attribution of a causal agent to disturbances resulting from an ice storm and from partial harvests in northern hardwood forests.This study also aimed to identify indices that are appropriate to assess the severity of these disturbances.The influence of the pre-disturbance conditions of the stands on the relationships between changes in surface reflectance and the severity of the disturbance was also investigated, and case studies were produced in which the severity of the two types of disturbance was mapped on samples of the study area.

Study area
The study area is located in southern Quebec, Canada, within woodlots owned by Domtar Corporation (Figure 2).Totaling more than 160 000 hectares, the forest estate stretches through two bioclimatic subdomains i.e., the American basswood (Tilia americana L.) and the sugar maple (Acer saccharum Marsh.)-yellowbirch (Betula alleghaniensis Britt.)subdomains (Saucier et al. 2009).In the ecological districts where the woodlots are located, the mean annual temperature ranges from 2.5-5 C and the growing season length varies between 160 and 190 days.These districts are characterized by mean annual precipitation between 1000 and 1100 mm (Gosselin 2005(Gosselin , 2007)).Forests on mesic sites typically consist of uneven-aged stands, dominated by northern hardwood species, mainly sugar and red maples (Acer rubrum L.), yellow birch and American beech (Fagus grandifolia Ehrh.), with a smaller component of paper birch (Betula papyrifera Marsh.), white ash (Fraxinus americana L.), balsam fir (Abies balsamea [L.] Mill.) and red spruce (Picea rubens Sarg.).According to projections made using a thirteen global circulation model ensemble in climate NA (Wang et al. 2016), the maximum temperature in both summer and winter over the study area is expected to increase by 4.8 C and 4.2 C, respectively, based on an average of projections for shared socioeconomic pathways (SSP) 1-3 and 5 (Riahi et al. 2017).

Disturbance data
The forest inventory dataset used in this study comes from an extensive network of permanent sample plots established in hardwood-dominated stands between 1984 and 2010 and revisited periodically.Plots established from 1997 onwards are circular, with a radius of 11.28 m (400 m 2 ) while plots established before 1997 had the same area but are rectangular.All trees with a diameter at breast height (DBH) greater than 90 mm were measured at each survey.Newly dead trees and newly recruited trees that had reached a DBH greater than 90 mm were also recorded at each visit.From all plots of the network, 140 were selected to be included in the current study.A partial harvest had occurred in 46 of these plots.Silvicultural treatments included shelterwood, selection and salvage cuts ranging from 16% to 50% of basal area removal, as assessed from the post-harvest field inventory.Partial harvest operations were conducted in these plots in one of fourteen years distributed in the period from 2001 to 2020.The exact dates were determined using the execution reports submitted after the operations.A minimum of one and a maximum of eight sample plots were harvested in each of the fourteen years of observation.
In addition, 47 plots that were damaged during a major ice storm event were selected.From 5 January 1998 to 9 January 1998, an unprecedented ice storm, depositing as much as 100 mm of ice in the worst affected areas, hit eastern Ontario, southern Quebec, and the northeast United States, causing devastation in the northern hardwood forests (Chabot et al. 1998;Irland 1998;Proulx and Greene 2001).Damage to crowns was assessed visually in the field during the growing season following the ice storm, by estimating the percentage of branches that were broken on each individual tree within the plot.The canopy damage at the plot level was obtained by averaging the damage to the crowns of each tree, weighted by their respective basal area.The plot-averaged damage to the canopy ranged from 5% to 67%.
Finally, 47 control plots were randomly selected among the remaining plots of the network located in the same area, ensuring no disturbance was recorded in the surveys conducted before and after the assigned disturbance year.Fictive disturbance years were assigned to the control plots using a stratified approach creating a distribution of years similar to that of disturbed plots.A description of the pre-disturbance stand characteristics is presented in Table 1.

Development of the Landsat Surface Reflectance composites
Landsat TM, ETM þ and OLI Tier 1 Surface Reflectance scenes with less than 70% cloud cover, geometrically and atmospherically corrected, were used to create median composites for various years over each of the 140 selected plots.Only scenes acquired during the growing season in this region (i.e., June 1-August 30) were used.Pixels containing clouds or cloud shadows were masked and an additional visual inspection of the scenes was done to ensure the absence of undetected cloud shadows or thin clouds over the plots.A total of 246 scenes from Worldwide Reference System Paths 12 and 13, Rows 28 and 29 were used to create the composites, the acquisition dates ranging from 1996 to 2021 (Table 2).
A series of composite images were built separately for each plot, at different intervals depending on the date of the disturbance.Two distinct annual median composites were first produced for each plot from the scenes acquired in the two growing seasons preceding the disturbance.These composites were then averaged to produce a single pre-disturbance composite referred to as the reference composite.Scenes acquired in the growing season immediately following the disturbance were used to produce an additional composite, hereafter referred to as the disturbance composite.The ice storm, as well as most of the harvest operations, occurred between the second and third growing seasons because the harvest occurred in winter.For plots harvested during a growing season, scenes from the growing season corresponding to the harvest were attributed either to the reference or the disturbance composites depending on whether the image was captured before or after the date of harvest.A last composite was produced from the scenes acquired in the second growing season following the disturbance, which will be referred to as the post-disturbance composite.Between 1 and 18 scenes were used to build the reference, disturbance and post-disturbance composites depending on the plot and type of disturbance (Table 3).The Google Earth Engine platform (GEE, Gorelick et al. 2017) was used for scene extraction and compositing.

Calculation of vegetation indices and change metrics
Five vegetation indices were calculated from the median composites (Table 4).The NBR was originally developed to assess burn severity after wildfires (Key and Benson 2006), but has also proven useful to discriminate between different harvest intensities and detect a large variety of disturbances (Jarron et al. 2017;Kennedy et al. 2010).Two indices derived from the Tasseled Cap transformation, the TCW and TCA, were also included because of their known association with the density of the forest cover and their potential use in detecting forest removal and other disturbances (Franklin et al. 2000;Healey et al. 2006;Jin and Sader 2005).An index entirely based on the visible bands, the Green-Red Vegetation Index (GRVI) was also included in the study because of its sensitivity to phenology and subtle changes in the canopy of broadleaved forests (Motohka et al. 2010;Muraoka et al. 2013).The NDVI was also included as a reference, because of its extensive early studies on disturbance detection.
The average value of the vegetation indices was extracted from the composites over a 3 Â 3-pixel window enclosing the location of the plots, resulting in a short time series of three values per vegetation index per plot, representing the reference, disturbance and post-disturbance conditions, respectively.These values were used to compute two change metrics for each  vegetation index.The change in the index value due to the disturbance (Equation 1), hereafter referred to as the delta (d) was calculated as the difference between the value of the vegetation index in the growing season immediately following the disturbance (VI t ) and the reference index value (VI pre ).The recovery of the index value one year after disturbance (Equation 2), referred to as the recovery, was calculated by subtracting the value of the vegetation index taken from the disturbance composite (VI t ) to its value in the second growing season following the disturbance, taken from the post-disturbance composite (VI tþ1 ).
Recovery was included in the current study because the persistence of the changes in the canopy spectral reflectance has proven useful in discriminating between disturbance agents (Coops et al. 2020).

Modeling
Attribution of the disturbance agent Multinomial logistic regression (MLR) was used to investigate how the Landsat-derived metrics could allow discriminating between the effects of the two disturbances and differentiate stands that were disturbed from those that remained undisturbed.MLR is an extension of binary logistic regression used when the categorical dependent variable has more than two categories.The outcome from MLR is a probability associated with each class of the dependent variable (Jobson 2012).MLR has been used extensively to solve classification problems in environmental and remote sensing studies, for instance, to classify land cover types from satellite imagery (e.g., McRoberts 2011) soil types from topographic data (Debella-Gilo and Etzelm€ uller 2009) and crop disease status using a hyperspectral radiometer (Prabhakar et al. 2013).The assumptions of MLR do not include normality, linearity, or homoscedasticity, making it less restrictive than other classification methods used in remote sensing such as discriminant analysis.The results from MLR are also easy to interpret compared to the results of other machine learning approaches (Hogland et al. 2013).Penalized MLR was used in the current study because it allows removing the bias in the maximum likelihood estimates of the parameters when the number of outcome categories is large or, as in the current study, when the sample size is small (Bashir and Carter 2010;Firth 1993;Kosmidis and Firth 2011).In penalized MLR, a penalty function is applied to shrink the coefficients of the less contributive variables toward zero, reducing overfitting issues (de Jong et al.

2019).
A set of candidate models was built from the Landsat metrics calculated for the vegetation indices presented above.The reliability of penalized MLR remains affected by the number of events per variables (EPV), which represents the number of events in the smaller outcome group divided by the number of regression coefficients estimated.While no absolute recommendations have been formulated in the case of MLR, a conservative minimal number of EPV (i.e., 23) was maintained in the current study by including a maximum of two explanatory variables in the candidate models (Austin and Steyerberg 2017;de Jong et al. 2019).Only combinations of variables with a Pearson R coefficient under 0.7 were included in a given model.In the first group of candidate models, two delta metrics susceptible to be indicative of different changes in canopy condition were included, while a second group included both one delta and one recovery metric.Implementation of the models was done in the R programming environment using the 140 plots distributed in the three disturbance categories.The final model was identified after a modelselection procedure based on the Second-order Akaike Information Criterion (AICc) conducted using R package AICcmodavg (Mazerolle 2020).(Crist 1985;Powell et al. 2010) Assessment of the severity of the disturbance Examination of the data revealed a linear association between the dependent and independent variables.The assumptions of normality of the residuals and homoscedasticity were also fulfilled, it was deemed that ordinary least square regression was the most suitable and direct approach to create the disturbancespecific severity assessment models.Linear regression allowed modeling the severity of the disturbances as a continuous variable and the creation of easy-to-interpret models.A first group of candidate models was created for each disturbance type and included the delta metrics derived from either one or two vegetation indices, ensuring the Pearson R value remained under 0.7 when multiple predictors were included.To investigate the effect of the pre-disturbance stand condition on the relationship between disturbance severity and changes in the canopy spectral reflectance, we created a second group of candidate models that included the basal area of the plot prior to the disturbance as a covariate.Model selection based on AICc was conducted separately for both disturbance types but involved the same set of candidate models.
The number of available permanent sample plots was deemed insufficient to split the data into training and validation datasets.Therefore, all 47 plots that sustained damage from the 1998 ice storm were used in the development of the model dedicated to the assessment of the damage to the canopy caused by the ice storm, hereafter referred to as the ice storm model.Likewise, all 46 plots that were harvested were used in the development of the model dedicated to the assessment of the basal area removal after partial harvesting, hereafter referred to as the partial harvest model.The normality of the residuals was verified graphically and using the Shapiro-Wilk test, while the assumption of homoscedasticity was verified using the Breusch-Pagan Test.

Model application
To visualize the behavior of the disturbance-specific models and how they can serve to illustrate the variability in disturbance severity at both the landscape and stand levels, we mapped the predicted damage from disturbance events at specific locations within the study area.The selected subsets of our study area are located southwest of Lake Megantic in Southern Quebec (Figure 2(a)).This area was affected by the severe ice storm of 1998, and damage to the canopy ranging from light to severe was confirmed over the region by an aerial survey conducted in February 1998 (Chabot et al. 1998;Majcen et al. 1999).In the same region, we also selected six adjacent stands totaling 373.4 ha, in which partial harvest operations were conducted in 2020, after the end of the growing season (Figure 2(b)).
Landsat composites representing the reference and disturbance periods were created using the GEE implementation of the best available pixel (BAP) composite algorithm (White et al. 2014), allowing the creation of cloud-free, radiometrically and phenologically consistent surface reflectance composites, which are spatially contiguous over the targeted areas.As in the model-development step, two composites were first created for the two growing seasons preceding the disturbances.These composites were averaged to create a single composite representing the reference condition of the forest.The reference composite used to assess the severity of the damage caused by the ice storm of 1998, therefore, included pixels from the Landsat scenes acquired in the growing seasons of 1996 and 1997.The reference composite for the assessment of the basal area removal consecutively to the partial harvests included pixels from scenes acquired in the growing seasons of 2019 and 2020.Disturbance composites were created from scenes acquired in the growing season of 1998 for the ice storm-affected area, and 2021 for the harvested area.The vegetation indices and their associated delta metrics included as predictors in the ice storm and partial harvest models were then calculated from the reference and disturbance composites.
The ice storm model was applied to all broadleaved dominated stands of the area that had not been harvested in that same year.The sum of stands meeting these criteria had a combined area of 3974.1 ha in the first application subset.Using data from the forest inventory layer and a one-meter digital elevation model, we analyzed how the damage to the canopy predicted using the ice storm model varied as a function of biophysical characteristics such as the stand age, species composition, elevation and slope aspect.Analyses of variance and multiple comparisons with the Tukey honestly significant difference (HSD) test were used to investigate for statistically significant differences.The partial harvest model was applied to the six stands harvested in 2020, where the basal area removal levels ranged from 24% to 39% depending on the stand.Predictions from the model were averaged at the stand level, using polygons from the forest inventory layer, and compared to the post-harvest inventory data included in the geospatial layer made available by the landowner.

Attribution of the disturbance causal agent
Among the tested MLR candidate models for the attribution of the disturbance causal agent, the most accurate comprised the dNBR and dTCA as predictors (Table 5).The model yielded an overall accuracy of 82.9% (Table 6).The lowest producer's accuracy was obtained for undisturbed plots while the highest was obtained with the ice storm damage.A proportion of 22.0% of the undisturbed plots was misclassified as disturbed, while only 13.3% of the ice storm-damaged plots were misclassified either as undisturbed (11.1%) or harvested (2.2%).The producer's accuracy for the harvested plots lay in between, with 15.6% of misclassified plots, either as undisturbed (6.7%) or damaged by the ice storm (8.9%).The user's accuracy values were very similar for all three categories.A proportion of 17.0% of plots classified as undisturbed were actually disturbed plots.A same proportion of 17.0% of the plots classified as damaged by an ice storm was in reality harvested (8.5%) or undisturbed (8.5%) plots and 17.4% of the plots classified as harvested were, in reality, undisturbed plots (15.2%) or plots that had sustained some level of damage from the ice storm (2.2%).
The confusion matrix is presented in Table 7, and the individual effect of the predictors included in the final model is shown in Figure 3. Pixels exhibiting a high dNBR value, but only minor to no change in TCA were identified as damaged from an ice storm, while pixels exhibiting a moderate change in NBR but a large variation in TCA were identified as harvested.Pixels for which both dNBR and dTCA values were small were identified as undisturbed.

Ice storm
The best-fit model among the tested candidates for assessing the severity of the damage to the canopy after the ice storm included two delta metrics, i.e., the difference in Tasseled Cap Wetness (dTCW) and the difference in the Green-Red Vegetation Index (dGRVI).The selected model also included the   basal area of the plot prior to the disturbance as a covariate (Table 8).
The model reached an r-squared value of 0.59, with a root mean square error of 10.28.Estimates of the model parameters are shown in Table 9, while model predictions and unconditional 95% confidence intervals for the best-fit model parameters are shown in Figure 4.There was a positive relationship between the Landsat-derived metrics and the predicted damage to the canopy, and a higher basal area prior to the disturbance tended to be associated with lower damage for a given change in TCW and GRVI.High dTCW values combined with high dGRVI values and a low pre-disturbance basal area yielded the highest predicted damage to the canopy.The actual versus predicted damage to the canopy for the 47 plots used in model development is presented in Figure 5a.
Model applicationarea damaged by the ice storm of 1998 Although the basal area of the stand prior to the disturbance was among the predictors retained in the ice storm model, this information was not available in the historical forest inventory layers available for the study site.We, therefore, applied the predictions of the third-ranked model, which included only the two delta metrics from the best-fit model, and had a relatively small difference in AICc (Di ¼ 3.52) compared to the first-ranked model (Figure 6).The median of the predicted damage to the canopy over the area of application of our model was 26.3% and 79.9% of pixels (3175.7 ha) were deemed to have sustained over 10% of damage to the canopy.The predicted damage was lower than 50% in 90% of the area.
Within the model application area, young evenaged stands exhibited a significantly lower (p < 0.0001) average level of damage to the canopy than other stand types.Stands that were harvested either by clearcutting or strip cutting in the fifteen years prior to the ice storm had a predicted average damage to the canopy of 15.3% compared to 25.1% for other stand types.Similarly, the analysis of variance conducted over the area revealed that broadleaved-dominated mixed stands sustained significantly (p < 0.001) less damage than pure broadleaved stands during this disturbance event, with predicted proportions of damage to the canopy of 23.4 and 26.5%, respectively.There was also a significant difference (p < 0.0001) between the levels of damage predicted in stands at lower altitudes ( 500 m, 23.1%) and stands located at higher altitudes where predictions were higher (>500 m, 29.4%).The slope aspect also had a significant (p < 0.001) effect on the level of damage as stands sustained higher damage on slopes facing East or South (29.1% and 28.4%, respectively) compared to stands either on North (24.1%) or West (24.1%) facing slopes.Partial harvest Among the tested candidates, the best-fit model for assessing the basal area removal from partial harvest included a single predictor i.e., the change in NBR value (dNBR, Table 10).The r-squared value of the final partial harvest model was 0.57, with a root mean square error of 5.76.Estimates of the model parameters are shown in Table 11, and model predictions and unconditional 95% confidence intervals for the model parameters are presented in Figure 7.
There was a positive relationship between dNBR and the harvested proportion of the basal area; the removal of a greater proportion of the basal area resulted in larger dNBR values.The actual versus predicted basal area removal for the 46 plots used in model development is presented in Figure 5(b).

Model application -Area harvested in 2020
The selected partial harvest model was applied over the harvest area composed of six different stands that were harvested after the 2020 growing season (Figure 8).
The actual levels of basal area removal ranged from 24.0% to 39.0%, while the predicted removal levels ranged from 24.0% to 36.3% (Table 12).The absolute prediction error ranged from 0.7% to 16.6% for five stands out of six, but reached 38.5% in the case of stand E, for a median error of 12.6%.

Attribution of the disturbance causal agent
The two disturbances included in this study led to contrasting changes in the canopy spectral reflectance.Our results show that a combination of Landsatderived vegetation indices can help identify the causal agent of a disturbance.The final MLR classification model included two indices that are each associated with different changes within the canopy.The dNBR is increasingly used in disturbance detection algorithms, such as LandTrendr (Kennedy et al. 2010) and C2C (Hermosilla et al. 2015b) to detect both standreplacing and non-stand-replacing disturbances.For the disturbance events investigated in the current study, the average dNBR value was greater for plots damaged by the ice storm than for harvested plots.This observation highlights the importance of correctly identifying the nature of a disturbance before assessing its severity.While the damage caused to crowns by an ice storm can be considerable, its impact on key ecological attributes, such as the amount of living biomass, is not comparable to that of a partial harvest where trees are removed from the stand.The larger variation in NBR observed in plots affected by the ice storm may be explained by the fundamentally different nature of the two types of disturbance.In the case of partial harvest, operations are usually conducted through a limited number of logging trails in which the machinery operates.The basal area removal is mainly concentrated within these trails (Moreau et al. 2019), which can be occluded from airborne observations by the crowns of the residual trees.As the targeted basal area removal increases, more trees are going to be harvested   between the trails.The effect on the canopy spectral reflectance varies locally and may be relatively subtle at a 30 m resolution (Franklin et al. 2001).While the damage caused to the canopy after an ice storm is also variable, the freezing rain is distributed uniformly, at least at the plot scale, and hits the upper layer of the canopy first.Broken branches and the associated decrease in leaf area index (LAI) may therefore be captured more easily through airborne observations than the effect of partial harvesting.Despite this contrast between disturbance types, the range of dNBR values of plots that were affected by the ice storm (À0.01 dNBR 0.20) remained very similar to that of plots that were harvested (À0.01 dNBR 0.22).The inclusion of the dTCA in the model allowed a better differentiation of the respective effects of both disturbances.Indeed, the overlap in dTCA values between plots damaged by the ice storm (À0.04 dTCA 0.02) and harvested plots (À0.01 dTCA 0.14) was very small.The dTCA was significantly higher (p < 0.0001) in harvested stands compared to stands affected by the ice storm, while there were no significant differences (p ¼ 0.70) in dTCA between the latter and the undamaged plots.These results are in line with previous studies that used the TCA as an estimator of living biomass and forest  cover density (Ahmed et al. 2014;G omez et al. 2011;Powell et al. 2010).Although an ice storm causes damage to tree crowns, it does not affect the structure of a stand in the way a partial harvest does.

Disturbance-specific models
The Landsat-derived indices included in the final disturbance-specific models allowed assessing the severity of both types of disturbances in a continuous manner with a moderate proportion of the variance is accounted for.The partial harvest model relied on the dNBR as a single predictor.This index has been previously used to study partial harvesting in different ecosystems (Jarron et al. 2017;Tortini et al. 2019) and to monitor recovery following stand-replacing disturbances in boreal forests (Pickell et al. 2016;White et al. 2018).Our results suggest this index is appropriate to measure the intensity of partial harvests in northern hardwood forests.The NBR is a normalized difference of the Landsat near-infrared and second shortwave-infrared bands.Photosynthetically active vegetation shows a high reflectance in the near-infrared region compared to soil or non-photosynthetic vegetation (Clevers 1988).Soil and non-photosynthetic vegetation such as stems and branches show higher reflectance in the shortwave infrared regions compared to photosynthetically active vegetation, due to fewer water absorption features (Asner 1998).Partial harvest reduces the number of living trees and exposes bare soil and dry wooden debris, increasing the spectral reflectance in both the near-infrared and shortwave infrared regions of the spectrum, resulting in lower NBR as harvest proportion increases.It was challenging to assess the intensity of harvests with less than 30% of the basal area removed.At such low levels of harvest, the intensity was generally overestimated.This may be attributable to the presence of harvest trails in which all trees are removed to allow the machinery to circulate.In low-intensity partial harvests, the removal of trees within trails can be sufficient to reach the targeted basal area removal (Moreau et al. 2019).The canopy opening in logging trails may cause a substantial reduction in NBR values compared to the pre-harvest conditions, even at low levels of harvest.The additional removal of trees between the trails in harvests of higher intensity may produce more subtle variations in the index value.
The first-ranked model for the assessment of canopy damage due to the ice storm included two Landsat-derived indices i.e., the dTCW and dGRVI, in addition to the pre-disturbance basal area of the stand.The TCW was found to be a good indicator of the structural attributes and development stage in closed canopy forests (Cohen et al. 1995;Czerwinski et al. 2014;Jin and Sader 2005), and is also used in disturbance detection algorithms (e.g., LandTrendr, Kennedy et al. 2010).When exploring the relationships between the TCW and the structure of various forest stands, Hansen et al. (2001) found that it was strongly correlated to attributes such as crown closure, average crown diameter and structural complexity, all of which are likely to be altered in a stand damaged by an ice storm.
The GRVI has been used as an indicator of seasonal changes in foliage in broadleaved forests and is considered an indicator of the photosynthetic rate of the canopy (Motohka et al. 2010;Muraoka et al. 2013;Yin et al. 2022).Its relevance for assessing damage to crowns is in line with the effects of an ice storm on the spectral reflectance of the canopy.Forest canopies impacted by ice storms experience a reduction in LAI that can last for up to two years after the event, because of the physical damage caused to stems and branches (Rhoads et al. 2002;Zarnovican 2002).Photosynthetically active vegetation strongly absorbs visible light and shows particularly strong absorption features in the red region (Clevers 1988;Federer and Tanner 1966).Damage to the canopy results in a decrease in these strong absorption features, captured in the red spectral reflectance band, and as a consequent lowering of GRVI values.The relationship between changes in LAI and spectral reflectance of forest canopies in the red region has been observed in several studies involving both satellite and aerial imagery, some of which were conducted consecutively to ice storms (King et al. 2005;Nemani et al. 1993;Pellikka et al. 2000;Rautiainen 2005).

Influence of the pre-disturbance basal area
The basal area of the stands prior to the disturbance influenced the relationships between the severity of the disturbance and the Landsat-derived metrics only in the case of ice storm damage.For the same level of damage to the canopy measured on the ground, dTCW and dGRVI values were higher in stands with a higher pre-disturbance basal area.The steeper relationship between the delta indices and the damage to crowns could be attributable to a more complex vertical structure in stands of high basal area.During an ice storm, the larger trees, which form the top layer of the canopy, are generally the most damaged (King et al. 2005;Rhoads et al. 2002).Therefore, the nature of satellite-based vegetation indices, which are mainly focused on the upper layer in closed-canopy forests, could cause an overestimation of the damage, while trees in the lower layers, less visible from airborne observation, are likely to have sustained less damage.

Model application
Using disturbance-specific models to map the intensity of ice storm damage and partial harvesting over a subset of the study area revealed within-stand variability in both the damage to the canopy and harvest intensity.Consecutively to an ice storm, different levels of damage to the canopy may require specific silvicultural actions to restore the health and vigor of the stands and to account for the changes in the light regime caused by the disturbance.In severely damaged stands, salvage harvesting operations may be required to prevent losses due to mortality or degradation (Boulet et al. 2000;Bragg et al. 2003).Information on damage to the canopy at this resolution would therefore be useful in planning such interventions without the need to deploy a costly, time-consuming field sampling assessment.The application of the ice storm model over a large area allowed concluding significant correlations between biophysical variables and the level of damage to the canopy.In the case study dedicated to the ice storm event, young even-aged stands were less affected than older stands, which is in line with previous findings concluding that even-aged stands in northern hardwood forests are at low risk of damage from ice storms until they reach an age of 15-20 years (Pellikka et al. 2000;Rhoads et al. 2002).A comparison of the level of damage sustained by pure broadleaved stands compared to stands composed of a proportion of conifers revealed significant differences between the two groups, which is also consistent with findings from previous studies on the subject (Millward and Kraft 2004;Van Dyke 1999).The higher level of damage to the canopy for stands located at a higher elevation and the differences in damage between stands located on slopes facing different directions also agrees with what was found by other authors (Isaacs et al. 2014;King et al. 2005;Rhoads et al. 2002).
Regarding the application of the partial harvest model, the predicted basal area removal level was in general consistent with data from field inventories conducted consecutively to the harvest operations, with a median absolute error of 12.6%.Overall, the proportion of the basal area harvested tended to be slightly underestimated, although this was mainly the case for one stand in which the difference between the predicted and actual level of the harvest was quite large.An indication that the model allowed a representative visualization of the variability of basal area removal within and between stands comes from the fact that Figure 8 showed a clear distinction between the silvicultural treatments that were applied in the harvest area.The higher level of variability of harvest levels in stand A, for example, was the result of a form of salvage logging, which aimed to recover unhealthy stems in a depleted stand.This resulted in higher spatial variability in basal area removal compared to other stands, which were all subjected to selection cuts characterized by a more uniform basal area removal to favor shade-tolerant species (Leak et al. 2014;Nyland 1998).

Limitations and needs for future research
The study of non-stand replacing disturbances through remote sensing remains at its early stages, especially in broadleaved forests where most of the previous work has focused on a single disturbance agent.The current study proposes a hurdle approach in which disturbances are first detected and attributed to a causal agent using MLR before disturbancespecific linear models are used to estimate the severity of the disturbance as a continuous variable.Although the decision to model severity as a continuous variable has the advantage of revealing fine spatial variability in the severity of the disturbances, the proposed disturbance-specific models only explained a moderate proportion of the variance.Addressing some of the limitations raised in this section could help increase our ability to estimate more accurately the level of change caused to forest ecosystems by non-stand replacing disturbances and to attribute their causal agents.
The range in severity of the disturbance events included in the current study was dictated by the occurrence of such events during the study period and the availability of data from the permanent sample plots.Conclusions on the accuracy of the approach used to attribute the causal agent to the disturbances are therefore only applicable considering the specific distribution of severity among the disturbances included in the sample.Likewise, the establishment of a threshold for detectability of damage to the canopy or basal area removal using the MLR model falls beyond the scope of the current study.Additional work is needed to identify such a threshold and how it may vary between different types of disturbances.
The relevance of a severity threshold should also be investigated in the case of disturbance-specific models, especially for partial harvests, for which the intensity was challenging to assess accurately when the BA removal was low.
While only two types of non-stand replacing disturbances were included in the current study, northern hardwood forests are subjected to other types of small-to moderate-scale disturbances from various causal agents such as insects and wind.The detectability of these disturbances and the possibility to distinguish their effect on the spectral reflectance of the canopy using the limited set of vegetation indices tested in this study needs to be investigated.The success of such an undertaking is likely influenced by the time of the year when the disturbance occurs.In the context of the present study, this is especially relevant for the partial harvest.Operations conducted during the winter on a snow cover may preserve the understory vegetation, which will benefit from the opening of the canopy and proliferate during the following growing season.Changes to the spectral reflectance of the canopy consecutively to such harvests may be more subtle than for harvests conducted within the growing season, for which the bare soil is more likely to be exposed by the circulation of the machinery (Wolf et al. 2008).
Our results have also shown that, in the case of an ice storm, the condition of the stands prior to the disturbance must be considered in the estimation of the damage caused to the canopy.There is a need to further investigate how this affects the accuracy of methods used to study non-stand replacing disturbances using satellite imagery and to better understand how it may vary among causal agents and forest types.To that extent, the retrospective estimation of forest structural attributes including basal area using Landsat time series in combination with other remote sensing tools (e.g., Matasci et al. 2018), appears promising, since it could provide valuable information on initial stand conditions when there is a lack of historical field inventory data.

Conclusions
Results from the current study confirm the feasibility and relevance of using Landsat-derived vegetation indices to discriminate between two different types of non-stand replacing disturbances in northern hardwood forests and assess their severity.Using a combination of vegetation indices allowed capturing differences in the effects of two non-stand replacing disturbances on the spectral reflectance of the canopy.Our results confirmed the importance of identifying the causal agent before assessing the severity of a disturbance.For instance, the range of dNBR, one of the most commonly used metrics in disturbance detection, was very similar for plots that sustained damage after an ice storm and plots that were harvested.Yet, the effects of these two disturbances on the ecosystem dynamics are contrasting, both in their persistence and in severity.The different effects of ice storms and partial harvests on the spectral reflectance of the canopy have translated into the inclusion of a distinct set of predictors in the disturbance-specific severity assessment models.The application of these models over a subset of the study allowed for generating information on the spatial variability of the disturbances at a finer resolution than what is generally available to forest practitioners.This information could be useful in the planning and application of adapted silvicultural practices.

Figure 1 .
Figure 1.(a) Hardwood stand with soil exposed in the logging trails after a partial harvest.(b) Hardwood stands in which crowns were severely damaged consecutively by an ice storm.

Figure 2 .
Figure 2. Location of the study area and of the two subsets used for model application.(a) Area of broadleaved-dominated stands on which the severity of the damage resulting from the ice storm was mapped.(b) Harvest area composed of six stands harvested in 2020 on which the basal area removal consecutively to the partial harvests was mapped.

Figure 3 .
Figure 3.The individual effect of the predictors are included in the final disturbance classification model.

Figure 4 .
Figure 4. Model predictions and unconditional 95% confidence intervals for the best-fit model parameters for ice storm damage.

Figure 5 .
Figure 5. (a) Actual versus predicted damage to the canopy of the 47 plots included in the development of the ice storm model (b) Actual versus predicted basal area removal of the 46 plots included in the development of the partial harvest model.The solid line represents the 1:1 relationship.

Figure 6 .
Figure 6.Predicted damage to canopy (%) consequently to the ice storm of 1998 over a subset of the study area.

Figure 7 .
Figure 7. Model predictions and unconditional 95% confidence intervals for the best-fit model parameter for partial harvest.

Figure 8 .
Figure 8. Predicted basal area removal (%) over the area harvested in 2020 within the study site.

Table 1 .
Pre-disturbance characteristics of the plots included in the sample.

Table 2 .
The number of scenes from each Landsat satellite used to build the composites.

Table 3 .
The number of Landsat scenes used to build the composites for plots in each disturbance category and timing relative to the disturbance event.

Table 4 .
Vegetation indices, equations and references used in the current study.

Table 7 .
Confusion matrix after applying the multinomial classification model.

Table 5 .
Ranking of the candidate models for the assessment of ice storm damage.
AICc is the Akaike Information Criteria, corrected for small sample sizes, with deltaAICc (Di), model likelihood (Mk) and AICc weight (Wt i ).Values in bold indicate the model and associated explanatory variables that were retained after the model selection procedure

Table 6 .
Accuracy metrics and Kappa statistics for the selected multinomial classification model.

Table 8 .
Ranking of the candidate models for the assessment of ice storm damage.

Table 9 .
Estimates and standard error (SE) for the best-fit ice storm model parameters.

Table 10 .
Ranking of the candidate models for the assessment of basal area removal.AICc is the Akaike Information Criteria, corrected for small sample sizes, with deltaAICc (Di), model likelihood (Mk) and AICc weight (Wt i ).Values in bold indicate the model and associated explanatory variables that were retained after the model selection procedure.

Table 11 .
Estimates and standard error (SE) for the best-fit partial harvest model parameters.

Table 12 .
Actual and predicted basal area removal for the stands within the cut block (area) harvested in 2020.