## 2 Materials and Methods
### Analysis
#### 2.3.2 Lidar-Derived Canopy Structure  
To ensure alignment of bi-temporal lidar data, point cloud elevations for the 2000 dataset were block-adjusted to the corresponding 2016 and 2018 datasets by comparing point clouds over unchanging open terrain features. This corrected minor elevation biases (0.3-0.7 m) likely caused by GNSS trajectory error in the earlier generation lidar sensor (typically included in most lidar data processing steps, e.g., [PERSON] et al. [103]), assuming higher vertical accuracy in the newer datasets. Last and single return points  
Figure 2: The analytical workflow. L2 (Level 2); MNDWI (Modified Normalized Difference Water Index); NDVI (Normalized Difference Vegetation Index); EVI (Enhanced Vegetation Index); LST (Land Surface Temperature); DEM (Digital Elevation Model); DSM (Digital Surface Model); CHM (Canopy Height Model); DUC (Ducks Unlimited Canada). Bi-temporal CHM changes refer to observed differences in canopy height (>1 m) between 2000 and 2016/18, while ecotnal change refers to spatial extent variations in the 90 th percentile CHM change for the same period (see text). were then classified into ground and non-ground in TerraScan, with first and intermediate returns assigned to the non-ground class. Ground points were interpolated to a 5 m digital elevation model (DEM) using Tri-angulation with Linear Interpolation. All remaining data points were used to create a 5 m digital surface model (DSM) using a local maximum height filter in Surfer (Golden Software Inc., Golden, CO, USA). A 5 m resolution CHM was constructed by subtracting the DEM from the DSM and then compared through time, hereafter referred to as CHM5. Low CHM heights (<1 m) were ignored in subsequent analyses to mitigate influences of seasonal vegetation phenology variations. To compare with 30 m Landsat VIs, CHM5 was aggregated to 30 m resolution (CHM30) using a block mean in ArcGIS, where each 30 m pixel represents the average of the underlying 36 (5 m) pixels. This provides an index of both height and canopy cover. Changes in canopy heights >1 m were then calculated by subtracting the year 2000 CHM30 from 2016 and 2018 CHM30. To investigate spatial correspondence between ecotone expansion or contraction and any associated VI trends over the study period, tall canopy height grid cells were targeted and smoothed to emphasize locally contiguous regions of shrub and tree, or woody biomass change. Specifically, the 90 th percentile of canopy height (P90 CHM) was derived from the CHM5 using a 9 \(\times\) 9 (45 m \(\times\) 45 m) grid cell moving window. To determine the optimal window size, four window sizes from 5 \(\times\) 5 to 11 \(\times\) 11 were tested and 9 \(\times\) 9 was chosen based on visual patterns of canopy continuity and similarity compared with the 30 m resolution Landsat VIs (Figure 3). P90 was used as this correlates well with field-measured tree height at the plot-level [78]. Figure 3: Schematic of ecotonal expansion and contraction using 90 th percentile canopy height (P90 CHM). Change detection shows forest expansion (non-forest to forest) and contraction (forest to non-forest) between 2000 and 2016/18. To define mature forest ecotone boundaries for each time period, a height threshold approach was adopted to differentiate between taller forest vegetation and all other classes of herbs, shrubs or juvenile trees. According to [PERSON] et al. (2014) and [PERSON] (2018), willlow shrubs--the tallest non-forest landcover type in the PAD--can reach up to 5 m tall, while mature forest species tend to be many meters taller. Therefore, P90 CHM values were classified into short (<5 m) and tall (>5 m) vegetation. For the purpose of our analysis, the shared tall/short class boundary was assumed to be a proxy for the forest ecotone boundary. A categorical change detection analysis was then performed to illustrate the changes occurring in the forest and non-forest class either side of the boundary, as well as allowing a spatial assessment of the change in the ecotone boundary position through time (Figure 3).