Quantifying structural diversity to better estimate change at mountain forest margins

Global environmental changes are driving shifts in forest distribution across the globe with significant implications for biodiversity and ecosystem function. At the upper elevational limit of forest distribution, patterns of forest advance and stasis can be highly spatially variable. Reliable estimations of forest distribution shifts require assessments of forest change to account for variation in treeline advance across entire mountain ranges. Multispectral satellite remote sensing is well suited to this purpose and is particularly valuable in regions where the scope of field campaigns is restricted. However, there is little understanding of how much information about forest structure at the mountain treeline can be derived from multispectral remote sensing data. Here we combine field data from a structurally diverse treeline ecotone in the Central Mountain Range, Taiwan, with data from four multispectral satellite sensors (GeoEye, SPOT-7, Sentinel-2 and Landsat-8) to identify spectral features that best explain variation in vegetation structure at the mountain treeline and the effect of sensor spatial resolution on the characterisation of structural variation. The green, red and short-wave infrared spectral bands and vegetation indices based on green and short-wave infrared bands offer the best characterisation of forest structure with R2 values reported up to 0.723. There is very little quantitative difference in the ability of the sensors tested here to discriminate between discrete descriptors of vegetation structure (difference of R2MF within 0.09). While Landsat-8 is less well suited to defining above-ground woody biomass (R2 0.12–0.29 lower than the alternative sensors), there is little difference between the relationships defined for GeoEye, SPOT-7 and Sentinel-2 data (difference in R2 < 0.03). Discrete classifications are best suited to the identification of forest structures indicative of treeline advance or stasis, using a simplified class designation to separate areas of old growth forest, forest advance and grassland habitats. Consequently, our results present a major opportunity to improve quantification of forest range shifts across mountain systems and to estimate the impacts of forest advance on biodiversity and ecosystem function.


Introduction
Rapid changes in global climate and land-use are driving shifts in forest distribution (Améztegui et al., 2016(Améztegui et al., , 2010Harsch et al., 2009). Mountain ecosystems are expected to experience higher than average temperature increases (Dirnböck et al., 2011;IPCC, 2013;Pepin et al., 2015) and are often subject to land abandonment as agricultural practices change (Haddaway et al., 2014;MacDonald et al., 2000). Naturally occurring elevational treelines (limits of forest distribution) are predominantly climatically determined (Körner and Paulsen, 2004), while the exact position of the treeline and response to environmental change varies due to topographic or geological controls (Butler et al., 2007;Malanson et al., 2011) as well as anthropogenic land-use (Améztegui et al., 2016(Améztegui et al., , 2010. Consequently, shifts in mountain forest distribution have been used as indicators of the impacts of global environmental change (Martin and Bellingham, 2016). Increased forest area and tree growth rates in mountain areas are expected to alter ecosystem service provision, most notably increasing the carbon storage potential of montane forests (Devi et al., 2008;Peng et al., 2009;White et al., 2000). However, forest expansion is considered one of the most significant threats to grassland biodiversity world-wide (Bond and Parr, 2010) and is of concern in mountain ecosystems where disproportionately high numbers of endemic and rare species are found (Steinbauer et al., 2016).
There is limited understanding of how shifts in forest distribution will impact biodiversity and ecosystem function across entire mountain https://doi.org/10.1016/j.rse.2019.01.027 Received 19 September 2018; Received in revised form 17 January 2019; Accepted 21 January 2019 ranges due to variation in the response of mountain forests to environmental change both within-species and between geographic areas . A meta-analysis of forest responses at the upper elevational and latitudinal treelines, found that of 166 sites investigated, 52.4% show upward or poleward migration of forest, 46.4% show no change and 1.2% show downslope movement (Harsch et al., 2009). In areas where change in treeline elevation is not exhibited, increased tree density below the upper limit of forest distribution and across-slope movement have been observed (e.g. Bharti et al., 2012;Klasner and Fagre, 2002). Accurately identifying geographic variation in mountain forest response to environmental change is, therefore, essential to improve the understanding of drivers of forest change and to enable assessments of the impacts of changing treeline position and structure on biodiversity and ecosystem function.
Remote sensing provides an opportunity to expand the scope of field surveys which are often restricted to localised, easily accessible mountain areas due to high time and financial costs. However, characterising variation in forest structure from remote sensing data presents significant challenges for accurately quantifying variation in forest response to environmental change. In some areas, a sharp, welldefined boundary between the forest and the grassland exists. However, mountain treelines are often represented as an ecotone, with a gradual transition between forest and grassland habitats. Despite the prevalence of gradual forest changes globally, there is no optimal method for characterising structural variation in mountain forest-grassland transitions that lack clear boundaries between vegetation classes (Fortin et al., 2000;Hill et al., 2007).
The type of sensor and platform used to acquire remotely sensed data will impact the degree of forest structural information that can be identified and the geographic extent of investigations at mountain treelines. Airborne Laser Scanning (ALS) data are an attractive remote sensing data source for detecting vegetation boundaries across treeline ecotones because of the ability to determine 3-dimensional vegetation structure (Bolton et al., 2018;Coops et al., 2013;Ørka et al., 2012). ALS data have been used to describe vegetation structure within the treeline ecotone (Coops et al., 2013), have been integrated with multispectral satellite imagery to produce maps of vegetation cover types over large areas (Ørka et al., 2012) and have helped improve the interpretation of spectral trends identified from the Landsat data archive (Bolton et al., 2018). Despite the benefit of capturing 3-dimensional information on vegetation structure and the possibility of integrating ALS data with other remote sensing datasets, ALS data are not widely available in many mountainous areas and acquisition of new data sets can be prohibitively expensive. Consequently, there are relatively few published studies using ALS data in mountain treeline ecotones (e.g. Bolton et al., 2018;Coops et al., 2013;Naesset and Nelson, 2007;Ørka et al., 2012).
Synthetic Aperture Radar (SAR) data are sensitive to vegetation structure and, combined with data available from satellite-borne platforms, are attractive for identifying variation in vegetation structure at the treeline ecotone across large areas. Despite the rapid expansion of SAR data availability and the reducing cost of data acquisition, with data from sensors such as Sentinel-1 freely available and high-resolution SAR data available from commercial providers, the use of SAR data in mountain ranges has been restricted due to challenges associated with image processing in mountainous regions. The use of a directional signal in areas with complex and steep terrain often results in geometric distortion of the land surface and occultation due to layover and radar shadowing (Sinha et al., 2015). The capability of Synthetic Aperture Radar (SAR) to penetrate cloud presents obvious benefits for characterising forest structure at mountain treelines. However, there remain significant difficulties in obtaining and processing SAR data with suitable geometric and radiometric properties (Shimada and Ohtaki, 2010) that could be used for large area assessments of forest distribution shifts in mountain ranges.
The most common source of remote sensing data used in the assessment of mountain treeline change to date has been aerial photography or multispectral satellite remote sensing data (Morley et al., 2018). Many studies examine change in forest distribution by classifying multiple remotely sensed multispectral images into forest/nonforest classes, identifying changes in maximum elevation and forest extent over time (e.g. Dinca et al., 2017;Luo and Dai, 2013;Mihai et al., 2017). The simple forest/non-forest definition is an efficient descriptor of the forest-grassland transition and can provide an accurate indicator of change in forest extent if assessed in images from multiple dates. The definition does not, however, capture sufficient information about forest structure to improve the characterisation or understanding of variation in forest response to environmental change (Table 1).
Defining intermediate classes between areas of old-growth forest and treeless habitats improves the representation of structural variation contained within the treeline ecotone. Grouping forest margins into areas that share similar structural characteristics, such as tree canopy cover, density, size and growth form, allows classes to be identified that have reasonably homogeneous within-class forest structure while emphasising between-class variation. Underlying biotic and abiotic processes determine forest structural classes at the treeline Harsch and Bader, 2011) and the impact of forest distribution change on biodiversity and ecosystem function will depend on the forest structure (Greenwood et al., 2016;Tomback et al., 2016). Consequently, using structural classes to represent heterogeneity in the treeline ecotone allows us to characterise variation in changes in forest extent and structure in a manner that improves our understanding of shifts in forest-grassland transitions and their implications (Table 1). However, this level of structural detail is uncommon in studies utilising remote sensing to examine mountain treelines (e.g. Allen and Walsh, 1996;Klasner and Fagre, 2002;Resler et al., 2004). This deficiency exists despite structural classes being sensible ecological units and being efficient to survey. There is also the possibility to identify classes by image classification or by manual interpretation of aerial photography or satellite images with a spatial resolution of 2 m or better (Table 1; Allen and Walsh, 1996).
Defining structural classes requires boundaries to be imposed onto the mountain treeline ecotone. The decision of where to define boundaries along a continuum of differing tree density, size and spatial arrangement that vary over time is non-trivial (Arnot et al., 2004). Table 1 Relative merits of different definitions of vegetation structure at the treeline ecotone for characterising variation in treeline response to environmental change and potential ecological interpretations.
Continuous variables can be used to map the mountain forest transition and offer an attractive alternative to structural classes due to the ability to represent vegetation heterogeneity on a continuous scale, avoiding the use of subjective boundaries (DeFries et al., 2000;Hill et al., 2007). Indeed, classifications of continuous forest descriptors have been used to improve the structural representation of the treeline ecotone (Hill et al., 2007;Král, 2009) and to identify changes in vegetation abundance over time . However, canopy cover, the most commonly used descriptor, is not always appropriate for monitoring change in mountain treelines because of an inability to distinguish differences in tree size class (Morley et al., 2018).
While multi-spectral sensors are the most commonly used source of remote sensing data used in studies of the mountain treeline; there are uncertainties in the ability of multi-spectral sensors to resolve structural variation in mountain treelines. This uncertainty had led to a poor understanding of which spectral properties best characterise structural variation within the treeline ecotone (Morley et al., 2018). Vegetation indices are used to transform two or more spectral bands into indices that emphasise key biophysical characteristics of vegetated ecosystems. The Normalised Difference Vegetation Index (NDVI) correlates with Leaf Area Index (Wang et al., 2005) and fractional vegetation cover (Carlson and Ripley, 1997) and has been used in estimates of tree canopy cover at the mountain treeline (Hill et al., 2007). Green-based indices, such as the Green Normalised Difference Vegetation Index (GNDVI) or the Green-Red Vegetation Index (GRVI), show an increased sensitivity to chlorophyll-a concentration, and consequently, have been suggested to improve the characterisation of subtle differences among ecosystem types (Gitelson et al., 1996;Motohka et al., 2010). The characterisation of forest structure has also been improved by using vegetation indices based on short-wave infrared due to an increased sensitivity to foliar moisture and vegetation density (Schroeder et al., 2011). Indices that make use of shortwave infrared bands have been used to monitor vegetation regrowth following disturbance events, with particular emphasis placed on monitoring post-fire recovery. While indices such as the Normalised Burn Ratio Index (NBRI) were first conceived for monitoring vegetation regrowth post-fire, the sensitivity to foliar moisture and vegetation density makes them potential candidates for characterising variation in vegetation structure in areas of ecological succession, such as across the treeline ecotone.
In addition to remotely sensed vegetation indices, textural features that describe the statistical distribution of pixel data within a defined neighbourhood have been shown to correlate with forest structural variables, such as tree density or average stem diameter (e.g. Meng et al., 2016;Ozdemir and Karnieli, 2011). Sensors with a finer spatial resolution allow for textural features to be calculated at the scale of the individual plots. However, consideration is required to determine if the increased number of textural parameters that can be defined from imagery of higher spatial resolution results in data that will be ecologically meaningful if used in image classification algorithms given the high degree of collinearity present in spectral remote sensing data. Identifying spectral features that show the strongest relationship with forest structure is important to maximise the amount of structural information that can be resolved in multispectral remote sensing data.
Delineation of structural variation at forest margins is required to improve our understanding of the underlying processes that govern variation in forest response to environmental change and to estimate the impacts of forest distribution shifts on biodiversity and ecosystem services. Here we focus on the use of multispectral satellite remote sensing data because it is the most accessible form of remotely sensed data for assessing mountain treeline change across large areas. However, the issue of how best to characterise variation in vegetation structure at the mountain treeline using multispectral satellite remote sensing data remains unresolved. To address this knowledge gap, this work aims to improve the characterisation of mountain treeline ecotones by i) determining which spectral features derived from multispectral satellite remote sensing best explain variation in vegetation structure at the mountain treeline and ii) quantifying the ability of sensors with different spatial resolutions to resolve variation in vegetation structure at the mountain treeline.

Study location
The Central Mountain Range of Taiwan has > 200 mountains over 3000 m a.s.l., the highest of which, Yushan (Jade Mountain), reaches 3952 m. Although Taiwan spans the Tropic of Cancer, the highest elevations experience temperate and alpine conditions. At the highest elevation of forest distribution, the canopy is dominated by four conifer species, primarily Abies kawakamii and Tsuga chinensis with areas of Pinus taiwanensis and Pinus armandii establishment. The adjacent grassland is dominated by the bamboo Yushania niitakayamensis which extends to the peaks with a low density of shrubby species, of which Juniperus spp. and Rhododendron spp. are the most common.
Climate is considered to be the primary regulatory factor of the treeline in the Central Mountain Range, with temperature and topographic sheltering identified as two fundamental controls on treeline structure, position and advance . Natural disturbances caused by small-scale fires and landslides that result in a localised reduction of the treeline and removal of substrate affect the treeline sporadically. However, routine disturbance events are considered to be of low impact at the landscape scale, with little evidence to support widespread anthropogenic disturbance or grazing by large herds of herbivores (domesticated or wild).

Field data
To identify limitations to an accurate characterisation of structural variation this work considers three definitions of vegetation structure at the mountain treeline that have been used across the ecological, biogeographical and remote sensing literature. The forest/non-forest definition is based on the FAO Global Forest Resources Assessment (2018) criteria of a forest with at least 10% canopy cover and trees higher than 5 m or able to reach these thresholds in situ. The FAO (2018) definition was chosen because the leading edge of forest expansion is often characterised by a few trees < 5 m in height. Consequently, the FAO's forest definition aligns with ecological and biogeographic studies investigating pattern-process responses of the treeline ecotone because it captures a greater area of the forest-grassland transition than forest definitions with a higher canopy cover threshold. Six structural classes were identified based on criteria proposed by Harsch and Bader (2011) and subsequently adapted by  for the A. kawakamii treeline in Taiwan. Areas of forest advance were first identified using repeat aerial photography and subsequently defining the structural characteristics of forest plots in field surveys ; Table 2). Structural classes are based on differences in stand density, average tree size and successional stage that the dominant canopy forming species belongs to and include classes that exhibit sharp boundaries as well as diffuse areas (Fig. 1). Complete species separation was not possible due to insufficient field data for species that are sparsely distributed at high elevation in the Central Mountain Range. Consequently, two forest successional stages are defined; the late successional stage is defined as a canopy dominated by A. kawakamii or T. chinensis and the early successional stage dominated by P. taiwanensis or P. armandii. Above-ground woody biomass is investigated as the continuous variable in this analysis because of the correlation with tree size and density from which the categorical groupings are defined ( Fig. 1) and its importance for the estimation of global carbon storage as an Essential Climate Variable (Bojinski et al., 2014).
Field data were collected in the Mt. Hehuan area of the Central Mountain Range. A purposive sampling strategy was used to ensure representation of all forest sub-classes present at the Hehuan treeline. A Table 2 Description of full structural classes based on successional stage and stand structure identified in the Mt. Hehuan region of the Central Mountain Range, Taiwan, and the number of sampling plots.

Vegetation class (Number of plots) Description
Late -Old growth forest (33) Canopy dominated by A. kawakamii or T. chinensis. Areas of the forest interior where the forest has persisted for many years and characterised by a few large trees. Late -Static treeline (12) Canopy dominated by A. kawakamii or T. chinensis. Forested areas at the forest-grassland boundary with trees representative of old growth and no signs of forest advance. Usually with a sharp boundary with the adjacent grassland. Late -Abrupt advancing treeline (24) Canopy dominated by A. kawakamii or T. chinensis. Areas of forest advance that have a high density of establishing trees usually over short distances and a sharp, well-defined boundary with the adjacent old growth forest and grassland. Late -Diffuse advancing treeline Canopy dominated by A. kawakamii or T. chinensis. Areas of forest expansion with a low density of establishing trees usually over long distances and a diffuse, poorly defined boundary with the adjacent grassland. Early -Diffuse advancing treeline (22) Canopy dominated by P. taiwanensis or P. armandii. Areas of forest expansion with a low density of establishing trees usually over long distances and a diffuse, poorly defined boundary with the adjacent grassland. Grassland Areas devoid of tree species but may include a low density of shrubs. Tree density (ha −1 ) Fig. 1. Definitions of structural classes are based on differences in tree size and tree density (median values for stand height and tree density are shown with standard error). The combination of tree size and density results in a significant difference in above-ground woody biomass between structural classes (ANOVA: F (5,148) = 55.96, p < 0.001), with the Early-Diffuse advancing and Late-Diffuse advancing treeline classes not separable due to a similar forest structure but defined by a different successional stage (median values for above-ground woody biomass are shown with standard error).

Grassland
total of 154 plots were sampled, split among the different vegetation classes (Table 2). Early successional species are only found in localised areas of low-density establishment, therefore, are only represented in a single structural class. Data from  were combined with data from an additional survey conducted by the authors in 2016. To retain consistency in plot size, transect data from  were split into 84 subplots measuring 20 × 20 m returning a sample area of 0.04 ha. The 70 plots surveyed in 2016 used a 10 m fixed radius design returning a sample area of 0.03 ha. Since quantity measures are area-standardised, this difference in plot size has no consequence for subsequent analyses. Field plot location were recorded using a handheld Garmin GPSMAP 62s (best accuracy ± 3 m). All trees were measured for Diameter at Breast Height (DBH) at 1.3 m and the height of all live saplings < 1.3 m in height was recorded in all plots. During the 2016 survey, a sample of live trees within each plot was also measured for height. Height was related to DBH using nonlinear least squares regression, thereby allowing estimation of height for any plots where it was not recorded (data not shown). Stand above-ground woody biomass was calculated from stand basal area and median stand height, accounting for differences in specific wood gravity between species. Sapling data were used to inform the designation of structural classes but were not used to calculate stand above-ground biomass values.

Earth observation data
To investigate the importance of sensor spatial resolution on the characterisation of treeline structural variation, data from four multispectral satellite-borne sensors are compared; 2 m pixel size GeoEye multispectral data captured in October 2012, 6 m pixel size SPOT-7 multispectral data captured in October 2016, 10 m and 20 m pixel size Sentinel-2 MSI data captured in October 2016 and 30 m pixel size Landsat-8 OLI data captured in January 2017. Sentinel-2 and Landsat-8 are delivered as orthorectified products and GeoEye and SPOT-7 images were orthorectified using a 30 m resolution SRTM DEM. The spectral bands were calibrated and converted to top-of-atmosphere reflectance in ENVI 5.3 using gain and offset values, accounting for solar irradiance, sun elevation and time of image acquisition. Atmospheric correction was not implemented as single date images are considered independently and pseudo-invariant features (roads and buildings) did not indicate differences in atmospheric conditions between individual images (Song et al., 2001; data not shown). All images were collected in the same season (Autumn-Winter) to avoid differences in vegetation phenology.
Where available, up to seven spectral bands (blue, green, red, near infrared (NIR), red-edge and two short-wave infrared (SWIR) bands) and four vegetation indices were considered (Table 3). The vegetation indices considered were the Normalised Difference Vegetation Index (NDVI), calculated as: where NIR and RED are the near infrared and red spectral bands respectively; the Green-Red Vegetation Index (GRVI), calculated as: where GREEN and RED are the green and red spectral bands respectively; the Green Normalised Difference Vegetation Index (GNDVI), calculated as: where NIR and GREEN are the near infrared and green spectral bands respectively; and Normalised Burn Ratio Index (NBRI), calculated as: where NIR and SWIRII are the near infrared and second short-wave infrared (approx. 2200 nm) spectral bands respectively. Five measures Table 3 Remote sensing spectral features used to investigate the relationship between the statistical distribution of spectral values and vegetation structure at the mountain treeline. The textural features of the centre (mean), dispersion (standard deviation, coefficient of variation) and shape (skewness, kurtosis) were calculated from the pixel values recorded within the boundary of field plots.  of statistical distribution were considered as textural features, the mean, two measures of dispersion (standard deviation and coefficient of variation) and two measures of shape (skewness, kurtosis) were calculated for each spectral band and vegetation index in all sample plots. Due to differences in spatial resolution of the different remote sensing images, we are unable to calculate dispersion and shape statistics for all sample plots. Consequently, the number of sample points that are used in the analysis of the dispersion and shape statistics varies; 154 plots for GeoEye and SPOT-7 and 101 for 10 m pixel size Sentinel-2 (Table 3). It was not possible to calculate dispersion and shape statistics at the plot scale from 20 m pixel size Sentinel-2 data or 30 m pixel size Landsat-8 data, consequently, only the mean spectral response is considered and no other descriptors of spectral response are investigated. Statistical descriptors of spectral response were calculated in R (R Core Team, 2017) using packages raster (Hijmans, 2016) and rgdal (Bivand et al., 2016).

Statistical analysis
Each spectral feature was regressed independently to assess the strength of the relationship between each band or vegetation index and the forest biophysical properties due to the high degree of correlation between spectral bands. The relationship between the forest/non-forest definition and spectral features was assessed using binomial logistic regression with a logit link function. Multinomial logistic regression was used to investigate the probability of separating different structural classes (Table 1) from the spectral data. The results from multinomial regression of the full structural classes indicated that some forest classes could not be separated and consequently, class simplification was carried out. The simplified class structure considered three vegetation classes in multinomial logistic regression: Old-growth forest (an amalgamation of the old-growth forest and static treeline classes), areas of forest advance (an amalgamation of the three classes of forest advance: Early-Diffuse, Late-Diffuse and Late-Abrupt advancing treeline classes) and the grassland class. Least squares regression was used to explore the relationship between above-ground woody biomass and spectral features. The inclusion of zero values of above-ground woody biomass from grassland plots caused heteroscedasticity in model residuals. Consequently, the analysis was conducted as a two-stage procedure, first considering the forest/non-forest definition through binomial logistic regression and subsequently conducting least squares regression on data points from the forest class with a log transformation on aboveground woody biomass.
Multiple regression was carried out for each of the four definitions of vegetation structure (Forest/non-forest, full and simplified structural classes and above-ground woody biomass) to ascertain if the characterisation of structural variation at the treeline could be significantly improved by using multiple spectral predictors. Multi-collinearity was tested for using variance inflation factors and, where present, spectral variables were removed to reduce the severity of multi-collinearity. Model simplification was carried out using partial F-tests to identify the minimum adequate model required to explain variation in the response. To identify potential strengths or limitations of the different forest definitions, the probability of class assignment or above-ground woody biomass was estimated for a subset of the Mt. Hehuan study area using the GRVI derived from Sentinel-2 imagery. All statistical analyses were carried out in R (R Core Team, 2017) using packages boot (Canty and Ripley, 2016) and nnet (Venables and Ripley, 2002), variance inflation was tested for using the car package (Fox and Weisberg, 2011). Statistical significance was considered at p < 0.05 and the coefficient of determination (R 2 ) used to gauge the strength of the relationship (the reported R 2 of the binomial and multinomial logistic regression is McFadden's pseudo R 2 -henceforth R 2 MF ). The reported R 2 for aboveground biomass are from least-squares regression of biomass against spectral properties in which zero values of biomass in grassland areas were excluded due to heteroscedasticity in the residuals.

Textural features
Of the textural features considered, the mean spectral response has the strongest relationship with the four treeline ecotone definitions investigated, particularly with the green and red spectral bands and the GRVI ( Table 4). The dispersion measures (standard deviation and coefficient of variation) of the near-infrared band from the GeoEye sensor show a significant relationship with each of the four definitions of vegetation structure (Table 4). However, the strength of the relationship between dispersion measures and forest-grassland definitions is considerably lower when data from either SPOT-7 or Sentinel-2 are considered (e.g. the strongest measure for forest/non-forest: GeoEye NIR Coef. of Var. R 2 MF 0.464, p < 0.01; SPOT NIR Coef. of Var. R 2 MF 0.200, p < 0.01; Sentinel-2 NIR Coef. of Var. R 2 MF 0.086, p < 0.01; Supplementary 1). While the strength of the relationship between definitions of vegetation structure and the dispersion and shape features derived from the spectral bands are typically better than dispersion measures derived from vegetation indices, the dispersion measures of the GNDVI measured from GeoEye data have a comparable, significant relationship with above-ground woody biomass (St. dev. R 2 = 0.452, p < 0.01; Coef. of Var. R 2 0.469, p < 0.01; Table 4).

Spectral bands
When considering the mean response of spectral bands from all the sensors investigated, the blue, green, red and shortwave infrared bands show the strongest relationship with the definitions of vegetation structure at the treeline ecotone (Highest R 2 (MF) for forest/non-forest: Sentinel-2 Red 0.642, p < 0.01; Full structural Classes: SPOT-7 & Sentinel-2 Blue 0.396, p < 0.01; Simplified structural Classes: Sentinel-2 Blue 0.531, p < 0.01; Above-ground woody biomass: GeoEye Green 0.704, p < 0.01; Table 5), while the near-infrared and red edge bands show a weak relationship (Table 5). The spatial resolution of the sensor has little effect on the strength of the relationship with either the forest/non-forest definition or full structural classes (difference of R 2 MF within 0.06 and 0.09 respectively). When the simplified structural classes are considered the difference in R 2 MF in the visible wavelengths is < 0.08, however, in the near-infrared this difference increases to 0.19 due to a reduced strength of the relationship between the mean response of the near-infrared band from Landsat-8 and the simplified structural classes. Similarly, the strength of the relationship between above-ground woody biomass and the visible and near-infrared bands from the Landsat-8 sensor are consistently 0.12-0.29 R 2 lower than the strongest relationship with the alternative sensors considered here. However, the difference in R 2 among the remaining three sensors in the visible range is < 0.03 and in the nearinfrared there is a difference in R 2 of 0.1 when above-ground woody biomass is used to describe the treeline ecotone (Table 5).

Vegetation indices
When vegetation indices are considered, the mean response of the GRVI and NBRI show the strongest relationships with each of the four forest-grassland transition definitions considered (Highest R 2 for forest/ non-forest: Landsat-8 GRVI 0.623, p < 0.01; Full structural Classes: GeoEye GRVI 0.329, p < 0.01; Simplified structural Classes: SPOT GRVI 0.425, p < 0.01; Above-ground woody biomass: GeoEye GRVI 0.586, p < 0.01; Table 6). The strength of the relationship between the mean response of the Green-Red vegetation index does not depend on the spatial resolution of the sensors compared in this study when categorical definitions are used to describe vegetation structure across the mountain treeline (difference of R 2 MF within 0.07 for forest/non-forest; 0.04 for full structural classes; 0.05 for simplified structural classes, Table 6). However, when above-ground woody biomass is used the P.J. Morley et al. Remote Sensing of Environment 223 (2019)     difference in R 2 between sensors tested here increases to 0.12 (Table 6).

Multiple regression
When considering the GeoEye sensor, the use of multiple spectral bands as predictor variables leads to a significant increase in variance explained across all four of the definitions of vegetation structure considered here (Table 7). However, when considering the alternative sensors tested in this study the benefit of including multiple spectral bands as predictors depends on the definition used to describe vegetation structure across the treeline ecotone. For example, for SPOT-7 and Sentinel-2 the strength of the relationship between spectral bands and both the full and simplified structural class definitions increases when using multiple regression (Table 7). However, when a simple forest/ non-forest or above-ground woody biomass definition is used to describe the vegetation structure at the treeline ecotone, linear models with R 2 above 0.6 can be derived from a single spectral band such as the green or red spectral bands (Table 7). When spectral bands from Landsat-8 are used, the minimum adequate model uses only a single spectral band due to multi-collinearity in the spectral data. However, when vegetation indices are derived from Landsat-8 spectral bands, the combination of GRVI and the NBRI in regression models significantly improves the strength of the relationship with the full structural classes, simplified structural classes and above-ground woody biomass (Table 7).

Data visualisation
The probability of class assignment and above-ground woody biomass were estimated to identify potential strengths or limitations of the different definitions of forest structure at the mountain treeline. Tables 5 and 6 indicate that several spectral bands or vegetation indices would be suitable for this purpose, with the green and red spectral bands and GRVI among the highest performing spectral variables considered. While the NBRI shows a similar strength of relationship to the GRVI, it is only possible to calculate the NBRI from two of the four sensors considered (Sentinel-2 and Landsat-8, Table 6). Therefore, the GRVI is used to compare definitions of forest structure at the mountain treeline because the GRVI is the best performing vegetation index that can be derived from all four of the sensors investigated in this study.
The estimated probability of the forest/non-forest binomial logistic regression shown using the GRVI derived from Sentinel-2 imagery highlights significant areas of the grassland that would be estimated as forest under a definition based on 10% canopy cover if a maximum probability classifier were implemented (Fig. 2). The increased probability of forest occurrence in areas that can be visually identified as grassland coincides with increased standard error compared to neighbouring forested areas (Fig. 2). Based on a maximum probability classification, the mean response from the GRVI would be expected to identify up to three of the six classes considered here: grassland, latediffuse advancing treeline and old growth forest (Fig. 3). This approach does not determine differences in community composition between the early-and late-successional stages in the diffuse advancing structure, areas of late-abrupt advance or differences between the static treeline and old growth forest (Fig. 3). Simplifying intermediate classes, by amalgamating the late-static class with the old growth forest as well as combining the three classes indicative of forest advance, results in a better distinction between old growth forest and grassland areas as well as a reduction in the area of grassland that would be incorrectly classified as forested under a forest/non-forest approach (Fig. 4). While the simplified class structure leads to better discrimination of areas of forest advance from areas of old-growth, this approach is unable to resolve heterogeneity in forest structure at the mountain treeline due to overlap in the spectral properties of forest structural classes (Fig. 3). The estimation of above-ground biomass indicates an improved ability to estimate differences in structural heterogeneity at the mountain treeline showing good correspondence to the true colour image. However, the characterisation of areas with biomass values above 25 t C ha −1 are likely to be inaccurate due to a saturation effect that occurs in the relationship between biomass and the spectral properties (Fig. 5). In addition, the spectral signature of some grassland areas still results in an elevated above-ground woody biomass estimation in areas not undergoing forest expansion. However, estimated biomass values of these grassland areas are reduced when compared against the neighbouring forested areas (Fig. 5).

Discussion
Here we show that the ability to identify variation in forest structure at the mountain treeline using multispectral satellite remote sensing data is best achieved when above-ground woody biomass is used to describe variation in vegetation structure (Tables 4-7). Furthermore, we show that a simplified class structure that considers areas of forest advance separately to old growth forest improves the discrimination of areas indicative of forest advance or stasis. The relationships defined here between four definitions of vegetation structure at the mountain treeline and spectral features highlight little quantitative difference between the remote sensing sensors tested here (difference of R 2 (MF) within 0.03 for forest/non-forest; 0.11 for full structural classes; 0.05 for simplified structural classes; 0.16 for above-ground woody biomass; Table 7). Consequently, effective use of multispectral satellite remote sensing data presents a major opportunity to improve the ecological understanding of range shifts in mountain forests and estimate their subsequent impacts to biodiversity and ecosystem function.
The relationship between treeline definition and spectral variables is strongest when using above-ground woody biomass or the forest/ non-forest definition to describe the forest-grassland transition (strongest R 2 for above-ground woody biomass: GeoEye Green mean & NIR Coef. of Var. 0.723, p < 0.01; strongest R 2 MF for forest/non-forest: SPOT-7 Green mean 0.645, p < 0.01; Table 7). However, both definitions show limitations in the ability to characterise areas indicative of forest advance or stasis across forest-grassland transitions in mountain ecosystems. When using a forest/non-forest definition, thresholds of canopy cover used to delineate a forest boundary in an ecotone are difficult to define because areas of grassland and areas with a low forest canopy cover can have similar spectral responses which in turn influences estimates of forest extent ( Fig. 2; Arnot et al., 2004;Hill et al., 2007). Song et al. (2014) found that varying the threshold of canopy cover from 20 to 30% resulted in considerable disagreement in forest cover estimates, resulting in a significant under-representation of diffuse forest expansion when a threshold of 30% canopy cover was used. Treeline ecotones are often characterised by areas with sparse and discontinuous tree cover and, consequently, assessments of change must be able to identify areas of diffuse forest expansion accurately. The sensitivity of the canopy cover threshold used to define the forest/ non-forest boundary highlighted above not only leads to high uncertainty in estimates of forest expansion but also understates the variety of responses of the mountain treeline, restricting the ecological interpretation of forest change in mountain treeline ecotones Broll, 2017, 2007).
The representation of the mountain treeline ecotone is improved when using above-ground woody biomass to characterise vegetation structure. However, using above-ground woody biomass to represent the treeline ecotone is likely to under-estimate biomass in areas of old growth forest because spectral reflectance has an asymptotic relationship with plant biomass, which leads to a saturation effect in dense vegetation (Fig. 5;Asner et al., 2003;Huete et al., 1997). In upland plantations, Puhr and Donoghue (2000) highlighted that the spectral signature of conifer trees converges with increasing size and as the canopy approaches closure. Consequently, once coniferous forest stands approach 13 m in height and the basal area exceeds 40 m 2 ha −1 identifying differences in forest structure is problematic, and predictions are likely to become unreliable (Puhr and Donoghue, 2000). The saturation of GRVI indicated in Fig. 5 coincides with the class average biomass values of the abrupt advancing treeline where trees establish in high density and reach canopy closure quickly (Fig. 1). Consequently, it may be possible to identify changes in above-ground woody biomass in areas undergoing forest expansion and so improve the characterisation of vegetation structure. However, characterisation of areas with biomass values above 25 t C ha −1 are likely to be inaccurate and thus limits the use of above-ground biomass as a single predictor of forest structure at mountain treeline ecotones.
Structural classes that define intermediate classes between areas of old-growth forest and treeless habitats have not been widely used in studies using remote sensing data to identify shifts in mountain forests. Consequently, there was uncertainty surrounding the degree of structural information that multispectral satellite remote sensing is able to resolve (Morley et al., 2018). The high spectral similarity between the late successional -old growth forest, late-static treeline and late-abrupt advancing treeline classes indicates a saturation in the spectral properties during the transition between the closed canopy, abrupt advancing treeline and old growth forest structures (Fig. 3). Amalgamating the static treeline and old growth forest classes leads to an improvement in the ability to identify areas of old-growth forest (Fig. 4). Similarly, amalgamating the early-diffuse, late-diffuse and late-abrupt advancing forest classes leads to an improved ability to separate areas at the leading edge of forest advance from areas of old-growth forest (Fig. 4). Consequently, we find that the use of simplified structural classes improves the characterisation of areas indicative of treeline advance or stasis. However, overlap in the spectral properties means that it is not possible to identify variation in forest structure within areas of forest advance using discrete classes.
The spectral similarity highlighted by the probability estimates of the forest/non-forest regression persists between areas of diffuse forest advance and some areas of the grassland (Fig. 3). At the leading edge of forest advance, the diffuse advancing treeline class is characterised by a few trees < 5 m in height. Consequently, the size and density of establishing trees are not sufficient to provide a significant difference in spectral reflectance in some grassland areas when using a single date image, leading to an over-estimation in the extent of diffuse treeline advance (Figs. 3 & 4). The spectral similarity between some grassland areas and the diffuse advancing treeline exists in areas where changes in forest extent and structure are occurring most rapidly. Consequently, the identification of areas of forest advance can be improved by comparing images over time (e.g. Dinca et al., 2017;Mihai et al., 2017;Bharti et al., 2012). The Landsat archive offers the most consistent source of multispectral satellite data with images dating back to the 1980s at 30 m pixel size. Concerns had been raised over the potential suitability of data from the Landsat archive to characterise vegetation heterogeneity due to the spatial resolution of the spectral data (Bharti et al., 2012;Buchanan et al., 2015;Chen et al., 2015). It is often perceived that imagery with a high spatial resolution will improve the characterisation of habitat heterogeneity because of the ability to identify small objects, e.g. individual trees. However, we show that spectral features derived from Landsat-8 data have a comparable strength of relationship to higher resolution imagery when simple structural classes are used, despite the ability to include multiple measures of spectral texture at the plot scale from imagery with a high spatial resolution (difference of R 2 MF between sensors tested within 0.03 for forest/non-forest and 0.05 for the simplified structural classes, Tables 7; Donoghue and Watt, 2006). Consequently, exploiting the long-term, open-access Landsat archives to identify changes in forest extent over time will improve estimates of forest distribution change at the leading edge of forest-grassland transitions (e.g. Dinca et al., 2017;Mihai et al., 2017).
While exploiting the Landsat archive is beneficial for identifying areas of treeline change, there are still difficulties in characterising variation in forest structure within areas of change using structural classes. Identifying change between simplified class assignment over time using archived Landsat data would allow estimation of the extent of forest change and stand age of advancing forest, however, this approach would not directly characterise forest structure. An alternative solution is to estimate above-ground woody biomass directly using imagery of higher resolution within areas identified as advancing forest. While data from the GeoEye sensor returns the highest correlation coefficient with above-ground woody biomass (R 2 = 0.723, p < 0.01), the difference in correlation coefficient using data from Sentinel-2 or SPOT-7 are within 0.04 (R 2 = 0.681, p < 0.01; R 2 = 0.682, p < 0.01 respectively). Consequently, there is an opportunity to make use of freely available Sentinel-2 data to estimate aboveground woody biomass at the stand scale, thereby allowing variation in forest structure to be characterised within areas undergoing forest advance at the mountain treeline. This two-stage approach would allow variation in vegetation structure at the mountain treeline to be estimated with higher confidence than using either a forest/non-forest classification, which over-estimates forest cover, or the direct estimation of above-ground woody biomass, which would give an unreliable estimate of biomass past a saturation threshold.
Structural classes are broad enough to be used at a global scale, and while not necessarily present in every mountain forest, are suitably flexible to be adapted to the local community and structural composition. Here the spectral similarity of conifer species meant it was not possible to distinguish between the early successionaldiffuse advancing treeline and the late successionaldiffuse advancing treeline. . Based on a maximum probability approach, three of the six vegetation classes would be estimated while the early-diffuse advance, late-abrupt advance and Late-static forest structural classes are unlikely to be identified.
However, Bharti et al. (2012) have shown the ability of multispectral Landsat imagery to separate coniferous and broadleaf species at other elevational treelines. Therefore, in other mountain areas, broader species differences could be incorporated to account for species-specific responses to environmental change providing suitable training data are available. When adapting structural classes, the spectral similarity between classes highlighted above emphasises the importance of independent accuracy assessments, the process by which image classification algorithms are trained and validated using subsets of the groundtruthing data set (see Castilla, 2016;Olofsson et al., 2014Olofsson et al., , 2013. While accuracy assessments are required of any study using remotely sensed data to estimate land-surface properties; they remain an element absent from many previous studies of mountain treelines (Morley et al., 2018). At the leading edge of mountain forest distribution, some areas respond to environmental change very quickly while other areas slowly or not at all (e.g. Harsch et al., 2009;Lloyd, 2005). Consequently, as the uptake of remote sensing technology in assessments of changing mountain forest distribution increases, there is a need to ensure that conclusions drawn about changes in forest extent and structure at forest margins are reliable when scaled up to assess entire mountain ranges.

Conclusion
Obtaining estimates of changes in forest distribution over large areas is challenging in mountain areas where steep terrain often restricts the geographical scope of field campaigns. Change assessments must account for variation in forest structure to make reliable estimations of the impacts of distribution shifts on biodiversity and ecosystem function. By comparing different satellite sensors against four definitions of vegetation structure at the mountain treeline that are widely used in the ecology, biogeography and remote sensing literature, we demonstrate that the identification of areas indicative of forest advance or stasis is best achieved using a simplified class structure while variation in structure within areas of forest advance is best characterised in multispectral satellite remote sensing using above-ground woody biomass to describe forest structure. There is very little difference in the ability of the sensors tested here to discriminate between categorical descriptors of vegetation structure, and while Landsat 8 is less well suited to defining above-ground woody biomass there is little difference between the relationships defined for GeoEye, SPOT-7 and Sentinel-2 data. The results presented here enable structural variation in mountain forest margins to be identified in multispectral satellite remote sensing, facilitating research in mountain areas where significant fieldwork is not possible. Consequently, the methods described in this paper will advance our understanding of the ecological mechanisms driving forest distribution shifts across mountain ranges and improve estimates of the impacts that changes in forest distribution will have on biodiversity and ecosystem function.