Evidence of vegetation greening at alpine treeline ecotones: three decades of Landsat spectral trends informed by lidar-derived vertical structure

Monitoring changes in vegetation at high-latitude and alpine treeline ecotones is critical for characterizing changes to carbon and energy budgets, plant species richness, and habitat suitability and is often considered a bellwether of a changing climate. Herein, we used transects of airborne laser scanning (ALS) data to identify alpine treeline ecotones in the Yukon Territory of Canada, and assessed changes in vegetation greenness using a time-series of Landsat imagery over a 30 year period from 1985 to 2015. Specifically, we calculated the enhanced vegetation index (EVI) from annual Landsat composites and assessed temporal trends within 500 m of detected forest-lines (i.e., transition point from continuous forest into treeline ecotones) using Theil–Sen’s nonparametric regression. Across 74 detected treeline ecotones, 27.5% of Landsat pixels displayed a significant positive trend in EVI and 5.6% of pixels displayed a significant negative trend (p < 0.05). By using ALS data to determine vegetation structural class, we found that non-treed pixels had the highest percentage of significant positive trends in vegetation greenness (40.8%), followed by shrubs (30.5%), with lower percentages in sparse forests (18.9%) and open/dense forests (13.3%). These results suggest herbaceous and shrub vegetation types are undergoing the most significant changes in greenness, likely due to increases in shrub cover and herbaceous biomass in areas associated with these alpine treeline ecotones. The limited increases in EVI in forests likely indicates that vegetation cover is changing less rapidly in forests than in shrub and herbaceous vegetation types. Moreover, EVI may not be capturing increased height growth in forests near the treeline. Combining ALS data and Landsat time-series data provides a useful approach to locate and characterize alpine treeline ecotones, and enables the direct assessment of which vegetation structural classes are experiencing the greatest greening trends, thereby providing new insights to ecosystem change.


Introduction
Changes in vegetation at both high-latitude and alpine tundra ecosystems have been well documented in recent years due to warmer temperatures, longer growing seasons, and longer snow-free periods (Harsch et al 2009, Wilson and Nilsson 2009, Danby et al 2011, Myers-Smith et al 2011, Villarreal et al 2012. Research has documented shifts in the limits to tree establishment, both in elevation and latitude. Harsch et al (2009) documented tree advancement at 52% of 166 globally distributed high-latitude and alpine tundra sites. Additionally, increases in shrub cover beyond the treeline have been observed (Cannone et al 2007, Dial et al 2007, Wilson and Nilsson 2009, Myers-Smith et al 2011, which can potentially lead to a decrease in alpine species richness (Wilson andNilsson 2009, Pajunen et al 2011), as shade intolerant herbaceous species are outcompeted by shrub species. Alternatively, in areas with minimal changes in shrub cover, herbaceous biomass and alpine species richness have been shown to increase with warmer temperatures and longer growing seasons (Danby et al 2011, Villarreal et al 2012. Changes in tundra ecosystems can have profound impacts on both carbon and energy budgets, as well as impacts on habitat suitability (Harsch and Bader 2011, Myers-Smith et al 2011, Epstein et al 2013. Specifically, the movement of trees and shrubs into ecosystems previously dominated by herbaceous vegetation will increase carbon storage in aboveground biomass. Carbon storage in aboveground biomass will also increase in areas not colonized by trees and shrubs if herbaceous biomass increases (Danby et al 2011). The replacement of herbaceous vegetation with trees and shrubs will also impact regional energy budgets, as greater tree and shrub cover could result in increased radiation absorption and regional warming (Harsch and Bader 2011). Further, as herbaceous plant species are a critical food source for a range of species in alpine and high-latitude environments (McIntire andHik 2005, Munro et al 2006), the reduction in herbaceous plant diversity as herbaceous species are replaced by trees and shrubs could reduce habitat suitability for some species in high-latitude and alpine environments. While monitoring changes to alpine and highlatitude environments is critical for understanding changes to carbon and energy budgets, as well as habitat suitability, monitoring these changes over large areas can be difficult due to the vastness and remoteness of many tundra environments, and the need for multi-temporal observations (Beck andGoetz 2011, Andrew et al 2012).
Measurements of vegetation productivity derived from satellite imagery, most notably the normalized difference vegetation index (NDVI), have been used extensively to monitor changes in tundra vegetation over both regional and continental scales (Beck and Goetz 2011, Fraser et al 2011, McManus et al 2012, Walker et al 2012, Ju and Masek 2016, Carlson et al 2017. Beck and Goetz (2011), for example, observed greening trends across the high-latitude tundra from 1982 to 2008 using NDVI at 8 km spatial resolution from the advanced very high resolution radiometer. While studies using coarse-resolution remote sensing products are critical to establish continental and global trends, the spatial coarseness of these products severely limits finer understanding of the changes occurring across the landscape. Alternatively, McManus et al (2012) used 30 m spatial resolution Landsat data in combination with land cover information to assess changes in tundra vegetation between 1986 and 2010 in Quebec, Canada. They found that NDVI increased prominently in shrub and graminoid land cover classes, suggesting increases in shrub cover and herbaceous biomass. Similarly, in an alpine context, Carlson et al (2017) used Landsat data in combination with a land cover map in the European Alps, and found the most significant increases in NDVI in rocky habitats, as well as grasslands and lowshrubs.
While spectral information from Landsat can be used to assess vegetation change, the lack of threedimensional information from passive optical data limits our ability to assess changes in vegetation structure with Landsat data alone (Goetz and Dubayah 2011). Alternatively, airborne laser scanning (ALS), an active remote sensing technology, enables three-dimensional vegetation structure to be characterized (Lim et al 2003), providing a tool to assess vegetation structure across alpine and high-latitude treelines (Rees 2007, Ørka et al 2012, Coops et al 2013. A range of structural attributes can be estimated from ALS data that are of interest to researchers studying changing tundra ecosystems, such as the proportion of tree cover and the proportion of shrub cover. For example, Ørka et al (2012) delineated forest, subalpine, and alpine vegetation classes in Norway by calculating measures of forest cover (% cover>5 m in height) and shrub cover (% cover>0.5 m in height) directly from ALS data. However, while ALS data provides a detailed snapshot of structure at one time period for the areas sampled, ALS data is not typically multi-temporal, limiting our ability to assess temporal changes to structure at treeline ecotones.
By combining structural information from ALS data with temporal information from Landsat, changes in vegetation condition through time can be linked to specific structural classes from ALS data, allowing changes in vegetation structure to be assessed through time. Here, our objective is to investigate 30 years of change at alpine treeline ecotones in the northwestern Canadian boreal using a combination of ALS data and Landsat time-series data. Specifically, we use ALS data to first locate treeline ecotones where forests transition to shrubs and alpine vegetation, and use Landsat timeseries to assess which vegetation structural classes are experiencing the most change in these ecotones. By characterizing the structural classes that are experiencing spectral changes during the Landsat record, we can better understand how Canada's extensive mountain treeline ecotones are being impacted by a changing climate, allowing impacts on carbon and energy budgets, as well as impacts on alpine habitats to be assessed over large regions.

Study area
The Canadian boreal spans over 550 million ha, of which 270 million ha is treed (Brandt et al 2013), and is dominated by cold-tolerant coniferous species, such as black spruce (Picea mariana), white spruce (Picea glauca), and balsam fir (Abies balsamea). The majority of the northern boreal is 'de facto' protected, as access to these forests is limited (Andrew et al 2012), resulting in forested ecosystems that are dominated by natural disturbance and recovery dynamics . This study focused on the northwest of the Canadian boreal, primarily in the Yukon and the Boreal Cordillera ecozone (figure 1). The Boreal Cordillera is dominated by mountainous terrain, with trees such as black and white spruce growing on lower mountain slopes and transitioning to shrubs, herbs, moss, and lichen above alpine treelines (Ecological Stratification Working Group 1995).

Data sources 2.2.1. ALS transects
In the summer of 2010, 34 transects of small-footprint ALS data were collected over the Canadian boreal by the Canadian Forest Service, in collaboration with the Canadian Consortium for Lidar Environmental Applications Research and the Applied Geomatics Research Group. The data were collected by an Optech ALTM 3100 discrete return sensor, with a pulse repetition frequency of 70 kHz and a fixed scan angle of 15°between the altitudes of 450 and 1900 m. The transects had an average pulse density of 2.8 returns m −2 , a swath width ranging from 400 to 800 m, and totaled ∼25 000 km in length . The classification of points into ground and non-ground returns, along with other preprocessing steps, was completed using customized software tools designed to deal with large transect files (Hopkinson et al 2011). From these ALS data, a canopy height model (CHM) was calculated at a 2 m spatial resolution across all transects.

Landsat best available pixel (BAP) composites
Landsat time-series data and information on forest disturbance were derived from gap-filled Landsat surface reflectance composites for 1984-2016 produced following the Composite 2 Change (C2C) approach (Hermosilla et al 2016. Specifically, BAP image composites were first produced from Landsat imagery by selecting observations for each pixel within a specific date range (August 1+/−30 days) based on the scoring functions defined by White et al (2014) and Hermosilla et al (2016), which rank the presence and distance to clouds and their shadows, the atmospheric quality, and the acquisition sensor. Next, these image composites were further refined by removing noisy observations (e.g., haze and smoke) and infilling data gaps using spectral trend analysis of pixel time-series (Hermosilla et al 2015). Implementation of these rules and subsequent temporal analysis allows for the production of seamless annual surface reflectance composites for all of Canada, as well as the detection and characterization of forest change events (refer to Hermosilla et al 2016 for more details). Availing upon the BAP image composites and associated attributed change information produced by C2C, an annual, change informed data cube of land cover has also been produced (Hermosilla et al 2018).

Identifying alpine treeline examples
In order to investigate changes at alpine treeline ecotones, treeline ecotones were first located using the 2010 ALS transects. To prevent confusion between treeless alpine environments and recently disturbed forests, areas that were harvested or burned between 1985 and 2015 were first masked using the C2C Landsat-derived disturbance dataset. Additionally, areas detected as burned between 1960 and 1985 in the Canadian National Fire Database (CNFDB, http:// cwfis.cfs.nrcan.gc.ca/ha/nfdb) were also masked, along with areas classified as water or wetland in a Landsat-derived land cover map . Spatial segments of the ALS transects were then selected by identifying portions of transects that were uninterrupted by disturbance, water, or wetland, and crossed alpine terrain. In total, 11 segments were selected across the Yukon, ranging in length from 8 to 80 km, with all samples including multiple forest/ alpine interfaces (figure 1).
In order to determine the location of treeline ecotones with ALS data, we first developed a structural definition of the ecotone. Based on a definition developed by Holtmeier and Broll (2005), Harsch and Bader (2011) define the treeline ecotone as 'an ecotone delimited at the upper end by the tree species limit, the uppermost elevation or latitude at which tree species occur as trees or krummholz regardless of height, and at the lower end by continuous forest >3 m tall'. While detecting the upper end of this ecotone would be difficult with ALS data (i.e., 'the uppermost elevation at which tree species occur as trees or krummholz regardless of height'), detecting the lower end of this ecotone is well suited to ALS data (i.e., the transition between continuous forest and the treeline ecotone). Therefore, our approach aimed to detect this interface between continuous forest and the treeline ecotone, which Bryn and Potthoff (2018) refer to as the 'forestline'. As 10% tree cover is commonly used to define a forest (FAO 2006), we defined the forest-line as the location at which tree cover above 3 m drops below 10% (figure 2). Here, we calculated tree cover above 3 m as the percentage of CHM cells above 3 m in height within each 30 m pixel. The CHM was derived from ALS data at 2 m spatial resolution, with each 2 by 2 m cell representing the height of the tallest vegetation in the cell. To search for a forest-line along the ALS segments, a moving window of 300 m (swath width) by 500 m (fight direction) was passed over the ALS segments, and a forest-line was identified as the location at which 90% of the 30 m cells dropped below 10% tree cover above 3 m. In order to identify a second forest-line along each ALS segment, the ALS segment first needed to pass back through a forested area (90% of the CHM cells in the 300 by 500 m area with tree cover >10%). The moving window was then run in the reverse direction along the transects to identify transitions from forest to non-forest that occurred in the opposite direction of the original search. Once transitions from continuous forest into the treeline ecotone were identified, a window of 500 m on either side of the forest-line was identified to capture the treeline ecotone and the forests immediately surrounding each treeline ecotone (figure 2). We recognize that a 500 m buffer may not capture the entire area of the treeline ecotone in all cases. However, as the transition from forest to alpine vegetation is often abrupt (Harsch and Bader 2011), the entirety of most treeline ecotones was likely captured by this buffer.

Assessing 30 year spectral trends in alpine ecotones
For each forest-line identified, spectral changes in the enhanced vegetation index (EVI) were investigated for each 30 m pixel within 500 m of a forest-line from 1985 to 2015. EVI was calculated for each Landsat composite as: NIR Red NIR 6 Red 7.5 Blue 1 .
The EVI was developed to reduce the influence of atmospheric anomalies and understorey vegetation that are a factor in the application of NDVI (Huete et al 2002), as well as to provide improved sensitivity in areas with high biomass (Jiang et al 2008). Additionally, Xiao et al (2005) found that EVI was more sensitive to seasonal dynamics in plant productivity than NDVI, suggesting that EVI will be more sensitive to changes in plant productivity in treeline ecotones through time. We examined trends in the EVI in the following structural classes defined using the 2010 ALS data: dense/open forest (>30% tree cover), sparse forest (10%-30% tree cover), shrub, or non-treed. Tree cover was determined as the percentage of CHM cells >3 m. Shrub was determined as pixels that had >10% cover above 0.5 m, but <10% cover above 3 m, while non-treed were defined as <10% cover above 0.5 m (

Results
In total, 74 alpine forest-lines were identified along the ALS transects, ranging in elevation from 850 m above sea level (asl) to 1600 m asl, with the majority of forestlines occurring between 1200 and 1600 m asl (figure 3). Within 500 m of the detected forest-lines, 22.2% of pixels were dense or open forest (5204 pixels), 19.4% were sparse forest (4559 pixels), 28.4% were shrubs (6661 pixels), and 30.0% were non-treed (7025 pixels) according to the ALS structural classification. Across all sampled alpine forest-lines, 27.5% of pixels displayed a significant positive trend in EVI and 5.6% of pixels displayed a significant negative trend (table 2). The non-treed class had the highest percentage of pixels with significant positive EVI trends (40.8%), followed by shrubs (30.5%), and lower percentages for both sparse (18.9%) and open/dense forest (13.3%, table 2, figure 4). Dense/open forests had the highest proportion of significant negative slopes (13.8%), followed by 5.7% for sparse forests, 3.4% for shrubs, and only 1.5% for non-treed. Figure 5 displays the proportion of significant trends near each of the 74 alpine forest-lines for forested (dense/open and sparse) and non-forested (shrub and non-treed) pixels. Significant positive trends were more consistent for non-forested pixels, as 48 of 74 forest-lines had significant positive trends in at least 20% of non-forested pixels. Alternatively, only 21 forest-lines had significant positive trends in at least 20% of forested pixels. While negative trends in EVI were more prevalent in the forested pixels, the majority of these significant negative trends occurred at a small subset of forest-lines (60% of the negative slopes Table 1. Structural classes used to investigate temporal trends in Landsat data. Structural class was assessed using ALS derived canopy height models (CHMs).

Structure class Description
Dense/Open forest >30% CHM cells above 3 m Sparse forest 10%-30% CHM cells above 3 m Shrub <10% CHM cells above 3 m, but >10% above 0.5 m Non-treed <10% CHM cells above 0.5 m occurred in eight out of 74 forest-line examples, concentrated in the southeast portion of the study area). Figure 6 displays the EVI slope for each pixel against the EVI value at the start of the time-series (i.e., median EVI from 1985 to 1989). The majority of pixels for shrubs and non-treed pixels had an EVI slope above zero (mean slope=0.0008 and 0.0012 for shrubs and non-treed, respectively), while forested pixels tended to cluster closer to zero (mean slope=−0.0001 and 0.0004 for dense/open and sparse, respectively). The majority of negative slopes corresponded to dense/open pixels with high initial EVI values.
To investigate the spatial relationship between positive EVI trends and the forest-line location, figure 7 displays the positive EVI trends as a function of distance to the forest-line, wherein only forest-lines displaying positive trends in >20% of the pixels are shown, representing 44% of the forest-lines. Negative distances correspond to the forested side of the forestline, while positive distances correspond to the nonforested side of the forest-line. The percentage of positive trends in vegetation greenness was higher on the non-forested side (54% of pixels with positive trends) than the forested side (32% of pixels with significant positive trends). On the forested side, the percentage of significant positive slopes was highest in pixels that were closest to the forest-line.

Discussion
Our results demonstrate the utility of combining ALS and Landsat time-series data to examine changes in alpine tundra ecosystems. Shrub and non-treed structural classes showed the most significant increases in greenness, suggesting an increase in shrub density and herbaceous vegetation cover in the alpine tundra ecotone. These findings are in agreement with other studies in the Yukon, as Danby et al (2011) reported an increase in alpine herbaceous vegetation biomass in the Ruby Range Mountains between 1968and 2008, and Myers-Smith (2011 observed an expansion of shrubs in an alpine tundra environment of southwest Yukon. Similarly, in Sweden, Hallinger et al (2010) documented increased shrub growth (both vertical and radial) in alpine tundras as well as shrub colonization at higher elevations.
Our results are in agreement with past remote sensing studies that measure vegetation greening at treeline ecotones. Carlson et al (2017) found a strong greening signal in rocky habitats, grasslands, and lowshrubs within the alpine of the Alps, while forested classes showed minimal greening. Similarly, McManus et al (2012) found the most evidence of greening in shrub and graminoid classes at high-latitudes in Quebec. While previous studies relied on land cover products to differentiate greening trends, this current study was able to locate the treeline ecotone directly using vegetation structure information from ALS data, and investigate how greening trends varied as a function of distance to the forest-line and as a function of structural class. Understanding which structural classes are experiencing greening trends is critical for understanding the implications of the ecosystem change. For example, greening of herbaceous alpine vegetation may signal an increase in food availability, as herbaceous plants are a critical food source for a range of species (McIntire andHik 2005, Munro et al 2006), while the encroachment of shrubs or trees may signal a decrease in herbaceous species richness, and therefore a potential decrease in available food (Wilson andNilsson 2009, Pajunen et al 2011).
While an advancement in treeline has been documented in the Yukon (Danby and Hik 2007) and Alaska (Dial et al 2007), as well as in Scandinavia (Bryn and Potthoff 2018), minimal evidence of greening in forested structural classes was observed in this analysis. As Danby and Hik (2007) found, most recruitment Table 2. Percentage of pixels with significant (p < 0.05) positive and negative Theil-Sen slopes in EVI by structural class. The value in parenthesis is the total number of pixels.  of trees at the upper limit of the forest occurred in the early to mid 20th century, indicating that it is possible that the highest elevation trees in this study were established prior to the start of our time-series in 1985, preventing strong evidence of greening in the treed classes. However, even if vertical changes in the forested classes occurred (e.g., increased growth rates), these changes may not be accompanied by strong increases in EVI, as changes in vegetation height will exhibit a smaller spectral change than changes in vegetation cover. Therefore, the lower levels of changes in EVI for the forested classes should not be interpreted as a lack of vegetation change, but rather, should be interpreted as less change in vegetation cover relative to the shrub and non-treed classes. Ultimately having multi-temporal acquisitions of ALS data at treeline ecotones could allow for a more informed assessment of the vertical forest structure (e.g., Nelson 2007, Hauglin et al 2018). Kambo and Danby (2018), working in sub-arctic treelines in the Yukon, explored the drivers of tree seedling establishment at alpine treelines. The authors found that the occurrence of tree seedlings across 640 sites was best explained by the proximity, height, and upslope orientation of shrubs. Shrub cover was found to influence the establishment of seedlings and thus the authors concluded that changes in shrub distributions with climate change will in turn influence the establishment of tree seedlings.
While the greening signal was strongest in the non-treed class, it is possible that dwarf shrubs contributed to this greening trend, as dwarf shrubs range in height from 0 to 0.4 m (Myers-Smith et al 2011), and would therefore not be captured by ALS estimates of shrub cover above 0.5 m in height. Dwarf shrubs can grow laterally and form a dense layer of vegetation along the ground surface (Myers-Smith et al 2011), which can negatively impact growing conditions for a range of herbaceous plants, and potentially reduce herbaceous species richness. Therefore, it is important to note that greening within the non-treed class does not necessarily indicate enhanced herbaceous vegetation growth, but could also indicate an increase in dwarf shrub cover.
As captured in figure 7, the detected forest-lines represented a clear divide between weaker greening signals in forests and the stronger greening signals in alpine vegetation, demonstrating the feasibility of identifying forest/alpine treeline ecotones with ALS data. Identifying the location of these ecotones is critical in order to understand the spatial relationship between greening trends and potential structural changes occurring on the landscape, and to monitor changes in treelines over time (Ørka et al 2012).
Changes in seasonality through the time-series could have implications on the greening trends observed. The Landsat compositing approached used in this study takes the best available pixel within 30 days of August 1st for each year. As the number of growing days is short in alpine environments due to late snow melt (Galen and Stanton 1995, Fu et al 2014, these alpine environments may not reach peak greenness by the start of the compositing period (early July). Therefore, the compositing approach is more likely to capture peak greenness if the growing season is lengthening due to warming temperatures and accompanied by an earlier snow melt. While changes in the nontreed and shrub classes are clear, these changes may represent an increase in growing season length in addition to increases in vegetation density and cover Klein 2011, Ernakovich et al 2014).
The timing of the ALS acquisition (2010) in relation to the time period of analysis  will also impact our findings. As the vegetation classes considered herein are dynamic, some forest pixels that exhibited positive EVI trends, for example, may have transitioned from shrubs in 1985 to forest by 2010. Alternatively, some herbaceous pixels with positive trends may have transitioned from herbaceous in 2010 to shrub by 2015. Therefore, it is important to recognize that greening trends may represent multiple structural changes within each class (e.g., increased herbaceous biomass versus transition from herbaceous to shrub).

Conclusions
Combining ALS data and Landsat time-series allows for a powerful assessment of changes in vegetation greenness near alpine and high-latitude treeline ecotones. When using spectral indices to assess trends in vegetation greenness, it is critical to identify the vegetation structural classes that are experiencing greening trends in order to understand the implications for ecosystem change and for the carbon cycle. Our results suggest that herbaceous and shrub classes are experiencing the greatest amount of change, signaling a potential increase in shrub cover and herbaceous biomass. However, it is important to note that vertical changes in forest classes may also be occurring, albeit at a slow rate, but these changes are not accompanied by the same significant increases in EVI found in the herbaceous and shrub classes. Future studies should utilize multi-temporal ALS acquisitions in order to directly assess the changes occurring in these forest classes near treeline ecotones. Herein we have demonstrated that single-date acquisitions of ALS data can provide a valuable source of information for locating treeline ecotones in alpine environments, and provide insights to the structural classes in which greening trends are being observed. The methods and logic developed here should be utilized in future studies to monitor changes across high-latitude treeline ecotones, in order to provide a stronger understanding of the vegetation classes most strongly influenced by a changing climate.

Acknowledgments
This research was undertaken as part of the Assessment of Wood Attributes using Remote sEnsing (AWARE) project (NSERC File: CRDPJ 462973-14, Grantee: N C Coops, FRM, UBC), in collaboration with Canadian Wood Fiber Centre (CWFC), FP-Innovations and Kruger Inc. Additionally, support was provided by the 'National Terrestrial Ecosystem Monitoring System (NTEMS): timely and detailed national cross-sector monitoring for Canada' project jointly funded by the Canadian Space Agency (CSA), Government Related Initiatives Program (GRIP), and the Canadian Forest Service (CFS) of Natural Resources Canada. This research was enabled in part by support provided by WestGrid (www.westgrid.ca) and Compute Canada (www.computecanada.ca). Piotr Tompalski (University of British Columbia) is thanked for providing valuable discussions and feedback on detecting treeline ecotones with ALS data. Chris Hopkinson (with the University of Lethbridge) is thanked for his transect project partnership, which was critical in obtaining the research data used in this study.