Assessing Legacy E ﬀ ects of Wildﬁres on the Crown Structure of Fire-Tolerant Eucalypt Trees Using Airborne LiDAR Data

: The ﬁre-tolerant eucalypt forests of south eastern Australia are assumed to fully recover from even the most intense ﬁres; however, surprisingly, very few studies have quantitatively assessed that recovery. The accurate assessment of horizontal and vertical attributes of tree crowns after ﬁre is essential to understand the ﬁre’s legacy e ﬀ ects on tree growth and on forest structure. In this study, we quantitatively assessed individual tree crowns 8.5 years after a 2009 wildﬁre that burnt extensive areas of eucalypt forest in temperate Australia. We used airborne LiDAR data validated with ﬁeld measurements to estimate multiple metrics that quantiﬁed the cover, density, and vertical distribution of individual-tree crowns in 51 plots of 0.05 ha in ﬁre-tolerant eucalypt forest across four wildﬁre severity types (unburnt, low, moderate, high). Signiﬁcant di ﬀ erences in the ﬁeld-assessed mean height of ﬁre scarring as a proportion of tree height and in the proportions of trees with epicormic (stem) resprouts were consistent with the gradation in ﬁre severity. Linear mixed-e ﬀ ects models indicated persistent e ﬀ ects of both moderate and high-severity wildﬁre on tree crown architecture. Trees at high-severity sites had signiﬁcantly less crown projection area and live crown width as a proportion of total crown width than those at unburnt and low-severity sites. Signiﬁcant di ﬀ erences in LiDAR -based metrics (crown cover, evenness, leaf area density proﬁles) indicated that tree crowns at moderate and high-severity sites were comparatively narrow and more evenly distributed down the tree stem. These conical-shaped crowns contrasted sharply with the rounded crowns of trees at unburnt and low-severity sites and likely inﬂuenced both tree productivity and the accuracy of biomass allometric equations for nearly a decade after the ﬁre. Our data provide a clear example of the utility of airborne LiDAR data for quantifying the impacts of disturbances at the scale of individual trees. Quantiﬁed e ﬀ ects of contrasting ﬁre severities on the structure of resprouter tree crowns provide a strong basis for interpreting post-ﬁre patterns in forest canopies and vegetation proﬁles in Light Detection and Ranging (LiDAR) and other remotely-sensed data at larger scales.


Introduction
Wildfire impacts on forest ecosystems worldwide are predicted to change under fire and climate regimes [1][2][3][4]. The incidence of uncontrolled large wildfires has increased in temperate regions in particular, resulting in changing effects on forest structure and sustainability [5]. Emerging evidence suggests that the frequency of wildfires in temperate south-eastern Australia has increased in recent years [6][7][8]. For example, during 2003-2014, more than 4.3 million ha of eucalypt forest in the south-eastern state of Victoria was burnt by wildfires-a much greater area than in previous decades [6]-including 350,000 ha that was burnt twice in those 11 years [8].
Impacts of fire at the forest-stand level can range from the partial consumption of surface fuels and understorey vegetation (low severity) to the burning of all vegetation strata including the crowns of the overstorey trees; i.e., high-severity fires, also known as "crown fires" [9,10]. Fires of high severity, by definition, impact tree crowns. The crown damage and loss associated with high-severity fires will impact tree carbon uptake, influencing stand productivity at least until the functional leaf area is recovered [11,12]. In addition to productivity, the impacts of high-severity fires on the height and architecture of tree can have practical implications [13], including the accurate estimation of canopy biomass and carbon stores based on tree allometry. Tree-based allometric equations assume a constant relationship of tree biomass or volume with measurable attributes such as diameter and height and do not account for either short-term decreases in tree foliage and branches after fire or more persistent effects on stem or branch architecture after severe fire [14]. The impacts of fire on tree crowns will likely vary depending on the pre-fire architecture [15] and the evenness of fire severity, as well as the evenness of resprouting in fire-tolerant crowns, leading to variable changes in both the horizontal and vertical distribution of foliage and branches. Despite the many potential effects of fire on forest canopy structure, recent small-scale studies in Western Australia and Tasmania [14,16] have highlighted how little is known about the lasting effects of fires on eucalypt crowns. While broad-scale studies of landscapes dominated by eucalypt forests indicate the recovery of remotely-sensed spectral responses within about six years of wildfire [17], the contribution of eucalypt crowns versus understoreys to these spectral data, and in particular the time taken for fire-tolerant eucalypts to re-establish a full crown after severe fire, remain under-examined [18].
Understanding the impacts of fire on tree crown structure and recovery requires the accurate assessment of both horizontal (x and y dimension) and vertical (z dimension) crown attributes. Crown attributes include the structure of foliage, stem, twigs and branches of the individual tree, which can be characterised and assessed either directly from field-based measurements or indirectly from passive or active remote sensing techniques [19] such as LiDAR (Light Detection and Ranging). Forest crown fires can influence key attributes such as crown base height, crown depth, crown width and crown projection area (CPA) [20,21]. Estimating these and more complex crown attributes (e.g., canopy volume, foliage density profiles) in the field is labor-intensive, time-consuming and can be prone to error [22]. In comparison, airborne LiDAR data can provide estimates of canopy structural parameters-i.e., canopy height and canopy cover [23,24]-at a range of scales from the individual tree [25,26] to the plot or stand [27,28] to entire landscapes [29][30][31]. Discrete return small-footprint LiDAR has been successfully used to estimate tree heights [25,32], to characterise the canopy structure of several forest types [33][34][35], to measure and map canopy gaps [22,36] and to assess canopy height changes over time [37].
Various studies indicate the potential of precise estimates of horizontal and vertical canopy structure using airborne LiDAR data with a 10 to 30 cm footprint [27,38,39]; however, there are few studies that have used these data to assess the spatial and temporal dynamics of the forest canopy structure and its recovery after natural disturbances [17,22,36]. Most LiDAR-based studies have examined the dominant tree species either in Northern-Hemisphere boreal forests [40,41] or in rainforests [42]. In addition, several studies have examined the impacts of fire severity on fire-tolerant and fire-sensitive forests in the Northern Hemisphere, specifically addressing canopy structural attributes [43][44][45] and dynamics at landscape scales. However, while a handful of studies have used airborne LiDAR data to describe the horizontal [46,47] and vertical structure [29,31,48] of Australian forest canopies at stand and landscape-scales, no studies to our knowledge have examined the effects of high-fire severity wildfires on the individual tree crowns of temperate eucalypts, which dominate the fire-prone forests of southern Australia.
Fire-tolerant eucalypts are assumed to fully recover from even the most intense fires through epicormic (stem and crown) and/or basal resprouting [49]; however, surprisingly, very few studies have quantitatively assessed that recovery. This has implications for understanding the post-fire recovery of forest productivity [50], particularly because the fire-tolerant eucalypt forests of temperate Australia are often assumed to be resilient stores of carbon under changing fire regimes [51]. In this paper, we assess the tree-crown architecture of fire-tolerant eucalypts 7 to 8.5 years after they were burnt by a high-severity wildfire that started under extreme weather conditions on "Black Saturday", 7 February 2009 [52]. We assess individual trees in the field across replicated sites of four fire severities (unburnt, low, moderate, high), and also use high-density airborne LiDAR data to develop tree-based metrics that quantitatively describe the crown architecture (e.g., cover, density, leaf area density, evenness, clumping) as it is influenced by fire. Based on knowledge of fire crown scorch and epicormic recovery, and including the development of unique tree-based LiDAR crown metrics, our specific objectives were to 1) examine the persistence of changes in crown horizontal structure related to the consumption of small branches, and 2) quantitatively describe changes in the crown vertical structure associated with epicormic growth down the stem. Overall, we aimed to quantitatively assess any legacy effects of wildfires on fire-tolerant tree crowns so that we can better understand the longevity of wildfire impacts and the post-fire patterns in forest canopies and vegetation profiles at larger scales.

Study Area
Our study area is in the Central Highlands region of south-eastern Victoria, Australia, approximately 100 km north-east of Melbourne (centred on the township of Marysville: 37.510 S, 145.750 E; Figure 1). The long-term mean annual rainfall is 1343 mm, with the majority (464 mm) falling in winter (compared with 209 mm in summer). The climate of the study area is temperate with mean daily minimum temperatures ranging from 8 • C in July (winter) to 23 • C in February (summer), and mean daily maximum temperatures from 11 • C (July) to 28 • C (February; Marysville weather station) [53]. The landscapes of the study area have a variable topography and complex terrain (elevation: 352-762 m; slope: 2 • -30 • ), with many of the study's sites in remote locations. We focused on the herb-rich Foothill Forest [54,55] which is an "open forest" (with tree heights of 10 to 30 m and a projective foliage cover of 30%-70% [56]). This forest is typically dominated by a mix of fire-tolerant eucalypt species that survive wildfires with both epicormic and basal resprouting, such as Eucalyptus obliqua L'Hér (messmate stringybark), E. radiata Sieber ex DC. (narrow-leaf peppermint), E. dives Schauer (broad-leaved peppermint), E. viminalis Labil. (manna gum), and E. cypellocarpa L.A.S. Johnson (mountain grey gum).

Figure 1.
Location of the study sites (as coloured dots), extent of 2016 LiDAR data, and the study area (defined by the overlap of the 2009 wildfire, the LiDAR data, and the distribution of the herb-rich Foothill Forest).

Study Design
Our study design involved clusters of plots at each of the 17 sites to improve sampling efficiencies in the logistically challenging (steep and isolated) terrain. These 17 sites represented the most accessible combinations of the overlap of the wildfire complex with our selected forest type and the extent of the LiDAR data. They were four replicates each of four wildfire severity types-i.e., "unburnt" (UB), "low severity" (LS), "moderate severity" (MS) and "High severity" (HS)-plus an additional MS site to ensure more comprehensive coverage of the study area. All sites were within or around the boundaries of the Kilmore East-Murrindindi fire complex, which burnt over 400,000 ha in a three-week period starting on "Black Saturday", 7 February 2009 ( Figure 1). Wildfire severity was classified at a 1:25,000 scale in 2009 by the state government environmental department using optical pre and post-fire remote sensing data and extensive ground-truthing in the two months after the wildfire [9]. Study sites were all located within the herb-rich Foothill Forest, and were carefully selected to ensure they fully occurred within one fire severity type as described in Department of Sustainability and Environment report (DSE) [9], which was verified either from field data from a previous study [55] or by remote-sensing estimates of the Composite Burn Index (CBI) [57]. Selected unburnt sites had no signs of wildfire in the 2009 assessment, had not been burnt by any fire for more than 30 years prior to the 2009 fire, and were as close in distance as possible to the other severity sites around the wildfire boundary. Low severity sites were burnt by low-intensity surface fire as a mosaic of mostly unscorched crowns or with light crown scorching (<35%, severity class 4 & 5a). Moderate severity was a mixed mosaic of crown burning and scorching ranging from 30%-65% (severity class 3), and high severity was an intense fire that consumed the entire understorey and burnt or scorched 60%-100% (severity class 1 & 2) of the canopy [9]. Our overall study design, including the links between field and LiDAR data, is summarised in Supplementary Figure S1.

Study Design
Our study design involved clusters of plots at each of the 17 sites to improve sampling efficiencies in the logistically challenging (steep and isolated) terrain. These 17 sites represented the most accessible combinations of the overlap of the wildfire complex with our selected forest type and the extent of the LiDAR data. They were four replicates each of four wildfire severity types-i.e., "unburnt" (UB), "low severity" (LS), "moderate severity" (MS) and "High severity" (HS)-plus an additional MS site to ensure more comprehensive coverage of the study area. All sites were within or around the boundaries of the Kilmore East-Murrindindi fire complex, which burnt over 400,000 ha in a three-week period starting on "Black Saturday", 7 February 2009 ( Figure 1). Wildfire severity was classified at a 1:25,000 scale in 2009 by the state government environmental department using optical pre and post-fire remote sensing data and extensive ground-truthing in the two months after the wildfire [9]. Study sites were all located within the herb-rich Foothill Forest, and were carefully selected to ensure they fully occurred within one fire severity type as described in Department of Sustainability and Environment report (DSE) [9], which was verified either from field data from a previous study [55] or by remote-sensing estimates of the Composite Burn Index (CBI) [57]. Selected unburnt sites had no signs of wildfire in the 2009 assessment, had not been burnt by any fire for more than 30 years prior to the 2009 fire, and were as close in distance as possible to the other severity sites around the wildfire boundary. Low severity sites were burnt by low-intensity surface fire as a mosaic of mostly unscorched crowns or with light crown scorching (<35%, severity class 4 & 5a). Moderate severity was a mixed mosaic of crown burning and scorching ranging from 30%-65% (severity class 3), and high severity was an intense fire that consumed the entire understorey and burnt or scorched 60%-100% (severity class 1 & 2) of the canopy [9]. Our overall study design, including the links between field and LiDAR data, is summarised in Supplementary Figure S1.

Field Assessments
We assessed our study sites in the field from September to November 2017. Each of our 17 sites consisted of a cluster of three circular plots (12.5 m nominal radius, 0.05 ha; 51 plots in total), which were located 50 m apart (centre to centre) to reduce within-site spatial correlations (Supplementary Figure S2). The locations of our plot centres (as well as the within-plot location of selected trees; see below) were recorded using a Trimble Geo 7x (Sunnyvale, California, United States) with post corrections by the GPS provider (Ultimate Positioning Group, Melbourne, Australia) based on receiver data from a base station located in Marysville, Victoria [58]. We left the differential GPS in place at the plot centre for at least 20 min to maximise the potential for an accurate reading [59]; for logistical reasons, this time was reduced to 2 min when recording the location of selected trees. Given that the accuracy of a differential GPS can be reduced by continuous tree cover [59], we also recorded tree locations relative to the plot centre using a TruPulse360B Laser Rangefinder (Measurably Superior, Colorado, USA).
We selected six to eight live trees per plot for detailed characterisation, hereafter known as candidate trees (with a total of 342 trees; Supplementary Table S1). We divided each plot into four quadrants to minimise the potential for overlap of the candidate-tree crowns and then selected one to two live trees per quadrant using the following criteria (in order of importance): (1) dominant crown position (not obviously over-topped by neighbouring trees); (2) largest available diameter in that quadrant (nominally >20 cm DBHOB (overbark diameter at breast height), 1.3 m); and (3) within 12.5 m of the plot centre (or up to 15.8 m, dependent on the closest suitable tree per quadrant). We assessed the following attributes of the live trees (as illustrated in Figure 2): (1) species; (2) DBHOB; (3) top live height (to highest live leaf); (4) total height (to highest live leaf or dead branch); (5) crown base height (to the origin of the lowest substantial branch or fork); (6) crown width (north-south and east-west for both the bulk of the live crown and the total crown if it included dead protruding branches); (7) presence or absence of stem epicormic resprouts; (8) presence or absence of basal resprouts; and (9) the height of fire scarring as a proportion of total tree height. All tree heights were measured using a Haglöf Vertex IV Ultrasonic Hypsometer (Haglöf, Sweden), and crown widths were measured with the TruPulse360B Laser Rangefinder. We also photographed all candidate trees from the base (panoramic and simple portrait, Samsung S8+ smart phone, Samsung Electronics, China) as a visual record to check with field data and compare with LiDAR data as required.
To assist with the interpretation of LiDAR data, we also measured all trees greater than 20 cm DBHOB in each plot for (1) status (live, dead), (2) DBHOB and (3) location (relative to the plot centre, using a TruPulse360B Laser Rangefinder). In addition, we measured the heights of a selection of trees representing the full range of diameters (typically, an additional five tree heights per plot).

Figure 2.
Field measurements of unburnt and severely burnt trees, including tree height metrics, crown-width metrics (bulk crown width of the live crown, total crown width to outermost extremities, which could be dead) and the presence/absence of epicormic (stem) and basal resprouts.

LiDAR Data Acquisition and Processing
Our study used small-footprint airborne LiDAR datasets acquired from 9 January to 6 August 2016 (about one to two years before field data were collected), which encompassed our 17 field sites ( Figure 1). The datasets were acquired from fixed-wing aircraft and were characterised by a high point density (average 24.34 m 2 ) and horizontal and vertical accuracies of ≤20 cm (additional details in Supplementary Table S2). The LiDAR point cloud data were provided by the state environment department (Department of Environment, Land, Water and Planning (DELWP)) in the form of classified LiDAR returns containing ground, vegetation and other classes. Data were quality checked using a visual inspection of preliminary vertical profiles within our 51 field plots in "3D viewer" using ArcMap 10.5.1 software [60] to verify that the point cloud corresponded with the field-assessed forest structure. We also confirmed that the point cloud data were in the required projection system and included the expected number of returns, proportion of returns by classes, and the data header ("lasinfo" function; '"astools" software (rapidlasso GmbH, Germany)) [61,62], and that there were no height displacements within flight overlap areas [62].
Prior to the calculation of individual tree-crown metrics, point cloud data were heightnormalised using the Lastools "lasheight" function after the development of an accurate Digital Terrain Model (DTM) based on the ground-class points. After normalization, a Z-value histogram for each field plot was used to filter out data noise above the maximum tree height (45 m, from field measurements) and near-ground points below 0.6 m, which is recommended as a minimum terrain value for optimising a DTM [63]. A pit-free Canopy Height Model (CHM) was generated from the Field measurements of unburnt and severely burnt trees, including tree height metrics, crown-width metrics (bulk crown width of the live crown, total crown width to outermost extremities, which could be dead) and the presence/absence of epicormic (stem) and basal resprouts.

LiDAR Data Acquisition and Processing
Our study used small-footprint airborne LiDAR datasets acquired from 9 January to 6 August 2016 (about one to two years before field data were collected), which encompassed our 17 field sites ( Figure 1). The datasets were acquired from fixed-wing aircraft and were characterised by a high point density (average 24.34 m 2 ) and horizontal and vertical accuracies of ≤20 cm (additional details in Supplementary Table S2). The LiDAR point cloud data were provided by the state environment department (Department of Environment, Land, Water and Planning (DELWP)) in the form of classified LiDAR returns containing ground, vegetation and other classes. Data were quality checked using a visual inspection of preliminary vertical profiles within our 51 field plots in "3D viewer" using ArcMap 10.5.1 software [60] to verify that the point cloud corresponded with the field-assessed forest structure. We also confirmed that the point cloud data were in the required projection system and included the expected number of returns, proportion of returns by classes, and the data header ("lasinfo" function; '"astools" software (rapidlasso GmbH, Germany)) [61,62], and that there were no height displacements within flight overlap areas [62].
Prior to the calculation of individual tree-crown metrics, point cloud data were height-normalised using the Lastools "lasheight" function after the development of an accurate Digital Terrain Model (DTM) based on the ground-class points. After normalization, a Z-value histogram for each field plot was used to filter out data noise above the maximum tree height (45 m, from field measurements) and near-ground points below 0.6 m, which is recommended as a minimum terrain value for optimising a DTM [63]. A pit-free Canopy Height Model (CHM) was generated from the normalized point cloud data using the default settings of the algorithm developed by Khosravipour, et al. [64] included in the Lastools software. A range of grid resolutions (0.25, 0.5. 0.75, 1.0, and 2.0 m) were used to initially generate the CHM, and a final resolution of 0.5 m was chosen for subsequent steps because, based on visual comparisons with known locations of trees in the field, it provided sufficient detail to delineate the individual-tree crown boundaries and to identify tree positions.

Individual Tree Detection and Manual Crown Delineation
We used several manual (rather than automated) steps involving multiple lines of evidence to identify individual tree crowns in the normalised point cloud and pit-free CHM, which then formed the basis for identifying tree locations. The automatic identification of field-assessed individual trees in the LiDAR point clouds was not possible because crowns frequently overlap in these natural forests with a complex structure (see Figure 3) and because the differential GPS is less reliable under continuous tree canopies and in remote and topographically variable terrain [59].

Individual-Tree Crown Metrics (Field and LiDAR Data)
Field-based horizontal and vertical crown metrics were calculated for the 342 field-measured candidate trees (Table 1). Field metrics were analysed separately for two tree-size classes based on DBHOB: 20 to 50 cm (representing mature trees) and >50 cm (representing the largest, potentially over-mature trees) [55].
Horizontal and vertical LiDAR metrics (Table 1) were calculated for the 196 trees with nonoverlapping crown boundaries (90 candidate trees plus 106 non-field-assessed trees). The individualtree metrics were calculated using LiDAR data that were clipped to the manually delineated crown boundaries (CPA; Table 1). The LiDAR metrics were mostly quantified from the vertical profile of LiDAR returns classified as canopy hits per height bin [68][69][70] between the tree's top height and 7.5 m, which was the maximum height of understorey plants measured in the field (Figure 4). LiDAR We used ArcMap 10.5.1 (Redlands, CA, USA) [60] to define individual-tree locations, which were based on both the post-processed (differentially corrected) tree coordinates as well as the field-based measurements of the distance and azimuth from the plot centre; we then overlaid the lidar point cloud and drew an indicative crown boundary based on the field-measured tree live-crown width. We used these indicative crown boundaries to identify the local height maxima within the pit-free CHM and LiDAR point cloud. The local height maxima were evaluated as the potential top height of an individual tree using four fixed-sized windows and one variable-size window. The fixed-sized methods were based on (i) ArcMap's Focal Statistics tool; (ii) the Tree Identification and Delineation Algorithm (TIDA) [65]; (iii) the tree_detection function in the lidR package [66] using the CHM; and (iv) the tree_detection function in the lidR package using the point cloud data. For the variable-size method, we used the CHM and the CanopyMaxima function in the Fusion software [67], which indicated tree-top coordinates based on the coefficients of the polynomial relationship between height and crown width of the trees as they were measured in the field. We then finalised a tree's most likely location (as an adjusted shapefile) based on the following steps: (i) a provisional tree top was determined based on close agreement (within 1 m) of at least three of the five height maxima; (ii) the tree top was adjusted based on the LiDAR profile view, indicating the highest point of the tree closest to the provisional tree top; (iii) the provisional tree top was verified by cross-checking the tree height with the field-measured height; and (iv) the tree form was also verified by cross-checking the LiDAR profile with additional field measures (e.g., crown base height and crown width) and the panoramic photo of the tree taken during field measurements. This process led to the confident identification of 90 candidate trees (i.e., those assessed in the field) that had negligible canopy overlap in the LiDAR data. In addition, a further 106 "non-candidate" trees with minimal crown overlap were confidently identified in the LiDAR data either based on their field-mapped TruPulse locations (i.e., without supporting field-based tree height and photographic information), or from the LiDAR profile view just beyond the plot boundaries (and still clearly within the same forest type and wildfire severity types; Supplementary Table S1). The full process led to the final delineation of 196 individual tree crown boundaries in the CHM and LiDAR point cloud data.

Individual-Tree Crown Metrics (Field and LiDAR Data)
Field-based horizontal and vertical crown metrics were calculated for the 342 field-measured candidate trees (Table 1). Field metrics were analysed separately for two tree-size classes based on DBHOB: 20 to 50 cm (representing mature trees) and >50 cm (representing the largest, potentially over-mature trees) [55].
Horizontal and vertical LiDAR metrics (Table 1) were calculated for the 196 trees with non-overlapping crown boundaries (90 candidate trees plus 106 non-field-assessed trees). The individual-tree metrics were calculated using LiDAR data that were clipped to the manually delineated crown boundaries (CPA; Table 1). The LiDAR metrics were mostly quantified from the vertical profile of LiDAR returns classified as canopy hits per height bin [68][69][70] between the tree's top height and 7.5 m, which was the maximum height of understorey plants measured in the field (Figure 4). LiDAR metrics were calculated using the lidR package in R, or the Lastools software, while ArcMap 10.5.1 and the Fusion software [67] were also used to prepare data for metric calculation.
We also calculated all LiDAR metrics using a fixed (unburnt) CPA per tree diameter as tree CPA was visibly influenced by severe wildfire. Fixed CPAs were estimated from a regression of field-measured DBHOB, top live height and crown width for 79 unburnt candidate trees (r = 0.75; P value <0.001; R2 adjusted = 0.60) to provide a representation of pre-fire crown area by tree size. These fixed-CPA LiDAR metrics were only calculated for those LiDAR-delineated trees that had also been measured for DBHOB in the field (90 trees in total). Table 1. Horizontal and vertical crown metrics assessed for A) field-assessed trees (n = 342) and for B) trees delineated in LiDAR data (n = 196) across four fire-severity types (unburnt, low, moderate, high).

A. Field metrics
Crown projection area CPA (m 2 ) Vertical projection of a tree's crown live area on the horizontal plane estimated from east-west and north-south bulk crown widths using A = πab, where a and b are half the widths (Figure 2).
Live crown width LCW (%) Bulk crown width as a proportion of total crown width (Figure 2), both as the mean of east-west and north-south measurements.
Top live height TLH (m) Measured height to the uppermost live leaf of a tree ( Figure 2); top live height and total tree height might not be equal if the tree's uppermost point is dead.

Crown cover CC (%)
Estimated as the number of first returns above 7.5 m height (i.e., the top height of the understorey) divided by the number of all first returns ( Figure 4); a higher CC indicates a higher projected cover within the crown boundary on the horizontal plane [71][72][73].

Crown density CD
Estimated from the number of all lidar points above 7.5 m height divided by the number of all returns ( Figure 4); a higher CD represents a higher density of plant material within the crown boundary above 7.5 m; crown density can be used as a measure of the amount of leaf material available for maintaining tree productivity and vigour [72][73][74]. Clumping index CI (0-1) Calculated as effective leaf area index (LAIeff) [74] over leaf area index (LAI) [78]; unity represents completely random distribution of crown returns down the profile to 7.5 m height, whereas values closer to zero indicate more aggregated leaf distribution [79][80][81].
Mean leaf area density LAD_mean (0 −1) Mean leaf-area density (LAD) of the vertical profile of each 1 m bin between 7.5 m and the total tree height; values close to 1 indicate a denser crown biomass. The leaf area density profile is calculated by dividing LAI from dz, where LAI is estimated from gap fractions at each height interval with a given thickness value of dz [30,66].
Maximum leaf area density LAD_max (0 −1) Maximum LAD of the vertical profile (as above); a higher value could indicate a fuller crown [30,66].

Coefficient of variation of LAD LAD_cv (index)
Ratio of the standard deviation of the LAD from the crown profile over the LAD mean, as one measure of the evenness of crown returns from the profile (lower value indicating less variability and greater evenness of density).
Percentile height of maximum LAD HtMaxLAD (%) Height of LAD_max as a proportion of total tree height; a higher value indicates more foliage concentrated near the top of the crown [30,66].
Percentile height of minimum LAD HtMinLAD Height of the minimum LAD as a proportion of total tree height; a lower value indicates the minimum leaf density occurs towards the bottom of the crown [30,66].  We also calculated all LiDAR metrics using a fixed (unburnt) CPA per tree diameter as tree CPA was visibly influenced by severe wildfire. Fixed CPAs were estimated from a regression of fieldmeasured DBHOB, top live height and crown width for 79 unburnt candidate trees (r = 0.75; P value <0.001; R2 adjusted = 0.60) to provide a representation of pre-fire crown area by tree size. These fixed-CPA LiDAR metrics were only calculated for those LiDAR-delineated trees that had also been measured for DBHOB in the field (90 trees in total).

Statistical Analyses
We tested the effects of wildfire severity types (unburnt, low, medium, high) on individual-tree field and LiDAR-based crown metrics using generalised linear mixed effects models (LMEM) in the Lme4 package [82] in R version 3.5.0 [83] with a plot nested within the site as a random effect. Prior to the analysis, assumptions of normality and variance homogeneity were checked, and transformed dependent variables were used as necessary (natural logarithm for height and width variables, arcsine for percentages) [84]. Crown structure can vary with DBHOB and height; therefore, we tested those variables as covariates in our LMEM of field-assessed trees.

Statistical Analyses
We tested the effects of wildfire severity types (unburnt, low, medium, high) on individual-tree field and LiDAR-based crown metrics using generalised linear mixed effects models (LMEM) in the Lme4 package [82] in R version 3.5.0 [83] with a plot nested within the site as a random effect. Prior to the analysis, assumptions of normality and variance homogeneity were checked, and transformed dependent variables were used as necessary (natural logarithm for height and width variables, arcsine for percentages) [84]. Crown structure can vary with DBHOB and height; therefore, we tested those variables as covariates in our LMEM of field-assessed trees.
Models were run for 342 field-assessed trees and for the 196 LiDAR-assessed trees, defining all the crown metrics as a response variable and fire severity types and plot ID as predictor variable. We used 'lmerTest' [85] to estimate LMEM parameters and the "car" package [86] for the Wald chi-square tests, which provided the p-value for the model. Where LMEM models indicated significant differences among fire-severity types, post hoc pair-wise comparisons were made and adjusted for multiple comparisons using Fisher's Least Significant Difference (LSD) test in the "emmeans" R package [87].
Consistent with the definition of the fire-severity types, there was no evidence of fire scarring on trees at unburnt sites, whereas all trees at high-severity sites were burnt to their total tree height compared with the means of 64% of tree height at moderate-severity sites and 25% at low-severity sites ( Table 2). In addition, the mean percentage of trees with stem epicormic and/or basal resprouts increased in the order unburnt < low severity < moderate severity < high severity ( Table 2). In terms of tree size, the mean diameters (DBHOB) of the largest trees (>50 cm DBHOB) were not significantly different among fire-severity types; however, the mean DBHOB values of the smaller mature trees (DBHOB 20 to 50 cm) were significantly lower at moderate than high severity sites ( Figure 5). In addition, the mean top live heights of trees in both the 20-50 cm and >50 cm DBHOB classes were significantly greater at high than at both low and moderate severity sites. Nonetheless, effects of DBHOB and top live height as covariates in LMEM models were consistently non-significant and were subsequently excluded. Table 2. Characterization of 342 field-measured trees by fire-severity types (means + SE). P-values indicate the significance of the fire-severity effect in a linear mixed effects model, and different superscript letters indicate significant differences among fire-severity types (Fisher's LSD Test).  The live (bulk) crown width of field-measured trees at both moderate and high severity sites was significantly less than that at unburnt and low severity sites ( Figure 6). This resulted in a significantly lower crown projection area (CPA) and percentage of live crown width for both mature and larger trees at high severity sites ( Figure 6). In contrast to the horizontal metrics, vertical metrics of field-measured trees were not clearly associated with fire-severity types, as noted above for top live height.

LiDAR-Based Metrics
LiDAR-derived estimates of crown projection area (CPA) and total tree height were highly correlated with field-based measurements (r2 = 0.72 and 0.77 respectively), although the LiDARderived CPA was consistently lower than field-measured CPA for 10 trees above 100 m 2 (Supplementary Figure S3).
The mean crown cover (CC) of trees delineated in the LiDAR data was significantly lower at moderate and high severity sites than unburnt and low severity sites (Figure 7). In contrast, the mean density within individually delineated CPAs was not clearly influenced by fire severity, while (vertical) crown evenness significantly increased from unburnt and low severity sites to moderate severity sites, and again to high severity sites (Figure 7). Trends in these LiDAR-derived metrics among fire-severity types were similar when calculated within fixed-area CPAs that estimated the size of unburnt crowns, with the exception of the clumping index, which was not clearly influenced by fire severity when calculated for individually-delineated CPAs (Figure 7) but significantly increased with fire severity when estimated for fixed-area CPAs (Supplementary Figure S4).

LiDAR-Based Metrics
LiDAR-derived estimates of crown projection area (CPA) and total tree height were highly correlated with field-based measurements (r2 = 0.72 and 0.77 respectively), although the LiDAR-derived CPA was consistently lower than field-measured CPA for 10 trees above 100 m 2 (Supplementary Figure S3).
The mean crown cover (CC) of trees delineated in the LiDAR data was significantly lower at moderate and high severity sites than unburnt and low severity sites (Figure 7). In contrast, the mean density within individually delineated CPAs was not clearly influenced by fire severity, while (vertical) crown evenness significantly increased from unburnt and low severity sites to moderate severity sites, and again to high severity sites (Figure 7). Trends in these LiDAR-derived metrics among fire-severity types were similar when calculated within fixed-area CPAs that estimated the size of unburnt crowns, with the exception of the clumping index, which was not clearly influenced by fire severity when calculated for individually-delineated CPAs (Figure 7) but significantly increased with fire severity when estimated for fixed-area CPAs (Supplementary Figure S4). Mean leaf area density (LAD) profiles indicated a more even distribution of tree crowns down the stems at high and moderate severity sites than at low severity and unburnt sites (Figure 8). This change in crown structure from full crowns at the top of the stem in the unburnt condition to a vertically elongated crown after moderate and particularly high severity fire was reflected in the LAD metrics. For example, the mean maximum LAD and the coefficient of variation of LAD were significantly decreased at high severity sites and-to a lesser degree-at moderate severity sites (Figure 8). In addition, the percentile height of maximum LAD at high and moderate severity sites was consistently lower than low severity and unburnt sites, whereas the percentile height of minimum LAD tended to be higher (Figure 8). Similar trends and significant effects in LAD metrics were also evident when calculated for CPAs that were standardized by DHBOB to estimate the prefire crown area (Supplementary Figure S5). Mean leaf area density (LAD) profiles indicated a more even distribution of tree crowns down the stems at high and moderate severity sites than at low severity and unburnt sites (Figure 8). This change in crown structure from full crowns at the top of the stem in the unburnt condition to a vertically elongated crown after moderate and particularly high severity fire was reflected in the LAD metrics. For example, the mean maximum LAD and the coefficient of variation of LAD were significantly decreased at high severity sites and-to a lesser degree-at moderate severity sites (Figure 8). In addition, the percentile height of maximum LAD at high and moderate severity sites was consistently lower than low severity and unburnt sites, whereas the percentile height of minimum LAD tended to be higher (Figure 8). Similar trends and significant effects in LAD metrics were also evident when calculated for CPAs that were standardized by DHBOB to estimate the pre-fire crown area (Supplementary Figure S5).

Discussion
Studies of the effects of fire on forest canopy structure have mainly focused on stand and landscape scales with few exceptions [88]. Our study provides LiDAR-derived metrics which are uniquely applied at the individual-tree level to quantitatively describe fire-affected tree crowns that were supported by field measurements. Our data provide clear evidence of the legacy effects of a landscape-scale wildfire on the horizontal and vertical architecture of individual eucalypt crowns in a fire-tolerant forest that have persisted for nearly a decade.
Persistent changes were apparent in the horizontal crown structure of the fire-tolerant eucalypts 8.5 years after a wildfire. The consumption of small branches in burnt crowns led to significant decreases in crown width and the projected area of both medium and large-sized trees at moderate and high severity sites. There is a lack of comparative tree-scale studies; however, our results are consistent with high-severity wildfire leading to persistently more open forest canopies at stand and landscape scales at, for example, 5 to 7 years post-fire in mixed-conifer forests of the Sierra Nevada

Discussion
Studies of the effects of fire on forest canopy structure have mainly focused on stand and landscape scales with few exceptions [88]. Our study provides LiDAR-derived metrics which are uniquely applied at the individual-tree level to quantitatively describe fire-affected tree crowns that were supported by field measurements. Our data provide clear evidence of the legacy effects of a landscape-scale wildfire on the horizontal and vertical architecture of individual eucalypt crowns in a fire-tolerant forest that have persisted for nearly a decade.
Persistent changes were apparent in the horizontal crown structure of the fire-tolerant eucalypts 8.5 years after a wildfire. The consumption of small branches in burnt crowns led to significant decreases in crown width and the projected area of both medium and large-sized trees at moderate and high severity sites. There is a lack of comparative tree-scale studies; however, our results are consistent with high-severity wildfire leading to persistently more open forest canopies at stand and landscape scales at, for example, 5 to 7 years post-fire in mixed-conifer forests of the Sierra Nevada [89]; 26 years in conifer forests of Yosemite National Park, USA [45,90]; and up to 124 years in a coastal temperate rainforest of British Columbia, Canada [42].
The vertical crown structure of fire-tolerant eucalypts was also influenced by severe wildfire. Crowns of trees burnt by high-severity wildfire remained cylindrical in shape due to post-fire epicormic resprouting down the main stem and had not recovered the full, rounded crowns that characterized unburnt trees. These post-fire elongated crowns were evident in our tree-based LAD profiles [91], and in our evenness index (EI), which averaged 0.75 in moderate severity plots, and 0.85 in high severity plots, comparatively higher than an EI of 0.71 in mixed loblolly pine forest three years after a fire in Florida, USA [77]. Despite changes in the vertical crown structure, the top live height of our field-assessed trees was not clearly influenced by fire severity, which is in contrast to the significantly lower heights of fire-tolerant eucalypt trees in Western Australia 78 years after high-intensity fire [14]. This difference could be due to many factors including variation in the pre-fire top heights among fire-severity sites in our study (data not available), and/or the much greater diameter of the surviving eucalypts in the Western Australian study (>1 m cf.~40 cm in this study), particularly since stem breakages are common in large eucalypts [92].
Crowns of severely burnt trees held a similar amount of foliage to low severity and unburnt trees despite the change in shape. This is consistent with the rapid recovery of LAI after 3 to 6 years of fire at landscape scales in Australian eucalypt forests [17], although it is important to note that our finer-scale data provide more information regarding how canopy leaves are arranged in the post-fire period. Since eucalypt leaves have a lifespan of between 1.5 and 4 years [93][94][95], our data at nearly 10 years post-fire indicate that epicormic recovery can leave a persistent impression on the crown architecture of fire-tolerant eucalypts; that is, a similar amount of foliage but concentrated closely around the main branches and down the stem. The reasons for this persistent "epicormic structure" remain unclear but could be related to the severity and extent of the Black Saturday fire and its occurrence at the end of a prolonged drought; that is, high fire severity combined with prevailing dry conditions increased tree mortality [55], which-coupled with the contraction of surviving crowns-led to increased light infiltration to the lower strata, fostering the persistence of lower leaves that might otherwise be shaded out by full canopy recovery. The effects of these changed crown shapes on forest productivity remains unclear. Leaves lower down the stem might have less access to direct sunlight [96] and thereby have lower photosynthetic efficiency than leaves higher up in the canopy. This could mean that high-severity sites are less productive for many years after fire, although it should be noted that little is known about how the net ecosystem productivity varies with time after fire in these fire-tolerant forests [97]. Other potential changes associated with an opening up of the canopy in fire-tolerant forests could influence the likelihood of subsequent fires; that is, a more open canopy would presumably increase the heat and light exposure of the understorey, changing the moisture environment of fine fuels and influencing flammability [98].
Persistent changes in post-fire eucalypt crowns could also have implications for the accuracy of tree carbon stock estimations, which are often based on allometric equations that assume a scaled relationship with tree diameter [99]. Eucalypt allometric equations have traditionally been based on unburnt trees, whereas the persistent effects in our study suggest that new allometric relationships, at least for the crown component of tree biomass, might be needed for severely burnt trees that are yet to fully recover. This, in turn, will improve the estimation and validation of tree carbon stocks at larger scales, which are often based on plot-level aggregates of individual tree estimates.
Our LiDAR-based crown metrics proved effective in describing the complex crown architecture of evergreen eucalypt trees and in quantifying the long-lasting effects of severe fire at the individual-tree level. Our horizontal and vertical crown LiDAR crown metrics will inform the interpretation of canopy patterns in LiDAR data at landscape scales; for example, we anticipate that the lower cover of trees burned by moderate and high-severity fire will contribute to greater canopy gaps and more openness at stand and landscape scales. In addition, our tree-level characterization can inform the interpretation of remotely sensed optical data since a stronger characterization of the tree layer will assist with separating the canopy from the understorey in the remote sensing analyses of post-fire recovery [17].
In turn, this will assist with a stronger characterization of forest recovery-and associated services such as carbon sequestration and water provision-at larger scales.

Conclusions
Changes in the tree crown structure after severe wildfire are crucial to understanding the potentially persistent effects of fire on tree productivity and associated processes such as water use and carbon sequestration. This study confirmed that severely burnt trees had not recovered their full, rounded crown structure 8.5 years after a fire and instead were characterized by elongated, cylindrical crowns, reflecting the persistence of prolific epicormic resprouting down the stem in the post-fire environment. Our detailed tree characterization, based on unique tree-based LiDAR metrics, provides a robust assessment for comparisons with future LiDAR data to characterize recovery in the coming years. The full recovery of tree crowns at 8.5 years in our study could have been slowed by the severity of the Black Saturday fire coupled with prevailing dry conditions. Nonetheless, climate projections for south-eastern Australia indicate drier, hotter conditions and more frequent severe fire weather [100]. Our study could thus portend an increase in the incidence of persistent legacy effects on tree crowns, leading to more open canopies and flow-on effects to the fuel structure, microclimate and function of these fire-tolerant eucalypt forests. The implications of these legacy effects on, for example, tree productivity and the accuracy of biomass allometric equations warrant further study. Trees are keystone forest structures, and so an improved assessment of changes in tree structure will improve our understanding of forest landscape processes. Our unique tree-based LiDAR crown metrics demonstrate the utility of airborne LiDAR data in assessing post-fire severity effects on structurally complex tree crowns. Quantifying the effects of contrasting fire severities using this approach fully utilizes the ability of LiDAR data to provide a strong basis for interpreting post-fire patterns in forest canopies and vegetation profiles and other remotely-sensed data at larger scales.
Supplementary Materials: The following are available online at http://www.mdpi.com/2072-4292/11/20/2433/s1, Figure S1: Overarching schematic of study design, Figure S2: Sample plot layout for field-based assessments illustrating the arrangement of three plots per site (minimum distance of 50 m between the plot centres, orientated at 30 and 90 degrees from the centre of plot "A"), Figure S3: Comparison of LiDAR-derived and field-measured crown projection areas (CPA, left) and total tree height of candidate trees (n = 90) by fire-severity type. The dashed line is the 1 to 1 line. Figure S4: Comparison of LiDAR-derived crown metrics-crown cover, crown density, evenness index and clumping index-among fire-severity types for field-based trees that were also delineated in the LiDAR data. Metrics were calculated within standardised CPAs that were estimated from relationships with the DBHOB of unburnt trees, and are based on 11 to 24 trees per severity types (UB-21, L-22, M-24, H-23), Figure S5: Comparison of LiDAR-derived metrics extracted from leaf area density (LAD) profiles-maximum LAD (LADmax), mean LAD (LADmean), coefficient of variation of LAD (LADcv), percentile height of maximum LAD (HtmaxLAD) and percentile height of minimum LAD (HtminLAD)-among fire-severity types for field-based trees that were also delineated in the LiDAR data. Metrics were calculated within standardised CPAs that were estimated from relationships with DBHOB of unburnt trees. The mean LAD profile (blue line), percentile height of the minimum LAD (triangles), and percentile height of the maximum LAD are based on 41 to 60 trees per severity type (UB-41, L-43, M-52, H-60). The values in the figure are the mean metrics, and superscript letters indicate pairwise comparison at p < 0.027 (LADmax), p = 0.01 (LADcv), p = 0.557 (LADmean), p = 0.137 (HtmaxLAD), and p < 0.001 (HtminLAD). P-values are calculated from LMEM, and different letters indicate significant differences between fire-severity types as indicated by post hoc tests (Fisher's LSD test), Table S1: Numbers of individual trees assessed in the field (candidate trees) and delineated in the LiDAR data by wildfire severity types, Table S2: Summary of LiDAR acquisition and sensor specifications Author Contributions: Y.K.K. carried out the reference study design and field work, dataset design, the LiDAR data processing, analyzed the data and drafted the original text of the manuscript. L.T.B. contributed guidance on research methodologies and forest fire, contributed to the study design and assisted in the field work. She also supervised the steps of the reference plot design, result interpretation and the discussion and conclusion, as well as reviewing the manuscript. T.D.P. assisted in the fire regime understanding and the conceptualization and supervised the modelling as well as supporting in the use of R software. C.A. contributed to minor revisions, feedback and in supervising the ecological applications as well as supporting in the use of R. All co-authors provided comments and revised the manuscript.
Funding: This research received no external funding.
Water and Planning through the Integrated Forest Ecosystem Research program. We would like to thank Julio Cesar Najera-Umana, Ben Smith and Cordula Gutekunst for their diligent fieldwork. Similarly, we would also like to thank Paul Donald Bentley and Duncan Smith for his support in LiDAR data analysis using R software. Advice on statistical analyses was provided by Cameron Patrick from the Statistical Consulting Centre at The University of Melbourne, and GPS post-correction support was provided by Ultimate Positioning Group. Finally, we would like to thank the anonymous reviewers who helped to greatly improve the manuscript by their constructive comments.

Conflicts of Interest:
The authors declare no conflict of interest.