Decomposing the Interactions between Fire Severity and Canopy Fuel Structure Using Multi-Temporal, Active, and Passive Remote Sensing Approaches

: Within the realms of both wildland and prescribed ﬁre, an understanding of how ﬁre severity and forest structure interact is critical for improving fuels treatment e ﬀ ectiveness, quantifying the ramiﬁcations of wildﬁres, and improving ﬁre behavior modeling. We integrated high resolution estimates of ﬁre severity with multi-temporal airborne laser scanning data to examine the role that various fuel loading, canopy shape, and other variables had on predicting ﬁre severity for a complex of prescribed ﬁres and one wildﬁre and how three-dimensional fuels changed as a result of these ﬁres. Fuel loading characteristics were widely variable, and ﬁres were ignited using a several techniques (heading, ﬂanking, and backing), leading to a large amount of variability in ﬁre behavior and subsequent ﬁre e ﬀ ects. Through our analysis, we found that ﬁre severity was linked explicitly to pre-ﬁre fuel loading and structure, particularly in the three-dimensional distribution of fuels. Fire severity was also correlated with post-ﬁre fuel loading, forest structural heterogeneity, and shifted the diversity and abundance of canopy classes within the landscape. This work demonstrates that the vertical distribution of fuel is an important factor and that subtle di ﬀ erence has deﬁned e ﬀ ects on ﬁre behavior and severity. CHP class 12 had a positive inﬂuence in the model and had a bi-modal distribution of fuels that were nearly equal in the canopy and mid-canopy strata. Our results indicate that the vertical distribution of fuel is an important factor and that subtle di ﬀ erences have deﬁned e ﬀ ects on ﬁre behavior.


Introduction
Canopy structure is the vertical and horizontal distribution of plant material in a forested ecosystem and is a driver of many ecosystem functions and processes. Canopy structure influences the distribution of photosynthetic material and its efficiency [1], microclimate [2], stand species dynamics [3], energy balances [4], and species habitat and behavior [5]. Similarly, forest fluid dynamics are closely tied to forest structure, including atmospheric flow within forests [6] as well as the redistribution of intercepted precipitation [7].
In the case of wildland fire, canopy structure influences the fire dynamics directly as fuel and indirectly through its influence on other variables in the fire environment. For example, three-dimensional fuel structure has been shown to have an influence on fuel moisture regimes, with the forest floor being consistently more moist under dense canopies [8,9]. Structure also influences the fire environment, or factors affecting the dynamics of the fire, by influencing wind and energy flow through drag [10] and energy absorption [11]. Likewise, models of smoke dispersion through variable canopy architecture allow the prediction of potential smoke dispersion patterns and thus increase knowledge of the degree to which canopy architecture influences forest level air parcel movement [12,13]. As an assemblage Fire 2020, 3, 7 2 of 21 of fuel elements, the canopy directly contributes to the dynamics of wildland fires, including rate of spread and energy release [14,15]. Canopy fuels are thus not only vital to understanding fire spread and behavior [16], they are also inputs in both empirical [17] and computational fluid dynamics (CFD; [18]) models of wildland fire behavior.
Fire severity is a measure of the change in the properties of the surface and vegetation resulting from a fire event [19]. Fire severity varies greatly within and among burns in both vertical and horizontal dimensions, resulting in a complex landscape-wide mosaic of forest structure and composition [20]. Integrating measurements of forest structure from before and after a fire incorporates the full spatial dimensionality of terrestrial fire effects (note that some first and second order fire effects may occur below ground) that are directly relatable to management needs and ecological function in forests.
Field-based estimates of fire severity such as the Composite Burn Index (CBI; [21]) are often qualitative, relying on ocular estimates of "change" in various fuel strata, e.g., forest floor, shrubs, and canopies. A major confounding factor in CBI as a metric is that it does not consider pre-fire data, which can lead to the misattribution of an effect that was not actually due to the fire [22]. For example, McCarley et al. [23] illustrated the importance of considering pre-fire canopy structures with an example where the interaction of the effects of mountain pine beetle and severe wildfire reduced canopy cover more than severe wildfire alone. In this case, without the benefit of pre-fire data, CBI would have attributed these effects to the fire. Additionally, while CBI can be useful for estimating fire effects statistics for units on a landscape, they are poorly resolved spatially, can suffer from subjectivity, and have limitations for use in fire modeling or examining landscape pattern [24].
Fire severity can be effectively indexed using remote sensed spectral reflectance data at large scales and increasingly high resolutions [25][26][27]. The advantages of remote sensing analyses include the potential to map the entire spatial extent of a fire or region, and the potential to incorporate both before-fire and after-fire information. Furthermore, the increased availability of multi-temporal datasets has led to an increased understanding of spatial patterning and recovery through disturbance cycles. In many cases, studies have focused on post-fire landscape matrices by examining two-dimensional estimations of fire-severity, mostly derived from Landsat datasets [28,29]. These data are useful for operational purposes, chiefly to assess the spatial need for post-fire mitigation activities [24]. However, these spectral reflectance-derived indices are only useful for developing a two-dimensional understanding of severity analogs and are not suited for measuring forest structural variables in three-dimensions. Problematically, these estimates only provide an index of fire severity, as they are not a direct measurement of the changes in the fuel strata [30].
Characterization of forest structure, and the changes following fire, would therefore greatly enhance the understanding of fire severity. However, the field-based estimation of forest structural attributes can only be quantified in the field at very limited scales. First, field assessments of post-fire environments can be dangerous and difficult to access, limiting the potential spatial resolution and investment toward field collected data. Additionally, direct measures of structural attributes are challenging because of the difficulty in accessing the canopy. Thus, canopy structure attributes are typically inferred by models using measurements of stem height, stem diameter, and crown base height as inputs for species-specific models of canopy fuels distributions [31]. This approach is problematic in the context of fire effects on canopies. Pre-determined distributions typically used to populate canopy fuels can be insensitive to the heterogeneity typical of canopy structure. More importantly, in the case of non-stand replacing fires, stem height, stem diameter, and crown base height do not change appreciably even after changes in canopy architecture have occurred. This is particularly problematic in stands composed of species that resprout epicormically, where foliar elements of the canopy are completely consumed and then can grow back from the same stem. However, these plot-based estimates do provide vital input parameters for the development of models that scale measurements using remote sensing datasets.
Forest structural attributes can now be quantified using Light Detection and Ranging (LiDAR) sensors, which make three-dimensional measurements of the forest that can be analyzed in ways that are directly applicable to wildland fire science [32][33][34]. Using LiDAR, it is now possible to represent forest structure at the plot scale using Terrestrial Laser Scanning (TLS), at the landscape scale using Airborne Laser Scanning (ALS), and at the global scale using space-born laser sensors. Because ALS data can be collected at a scale and resolution that can encompass large fire perimeters, it is particularly well suited for providing structural fire effects data. Analysis approaches that leverage the full potential of these datasets can offer new insight into the alteration of forest function and processes due to disturbances. Voxel-based approaches for binning ALS returns into discrete height categories can preserve the structural information present in the ALS dataset that is otherwise lost in approaches that only make use of distributional statistics [35]. This information presents a direct representation of the horizontal and vertical distribution of fuels. Given that these attributes of structure are vital for improving our understanding of many aspects of forest ecology, particularly the behavior and effects of wildland fire, it is important to develop and demonstrate new analytical methods to extract more informative structural information about forest vegetation matrices. Additionally, pre-and post-fire LiDAR data has been shown to accurately estimate forest structural changes across landscapes in several cases. Exploring fire severity in three continuous dimensions via remotely sensed forest structure datasets as a change from a pre-fire condition, has received limited previous attention; however, such an approach is posed to facilitate robust analyses of pyrogenic landscapes in ways that expand our understanding of the dynamics of fire effects on forest structure in the context of both wildland and prescribed fire, as well as the ecological nuances of fire regimes.
The objectives of this study were: (1) to investigate whether pre-fire canopy structure influenced fire severity; (2) to examine how vertical structural heterogeneity was modified across a gradient of fire severities; and (3) to describe landscape-scale changes in spatial patterning of forest structure resulting from fire. To work towards these objectives, we decomposed the influence of pre-fire forest structure characteristics on fire severity and subsequent changes in post-fire canopy structure using a combination of active and passive remote sensing datasets. We simplified pre-and post-ALS datasets using a classification strategy that more effectively preserves information about three-dimensional forest structure heterogeneity than current commonly used approaches. Next, we used a previously published map of burn severity to stratify the classified ALS data and allow us to identify trends and patterns in these datasets. The intersection of the pre-fire fuel structure, fire severity, and post-fire fuel structure data products allowed us to directly address the objectives set forth above.

Study Area
The New Jersey Pinelands National Reserve (PNR) is a largely contiguous, fire-prone, landscape situated within 60 km of both New York City and Philadelphia. Despite the area's urban surroundings, the 400,000 ha PNR has remained densely forested due to its conservation status as a UNESCO world heritage site and its vast tracts, which are owned by the public and conservation groups. These forests are mostly comprised of Pinus rigida Mill. (pitch pine) or Quercus spp. overstories that are persistent across the region with a gradient of composition ranging from pure pine to pure oak stands. It should be noted that pitch pine resprouts epicormically, and that even in severe fires, there can be limited overstory mortality of larger stems [36]. The mid-story of these forests is made up of a variety of shrub oak species (e.g., Quercus ilicifolia Wangenh., Quercus marilandica Münchh.), mountain laurel (Kalmia latifolia L.), and suppressed overstory species as well. The understory of these stands is comprised mostly of Ericacious shrubs (e.g., Gaylussacia spp. and Vaccinium spp.), sheep laurel (Kalmia angustifolia L.), and a layer of fast-decomposing leaf litter.
Historically, the PNR has maintained a high frequency of wildfire, which peaked at the turn of the 20th century, and has subsequently been steadily decreasing [37]. Currently an average of approximately 3400 ha yr −1 is burnt. These wildfires are typically fast-moving, high intensity fires, which occasionally grow to reach over 10,000 ha. In addition to wildfires, there is a history of prescribed fire across a large proportion of this landscape, beginning in the modern era in the late 1930s. The current treatment regime averages approximately 6300 ha yr −1 of dormant season fire [38]. Typically, these fires are low-intensity and do not consume more than forest floor fuels and some portion of ladder fuels. However, occasional high intensity patches help to create a mosaic of post-fire conditions [14].
This study focuses on an area within the PNR that was burned in several prescribed fires and a single wildfire, with a range of meteorological conditions and firing techniques over a period of several days in February and March of 2016. The variation in meteorological conditions, fuels, and firing techniques for these burns resulted in a mosaic of fire severities over a total of 1112 ha [27]. The extents of the 16 burn units are depicted in Figure 1. Over the past 100 years, the entirety of this area was burned in at least one wildfire in 1929, 1943, 1946, 1964, or 1982, or in a prescribed fire in 2000 [39]. The study area represents a mixture of upland and lowland pitch pine-dominated forests, with some lowland Atlantic white cedar (Chamaecyparis thyoides (L.) Britton, Sterns, and Poggenb.) stands, and red maple (Acer rubrum L.) swamps which were relatively untouched by the wildfire or prescribed fires of 2016. In the north of the unit, there is an area of a genetically unique population of dwarf pitch pines, known as the Spring Hill Plains [40]. approximately 3400 ha yr −1 is burnt. These wildfires are typically fast-moving, high intensity fires, which occasionally grow to reach over 10,000 ha. In addition to wildfires, there is a history of prescribed fire across a large proportion of this landscape, beginning in the modern era in the late 1930s. The current treatment regime averages approximately 6300 ha yr −1 of dormant season fire [38]. Typically, these fires are low-intensity and do not consume more than forest floor fuels and some portion of ladder fuels. However, occasional high intensity patches help to create a mosaic of postfire conditions [14]. This study focuses on an area within the PNR that was burned in several prescribed fires and a single wildfire, with a range of meteorological conditions and firing techniques over a period of several days in February and March of 2016. The variation in meteorological conditions, fuels, and firing techniques for these burns resulted in a mosaic of fire severities over a total of 1112 ha [27]. The extents of the 16 burn units are depicted in Figure 1. Over the past 100 years, the entirety of this area was burned in at least one wildfire in 1929, 1943, 1946, 1964, or 1982, or in a prescribed fire in 2000 [39]. The study area represents a mixture of upland and lowland pitch pine-dominated forests, with some lowland Atlantic white cedar (Chamaecyparis thyoides (L.) Britton, Sterns, and Poggenb.) stands, and red maple (Acer rubrum L.) swamps which were relatively untouched by the wildfire or prescribed fires of 2016. In the north of the unit, there is an area of a genetically unique population of dwarf pitch pines, known as the Spring Hill Plains [40].

Spectral Reflectance Data
In a previous paper [27], we demonstrated that the properties of the WorldView-3 (WV-3) sensor allow for the estimation of burn severity from spectral reflectance in much the same way as burn severity is calculated from Landsat data [41], but at much higher spatial and spectral resolutions. This finer level (7.5 m) of detail minimizes spectral mixing and allows for better exploitation of the data when combined with higher-resolution products, such as from ALS datasets.
WV-3 data were collected before and after the fire event to provide high resolution, spatially explicit, estimates of fire severity within the area of interest (AOI). Imagery was collected on 8 December 2016 (pre-fire) and 13 April 2016 (post-fire) during leaf-off conditions. Here, we use the resulting map of categorized burn severity, which had an R 2 = 0.84 based on 156 field plots, to stratify the AOI into 6 severity strata ranging from 1 (unburned) to 6 (high severity) [27]. These results were used to spatially categorize the pre-and post-fire ALS data products that follow.

Airborne Laser Scanning Data
Two ALS acquisitions were used in this study. The first collection was contracted by the United States Geological Survey and collected in leaf-off conditions on 12 April 2015 over a broad geographic area with a mean scan density of 9.0 ± 2.8 pulses m −2 . The second collection was contracted by the USDA Forest Service for the same AOI as the WV-3 dataset and was collected on 19 January 2017, also during leaf-off conditions, and with a mean scan density of 9.5 ± 3.2 pulses m −2 . It is important to note that the second acquisition was flown one growing season following the February and March 2016 fires.

Preprocessing
Both ALS datasets were processed using an identical methodology. The Toolbox for Forest LiDAR data Filtering and Forest Studies (Tiffs; [42]) was used to classify ground and above-ground returns, generate a 1 m horizontal resolution digital elevation model from returns classified as ground, and then to bin counts of above-ground returns into continuous 7.5 m horizontal resolution and 1 m vertical resolution rasters to match the resolution of the WV-3 data product. Canopy height profiles (CHP) were then derived from these rasters following the methodology and equations presented in Skowronski et al. [35], producing a stack of 25 rasters representing vertical CHPs at 1 m increments (e.g., 1-2 m, 2-3 m, etc.). A 3 × 3 cell focal mean was used to smooth the data, in order to minimize spatial co-registration errors between the two datasets.
CHPs were converted to canopy bulk density (CBD) profiles following published equations for the New Jersey Pinelands [35]. These profiles were reduced to several variables that spanned the extent of the AOI ( Table 1). The continuous rasters of these variables were intersected with the continuous WV-3 severity map to describe changes in canopy structure based on different estimates of fire severity.  CBD for height bin 24-25 m CBD max maximum CBD in the profile CBD ladder sum of the CBD profile 1 ≤ 3 m CBD mid sum of the CBD profile 3 ≤ 8 m CBD canopy sum of the CBD profile > 8 m CBD total sum of the CBD profile CBD ladder : CBD mid (CBD ladder -CBD mid )/(CBD ladder + CBD mid ) CBD ladder : CBD canopy (CBD ladder -CBD canopy )/(CBD ladder + CBD canopy ) CBD mid : CBD canopy (CBD mid -CBD canopy )/(CBD mid + CBD canopy )

Classification of ALS Data
In their raw form, the CBD profiles are difficult to understand because of the volume and variability of the three-dimensional distributions. As a way of simplifying the 25 vertical bin CBD profiles, we used the Iterative Self-Organizing Data Analysis Technique (ISODATA) unsupervised classification in ERDAS Imagine [43] to create 20 classes based on a convergence threshold of 0.95. In this way, each cell was assigned one of the 20 classified shapes and allowed us to resolve the vertical distributions into a single raster layer. We decided to use 20 shapes with the criteria that each shape represented at least 1% of the landscape either pre-or post-fire. A total of 20 classes was chosen to provide a manageable number of classes, but also be large enough to capture the range of structural patterns evident in the data. We attempted three options for this classification process: (1) classify pre-and post-CBD stacks using individual classification schemes, (2) classify the pre-CBD stack and use this classification matrix to then classify the post-data, and (3) classify the post-CBD stack and use this classification matrix to classify the pre-data. The first option was dismissed because the classification results did not allow for a 1:1 comparison, as they provided very different representations of the classes between the two datasets. We settled on using the post-fire classification matrix to classify both datasets because very distinct canopy shapes were identified resulting from the fires and these shapes did not appear in the pre-fire classification sets.

Regression Analysis
To quantify relationships between pre-fire ALS-derived fuel variables and the continuous estimates of fire severity we followed the methodology of McCarley et al. [34]. As a preliminary screening, to eliminate poorly correlated variables, we calculated correlations between each available variable ( Table 1) and fire severity and eliminated those that had an absolute correlation coefficient less than 0.5. Next, we developed three models to predict fire severity on a per-pixel basis using ordinary least squares regression (OLS) in Systat 13 analysis software (Systat Software, Inc., Version 13). The first of the models only took into account the ALS-derived fuel parameters which had made it through the screening described above, the second model added the ISODATA canopy shape classification as indicator variables, and the third model added the ignition type (head fire, flanking fire, and backing fire) as indicator variables which were provided as vector data by the burn boss of these prescribed fires. We also considered a fourth model that considered the day each pixel was burned to broadly account for meteorological and fuel moisture variation, but do not report the results here because it provided no improvement over the prior models.

Landscape Pattern Analysis
We used Fragstats to parse the effects of the varying fire intensity on patch sizes [44]. Fragstats produces a wide range of spatial metrics at the patch, class and landscape scale. In the context of this analysis, a patch is a spatially contiguous group of pixels with the same class label (CHP class number in our study). We focused on three metrics: mean patch size, patch density, and Simpson's diversity index. Mean patch size is the arithmetic average of the area of all patches, irrespective of class, within a specified region (pixels with a specified burn severity stratum). Patch density is the average number of patches per hectare. Simpson's diversity index is a measure of the probability that any randomly selected two pixels are from different classes [45]. Thus, the higher the index, the less any one class dominates the landscape [44].

Canopy Fuel Characteristics
The spatial results of the estimation of three selected canopy fuel variables derived from the CHPs, alongside pre-and post-burn WV-3 short wave infrared (SWIR) imagery, are presented in Figure 1. It is notable that the pre-fire arrangement of fuels does not necessarily match the patterning of the fire, which indicates that ignition pattern and environmental variables played a large role in the fire dynamics and resulting fire severity. There is, however, strong agreement between the spatial patterning in the WV-3 SWIR data product and the post-fire fuel loading variables.

Canopy Classification Results
The ISODATA unsupervised classification was run using the full scene of post-fire canopy height profile data. We then used these post-fire classification signatures to classify both the pre-and post-fire datasets. The classification resulted in 20 canopy classes that varied in amplitude, number of peaks, and height. As an example of the variability of canopy shape across the study area, we emphasize four representative classes in Figure 2; the full set of canopy height profile graphs and associated statistics are presented in Appendix A. Table 2 summarizes these shapes for the entire study area in the terms of the ALS-derived canopy bulk density variables, CBD total , CBD max , and CBD ladder . In general, as the class number increases, so does the CBD total . The trends of CBD max and CBD ladder are similar but visibly weaker, illustrating the variability of fuel loading through various portions of the canopy in these stands that may be attributed to disturbance history rather than simple successional trajectories.

Landscape Pattern Analysis
We used Fragstats to parse the effects of the varying fire intensity on patch sizes [44]. Fragstats produces a wide range of spatial metrics at the patch, class and landscape scale. In the context of this analysis, a patch is a spatially contiguous group of pixels with the same class label (CHP class number in our study). We focused on three metrics: mean patch size, patch density, and Simpson's diversity index. Mean patch size is the arithmetic average of the area of all patches, irrespective of class, within a specified region (pixels with a specified burn severity stratum). Patch density is the average number of patches per hectare. Simpson's diversity index is a measure of the probability that any randomly selected two pixels are from different classes [45]. Thus, the higher the index, the less any one class dominates the landscape [44].

Canopy Fuel Characteristics
The spatial results of the estimation of three selected canopy fuel variables derived from the CHPs, alongside pre-and post-burn WV-3 short wave infrared (SWIR) imagery, are presented in Figure 1. It is notable that the pre-fire arrangement of fuels does not necessarily match the patterning of the fire, which indicates that ignition pattern and environmental variables played a large role in the fire dynamics and resulting fire severity. There is, however, strong agreement between the spatial patterning in the WV-3 SWIR data product and the post-fire fuel loading variables.

Canopy Classification Results
The ISODATA unsupervised classification was run using the full scene of post-fire canopy height profile data. We then used these post-fire classification signatures to classify both the pre-and post-fire datasets. The classification resulted in 20 canopy classes that varied in amplitude, number of peaks, and height. As an example of the variability of canopy shape across the study area, we emphasize four representative classes in Figure 2; the full set of canopy height profile graphs and associated statistics are presented in Appendix A. Table 2 summarizes these shapes for the entire study area in the terms of the ALS-derived canopy bulk density variables, CBDtotal, CBDmax, and CBDladder. In general, as the class number increases, so does the CBDtotal. The trends of CBDmax and CBDladder are similar but visibly weaker, illustrating the variability of fuel loading through various portions of the canopy in these stands that may be attributed to disturbance history rather than simple successional trajectories.    Table 2. Statistics of the 20 Canopy classes, including the percentage of the class found in the study area before and after the fires. Cells with a proportion of the area greater or equal to 10% are denoted with an "*".

Regression Results
The results of the regression analysis are presented in Table 3. In the first model, fuel loading alone, accounted for 37% of the variation in fire severity. The addition of the classified canopy shape variables minimally improved the predictive power of the model, increasing the R 2 by 5 percentage points. The final model shows a much larger improvement with the addition of firing pattern, with the final model explaining 61% of the variability in the dataset (Table 3).  Table 4 presents a comparison of the relative importance of each of the input variables for the prediction of fire severity in the study area. Notably, the ratio of CBD mid to CBD canopy has the largest absolute standardized coefficient of all variables, which indicates its importance for, and positive correlation with, predicting fire severity. CBD canopy also carried more weight than the other fuel variables, though these other variables were significant in the model (p < 0.001; Table 4).
For the most part, the canopy shape variables did not prove to be particularly important in predicting the severity of fire in the study area. However, there were a few notable exceptions. For instance, classes 2 and 12 both had an absolute relative influence in the dataset of over 0.2. However, these classes had opposite effects, with class 2 exhibiting a negative influence and class 12 exhibiting a positive one. There is a marked contrast in the distribution of fuels in these canopy shapes, with class 2 having low fuels mostly in the understory and class 12 having a large amount of fuel in both the mid-story and canopy strata (Figure 2; Appendix A). It is also notable that both classes 4 and 16 had high fuel loads in the understory, but a strong negative correlation with fire severity.
The firing technique used by the fire practitioners also proved to be important in the prediction of fire severity. Heading fire had a relative influence of 0.24 and was positively correlated with fire severity, while backing fire had a relative influence of −0.34 that was negatively correlated with fire severity. Flanking fire was considered but was not significant enough to enter the model.

Canopy Fuel Changes
To examine change in canopy fuels as a result of the fires, we spatially stratified the burn area by the fire severity estimated from WV-3 in our previous paper [27]. When excluding severity class 1, fire severity was negatively associated with pre-fire CBD total (R 2 = 0.82), but positively correlated with both CBD max and CBD total (R 2 = 0.79 and 0.62), respectively (Figure 3). For the post-fire variables, when excluding severity class 1, there is a negative linear relationship with R 2 values of 0.75, 0.78, and 0.79 for CBD total , CBD max , and CBD ladder , respectively ( Figure 3).

Canopy Class Change
Following the classification routine, each pixel in the pre-and post-ALS dataset was assigned a canopy shape class. The spatial results of this assignment are illustrated for the pre-and post-fire datasets along with the WV-3 derived fire severity for clarity, in Figure 4. The distributions of pre-and post-fire cover classes are stratified by severity strata in Figure 5. Pre-fire, there is a broad diversity of CHP classes in lower severity strata (1-3), but less diversity and several dominant CHP classes (e.g., class 6, class 16, and class 20) in the higher severity strata (4-6; Figure 5a). Following the fires, there is still a high diversity of CHP classes in the lower severity strata, but the class distributions have shifted. The higher severity strata maintain a low diversity of CHP classes, but there is a clear shift in the dominant classes in these strata, with CHP class 2 becoming the dominant feature in these areas (Figure 5b).

Canopy Class Change
Following the classification routine, each pixel in the pre-and post-ALS dataset was assigned a canopy shape class. The spatial results of this assignment are illustrated for the pre-and post-fire datasets along with the WV-3 derived fire severity for clarity, in Figure 4.

Canopy Class Change
Following the classification routine, each pixel in the pre-and post-ALS dataset was assigned a canopy shape class. The spatial results of this assignment are illustrated for the pre-and post-fire datasets along with the WV-3 derived fire severity for clarity, in Figure 4. The distributions of pre-and post-fire cover classes are stratified by severity strata in Figure 5. Pre-fire, there is a broad diversity of CHP classes in lower severity strata (1-3), but less diversity and several dominant CHP classes (e.g., class 6, class 16, and class 20) in the higher severity strata (4-6; Figure 5a). Following the fires, there is still a high diversity of CHP classes in the lower severity strata, but the class distributions have shifted. The higher severity strata maintain a low diversity of CHP classes, but there is a clear shift in the dominant classes in these strata, with CHP class 2 becoming the dominant feature in these areas (Figure 5b). The distributions of pre-and post-fire cover classes are stratified by severity strata in Figure 5. Pre-fire, there is a broad diversity of CHP classes in lower severity strata (1-3), but less diversity and several dominant CHP classes (e.g., class 6, class 16, and class 20) in the higher severity strata (4-6; Figure 5a). Following the fires, there is still a high diversity of CHP classes in the lower severity strata, but the class distributions have shifted. The higher severity strata maintain a low diversity of CHP classes, but there is a clear shift in the dominant classes in these strata, with CHP class 2 becoming the dominant feature in these areas (Figure 5b).  Figure 6 visually summarizes the pre-to post-fire CHP transitions for the six burn severity strata. For burn severity stratum 1, the most common transitions are unchanged class numbers (lying on the diagonal), or slightly increased class numbers (i.e., falling below the diagonal). Considering that the CHP class numbers are generally positively correlated with increasing CBDtotal values, this suggests that low severity burns are associated with no change, or a slight increase in CBDtotal. In contrast, for burn severity stratum 2, although the no-change values on the diagonal remain common, most transitions are to lower class numbers, i.e., generally lower CBDtotal values. Not surprisingly, as the burn severity strata increase from 3 to 6, the unchanged classes become very rare, and post-fire classes become concentrated in just a few low-numbered (i.e., low CBDtotal) classes, particularly classes 1 to 3. It is striking however, that burn severity strata 5 and 6 are both dominated by a small number of pre-fire CHP classes, particularly classes 7, 16 and 20.
Combining the transition matrices of Figures 6 and CBDtotal, CBDmax and CBDladder statistics of Table 2, allows the generation of summary statistics of changes in CBD measures associated with each burn severity class for the study area (Table 5). Although the values of each of the CBD measures show a strong decreasing trend from burn severity classes 1 to 6, it is not linear. This is particularly evident for CBDladder, where there is a particularly large decrease between classes 5 and 6. Furthermore, in the case of CBDmax, burn severity class 6 is associated with a smaller decrease than class 5. Table 5. Summary CBDtotal, CBDmax and CBDladder value by burn severity class.  Figure 6 visually summarizes the pre-to post-fire CHP transitions for the six burn severity strata. For burn severity stratum 1, the most common transitions are unchanged class numbers (lying on the diagonal), or slightly increased class numbers (i.e., falling below the diagonal). Considering that the CHP class numbers are generally positively correlated with increasing CBD total values, this suggests that low severity burns are associated with no change, or a slight increase in CBD total . In contrast, for burn severity stratum 2, although the no-change values on the diagonal remain common, most transitions are to lower class numbers, i.e., generally lower CBD total values. Not surprisingly, as the burn severity strata increase from 3 to 6, the unchanged classes become very rare, and post-fire classes become concentrated in just a few low-numbered (i.e., low CBD total ) classes, particularly classes 1 to 3. It is striking however, that burn severity strata 5 and 6 are both dominated by a small number of pre-fire CHP classes, particularly classes 7, 16 and 20.

Change in CBDladder
Combining the transition matrices of Figure 6 and CBD total , CBD max and CBD ladder statistics of Table 2, allows the generation of summary statistics of changes in CBD measures associated with each burn severity class for the study area (Table 5). Although the values of each of the CBD measures show a strong decreasing trend from burn severity classes 1 to 6, it is not linear. This is particularly evident for CBD ladder , where there is a particularly large decrease between classes 5 and 6. Furthermore, in the case of CBD max , burn severity class 6 is associated with a smaller decrease than class 5. Table 5. Summary CBD total , CBD max and CBD ladder value by burn severity class.

Spatial Patterning
We examined contrasts in the spatial patterning of classified canopy height profiles across the range of burn severity strata. Table 6 presents several of the landscape metrics derived from the Fragstats analysis. There are several notable patterns that occur in these data. Generally, for the prefire data, the landscape is more heterogeneous in the lower severity classes with patch size and density decreasing in the parts of the landscape that eventually were high-severity burns. Following the fires, the landscape heterogeneity followed the same pattern and magnitude, but in severity strata 3-6, these patches became larger and notably more homogeneous than the pre-fire landscape. This is most obvious in the areas with the highest severity, where mean patch size more than doubled and the Simpson's Diversity Index fell by more than half (Table 6). In this case, higher severity fire Figure 6. Transition matrices of CHP classes pre-fire to post-fire. A gray tone indicates classes with a frequency of occurrence of 0%. Locations down the diagonal indicate areas that were unchanged after the fire. Note that each image was individually contrast-enhanced in order to show the patterns within each matrix, and thus the colors are not directly comparable between the matrices.

Spatial Patterning
We examined contrasts in the spatial patterning of classified canopy height profiles across the range of burn severity strata. Table 6 presents several of the landscape metrics derived from the Fragstats analysis. There are several notable patterns that occur in these data. Generally, for the pre-fire data, the landscape is more heterogeneous in the lower severity classes with patch size and density decreasing in the parts of the landscape that eventually were high-severity burns. Following the fires, the landscape heterogeneity followed the same pattern and magnitude, but in severity strata 3-6, these patches became larger and notably more homogeneous than the pre-fire landscape. This is most obvious in the areas with the highest severity, where mean patch size more than doubled and the Simpson's Diversity Index fell by more than half (Table 6). In this case, higher severity fire homogenized the landscape relative to the unburned and low-severity states. This patterning is evident in Figure 4a,c.

Discussion and Conclusions
The objectives of this study were focused on examining the interaction between forest structure and fire severity using a time series of both active (ALS) and passive (WV-3) remote sensing datasets. To simplify the three-dimensionality of the canopy across the landscape, we used an ISODATA classification to derive categorical classes of canopy shapes, thus statistically separating the data and allowing for a two-dimensional analysis. Several classes dominated the landscape pre-fire, but transitioned to various new classes, depending on the severity of fire. Following the fires, the higher burn severity regions were dominated by classes with a larger CBD ladder component and generally lower CBD max position in the canopy (lower crowns). There were obvious patterns in the transitions between classes that were dependent on the severity of the fire, with mostly unchanged CHP shapes in the low severity classes, a wide diversity of transitions in the mid-severity classes, and transitions to a low diversity of post-CHP classes for the high severity pixels. The spatial patterning of the landscape changed relative to fire severity, with increasing fire severity resulting in larger patches and a lower diversity of canopy classes.
Using OLS-regression analysis, we were able to partition the effects of fuel loading, canopy shape, and firing pattern variables. We found that certain variables were distinctively more important than others in these regressions. CBD canopy , CBD mid :CBD canopy , CHP class 12, and heading fire ignition were relatively strong positive correlations in the model, while several CHP classes (2, 4, and 16) and backing fire ignition had the strongest negative correlations with fire severity. The results of the OLS analysis illustrate the importance of canopy fuel loading and canopy shape on fire severity. Surprisingly, most individual fuel loading variables that we considered did not have a strong influence. However, the strength of the CBD mid :CBD canopy ratio in the model illustrates that the proportion of fuels in different canopy strata have a great influence on fire behavior. This result is interesting because it illustrates the need to understand where fuels are, not just the total amount or the maximum amount of fuels. This is further demonstrated in our analysis of CHP classes as a few of these classes had large relative importance, either positively or negatively influencing fire severity. Of the classes with the strongest negative influence, classes 4 and 16 are very similar in that they have very high fuel loads at the base of the canopy and very little overstory fuel. Because of the high density of the understory, it is possible that fuels may have been wetter and not as receptive to fire or that high drag due to the high fuel loads may have negatively influenced fire behavior. Conversely, CHP class 12 had a positive influence in the model and had a bi-modal distribution of fuels that were nearly equal in the canopy and mid-canopy strata. Our results indicate that the vertical distribution of fuel is an important factor and that subtle differences have defined effects on fire behavior.
Factors like topography, fuel moisture status, meteorological conditions, and ignition technique obviously have a strong influence on fire behavior and resulting severity. We chose to ignore topographical influences in this study because there is very little topographic variation in the area of interest, though we may have missed the effect of lowland areas in the study which would have had higher fuel moisture contents. Following Wimberly et al. [46], we used the date of burn as a way to attempt to broadly account for meteorological and fuel status variation. However, likely because the burns happened over a short span of time, this variable did not increase the strength of the model and is therefore not reported here. Ignition pattern, as expected, had a very strong influence and greatly improved our prediction of fire severity. Not surprisingly, heading fire had a positive correlation and backing fire had a negative correlation on fire severity, while flanking fire did not have enough of an influence to enter the model.
The classification of ALS data into CHP profiles provides a way to simplify the representation of the vertical canopy structure, much the way that landcover is classified into representative groups using spectral reflectance datasets. Recently, Hoff et al. [47] used a classification approach to classify the landscape using various ALS-derived variables and aligned these with Landsat-derived burn severity data [48] for their area of interest. Using a classification scheme allowed us to group similar canopy profiles and to display them two dimensionally, thus allowing us to make use of more standard analysis techniques for two-dimensional datasets while still maintaining the ability to make inferences about the data in the context of vertical forest structure. Mutlu et al. [49] made use of ALS data for the classification of fuels on a landscape into seven of the Anderson fuel models [50]. While this work represented a step towards using ALS data as auxiliary data for developing maps of fuels for input into fire behavior models, such as FARSITE [17], the technique did not explicitly account for vertical fuel structure or provide a way to directly examine physical changes in canopy structure. As fire modeling continues to move toward a CFD approach that integrates physical estimates of structure [51], more detailed physical representations of forest canopy structure in three-dimensions are necessary. By representing forest structure using a classification scheme, our results allow for an examination of the transitions from pre-fire structure to the resulting post-fire structure. Looking closely at individual CHP class to class transitions is valuable because it could potentially allow prescribed fire managers to modify fire behavior with intentions of converting from one stand structure to another.
Our results indicate that varying levels of fire severity affect forest heterogeneity in much different ways. For instance, high severity fire has a much greater influence on Simpson's Diversity Index than low intensity fire. Additionally, looking broadly at the transition matrices illustrates that diversity decreases not only in the arrangement of classes across the landscape with severity, but also affects the number of classes that appear following the fire. So, both horizontal and vertical diversity are inversely related to higher fire severities.
Our study also contributes to the debate regarding the limitations of spectral indices of burn severity, such as the WV-3 data we used [19,22,23]. Figures 5 and 6 demonstrate that 2-D spectral indices are inherently limited, and that a spectral burn severity index value (such as our severity class 3) encompasses a wide range of structural changes to the forest. Indeed, any single measure of burn severity, even those derived from multi-temporal ALS data, such as the changes in CBD total , CBD max , and CBD ladder shown in Table 5, obscure distinctions in the changes and the resulting physical structure of the forest that may be important for resource managers and scientists alike. Therefore, we recommend that 3-D characterization of burn severity should be a major focus of future remote sensing fire research.
There are several limitations to this study. Field validation and verification is always challenging when considering forest structural characteristics that move beyond simple height or cover metrics. Because of the opportunistic nature of our study, we used previously developed models that estimated CBD for each bin [52]. These models were developed from plot data within 10 miles of our site and in the same forest type, so they were deemed to be suitable. It is notable though, that allometric estimates of canopy bulk density profiles, which typically rely on diameter and height measurements of individual stems, would not be effective in estimating the canopy bulk density of stems that had been fully or partially consumed by fire, thus making our technique superior as it directly estimated these changes. Ideally, fuel plots would have been installed prior to the fires and allometry used to estimate canopy bulk density [53], with models derived to estimate these variables to the ALS data [54]. In this study, though we use models and method verification from published studies [55], we are reliant on the ALS data to provide estimates of overstory consumption. There are also limitations directly associated with the parameters of the ALS data collections. Ideally, the ALS data would have been conducted immediately prior to and following the burn in order to best characterize the immediate impacts on the canopy. Additionally, there was likely some incongruity between the post-processed data products as flight altitude, sensor, sampling density, and other collection specifications were not identical. This is not a trivial issue for transfer of our methodologies to other regions because of the difficulty obtaining robust pre-and post-ALS surveys, especially in the context of a wildfire.
Further investigation of forest canopy structure, as it applies to wildland fire, is necessary in order to better provide the necessary inputs for CFD modeling, basic studies of heat transfer through fuel elements, and an understanding of combustion across landscapes. These studies of canopy structure are critical for improving our understanding of momentum fluxes through forest canopies and into and out of combustion zones. This work will also be important for the development of tools for pre-suppression planning for the mitigation of wildland fire, informing management objectives on the implementation of prescribed fire and for the evaluation of wildland effects and prescribed fire effectiveness. While ALS data has proven to be a good tool for landscape scale studies, these same CHP ideas may also transfer to both space-born LiDAR systems such as GEDI [55] and small-scale TLS data collections [56]. The power of these measurements and their application will be amplified as more spatially coincident collections are made which will allow for the examination of disturbance interactions and long-term trajectories across landscapes, as we now see coming to fruition with the Landsat data archive.
Our study aimed to develop our knowledge of the interactions between pre-fire fuel loading, fire severity, and post-fire fuel loading. We found that the pre-fire vertical distribution of fuels is an important factor in predicting fire severity across the landscape, most notably the ratio of mid-story to overstory fuels. Additionally, our classification approach demonstrated obvious patterns of transition CHP classes that were dependent on the severity of the fire, with mostly unchanged shapes in the low severity classes, a wide diversity of transitions in the mid-severity classes, and a low diversity of high frequency transitions for the high severity pixels. Finally, we found that the overall structural heterogeneity in our study area, both horizontally and vertically, was reduced in a way that was directly proportional to fire severity and resulted in a new set of conditions.