Rapid shrub expansion in a subarctic mountain basin revealed by repeat airborne LiDAR

As a consequence of increasing temperatures, a rapid increase in shrub vegetation has occurred throughout the circumpolar North and is expected to continue. Rates of shrub expansion are highly variable, both at the regional scale and within local study areas. This study uses repeat airborne LiDAR and field surveys to measure changes in shrub vegetation cover along with landscape-scale variations in a well-studied subarctic headwater catchment in Yukon Territory, Canada. Airborne LiDAR surveys were conducted in August 2007 and 2018, whereas vegetation surveys were conducted in summer 2019. Machine learning classification algorithms were used to predict shrub presence/absence in 2018 based on rasterized LiDAR metrics, with the best-performing model applied to the 2007 LiDAR to create binary shrub cover layers to compare between survey years. Results show a 63.3% total increase in detectable shrub cover >= 0.45 m in height between 2007 and 2018, with an average yearly expansion of 5.8%. These changes were compared across terrain derivatives to quantify the influence of topography on shrub expansion. Terrain comparisons show that shrubs are located in and are preferentially expanding into lower and flatter areas near stream networks, at lower slope positions and with a higher potential for topographic wetness. Overall, the findings from this research reinforce the documented increase in pan-Arctic shrub vegetation in recent years, quantify the variation in shrub expansion over terrain derivatives at the landscape scale, and demonstrate the feasibility of using LiDAR to compare changes in shrub properties over time.


Introduction
High-latitude ecosystems have experienced substantial warming over the past 40 years, which is expected to continue (Hinzman et al 2005, Overpeck et al 1997, Tape et al 2006. Consequently, an increase in vegetation growth has occurred throughout the circumpolar North as documented in remote sensing and plot-level studies (Epstein et al 2013, Myers-Smith et al 2011, Sturm et al 2001a, Tape et al 2006. A major component of this change is shrub expansion (shrubbing or shrubbification) in arctic and subarctic ecotones (Epstein et al 2013, Tape et al 2012. The rates of shrub expansion are not uniform, both at the ecotone scale and within individual catchments (Ackerman et al 2017, Myers-Smith et al 2020, Naito & Cairns, 2011. These changes are highly variable depending on plant species, topographic position, hydrology, soils and other ecosystem properties (Davis et al 2020, Elmendorf et al 2012, Epstein et al 2013. Documenting changes in northern vegetation is critical due to plants first-order control on water, energy and nutrient balances (Ackerman et al 2017, Epstein et al 2013. For example, shrub rainfall interception loss is a major component of the water balance in tundra environments (Zwieback et al 2019). Increases in shrub cover in northern basins lead to higher snow accumulation, largely from exposed shrubs increasing aerodynamic roughness , Pohl et al 2007. As a result of complex interactions between radiation, turbulent transfer, and sensible heat fluxes, snowmelt rates are generally enhanced under shrub canopies compared to sparsely vegetated tundra with strong influences on catchment hydrology and water chemistry . The relative influence of shrubs on local ground shading and increased active layer thickness due to snowpack insolation versus their effects on large-scale climate warming feedbacks lead to complex potential permafrost responses (Blok et al 2010, Lawrence & Swenson 2011, Wilcox et al 2019. The cycling of nitrogen, phosphorus, carbon, and methane in tundra ecosystems are all affected by interactions between the abiotic and biotic influences of shrub canopies (Myers-Smith et al 2011). In addition, both expansion and densification of shrubs decreases surface reflectance and increases absorption of solar radiation, leading to a positive global warming feedback (Ménard et al 2014). Understanding how vegetation development varies under different conditions is critical for predicting the future of northern watersheds under a rapidly changing climate.
Remotely sensed data can be used along with field surveys to map vegetation over large areas and evaluate how land cover and other properties change through time (Langley et al 2001). Although coarse-to mediumresolution satellite imagery such as Landsat (30m resolution) and MODIS (250m) is useful due to its widespread coverage of the pan-Arctic, these sensors are unable to capture the heterogeneity of finer-scale vegetation change within pixels (Naito & Cairns 2011, Tape et al 2012. Optical imagery is also limited by its inability to sense height or density of vegetation and lack of cloud cover penetration (Millard & Richardson, 2013).
Airborne laser scanning, commonly referred to as LiDAR (Light Detection and Ranging) is an active remote sensing method used to characterize objects at or near the Earth's surface (Farid et al 2008, Hopkinson et al, 2006, Wehr & Lohr, 1999. LiDAR has become one of the most efficient remote sensing techniques for acquiring detailed and accurate 3D data on landscape topography and vegetation, and can make up for several shortcomings of traditional satellite remote sensing (Wu et al 2016). Although airborne LiDAR can reliably estimate forestry metrics such as height, canopy cover, and biomass, obtaining accurate estimates for lowerstature vegetation such as shrubs is more challenging due to the low height, relatively open structure, and low foliage density (Estornell et al 2011, Greaves et al 2016, Hopkinson et al 2005, Streutker & Glenn 2006. Multitemporal LiDAR has been used to evaluate temporal changes in forestry metrics, which can be logistically difficult as the structure of point clouds is dependent on the LiDAR instrument used, sensor configurations, and flight details (Bater et al 2011, Roussel et al 2017. Although there is limited research on using airborne LiDAR to evaluate temporal changes in lower-stature vegetation, Streutker & Glenn (2006) and Estornell et al. (2011) suggested that threshold-based methods could feasibly be used to compare changes in shrub cover between surveys.
This study uses a combination of repeat airborne LiDAR and field surveys to measure changes in shrub vegetation over 11 years within a well-studied subarctic mountain basin. The potential for LiDAR data to accurately capture changes in shrub height is explored along with quantifying changes in detectable shrub cover using supervised machine-learning classifications. Terrain derivatives are used for comparisons of these changes over different landscape positions. The wide range of elevation, topography, and vegetation within the study area facilitate extension of these results to other northern alpine environments. Considering the rapid changes observed in circumpolar systems, results from this study will: (1) quantify how shrub vegetation cover has changed over time in an alpine subarctic ecosystem, and (2) link these changes to ecotone and physiographic variables such as elevation, aspect, slope position, topographic wetness, and proximity to riparian areas.

Study area
Granger Basin (GB) is a ∼7.6 km 2 subarctic headwater catchment located in the west-central region of the Wolf Creek Research Basin (WCRB) (McCartney et al 2006, Carey et al 2013. The WCRB is a ∼180 km 2 long-term research facility located ∼15 km south of Whitehorse, Yukon Territory (figure 1) with over 25 years of comprehensive hydrometeorological data (Janowicz, 1999, Rasouli et al 2019. WCRB is located in the traditional territories of the Kwanlin Dün, Ta'an Kwach'an Council and Carcross/Tagish First Nations. GB straddles the subalpine and alpine tundra ecozones of WCRB, with an elevation range of 1356 to 2080 m.a.s.l. The climate of GB is characterized as continental subarctic, with mean annual precipitation of ∼400 mm of and a mean annual temperature of −3°C (Carey et al 2013). The lower basin contains a mix of low-lying grasses, willow shrubs, and dwarf birch shrubs, with taller shrubs (>2 m) along the riparian corridor and few isolated patches of white spruce (Piovano et al 2019). At upper elevations, land cover is dominated by bare rock, short tundra mosses, and grasses (Rasouli et al 2019).

LiDAR acquisition and pre-processing
Airborne LiDAR surveys of WCRB were conducted in August 2007 and August 2018. The 2007 data collected using an ALTM 3100 discrete-return sensor were projected using the NAD 1983 UTM Zone 8N coordinate system, with elevations referenced to the CGVD28 geoid. The approximate survey resolution was 0.67 returns m −2 . The 2018 survey was conducted using a Riegl Q-780 full-waveform unit, delivered in the WGS 1984 UTM Zone 8N XY coordinate system with elevations referenced to the WGS84 ellipsoid, and an average return density of 11.9/m 2 . Details on the LiDAR surveys can be found in table S1 (available online at stacks.iop.org/ERC/3/071001/mmedia).
The 2018 point cloud was transformed to match the horizontal and vertical coordinate systems of the 2007 data, as they were surveyed using different spatial references. As lower sampling rates make LiDAR pulses more likely to miss the highest points of vegetation and underestimate heights compared to field measurements (Zhao et al 2018), the resolution differences between datasets had to be accounted for. The 2018 survey was thinned to give the area covering GB an average density of 0.67 returns m −2 to match with 2007 by retaining only a specified random fraction of the original returns.
Due to the inherent difficulties associated with airborne LiDAR measurements of low-stature vegetation, accurate ground segmentation is critical (Estornell et al 2011, Hopkinson et al 2005. The adaptive TIN densification algorithm in the LAStools software package (Isenburg, 2019) performed well in terms of removing vegetation from the bare-earth surface while preserving topographic detail. After classifying each dataset into ground-and non-ground returns, the height above ground (HAG) of the resulting point clouds was computed using the lasheight function.

Field data collection
To determine the accuracy of the LiDAR, 371 differential GPS (dGPS) points were collected using a Hemisphere S320 GNSS RTK system on hard, flat surfaces such as roads and established trails. 164 dGPS points collected during the vegetation surveys and shrubline mapping were used to evaluate the LiDAR elevation accuracy on terrain within GB of varying topographic complexity. 29 transects were distributed across GB to evaluate vegetation properties. Each 40 m transect consisted of 5 circular plots with a 1 m radius, spaced out every 10 m (figure 2). Plot centroids were recorded using the dGPS and photographed before taking measurements. Vegetation height was measured at 9 regularly-spaced points throughout each plot, with the highest species at each point recorded. The percent coverage of each species present within the plot was estimated by each surveyor and averaged to a final value.
To evaluate differences between true shrub cover and potential LiDAR artifacts, the shrubline of GB was delineated using the dGPS. With the highest shrubs recorded, any non-ground returns above them could be safely assumed to be either large rocks/boulders or LiDAR positional error.

LiDAR shrub height estimation
The dGPS coordinates of each plot centroid were used to create point shapefiles, which were then buffered to a 1 m radius to represent the true plot extent. Due to the high levels of vegetation surface underestimation compared to field heights typical of shrubs (Hopkinson et al 2005), both the average and maximum of the 9 fieldmeasured heights within each plot were computed and appended to the shapefiles. Shrub species were separated into willow, dwarf birch, and dead shrub categories based on dominant percentage of plot area. The resulting shapefiles representing 145 survey plots were used to compare field measurements within them to LiDAR heights.
LiDAR returns from the height-normalized point clouds were extracted for each field survey plot radius to compare to the average and maximum field-measured vegetation heights. LiDAR return heights were extracted from the 2007, original 2018, and the 2018 point cloud after thinning it to match the 2007 survey resolution. These values were used to evaluate how well the LiDAR could estimate heights measured within each 1 m plot.

LiDAR shrub cover estimation
Along with vegetation heights, the ability of the 2018 LiDAR to estimate shrub cover compared to field data was explored. Detectable shrub presence for the coverage models was defined using a cutoff due to variable shrub cover within the survey plots and operational lower limits for determining vegetation height (Streutker & Glenn, 2006). Only plots where the average field-measured height >=45 cm (∼RMSE of 2007 LiDAR last-returns on complex terrain) were retained as shrub presence points (n=114). An additional 114 absence points were created by generating random points within GB and interpreting which did not have shrub cover through LiDAR derivatives and WorldView-2 imagery.
Binary supervised classifications of shrub presence/absence were conducted using rasterized LiDAR metrics from the thinned 2018 point cloud through the Random Forests (RF) and Support Vector Machines (SVM) machine learning algorithms within the randomForest and e1071 R packages, respectively (Breiman, 2001, Liaw and Weiner 2002, Meyer et al 2019, R Core Team, 2019. Model inputs were gridded into rasters at 2 m spatial resolution, with values taken from a 3×3 pixel moving window. LiDAR derivatives used as shrub cover predictors were the ratio of returns above 45 cm to all returns (abbreviated to dns), elevation of the interpolated bare-ground surface (DEM), and maximum intensity (maxINT).
The effect of training and validation set selection on classification accuracy was evaluated by running 25 separate classifications for each method using different random seeds (Millard & Richardson 2015). Presence and absence points were separated into a 70% training set (n=160) to train the predictive models and a 30% validation set (n=68) to assess classification accuracies. Overall independent accuracies and Cohen's kappa coefficient were computed according to an error matrix (Congalton & Green 2008, Leutner et al 2019. Class stability over each of the iterations was assessed to see which areas of GB were most frequently classified as shrubs. The classification with the best overall accuracy and class stability using the thinned 2018 point cloud was applied to the 2007 LiDAR to create raster layers representing shrub cover for both years. These were used to create new change layers representing areas of shrub gain, shrub loss, consistent shrub cover, and consistent non-shrub. Due to the presence of large rock fields in the upper GB, a small percentage of shrub-classified pixels (0.14% in 2007, 0.53% in 2018) were erroneously located above the known shrubline (see section 2.3). These areas were masked out of the final shrub layers, along with pixels located within known water bodies.

Terrain comparisons
As the locations and rates of shrub expansion vary (Naito & Cairns, 2011, Tape et al 2006, Tremblay et al 2012, terrain derivatives were created from a bare-earth DEM using SAGA GIS software (Conrad et al 2015) so that changes in shrub cover could be compared across different landscape properties. The SAGA Wetness Index (SWI) is a modified version of the classic topographic wetness index, based on the ratio of catchment area to local slope gradient defined as ln(a/tanβ) (Bohner et al 2002, Sørensen et al 2006).
Terrain raster values that fell within each of the shrub change classes as well as the entire GB were extracted for statistical analysis. One-sided pairwise Wilcoxon rank-sum tests were used to determine whether the frequency distributions of pixels in increasing or consistent shrub classes were statistically different from areas of shrub loss, consistent non-shrub, or the total GB watershed area. The terrain derivatives used to evaluate shrub expansion and the alternative hypotheses used for the Wilcoxon tests are shown in table 1.
To illustrate the magnitude of shrub change variation across landscape properties, values from each terrain derivative were binned into 5 classes roughly based on the Jenks Natural Breaks classifier (de Smith et al, 2018, Tremblay et al 2012. Slope aspect was included as an additional categorical variable after binning values into cardinal directions. The percentage of total pixels in each shrub change category falling into individual terrain value classes was tabulated along with the LiDAR survey extents within GB as a whole.

LiDAR validation and shrub height estimation
Both LiDAR surveys were within their documented positional accuracies on hard and flat surfaces (table S2). Errors increased on the complex vegetated terrain within GB. This increased RMSE of 44 cm for last-returns in the 2007 dataset was used as a cutoff for low noise in predictor variables for shrub cover classifications.
The original high-density 2018 point cloud performed well estimating maximum shrub heights measured within the field survey plots (table S3), with an average underestimation of 10% of total vegetation height and an adjusted R 2 of 0.80. Field-measured plot heights had an average overestimation of 48% and lower R 2 of 0.68 (figure 3). All direct field vegetation height measurements are found in table S4. Maximum LiDAR return heights within field survey plots are in table S5.
The relationships between field-measured plot heights and LiDAR return heights are weaker after thinning the point cloud to match the 2007 return density of 0.67 returns m −2 (table S6). 60 plots lacked a LiDAR return above the 0.45 m cutoff after thinning the 2018 LiDAR (41% of total plots) compared to 34 plots in 2007. The average field-measured vegetation height within these plots is 0.55 m compared to 0.86 m for all plots. Maximum field-measured shrub plot heights were underestimated by 62.6% on average with an R 2 of 0.55 when compared to LiDAR return heights from the thinned 2018 point cloud.
The thinned point cloud performs better when comparing to the average heights measured within the plots instead of the maximum heights (figure 4), as noted in Hopkinson et al. (2005). Average field vegetation heights were underestimated by 38.7%, with an increased R 2 of 0.64. LiDAR heights are more distributed around the 1:1 line when compared to average heights, with a linear regression slope coefficient of 1.04 (figure 4).

Shrub cover estimation
When running the 25 seed iterations, RF provided a marginally higher average overall accuracy (93.1% versus 92.7%) and kappa coefficient (0.86 versus 0.85), with a smaller range between the minimum and maximum area classified as shrub (0.51 km 2 versus 0.63 km 2 ) than SVM. The overall independent accuracy, Cohen's kappa, and stability of shrub-classified area (table S7) were considered when selecting the final classification.
After the removal of pixels above the shrubline and within lake extents, the best-performing RF classifications with 97.1% overall accuracy in 2018 gave 1.60 km 2 of total shrubbed area with heights >= 0.45 m South-facing slopes are preferential Eastness (dimensionless) East-facing slopes are preferential Euclidean distance to stream network (m) Proximity to streams is preferential Overland flow distance to stream network (m) Proximity to streams is preferential Vertical distance to stream network (m) Proximity to streams is preferential SAGA wetness index (dimensionless) Higher topographic wetness is preferential in 2018 compared to 0.98 km 2 in 2007; an overall increase of 63.3% (figure 5). In terms of areal coverage, 0.77 km 2 of GB showed an increase in shrub cover with 0.15 km 2 of loss. Consistently shrub-covered pixels (figure 6) had an area of 0.83 km 2 , while consistently shrub-free areas made up 5.41 km 2 . All accuracies and shrub-classified areas predicted using the thinned 2018 LiDAR varied depending on the random seed for splitting points into the training and validation sets. Over the 25 classifications for each   algorithm, overall accuracies varied by 10.29% with RF and 11.76% with SVM, while the amount of shrubclassified area varied by 0.51 km 2 and 0.63 km 2 respectively (table S7). Additional measures of uncertainty were estimated using the top 5 highest performing classifications according to overall independent accuracy. The average amount of shrub expansion was 0.66 km 2 with a coefficient of variation of 0.08, resulting in an average percent increase of 67.8% and corresponding CV of 0.17.

Landscape-scale variation
The Wilcoxon rank-sum tests (table 2) indicate that for all terrain derivatives, consistent and expanding shrub pixels are preferentially located in areas that fulfill the alternative hypotheses when compared to the other change classes. For the slope gradient, northness, stream distance, and relative slope position, consistent shrub pixels had lower value distributions than shrub gain pixels, which were then lower than each other class. Elevation values in pixels exhibiting shrub gain were lower than in all other change classes. Shrub loss pixels also had more preferential terrain derivative values than areas of consistent non-shrub or the entirety of GB in all instances. For example, SWI value distributions in the shrub loss category were lower than in the consistent shrub or gain pixels, but greater than in the consistent non-shrub or the overall GB.  Greater Gain pixels>consistent non-shrub Consistent shrub pixels>gain, loss, consistent non-shrub Stream distance (overland flow) Less Gain pixels<loss, consistent non-shrub Consistent shrub pixels<gain, loss, consistent non-shrub SAGA wetness index Greater Gain pixels>loss, consistent non-shrub Consistent shrub pixels>gain, loss, consistent non-shrub * indicates p<0.05, all other p-values<0.01 As with the Wilcoxon tests, comparisons across individual terrain value classes (table s8) showed that consistent shrub pixels have the largest differences compared to GB overall, followed by areas of increasing and then decreasing shrub presence. For example, while only 31.9% of GB is within 500 m of the stream network, 63.7% and 66.6% of the consistently shrub-covered and increasing shrub areas respectively are within this distance class. Conversely, consistent and expanding shrub areas are preferentially located in areas of high topographic wetness, with 37.7% and 24.8% of these pixels having SWI values>5 compared to only 12.7% of the overall GB. While table 2 shows that value distributions of consistent and expanding shrub pixels have statistically higher eastness and lower northness than other classes, the differences were minimal when aspect was represented as cardinal direction (table s8). These results indicate that shrubs in GB are found in and are preferentially expanding into lower and flatter areas near stream networks, with lower slope positions and a higher potential for topographic wetness. Areas of consistent shrub cover, shrub gain, and shrub loss are visualized along with elevation, slope gradient, and SWI in figure 7.

Discussion
The relationships between original 2018 LiDAR and field-measured vegetation heights are stronger than in similar studies estimating shrub heights (Estornell et al 2011, Hopkinson et al 2005, Streutker & Glenn, 2006, largely as a result of the high return densities (11.9 returns/m 2 ). However, after thinning the point cloud to allow direct comparison with the 2007 survey, we were unable to calculate exact shrub height changes. However, LiDAR with similar densities to the 2018 survey should be sufficient for detecting changes in height in the future.
The average yearly increase in detectable shrub cover >= 0.45 m in height was 5.8%, which is greater than previous works on the expansion of shrub vegetation in northern regions (Naito & Cairns, 2011, Tape et al 2006, Tape et al 2012, Tremblay et al 2012, Davis et al 2020. In the Brooks Range and Northern Slope of Alaska, maximum yearly increases of 3.5% and 2.5% were found within terraces and floodplain shrub patches between the 1950s and 2000s (Naito & Cairns 2011, Tape et al 2006. Modelled shrub stem elongation rates in northern Labrador have been found to rapidly increase in recent years, from 1 cm year −1 in 2005 to 13.25 cm year −1 in 2014 (Davis et al 2020).
Classification accuracies and shrub-covered areas varied depending on training and validation set selection (table S7). The differences in accuracy between classifications were much more sensitive to which points were used for training and validation than the algorithm. As the average heights measured in the shrub presence points ranged from 0.45 m to 3.22 m, a large range of predictor variable values are captured within the model training and validation sets. The variations in area could be due to the different characteristics of the points used in the classifiers depending on the random seed.

Landscape-scale variation
The preferential location of expanding and consistent shrubs in areas of lower elevation, lower slope gradient, higher topographic wetness, and higher stream proximity are all consistent with previous findings and supported by established physical processes (Tape et al 2006, Naito & Cairns, 2011, Tremblay et al 2012, Davis et al 2020. As slopes in subarctic environments represent a gradient of hydrology, climate, nutrient and moisture availability, their influence on environmental phenomena result in vegetation patterns that are strongly correlated with relief (Walker et al. 2003). The association of high SWI values with shrub cover can be explained by its determination of nutrient resources influencing the distribution of woody plants (Wu and Archer 2005, Walker et al. 2003, Naito & Cairns 2011. Riparian areas are similarly higher-resource environments, which would explain the preferential location of consistent and expanding shrubs closer to streams (Myers-Smith et al 2011, Naito & Cairns, 2011, Tape et al 2012. The large elevation range within Granger Basin creates a gradient of climatic and hydrometeorological conditions (McCartney et al 2006), which would contribute to the contrast in shrub changes across elevation classes.

Implications for Granger Basin
The rapid rate of shrub expansion within GB has a wide range of implications for hydrological processes. Rainfall interception losses are a major component of the water budget in arctic shrub ecosystems, with effective rainfall reduced by up to 30% under birch canopies (Zwieback et al 2019). The impacts of increased shrub cover on the exchange of energy among atmosphere, vegetation, and soils are varied due to interactions between shrub canopies and snow cover (Myers-Smith et al 2011). Due to the increased aerodynamic roughness of exposed shrubs, a greater amount of canopy coverage will lead to increased snow storage . The insulative effects of these compacted snow layers can lead to winter soil temperatures under shrub canopies being up to 30°C greater than air temperatures (Sturm et al 2001b). Although the specific influence of shrub expansion on the response of the discontinuous permafrost within Granger Basin to warming are complex, an increase in thawing would impact groundwater and hydrological responses to change (Bonfils et al 2012, Carey et al 2013. While some works suggest that shrubs may preserve permafrost through ground shading, low birch shrubs have been found to increase frost table depth compared to bare tundra, likely due to the advancement of snowmelt timing (Blok et al 2010, Wilcox et al 2019. Shifts in Arctic and subarctic vegetation under future warming are projected to have an overall positive feedback to climate due to changes in the energy balance (Pearson et al 2013). Although increases in shrub biomass would sequester atmospheric carbon to a certain extent, the decreases in albedo and resulting greater summer sensible heat flux due to increased shrub cover will ultimately cause a positive atmospheric heating effect (Myers-Smith et al 2011, Lafleur & Humphreys 2018). This positive feedback on warming will have further cumulative effects on the energy balance and future climate of GB, as well as other similar northern catchments (Ménard et al 2014).

Areas for future research
Future studies should use coincident field surveys along with each LiDAR acquisition to ensure that vegetation metrics for each year have the best possible temporal correspondence. The LiDAR instruments used, acquisition parameters, and flight conditions should be standardized between surveys. Greaves et al (2016) used canopy metrics derived from high-resolution LiDAR and multispectral imagery to predict shrub biomass with high accuracy, while lower-resolution LiDAR has been used to estimate fractional vegetation cover in forest ecozones (Hopkinson & Chasmer 2009).
With increased LiDAR survey resolutions, better standardization between acquisitions, and incorporation of other high-resolution remotely sensed data, future works could potentially distinguish between in-filling, increasing shrub size, and new colonization, as well as measure changes in vegetation properties other than cover

Conclusions
This study has used multi-temporal LiDAR along with field vegetation surveys to evaluate shrub expansion between 2007 and 2018 in a well-studied subarctic mountain basin, and quantitatively compared how changes in vegetation cover varied over a range of landscape properties. Our results show that the use of airborne LiDAR is feasible for comparing shrub changes over time, though the potential for accurately measuring changes in properties other than cover are limited by return density. Shrubs were found to be preferentially located in and expanding into lower and flatter areas near stream networks, at lower slope positions and with a higher potential for topographic wetness. A greater understanding of how such changes vary at the landscape scale over different physiographic variables is crucial for predicting the future of northern watersheds under a rapidly changing climate.