Comparison of Empirical and Physical Modelling for Estimation of Biochemical and Biophysical Vegetation Properties: Field Scale Analysis across an Arctic Bioclimatic Gradient

: To evaluate the potential of multi-angle hyperspectral sensors for monitoring vegetation variables in Arctic environments, empirical and physical modelling using ﬁeld data was implemented for the retrieval of leaf and canopy chlorophyll content (LCC, CCC) and plant area index (PAI) measured at four sites situated across a bioclimatic gradient in the Western Canadian Arctic. Field reﬂectance data were acquired with an ASD FieldSpec (305–1075 nm) and used to simulate CHRIS Mode1 spectra (411–997 nm). Multi-angle measurements were taken corresponding to CHRIS view zenith angles (VZA) ( − 55 ◦ , − 36 ◦ , 0 ◦ , + 36 ◦ , + 55 ◦ ). Empirical modelling compared parametric regression based on vegetation indices (VIs) to non-parametric Gaussian Processes Regression (GPR). In physical modelling, PROSAIL was inverted using numerical optimization and look-up table (LUT) approaches. Cross-validation of the empirical models ranked GPR as best, followed by simple ratio (SR) with optimally selected NIR and red wavelengths, and then ROSAVI using its published wavelengths (mean r 2cv = 0.62, 0.58, and 0.54, respectively across all sites, variables, and VZAs). However, the best predictive performance was achieved by SR followed by GPR and ROSAVI (NRMSE cv = 0.12, 0.16, 0.16, respectively). PROSAIL simulated the multi-angle top-of-canopy reﬂectance well with numerical optimization ( r 2 = ~0.99, RMSE = 0.004 ± 0.002), but best performing LUT models of LCC, CCC and PAI were poorer than the empirical approaches (mean r 2 = 0.48, mean NRMSE = 0.22). PROSAIL performed best at the high Arctic sparsely vegetated site ( r 2 = 0.57–0.86 for all parameters). Overall, the best performing VZA was − 55 ◦ for empirical modelling and 0 ◦ and ± 55 ◦ for physical modelling; however, these were not signiﬁcantly better than the other VZAs. Overall, this study demonstrates that, for Arctic vegetation, nadir narrowband reﬂectance data used to derive simple empirical VIs with optimally selected bands is a more e ﬃ cient approach for modelling chlorophyll and PAI than more complex empirical and physical approaches. methodology, B.E.K.; validation, B.E.K.; formal analysis, B.E.K.; investigation, B.E.K.; resources, B.E.K., D.J.K., and J.D.; data curation, B.E.K.; writing—original draft preparation, B.E.K.; writing—review and editing, B.E.K. and D.J.K.; visualization, B.E.K.; funding acquisition, D.J.K. and J.D.


Introduction
Estimation and mapping of vegetation biophysical and biochemical characteristics is one of the most common applications of optical remote sensing and is implemented at all scales from local to global. Quantifying biochemical constituents is of great interest because of their importance to foliar and ecosystem processes (e.g., photosynthesis and nutrient cycling), their use for identifying plant 1. compare multi-angle vegetation spectra and measurements of leaf and canopy chlorophyll (LCC, CCC, respectively) and plant area index (PAI) at various locations in the Western Canadian Arctic that represent a latitudinal climate gradient; 2.
compare parametric linear regression combined with common VIs to non-linear non-parametric machine learning (GPR) for estimation of LCC, PAI, and CCC at all view angles and sites; 3.
compare the empirical modelling results to PROSAIL models inverted by LM and LUT methods; 4.
assess the effect of resampling spectral resolution and scaling-up field measurements on model results.

Study Sites
The study sites included Herschel Island (Yukon), Banks Island (Northwest Territories), Richards Island (Northwest Territories), and the Richardson Mountains (Yukon/Northwest Territories) ( Figure 1). They are situated in the Northern and Southern Arctic and Taiga Cordillera Terrestrial Ecozones [72] within the Western Canadian Arctic. Satellite and historical records in the region, such as ground-based photographs and plot-based studies, provide evidence of enhanced greening [7,10,[73][74][75][76][77][78][79][80][81] and rapid warming [82][83][84]. Changes to vegetation across this region may indicate how and how quickly other Arctic terrestrial ecosystems will respond to similar warming conditions [10]. Therefore, this region provides an ideal setting for vegetation related monitoring studies and methodological investigations.  [85] and Northern Geodatabase [86]. Scale bar is approximate, given that scale is variable with conical projections. Photographs of ground conditions are shown for Banks Island (top), Richards Island (middle), and Richardson Mountains (bottom).

Field Sampling Design
All field sites were visited between July 10th and August 4th in the years between 2011 and 2014 during the peak of the growing season. Data were acquired in both 2013 and 2014 for the Richardson Mountains site.
Sampling for field measurements at each site was conducted following the design schematic shown in Figure 2. Between 18 and 33 plots of 90 m × 90 m were sampled at each site depending on accessibility and weather. This plot size was selected to represent approximately 3 × 3 pixels in CHRIS Mode 1 (M1) (34 m pixels at nadir) imagery [91] and to account for vegetation heterogeneity and adjacency effects within satellite imagery. Locations of plots were determined using a stratified sampling approach: spectrally distinct and uniformly vegetated areas were identified through visual analysis of moderate resolution imagery (Landsat 5 and 8) and through on-the-ground surveys. Within these zones, field plot locations were randomly selected away from edges and distinct gradients. The central location of each 30 m × 30 m square in Figure 2 was recorded with a handheld Trimble Juno GPS unit accurate to~3 m. Subplots of 1 m × 1 m were placed within each plot as shown in Figure 2 (upper right). Within each subplot, five locations were selected for multi-angle spectral reflectance measurements ( Figure 2, bottom right, dotted circles) and ten locations were randomly selected for destructive and non-destructive chlorophyll measurements ( Figure 2, bottom right, red x's). A downward hemispherical photo was taken in each subplot for estimation of PAI [92,93] and destructive vegetation samples were taken for lab measurement of PAI from the five red subplots in Figure 2.

Field Measurements of Leaf Chlorophyll Content (LCC)
Measurements of non-senescent LCC (µg/cm 2 ) for the dominant species were made at each study site using a Minolta SPAD-502 chlorophyll meter (Minolta Camera Co. Osaka, Japan). Leaves were sampled randomly throughout the canopy and were sampled regardless of phenological state (both green and senescent leaves were included in the samples). The emphasis, however, was on sampling leaves (and other vegetation components) that would be detected by the optical sensor (i.e., leaves obscured from 90% of the canopy were not the primary sampling target), and so in this case, leaves occurring at the bottom of the canopy were not sampled. Species with extremely small needles were not sampled, as the SPAD field of view was too large for proper measurement. This instrument is commonly used to characterize chlorophyll concentration in plant species [19][20][21][94][95][96][97][98][99]. SPAD values are a unitless measure of LCC derived from transmittance in the red (650 nm) and near-infrared (NIR) (940 nm) using a calibrated light source [100]. Anatomical/morphological differences between specimens can confound readings, so caution should be exercised [91]. A predefined linear calibration function, or consensus equation [95], was used to convert the SPAD values to actual measurements of chlorophyll in µg/cm 2 . Within each of the subplots, 10 SPAD readings (e.g., at red x's in Figure 2) were averaged to produce one LCC measurement per subplot. To determine how well the SPAD measurements represented actual LCC values, lab-based measurements were made using destructive samples collected at the Richardson Mountains site. The lab-based methods for LCC, as well as carotenoids, non-pigments, equivalent water thickness (EWT), and leaf mass per area followed standard lab protocols (e.g., [101][102][103][104]). A full description of the lab-based procedures used to estimate Remote Sens. 2020, 12, 3073 7 of 41 these variables can be found in Kennedy [91]. The averaged pigment and non-pigment values for each species were then weighted by the species mass proportions in the plots to provide averaged values across the plot canopy. Using analysis of variance (ANOVA), no significant differences were found between the converted SPAD values and lab-derived LCC datasets, indicating that the converted SPAD (i.e., instrument-based measurements) LCC values for all sites could be used with confidence.

Field Measurements of Leaf Chlorophyll Content (LCC)
Measurements of non-senescent LCC (µg/cm 2 ) for the dominant species were made at each study site using a Minolta SPAD-502 chlorophyll meter (Minolta Camera Co. Osaka, Japan). Leaves were sampled randomly throughout the canopy and were sampled regardless of phenological state (both green and senescent leaves were included in the samples). The emphasis, however, was on sampling leaves (and other vegetation components) that would be detected by the optical sensor (i.e., leaves obscured from 90% of the canopy were not the primary sampling target), and so in this case, leaves occurring at the bottom of the canopy were not sampled. Species with extremely small needles were not sampled, as the SPAD field of view was too large for proper measurement. This instrument is commonly used to characterize chlorophyll concentration in plant species [19][20][21][94][95][96][97][98][99]. SPAD values are a unitless measure of LCC derived from transmittance in the red (650 nm) and nearinfrared (NIR) (940 nm) using a calibrated light source [100]. Anatomical/morphological differences between specimens can confound readings, so caution should be exercised [91]. A predefined linear

Field Measurements of PAI and Calculation of Canopy Leaf Chlorophyll Content (CCC)
Downward hemispherical photos (DHPs) were taken with a 12.3-megapixel Nikon D300 camera equipped with a Sigma 8 mm F3.5 EX DG Circular Fisheye lens. They were taken from an approximate height of 1.5 m at all subplot locations (white squares in Figure 2) for a total of 21 photos per plot to capture canopy spatial variability. Due to the limited time for this remote fieldwork, photographs were taken under all sky conditions and times of day. Photos were batch processed in Photoshop (Adobe Inc.) to correct white balance and exposure, and to partially eliminate shadowing, which varied between photos depending on sky conditions. Processing was limited to view zenith angles (VZA) equal to and less than 45 • to limit the number of mixed pixels (at greater view angles) and limit the field of view to be just larger than the subplot size (i.e., 1 m 2 ).
CAN-EYE software [105,106] was then used to derive PAI and average leaf inclination angle (ALA) (see Kennedy [91] for more details). PAI was used for this study because it considers all green vegetation components within the canopy (e.g., stems and leave) [58,92,106], whereas LAI, defined as the total surface area of green leaves relative to a defined ground surface area, does not capture/represent all photosynthetic components of tundra vegetation. Before processing, a calibration procedure was completed [106]. For each set of 21 DHPs per plot, gap fraction was calculated through two-class (vegetation, gap) supervised classification. Image classification was done in angular sectors with a zenith (∆θ) resolution equal to 2.5 • and an azimuth (∆ϕ) resolution equal to 2.5 • . Gap fraction (P o, CAN−EYE (θ)) was computed for each zenithal ring by taking the average values for all photographs across the 72 azimuthal sectors, excluding the masked areas. PAI was computed from the gap fraction using the Poisson law (distribution) as given in Welles and Norman [107] and Demarez et al. [108]. PAI values were calculated/integrated across transects and are therefore provided at the plot level.
To evaluate field CAN-EYE PAI values against a reference, the same Richardson Mountains site vegetation samples used for lab-based LCC measurement were laid out on a 1.8 m 2 white sheet in the lab and photographed (n = 99). Vegetation species such as moss may confound estimation of PAI due to its non-leafy structure and cover characteristics. Moss was not included in the destructive samples unless dominant and visibly exposed -only the top green components were harvested if present (i.e., the green components visible to the camera). Classification of vegetation was then conducted, and the relative number of vegetation pixels was taken as 'destructive' PAI [91]. Comparison of these to field-measured PAI values showed strong relationships with both 'effective' PAI (r 2 = 0.86, RMSE = 0.15), which does not account for clumping and 'true' PAI (r 2 = 0.77, RMSE = 0.21), which accounts for clumping. This gave confidence that the field PAI measurements for all sites were well representative and CAN-EYE derived effective PAI was selected for use in subsequent analyses.
Canopy chlorophyll content (CCC, g/m 2 ) was then calculated by multiplying LCC by PAI. CCC has been used in other studies as a means of scaling up leaf-scale measurements and assessing canopy-level chlorophyll content [19,20,59,109,110]. It provides a more useful predictor of canopy chlorophyll than simply using LCC, since the spatial extent of vegetation cover is considered in the scaling. Table 1 shows the descriptive statistics for LCC, PAI, and CCC measured at each field site. All variables showed a north-south gradient, where lowest overall values for each occurred at the Banks Island field site [91].

Field Spectral Reflectance Measurements and CHRIS/PROBA Simulation
Spectral reflectance measurements were acquired in each subplot ( Figure 2) using an ASD FieldSpec handheld spectroradiometer (Analytical Spectral Devices, Inc. Boulder, USA) with a 10 • fore optic. It uses a 512-channel silicon photodiode array, recording radiance between 305 nm and 1075 nm in 1.6 nm nominal bandwidths. Reflectance was derived by comparing radiance values to those for a Spectralon panel (Labsphere Inc., North Sutton, New Hampshire, USA). Measurements were taken at nadir, ±36 • and ±55 • (Figure 3), the VZAs of the CHRIS/PROBA sensor/platform, where positive angles correspond to forward scatter and negative angles to backscatter (see Kennedy et al. [71], which integrates the findings of this paper in vegetation modelling using multi-angle CHRIS M1 data). A tripod with fine angle adjustment was used and angles were measured using a digital inclinometer (accuracy ±0.2 • ) (Beall Tool Company, Newark, Ohio, USA) mounted on the spectroradiometer. Measurement height at nadir was 1.5 m and was reduced for ±36 • and ±55 • in an arc, similar to a goniometer. The resulting ground field of view (GFOV) varied from~26 cm (nadir) to~50 cm (±55 • ) but all GFOVs were less than one-half the subplot size. Furthermore, measurements were taken on the CHRIS/PROBA descending pathway angle of~7.84 • from true north to best simulate its expected VZA geometry. Timing of data acquisition was planned to align with CHRIS acquisitions, but this was not always possible due to weather and variable CHRIS/PROBA acquisition timing. Five measurements taken at the same position in each subplot were averaged to a single subplot value. These 21 subplot values were then averaged to derive a single value for each 90 m × 90 m plot ( Figure 2). In addition, subplot and plot averages were convolved to the 62 spectral bands (400-1050 nm; 10 nm bandwidth) of CHRIS M1 data using the equation presented in [29,111]. A full description of the methodology used to simulate CHRIS data can be found in Kennedy [91]. All measured field spectra were filtered with the Savitzky-Golay [112] technique using a second order polynomial and 15 data points (i.e., wavelengths). This technique is widely used in spectroscopy [113] for noise reduction while preserving the higher order moments of the original spectrum. Overly noisy portions at both ends of the spectral range were discarded, resulting in spectra between 400 nm and 1050 nm for the ASD full spectrum data, and 411 nm to 997 nm for the simulated 62-band CHRIS M1 data.
Although a thorough analysis was not completed for the exact position/geometry of the sun in relation to the sensor, sampling at each site was completed at comparable seasons (mid-July) and daily times (~10:00 h to 17:00 h) to reduce inter-site sun-induced variability. The associated polar plots of the sun/sensor positions for the high-Arctic and low-Arctic sites can be seen in Figure 4 (sampling periods are shown in green; 24-hour sun positions are shown in red). Sun zenith angles (SZA) typically varied between~46 • to 66 • and sun azimuth angles (SAA) varied from~139 • to 262 • for all field sites, where greater SZAs occurred at lower latitudes (i.e., Richardson Mountains). Finally, ground spectra were taken to best simulate the typical forward and backscatter reflectance observed in the CHRIS imagery (given the constraints imposed by timing of fieldwork and the need to produce adequate sample sizes). Backscatter reflectance was measured at relative azimuth angles (RAA) of less than 75 • , whereas forward scatter reflectance was measured at RAAs greater than 90 • . Remote Sens. 2020, 12, x FOR PEER REVIEW 10 of 42 were taken on the CHRIS/PROBA descending pathway angle of ~7.84° from true north to best simulate its expected VZA geometry. Timing of data acquisition was planned to align with CHRIS acquisitions, but this was not always possible due to weather and variable CHRIS/PROBA acquisition timing. Five measurements taken at the same position in each subplot were averaged to a single subplot value. These 21 subplot values were then averaged to derive a single value for each 90 m × 90 m plot ( Figure 2). In addition, subplot and plot averages were convolved to the 62 spectral bands (400-1050 nm; ~10 nm bandwidth) of CHRIS M1 data using the equation presented in [29,111]. A full description of the methodology used to simulate CHRIS data can be found in Kennedy [91]. All measured field spectra were filtered with the Savitzky-Golay [112] technique using a second order polynomial and 15 data points (i.e., wavelengths). This technique is widely used in spectroscopy [113] for noise reduction while preserving the higher order moments of the original spectrum. Overly noisy portions at both ends of the spectral range were discarded, resulting in spectra between 400 nm and 1050 nm for the ASD full spectrum data, and 411 nm to 997 nm for the simulated 62-band CHRIS M1 data.
Although a thorough analysis was not completed for the exact position/geometry of the sun in relation to the sensor, sampling at each site was completed at comparable seasons (mid-July) and daily times (~10:00h to 17:00h) to reduce inter-site sun-induced variability. The associated polar plots of the sun/sensor positions for the high-Arctic and low-Arctic sites can be seen in Figure 4 (sampling periods are shown in green; 24-hour sun positions are shown in red). Sun zenith angles (SZA) typically varied between ~46° to 66° and sun azimuth angles (SAA) varied from ~139° to 262° for all field sites, where greater SZAs occurred at lower latitudes (i.e., Richardson Mountains). Finally, ground spectra were taken to best simulate the typical forward and backscatter reflectance observed in the CHRIS imagery (given the constraints imposed by timing of fieldwork and the need to produce adequate sample sizes). Backscatter reflectance was measured at relative azimuth angles (RAA) of less than 75°, whereas forward scatter reflectance was measured at RAAs greater than 90°.

Parametric Linear Regression: Vegetation Indices
Eleven 2-band VIs comprised of a red and NIR band (e.g., simple ratio (NIR/Red); NDVI (NIR-Red/NIR+Red), etc.) were evaluated by implementing regressions for all possible two-band combinations against LCC, PAI, and CCC. These analyses were conducted for the ASD full spectrum and simulated CHRIS M1 datasets at each VZA, at both the subplot and plot scales, and for all four study sites. In total,~2700 tests were performed that allowed identification of the optimal band combination(s) for each VI and the best VI models for the vegetation variables. In addition, 65 published narrowband VIs with potential for chlorophyll and PAI estimation were regressed against the vegetation variables over the same site and data conditions listed above for a total of 20,000 regressions. Models based on the full ASD spectra used the wavelengths as specified for the given VI in the literature, while for the simulated CHRIS M1 datasets the closest CHRIS bands were used. Kennedy [91] lists the 11 two-band VIs and the 65 narrowband VIs that were tested.

Non-Parametric Gaussian Processes Regression
Gaussian Processes Regression (GPR), a non-parametric machine learning process was implemented to model and predict LCC, CCC and PAI for the same datasets as in the parametric regressions. It was selected based on a review of the literature on such methods and particularly based on Verrelst et al. [24,32], who compared several methods and selected GPR for estimation of LAI using hyperspectral data [24,32]. GPR assumes that an a priori Gaussian process governs the set of all possible unobserved latent functions that are consistent with the observed/known data [24,54]. The likelihood of the latent functions combined with the known data observations allows for posterior probabilistic estimates of the predicted variable [24]. Here, GPR was used to identify the optimally predictive wavelengths and to determine if the reduced 62-band simulated CHRIS M1 data could perform as well as the full ASD spectra in modelling the vegetation variables. To find the optimal wavelengths in the full ASD reflectance dataset, the GPR algorithm was programmed to remove 100 bands in each iteration (GPR-100); trials using single band removal per iteration proved to be too lengthy. For the simulated CHRIS M1 spectral subset, one band at a time (GPR-1) was removed, since there were only 62 bands in total.
The GPR algorithm provided the following outputs: the goodness-of-fit validation statistics (e.g., r 2 ) as a function of number of bands plotted over the sequentially removed bands until one was left (validation statistics are provided for each step); the associated wavelengths for each GPR model step; a list of the most important predictor variables (i.e., spectral bands); and the sum of the hyperparameter values for all bands in the dataset, where lower values provide an indication of higher predictive power of a spectral band. Full descriptions of the GPR implementation and algorithm is provided elsewhere [24,54,91].

Physical Modelling: PROSAIL
The PROSAIL RTM [33], a coupling of the PROSPECT-5 leaf model [34,35] and the SAILH canopy model [36,37,114], was used to simulate top of canopy (TOC) spectral reflectance and retrieve LCC, PAI, and CCC. PROSPECT-5 calculates the directional-hemispherical reflectance and transmittance of leaves based on two classes of input variables: (1) the leaf structure parameter (N) which is the number of compact layers specifying the average number of air/cell wall interfaces within the mesophyll; and (2) the leaf biochemical content (total chlorophylls, C ab ; total carotenoids, C ar ; brown pigment content, C bp ; equivalent water thickness, C w ; dry matter content, C m ) [33]. The PROSPECT-5 model parameters are shown in Table 2 with their respective units. The newest version, PROSPECT-D [39], which includes anthocyanins, was not used as it was released after the completion of most of this study. The spectral leaf optical properties (reflectance and transmittance) calculated by PROSPECT-5 are inputs into the SAILH canopy reflectance model. SAIL simulates the bidirectional reflectance factor (BRF) of turbid-medium plant canopies [33]. Defining a canopy as a turbid medium represents it as a horizontally homogenous and semi-infinite layer comprised of small vegetation elements of a given geometry and density that act as absorbing and scattering particles [20,21,33]. Consequently, the model is best adopted for use in homogeneous vegetation canopies [20,21,36,115,116]. Apart from leaf reflectance and transmittance, SAILH requires eight input parameters, including multi-angular spectral measurements: (1) leaf area index (LAI or PAI as used/defined in this study); (2) mean leaf inclination angle, (ALA); (3) hot spot size parameter (S L ), defined as the ratio between the average size of the leaves and the canopy height [37]; (4) background/soil reflectance (P S ); (5) fraction of diffuse incoming solar radiation (SKYL); (6) sun zenith angle (SZA); (7) sensor view zenith angle (VZA); and (8) relative azimuth angle (RAA), which is defined as the relative angle between the sun azimuth and sensor azimuth angles.
Inversion of PROSAIL was accomplished with two common methods: (1) the Levenberg-Marquardt iterative numerical optimization algorithm (LM) [117], and (2) a look-up table (LUT) approach. Iterative optimization was chosen because it is widely used for model inversion [30,33,34,39,61,118] etc., does not require a training/calibration dataset [40,41,118], and has shown similar performance as LUT-based inversion approaches in comparison studies [119,120]. LM was selected as the iterative optimization technique because it has been shown to converge quickly [121] and had recently produced excellent results with PROSAIL inversion for the retrieval of chlorophyll from multi-angle spectroscopic data [61]. A LUT approach was selected as a comparative approach because it has been shown to be a robust, physically sound method for vegetation variables retrievals, especially when regularization strategies are implemented to account for ill-posedness [20,59]. With the implementation of appropriate cost functions, the addition of noise to account for sensor and environmental uncertainty, and the use of the mean of multiple best solutions to account for sensor uncertainty, LUT-based inversions have been shown to produce excellent results [32,59].
LM was implemented through the MPFIT package for ENVI [122]. Table 2 shows the upper and lower bounds for each PROSAIL parameter. VZA, SZA, and RAA, which are always known and unique values, were fixed for each of the model inversions. For parameters that were field-measured, the mean value was used initially, and the possible values were constrained by setting their upper and lower limits. Unknown parameter values (N, C bp , S L , P S , SKYL) were based upon accepted values found in the literature [20,21,116,[123][124][125][126][127][128][129][130][131]. LM was implemented to successively iterate the PROSAIL parameters (i.e., 5000 times) until the simulated spectral curves matched as closely as possible the measured spectral curves (i.e., the value of the cost function was minimized for each plot).
The LUT-based approach was implemented through the ARTMO toolbox [32,59,[132][133][134]. Ten look-up tables with sufficient dimensionality to represent the range of measured and theoretical leaf and canopy conditions (~7.8×10 5 ) were generated for the two sensor configurations (ASD and CHRIS M1) and the 5 VZAs (+55 • , +36 • , 0 • , −36 • , and −55 • ) resulting in a LUT with approximately 7.8 × 10 6 unique spectral/parameter values. The parameter bounds and step used for each of the LUT variables are shown in Table 2. To investigate the use of a random LUT subset for modelling efficiency, 10 smaller LUTs of 100,000 were selected from the original LUT using Latin hypercube sampling (LHS) [135] and tested in conjunction with the full LUT. The original wet and dry spectra samples contained within the SAILH code were replaced with examples of wet and dry soil from the field spectral dataset. Two classes of cost functions were evaluated for each of the vegetation variables (i.e., LCC, PAI, CCC), these included three divergence measures (RMSE, Neyman chi-square, and Bhattacharyya) and three M-estimators (Geman and McClure, least absolute error (LAE), and least-squares estimator (LSE)). Divergence methods are based on the minimizing the distance between two spectra (measured vs. simulated), whereas M-estimators (maximum-likelihood type) seek to find the relationship between independent and dependent variables [59]. Both classes of cost functions (and those chosen within each class) were evaluated based on their published performances [59]. In addition, two regularization strategies were applied to the LUT inversions to help mitigate model ill-posedness: the addition of noise and the use of multiple best solutions. (1) Gaussian noise ranging from 0% to 25% (tested in steps of 5%) was added to the simulated leaf-canopy spectra to account for spectral uncertainties (i.e., potential sensor and environmental noise). Rivera et al. [59], showed that the addition of <~20% noise to the inversion procedure produces favorable results when retrieving biophysical and biochemical variables.
(2) The mean of the multiple best solutions ranging from 0% to 25% (using 5% steps) were evaluated as studies have shown that the single best parameter combinations, as calculated from the optimal cost distance, do not necessarily lead to the best inversion accuracies [19,40,59].
The M1-simulated multi-angle PROSAIL LHS spectral datasets were tested using the two classes of VIs (i.e., multi-band and narrowband) to provide a theoretical comparison to the empirical modelling retrievals and the LUT inversions results. VI retrievals were tested using 10% Gaussian noise added to the variable parameters (i.e., LCC, PAI, and CCC) and to the simulated spectral data to account for (theoretical) sensor and measurement uncertainty. The mean of the best 10% of spectral matches was also used. Inversions were completed with the RMSE cost function. A parameter iteration test was also conducted using the measured Cab, LAI, SZA, VZA, and RAA variable ranges to provide a comparison to the measured multi-angle field spectra.
It should be noted that the purpose of this study was to develop an understanding of the reflectance characteristics of Arctic vegetation at scales of current and expected hyperspectral satellite remote sensors (e.g.,~30 m pixels or greater) and the potential for physically based modeling, mapping, and monitoring of vegetation properties in such environments. Therefore, the focus of the physical modelling approach was at the plot scale/extent where fine spatial scale variability has been accounted for through the sampling approach. At this spatial scale, the tundra vegetation within these sampled sites was considered sufficiently mixed to mimic a homogenous surface (i.e., a turbid medium) (e.g., Figure 1 shows photographs of the various field sites).

Model Assessment and Validation
Validation of the parametric and non-parametric regression models was based upon k-fold cross-validation [136]. It randomly divides the modelling dataset into equal-sized k sub-datasets (k is defined by the analyst) from which k-1 sub-datasets are chosen for training and a k sub-dataset is chosen for validation. The process is repeated k-times; validation is then performed using each of the k sub-datasets [32]. The final step takes the mean score from each validation step to produce a single estimate. In this study, a 10-fold (k = 10) cross-validation procedure was used. Models were assessed in terms of cross-validated model fit (r 2 cv ), cross-validated predictive performance (RMSE cv ), and cross-validated Normalized RMSE (NRMSE cv = RMSE cv /range of measured data) [19,20], which was useful for comparison across VZAs and sites with different vegetation variable ranges. Although direct validation provides a more conservative estimate of model performance, a cross-validation procedure was selected because all data are used for both training and validation [32]. For physical modelling, the inverted PROSAIL TOC spectral reflectance and estimated leaf and canopy parameter values were compared to measured field variables (spectra, LCC, PAI, and CCC) using r 2 , RMSE and NRMSE directly (not with k-fold validation). A cross-validation procedure was not applied to either of the PROSAIL inversion techniques as they both provide a direct solution based on the input/parameterization data.

Field Reflectance Measurements
Multi-angle reflectance measurements showed a standard response across all sites; i.e., backscattered reflectance at −55 • and −36 • was greater at all wavelengths than both nadir (0 • ) and forward scattered measurements at +55 • and +36 • . These differences were likely due to combinations of canopy gaps and shadow effects [137,138], and were more prominent with increasing wavelength. No clear trends related to VZA and wavelength were observed in the coefficient of variation of reflectance for any of the field sites. Figure 5 shows average spectra and their coefficient of variation for each VZA at Banks Island and Richardson Mountains, the most northerly and southerly sites, respectively. Spectra for the other sites were similar and represented intermediate conditions with respect to the climatic N-S gradient described previously. Some distinctive elements are evident in Figure 5 (and in the spectra for the other sites) as follows. Absorption in the red wavelengths was not pronounced due to presence of non-green vegetation components and some non-vegetation background (e.g., soil/rocks), especially at Banks Island where vegetation was sparser. The red-edge between the red and NIR, which is commonly used to estimate LCC and detect vegetation stress [29,52,139], decreased in slope with increasing latitude and this corresponded to a decrease in measured LCC and CCC with latitude. NIR reflectance was lower (<0.5) than for typical dense green vegetation with the lowest values occurring in the most northern field sites. Similar low NIR reflectance has been reported in other studies [29,31,43,44,52] and the north-south gradient was due to an increasing proportion of vascular plants per unit area at lower latitudes, thereby enhancing the multiple scattering of NIR reflectance [25,44]. Also similar to these other studies, NIR reflectance increased slightly with increasing wavelength (a typical NIR plateau was not observed).
The coefficient of variation of reflectance values tended to decrease with increasing wavelength across all study sites. NIR reflectance had the lowest variability at all field sites. At the more sparsely vegetated Banks Island site, a greater degree of soil background radiance was sampled, which increased reflectance towards the red spectrum, thereby reducing the coefficient of variation towards longer wavelengths in the visible spectrum. At lower latitude sites such as Richardson Mountains, as shown in Figure 5, the lower variability of the~550 nm region may be an indication of leaf/canopy carotenoid content (i.e., xanthophylls) in relation to leaf/canopy chlorophyll content, as it has been shown that chlorophyll degrades more rapidly than carotenoid content when plants are undergoing stress or senescence [140][141][142].
The PROSAIL iterations using the range of field-measured vegetation and sensor parameters (Cab, PAI, SZA, VZA, and RAA) demonstrated that field spectra were well withing variability of the simulated multi-angle spectra. Full and detailed results of the ground spectra (multi-angle reflectance and its coefficient of variation) and associated PROSAIL sensitivity analysis can be found in Kennedy [91].

Multi-band Vegetation Index Models
In general, many of the VIs performed similarly, suggesting that the wavelength regions used in a VI are more important than the VI formulation itself [91]. The two-band combinations that produced the highest r 2 cv values for each VI were noticeably clustered in the visible around 350 nm and 550-600 nm, in the red-edge from 700 nm to 750 nm, and in the NIR from 950 nm to 1075 nm, with noticeable gaps of low r 2 cv between these regions. These regions are as expected and similar to those found in Arctic field-based studies (e.g., [29,31]). Another important result is that r 2 cv for the plot scale was almost always stronger than for the subplot scale. For the simple ratio (SR, discussed in more detail below), the differences in r 2 cv shown in Table 3 were significant (p < 0.05). Arctic vegetation canopies are highly variable at the 1 m 2 scale. This was also true for the narrowband VI regressions and GPR models and it concurs with Laidler et al. [49], showing that measurements

Multi-Band Vegetation Index Models
In general, many of the VIs performed similarly, suggesting that the wavelength regions used in a VI are more important than the VI formulation itself [91]. The two-band combinations that produced the highest r 2 cv values for each VI were noticeably clustered in the visible around 350 nm and 550-600 nm, in the red-edge from 700 nm to 750 nm, and in the NIR from 950 nm to 1075 nm, with noticeable gaps of low r 2 cv between these regions. These regions are as expected and similar to those found in Arctic field-based studies (e.g., [29,31]). Another important result is that r 2 cv for the plot scale was almost always stronger than for the subplot scale. For the simple ratio (SR, discussed in more detail below), the differences in r 2 cv shown in Table 3 were significant (p ≤ 0.05). Arctic vegetation canopies are highly variable at the 1 m 2 scale. This was also true for the narrowband VI regressions and GPR models and it concurs with Laidler et al. [49], showing that measurements should be averaged over larger homogeneous areas when retrieving Arctic biophysical and biochemical variables. Hereafter, only plot level (90 m × 90 m) results are presented. Table 3. Two-band VIs ranked. Overall ranking of the top three VIs based on mean r 2 cv and incorporates tests from all variables (e.g., LCC, PAI, and CCC), angles, scales, and spectral modes.  Table 3 shows the top three results (r 2 cv ) of the~2700 regression tests for the 11 multi-band VIs for the two spectral datasets (ASD full spectrum/simulated CHRIS M1 bands) and scales (subplot/plot) averaged over all sites, vegetation variables, and view angles (see [91] for a full list of indices and results). The evaluation of the VIs was based on the mean goodness-of fit statistics (r 2 ) of the modeling results for all three vegetation variables combined (see [91]).

ASD (Full
Given that most of the VIs performed quite similarly, SR was selected for further analysis since it is the simplest formulation, it was least affected by view angle geometry and it produced comparable results for LCC, PAI, and CCC. PAI model performance was best overall while LCC was the worst (Table 4). CCC model performance was in between PAI and LCC since it is a mathematical combination of both variables. The narrower bands of the ASD data produced slightly better model fit and accuracy than the resampled CHRIS M1 data. The VZA that produced the best overall model fit for all three vegetation variables was +55 • (mean r 2 cv = 0.60) for the full ASD dataset; however, mean r 2 cv values over all vegetation variables for each angle were comparable (within ±0.04). Modelled results were inclusive of all SZAs and SAAs for their respective forward and backscatter measurements; averaging of results dampened sun-sensor induced variations seen at the individual sites. Like the VI ranking ( Table 3) the evaluation of the VZAs was based on the mean goodness-of fit statistics (r 2 ) of the modeling results for all three vegetation variables combined [91].
Banks Island data produced better models (e.g., LCC: r 2 cv = 0.78, RMSE cv = 3.45; PAI: r 2 cv = 0.94, RMSE cv = 0.09; CCC r 2 cv = 0.96, RMSE cv = 0.03) than the other sites (e.g., Richardson Mountains 2013, LCC: r 2 cv = 0.52, RMSE cv = 3.88; PAI: r 2 cv = 0.69, RMSE cv = 0.12; CCC r 2 cv = 0.50, RMSE cv = 0.07). The strong results for Banks Island may be due to the large spatial variation in PAI values, where plots ranged from almost bare to completely covered with vegetation. Near zero values in the biochemical-biophysical variables help to ground the regression line and produce greater linearity. It is also possible that canopy height) played a role in the prediction strength. It was minimal in the high Arctic as compared to the greater vertical canopy complexity/height of the lower latitude sites, which can obscure some canopy elements. Figure 6 shows the best model predictions for each vegetation variable. The site and VZA are given; all best models were with ASD full spectrum data at the plot scale.
Remote Sens. 2020, 12, x FOR PEER REVIEW 18 of 42 The narrower bands of the ASD data produced slightly better model fit and accuracy than the resampled CHRIS M1 data. The VZA that produced the best overall model fit for all three vegetation variables was +55° (mean r 2 cv = 0.60) for the full ASD dataset; however, mean r 2 cv values over all vegetation variables for each angle were comparable (within +/−0.04). Modelled results were inclusive of all SZAs and SAAs for their respective forward and backscatter measurements; averaging of results dampened sun-sensor induced variations seen at the individual sites. Like the VI ranking ( Table 3) the evaluation of the VZAs was based on the mean goodness-of fit statistics (r 2 ) of the modeling results for all three vegetation variables combined [91].
Banks Island data produced better models (e.g., LCC: r 2 cv = 0.78, RMSEcv = 3.45; PAI: r 2 cv = 0.94, RMSEcv = 0.09; CCC r 2 cv = 0.96, RMSEcv = 0.03) than the other sites (e.g., Richardson Mountains 2013, LCC: r 2 cv = 0.52, RMSEcv = 3.88; PAI: r 2 cv = 0.69, RMSEcv = 0.12; CCC r 2 cv = 0.50, RMSEcv = 0.07). The strong results for Banks Island may be due to the large spatial variation in PAI values, where plots ranged from almost bare to completely covered with vegetation. Near zero values in the biochemicalbiophysical variables help to ground the regression line and produce greater linearity. It is also possible that canopy height) played a role in the prediction strength. It was minimal in the high Arctic as compared to the greater vertical canopy complexity/height of the lower latitude sites, which can obscure some canopy elements. Figure 6 shows the best model predictions for each vegetation variable. The site and VZA are given; all best models were with ASD full spectrum data at the plot scale.

Predefined Narrowband VI models
The best overall narrowband VI was the Revised Optimized Soil Adjusted Vegetation Index (ROSAVI) [91,143]: where R represents the reflectance at the indicated wavelength. ROSAVI produced an average r 2 cv = 0.36 at the subplot scale and r 2 cv = 0.54 at the plot scale for both the full ASD and simulated CHRIS M1 datasets. In total,~20,000 (84 VIs × view angle × variable × spectral dataset × spatial scale) tests were performed (see [91] for a full list of indices and results). There were no significant differences between the top 20 (at least) narrowband VIs. The Blue Green Index (BGI) [144] produced almost identical results to ROSAVI. Other high performing VIs were also based on the Soil Adjusted Vegetation Index (SAVI) [145] and included OSAVI [146], TSAVI [147], and RMCARI [143], which is not surprising given the sparse nature and abundant non-green components of Arctic tundra vegetation. Even though the 10 nm bandwidth CHRIS bands often did not match the central wavelength of a given VI perfectly (e.g., 705/750 ROSAVI vs. 703/748 nm for CHRIS), both ASD and simulated CHRIS data produced almost identical model equations and performance for most VIs. This indicates that, contrary to the results for the SR analysis, bandwidth differences of 1-10 nm or a slight shift in central wavelength of the CHRIS bands from the specified VI wavelengths do not affect the vegetation models. All narrowband VI models showed similar moderate strength at the plot scale, thus the choice of ROSAVI for further analysis (below) was somewhat arbitrary but representative of the other VIs. Table 5 shows that model fit and prediction errors were generally slightly poorer than for SR (Table 4). However, simulated CHRIS M1 model performance was much closer to that for ASD full spectrum data, whereas it was markedly poorer for SR models. Model performance did not vary as much as for SR models across vegetation variables and VZAs for both spectral datasets. Contrary to the SR models, the VZA that produced the best overall model fit for all three vegetation variables was −55 • (r 2 cv = 0.57) for the full ASD dataset; r 2 cv values were again comparable, being within ±0.05 for all VZAs. However, similar to the SR analysis, Banks Island (2012) data produced the strongest models (e.g., for −55 • : LCC, r 2 cv = 0.75, RMSE cv = 5.61; PAI, r 2 cv = 0.90, RMSE cv = 0.10; CCC, r 2 cv = 0.94, RMSE cv = 0.04), compared to the poorest models, which were for Richardson Mountains (2013) data, with an average r 2 cv = 0.38.

Non-Parametric Gaussian Process Regression Models
GPR models (Table 6) generally performed better than VI-based regressions. In GPR iterative processing, as the least contributing bands were sequentially removed, r 2 cv varied more for off-nadir view angles (Figure 7, LCC models given for illustrative purposes), and there was a slight increase in r 2 cv as the number of variables was reduced (Figure 7a,c). Associated with this, RMSE tended to decrease as variables were removed (Figure 7b,d). These trends were also manifested in the simulated CHRIS M1 data models (mean r 2 cv = 0.62), which were stronger than models using all ASD narrow spectral bands (mean r 2 cv = 0.56), and this result concurs with Verrelst et al. [24]. For LCC (Figure 7c), model performance was stable until approximately Remote Sens. 2020, 12, 3073 20 of 41 10 bands were left, whereas for PAI and CCC, results were stable until there were as few as four bands left. For all vegetation variables; however, using only one or two bands produced the worst overall results, indicating that optimal models would require 5-10 bands for optimal performance. Remote Sens. 2020, 12, x FOR PEER REVIEW 21 of 42 These trends were also manifested in the simulated CHRIS M1 data models (mean r 2 cv = 0.62), which were stronger than models using all ASD narrow spectral bands (mean r 2 cv = 0.56), and this result concurs with Verrelst et al. [24]. For LCC (Figure 7c), model performance was stable until approximately 10 bands were left, whereas for PAI and CCC, results were stable until there were as few as four bands left. For all vegetation variables; however, using only one or two bands produced the worst overall results, indicating that optimal models would require 5-10 bands for optimal performance. With respect to VZA, −55° again produced the best overall results (mean r 2 cv = 0.68 for the simulated CHRIS M1 data), while +55° produced the poorest models (mean r 2 cv = 0.58). Among the study sites, Banks Island data again produced the best models (e.g., for −55°: LCC, r 2 cv = 0.70, RMSEcv = 5.14; PAI, r 2 cv = 0.96, RMSEcv = 0.10; CCC, r 2 cv = 0.92, RMSEcv = 0.05).
Finally, as for the iterative band optimization using SR, the GPR process ranked the bands by their frequency as the top band in the cross-validated regressions. The most important bands were consistently in the 450-500 nm (adjacent to peak green reflectance ~550 nm) and ~650-750 nm (chlorophyll absorption and red-edge) ranges for all three vegetation variables. This was expected, given the covariation of leaf and canopy chlorophyll with plant area index. With respect to VZA, −55 • again produced the best overall results (mean r 2 cv = 0.68 for the simulated CHRIS M1 data), while +55 • produced the poorest models (mean r 2 cv = 0.58). Among the study sites, Banks Island data again produced the best models (e.g., for −55 • : LCC, r 2 cv = 0.70, RMSE cv = 5.14; PAI, r 2 cv = 0.96, RMSE cv = 0.10; CCC, r 2 cv = 0.92, RMSE cv = 0.05).
Finally, as for the iterative band optimization using SR, the GPR process ranked the bands by their frequency as the top band in the cross-validated regressions. The most important bands were consistently in the 450-500 nm (adjacent to peak green reflectance~550 nm) and~650-750 nm (chlorophyll absorption and red-edge) ranges for all three vegetation variables. This was expected, given the covariation of leaf and canopy chlorophyll with plant area index.

Spectral Curve Fitting Validation
The LM inversion curve fitting algorithm was capable of reproducing TOC directional-reflectance for all value ranges of LCC, PAI, and CCC, (i.e., low to high values) and at all VZAs (mean r 2 = 0.99; mean RMSE = 0.4% ± 0.2%). These results may be partially due to spatial averaging of subplots into plot scale measurements; more noise was observed in subplot spectra and this propagated through the inversion procedure. Figure 8a,b shows two example spectra representing the range from dense to sparse vegetation cover. Inverted modelled parameters are shown for each curve. Local minima were not encountered in the modelling process, supporting the approach of providing the best possible estimates for initial model values [121]. Optimum curve fitting was achieved with more than 100 iterations (5000 iterations were used for all final modelling); poor results were achieved with 1-50 iterations. The use of too many fixed variables also caused the LM inversion to behave erratically; using value ranges may be of greater importance than fixing parameters for such iterative optimization.
shown for each curve. Local minima were not encountered in the modelling process, supporting the approach of providing the best possible estimates for initial model values [121]. Optimum curve fitting was achieved with more than 100 iterations (5000 iterations were used for all final modelling); poor results were achieved with 1-50 iterations. The use of too many fixed variables also caused the LM inversion to behave erratically; using value ranges may be of greater importance than fixing parameters for such iterative optimization.
The greatest observed error between the measured and modelled spectra occurred in the lower wavelength regions (e.g., <~500 nm) and generally occurred as an over-prediction (i.e., modelled PROSAIL TOC reflectance values were greater than measured TOC reflectance values). Figure 8c,d shows the differences between the modelled and measured reflectance values for various sites and VZAs (modelled values were subtracted from measured values to better show over and under prediction of PROSAIL). Reflectance at these wavelengths was either not well measured, related to the inherent noise associated with the spectroradiometer, or not well modelled by the PROSAIL RTM. Similarly, large error was observed in the longer wavelength regions (>~900 nm); however, in contrast to the shorter wavelength regions, these regions commonly produced under-predictions of TOC reflectance. Error patterns (both over and under predictions) were observed across the peak green reflectance (~550 nm), chlorophyll absorption (~680 nm), and red-edge regions (~700-750 nm) for all sites and VZAs. Although it could be stated that the degree of error for most spectral regions was well within ±0.02 (or 2%), it is possible that these three regions are also not well modelled by PROSAIL in Arctic environments due to the complexity of the canopy (i.e., number of species, pigment pools) and the degree of spatial variation to canopy structure across all sites (i.e., gap fractions, soil contributions). The LUT-based inversions produced spectral fits that were poorer than the LM approach, often possessing reflectance errors that were in excess of 0.02 (2%). In general, the greatest errors occurred in similar wavelength regions as the LM inversion for all variables, VZAs, and field sites. Because the The greatest observed error between the measured and modelled spectra occurred in the lower wavelength regions (e.g., <~500 nm) and generally occurred as an over-prediction (i.e., modelled PROSAIL TOC reflectance values were greater than measured TOC reflectance values). Figure 8c,d shows the differences between the modelled and measured reflectance values for various sites and VZAs (modelled values were subtracted from measured values to better show over and under prediction of PROSAIL). Reflectance at these wavelengths was either not well measured, related to the inherent noise associated with the spectroradiometer, or not well modelled by the PROSAIL RTM. Similarly, large error was observed in the longer wavelength regions (>~900 nm); however, in contrast to the shorter wavelength regions, these regions commonly produced under-predictions of TOC reflectance. Error patterns (both over and under predictions) were observed across the peak green reflectance (~550 nm), chlorophyll absorption (~680 nm), and red-edge regions (~700-750 nm) for all sites and VZAs. Although it could be stated that the degree of error for most spectral regions was well within ±0.02 (or 2%), it is possible that these three regions are also not well modelled by PROSAIL in Arctic environments due to the complexity of the canopy (i.e., number of species, pigment pools) and the degree of spatial variation to canopy structure across all sites (i.e., gap fractions, soil contributions).
The LUT-based inversions produced spectral fits that were poorer than the LM approach, often possessing reflectance errors that were in excess of 0.02 (2%). In general, the greatest errors occurred in similar wavelength regions as the LM inversion for all variables, VZAs, and field sites. Because the mean of multiple best solutions approach was used for the LUT inversions (i.e., larger spectral ranges are considered valid in the inversion process), it was expected that spectral fits would not necessarily be as good as the LM curve fitting algorithm as valid spectral matches (and associated parameters) were allowed to range from 0 to 20%.

Vegetation Modelling: PROSAIL LM Inversion
The performance of the PROSAIL LM inversion to estimate LCC, PAI, and CCC (Table 7) was not as good as the empirical modelling (ASD spectra: mean r 2 = 0.33; simulated CHRIS M1 spectra: mean r 2 = 0.28). CCC and LCC were better modelled than PAI; this contrasts with the empirical models and with other studies using PROSAIL that have reported better models for LAI than LCC (e.g., [19,20]).
Of the VZAs, nadir (0 • ) produced the best overall results (r 2 = 0.44), followed by +55 • (r 2 = 0.39) for the full spectrum ASD dataset. Similar to the empirical analyses, Banks Island data produced the strongest models (LCC: r 2 = 0.76, RMSE = 9.94; PAI: r 2 = 0.64, RMSE = 0.64; CCC: r 2 = 0.79, RMSE = 0.12, for ASD spectra averaged over all VZAs). The Banks Island model fit for LCC was as strong as for the empirical models but the RMSE was much higher using the LM inversion.

Vegetation Modelling: PROSAIL LUT Inversion
The LUT inversion produced superior results (Table 8) to the LM inversion and to the LHS LUT for estimating LCC, PAI, and CCC (ASD spectra: mean r 2 = 0.48; simulated CHRIS M1 spectra: mean r 2 = 0.46). LHS field tests were subsequently discarded from the analysis. When used in conjunction with regularization approaches, such as the addition of noise and the use of multiple best solutions, results were comparable to the field-based ROSAVI retrievals (e.g., r 2 difference of 0.06; NRMSE difference of 0.08). The best ASD inversions were achieved with the addition of~9% Gaussian noise and the mean of the top~10% best spectral matches, whereas the CHRIS M1 benefited from~12% Gaussian noise and the mean of the top~12% best spectral matches. When averaged across variables, VZAs, and sensors configurations, the best performing cost function was RMSE (NRMSE=0.21); however, best results were achieved with the Neyman chi-square for both PAI and CCC when averaged across VZAs and sensor configurations (NRMSE=0.28 and NRMSE=0.23, respectively). Of the VZAs, nadir (0 • ) and the extreme view zenith angles (±55 • ) produced the best overall correlation results (mean r 2 = 0.50), followed by ±36 • (mean r 2 = 0.45) for the full spectrum ASD dataset; however, the best predictions accuracies were achieved at +55 for both sensor configurations (NRMSE = 0.17) and the worst at −55 • (NRMSE = 0.31). Like all previous analyses, Banks Island produced the strongest models (LCC: r 2 = 0.57, RMSE = 8.59; PAI: r 2 = 0.86, RMSE = 0.59; CCC: r 2 = 0.74, RMSE = 0.12, for ASD spectra averaged over all VZAs) although a positive bias was observed in all models for all vegetation variables. Figure 9 shows the best LUT inversion results for PAI using data for the entire bioclimatic dataset (a) and for the Banks Island site (b) for the +55 • VZA. PAI produced the best overall model fits followed by CCC though LCC had the lowest overall error (e.g., NRMSE = 0.14 for the CHRIS M1 dataset). Of the VZAs, nadir (0°) and the extreme view zenith angles (+/−55°) produced the best overall correlation results (mean r 2 = 0.50), followed by +/−36° (mean r 2 = 0.45) for the full spectrum ASD dataset; however, the best predictions accuracies were achieved at +55 for both sensor configurations (NRMSE = 0.17) and the worst at −55° (NRMSE = 0.31). Like all previous analyses, Banks Island produced the strongest models (LCC: r 2 = 0.57, RMSE = 8.59; PAI: r 2 = 0.86, RMSE = 0.59; CCC: r 2 = 0.74, RMSE = 0.12, for ASD spectra averaged over all VZAs) although a positive bias was observed in all models for all vegetation variables. Figure 9 shows the best LUT inversion results for PAI using data for the entire bioclimatic dataset (a) and for the Banks Island site (b) for the +55° VZA. PAI produced the best overall model fits followed by CCC though LCC had the lowest overall error (e.g., NRMSE=0.14 for the CHRIS M1 dataset). Results of the PROSAIL/VI sensitivity analysis (Table 9) showed the following. (1) Forward scatter spectra produced the best overall model fits and had the lowest overall error (e.g., SR VI r 2 = Figure 9. The best LUT inversion predictions for PAI for all sites combined (a) and for Banks Island (b) for the ASD +55 dataset. Banks Island produced the best overall model fits compared to other field sites at lower latitudes and to the entire dataset.

Comparison of VI Retrievals Using Simulated Multi-Angle Spectra
Results of the PROSAIL/VI sensitivity analysis (Table 9) showed the following. (1) Forward scatter spectra produced the best overall model fits and had the lowest overall error (e.g., SR VI r 2 = 0.71, NRMSE=0.12 for +55 • VZA) relative to backscatter receivals (e.g., r 2 = 0.57, NRMSE = 0.15 for −55 • VZA). The LUT-based inversions for both the ASD and M1 datasets produced similar trends for retrieval accuracy as they relate to VZA (i.e., increasing from forward to backscatter); however, model fits did not follow the same pattern as the theoretical VI retrievals (i.e., decreasing from forward to backscatter) but instead showed no clear trend. The SR VI empirical retrievals showed that the best overall model was achieved at +55 • VZA for the ASD dataset, although marginal, the ROSAVI and GPR showed the opposite. (2) CCC was the most accurately retrieved vegetation variable followed by PAI and LCC respectively for both the SR and ROSAVI models with simulated data. Modelling with field data produced similar results (r 2 ) for both ROSAVI (both datasets) and GPR (for the M1 dataset), but accuracies (NRMSE) were not necessarily the best for CCC. In contrast, all retrievals (empirical and physical) using the full bioclimatic field dataset showed that PAI had the highest overall errors. (3) The best overall retrievals using the simulated LHS dataset were achieved with the SR VI (r 2 = 0.66, NRMSE = 0.13) as compared to the ROSAVI (r 2 = 0.59, NRMSE = 0.15); however, ROSAVI produced slightly better results for LCC. The empirical modelling using the field dataset produced similar results, where the multi-band SR was overall better than using predefined narrow bands in the ROSAVI formulation. Finally, the VI field-based retrievals produced similar modelling results to the simulated retrievals when using the mean of all vegetation variables and VZAs (e.g., difference of r 2 = 0.08 and NRMSE = 0.01 for SR; and a difference of r 2 = 0.05 and NRMSE = 0.02 for ROSAVI), thus demonstrating that the PROSAIL model provides a reasonably accurate spectral representation of Arctic tundra directional reflectance (given the variable parameter space used) at the plot scale for the visited field sites. Table 9. SR and ROSAVI model performance for plot scale data averaged over all sites and VZAs using the LHS PROSAIL simulated CHRIS M1 dataset. RMSE cv units: LCC (µg/cm 2 ); PAI (unitless); CCC (g/m 2 ). NRMSE cv = RMSE cv /range of measured data. The top angle-based results are highlighted and bolded.

Discussion
A field-based, multi-angle, multi-method spectroscopic analysis of Arctic vegetation was conducted in the Western Canadian Arctic. New contributions include: (1) characterization of leaf and canopy vegetation biochemistry at four different sites across a bioclimatic gradient from north central Yukon Territory to the northern end of Banks Island in the Beaufort Sea; (2) multi-angle reflectance analysis across the same environments; (3) comparison of parametric and non-parametric empirical modelling approaches for estimation of LCC, PAI and CCC; and (4) comparison of empirical models with the physically based PROSAIL model for estimation of the same variables. Models were evaluated at two spatial scales (plot and subplot), using full visible-NIR spectra and simulated CHRIS M1 spectral bands.

Comparison of Field Measurements across a Bioclimatic Gradient
There are very few studies that have quantified the pigments of Arctic vegetation at the field scale. Tieszen and Johnson [148] measured pigment structures for four different species in northern Alaska; chlorophyll values of 32 to 77 µg/cm 2 were reported. Early spectrophotometry studies such as this are invaluable; however, the extinction coefficients used in the spectrophotometry analysis were from Arnon [149] and have been deemed incorrect [102]. In our research, field-measured chlorophyll concentrations of the dominant species sampled were between 10.5 µg/cm 2 to 57.2 µg/cm 2 (37.5 µg/cm 2 average). (Table 1). This average value is similar to the 38.0 µg/cm 2 reported in Tiezsen and Johnson [148]; however, the range of values reported were quite different despite measuring similar species (e.g., cotton grass). The values reported here were derived using updated and accepted extinction coefficients [102]. Leaf chlorophyll exhibited a gradient from higher LCC at lower latitudes (e.g.,~40.7 µg/cm 2 ) to lower LCC at higher latitudes (e.g.,~32.0 µg/cm 2 ), and this gradient was also evident in CCC values. Such latitudinal Arctic LCC and CCC values have not been previously reported. The destructive lab/field-based chlorophyll measurements showed no statistical difference (p < 0.05) to plot-scale handheld chlorophyll measurements when converted using a consensus equation [95], thus indicating the potential for the generalized use of such handheld chlorophyll meters in the Arctic.
Although not the primary focus of the modelling retrievals in this research, measurement of carotenoid and EWT values was important for parameterization of the PROSAIL model. Carotenoid concentrations of the dominant species ranged from 3.3 µg/cm 2 to 23.0 µg/cm 2 (average = 14.6 µg/cm 2 ), while EWT values ranged from 0.01 g/cm 2 to 0.04 g/cm 2 (average = 0.03 g/cm 2 ). Subsequent vegetation-based biochemical research conducted in the Arctic should focus on comprehensive characterization of pigment and non-pigment vegetation biochemistry across species, seasons, and representative environments. Not only would such a dataset serve as baseline/reference data for retrieval studies but would allow for further refinement of RTMs that are increasingly incorporating additional leaf/canopy parameters into their code (e.g., PROSPECT-D) [39].

Multi-Angle Spectroscopic Analysis across a Bioclimatic Gradient
Multi-angle spectroscopic analysis showed that top-of-canopy reflectance possessed angular dependencies, with peak variability centered on the chlorophyll absorption region (~680 nm), and low variability evident at peak green reflectance (~550 nm) and regions related to leaf and canopy structure (≥~700 nm). The low reflectance variability observed at all sites around~550 nm may be due to consistent quantities of leaf/canopy carotenoid content found within and across Arctic vegetation species/assemblages. Whereas the higher variability of reflectance found at~680 nm may be attributable to a high degree of leaf/canopy chlorophyll differences found between species and sites. Since chlorophyll degrades rapidly in comparison to carotenoids when plants undergo stress or senescence [140][141][142], the large differences in variability between these two spectral regions (i.e., chlorophyll:carotenoid ratios) may point towards a means of better understanding seasonal changes across a variety of environmental conditions, species, and functional types that are found in the Arctic [5,150].
Negative VZAs generally produced the greatest reflectance as expected for backscatter spectra where shadows are minimized by the relatively small angular difference between the sensor and illumination source (i.e., the sun) (Figure 5a,b). Nadir reflectance was generally lowest across the full spectral range, providing evidence that tundra does not act as a specular or diffuse reflector (more related to volume scattering) and that, like other ecosystems, nadir measurements are highly affected by canopy gaps and shadows [5]. Increased gap fraction can decrease the quality and magnitude of the reflectance signal [5,137,138]. Spectra measured across all field sites and VZAs were not typical of dense green vegetation, particularly for the northern sites, but were similar to those of other Arctic studies (e.g., [29,31,43]) where visible soil/substrates and non-photosynthetic vegetation components caused a general progressive increase in reflectance across the visible, red-edge, and NIR regions. Red-edge reflectance in the southern sites exhibited a greater slope compared to the northern sites for the same reasons; indeed, LCC and PAI (and associated canopy gap fractions) decreased with increasing latitude. These results concur with previous studies; however, the comprehensive analysis presented here for multiple field sites on a latitudinal gradient, different spectral data sets and different spatial scales (plot and subplot) provided a strong basis for the empirical and physical modelling at both the field and satellite imaging scales [71].

Empirical Modelling: Comparison of Parametric and Non-Parametric Retrieval Methods
Parametric and non-parametric regression analyses showed that prediction of LCC, PAI and CCC was most consistent and least sensitive to angular effects when using the simple ratio (SR) with optimally selected spectral bands (Table 4). Of the 11 two-band VIs tested, NDVI produced similar results to SR (Table 3). It has also been investigated in the Arctic using a similar band selection procedure with moderate to good success but using nadir data only (e.g., [29]). Optimal wavelengths found in this study concurred with Liu et al. [29] and were, as expected, in the red absorption region and the red-edge/NIR region influenced by vegetation structure. However, the large range of optimal wavelengths (~615-700 nm and~710-880 nm, respectively) is directly related to reflectance anisotropy, which has not previously been reported for Arctic tundra ecosystems with respect to vegetation parameter retrievals [69].
Non-parametric GPR produced the overall best r 2 cv (Table 6), and this was with the simulated CHRIS M1 dataset, but it had slightly poorer mean predictive capability than SR. The use of the full spectrum data led to poorer GPR performance as also reported by Verrelst et al. [24]. Multi-collinearity and noisy bands should be reduced/removed using band selection based on expert knowledge or quantitative techniques (e.g., principal component analysis). In this study, the simulated CHRIS M1 dataset of 62, 10 nm bands within the range of 411 nm to 997 nm accomplished both reduction of multi-collinearity and removal of noisy bands at both ends of the ASD spectral range, resulting in enhanced results over full spectrum narrowband data. This concurs with other studies for crops (e.g., [24,151]) but this is the first time it has been demonstrated in Arctic vegetation.
Many of the pre-defined narrowband VIs performed almost as well as the optimized SR and GPR. Of the 65 VIs, four SAVI-based VIs were among the top ten (ROSAVI, OSAVI, TSAVI, and RMCARI) [91]. ROSAVI performed the best overall irrespective of VZA and vegetation variable (Table 5); its wavelengths (705 nm and 750 nm) are within the range of the optimal spectral bands selected for SR. Soil-based VIs (e.g., SAVI variants) have been extensively tested in agriculture contexts but have been applied less in vegetation modelling/retrieval testing in natural ecosystems (e.g., [143]). Furthermore, SAVI-based VIs have been tested less in Arctic environments compared to other traditional VIs such as NDVI (e.g., [29,31]). The results here show that they are well suited to the sparse canopies and overall additive influence of soils and other non-green vegetation components of Arctic tundra.
The best estimates of LCC, PAI, and CCC were achieved at the plot scale for all empirical modelling techniques. Subplot (1 m 2 ) measurements were more spatially variable; reduction of such variability is important when scaling up to moderate resolution satellite scales. The best estimates for all techniques were also achieved in the high Arctic, concurring with Liu et al. [29] who also produced good models in the high Arctic using narrowband VIs. This suggests that low-statured two-dimensional canopies produce better models because vegetation elements (i.e., leaves) are not obscured as in more vertically complex canopies.
With respect to VZA, optimal r 2 cv values were generally achieved at −55 • with empirical approaches, although r 2 cv did not vary much with VZA. Relative model prediction error was lowest for 0 • but also did not vary much with VZA. Other studies have found that off-nadir observations are more useful since less background is visible within the field of view [152] and backscatter has been found to be better than forward scatter [153][154][155]. However, this does not always hold true. For example, Stagakis et al. [5] found retrieval of chlorophyll, carotenoids, and LAI to be best at +55 • within a Mediterranean ecosystem due to reduced soil background effects; i.e., for that ecosystem soil background influences were deemed to be more important than shadow effects. This reasoning may also hold in the Arctic due to the low stature and density of the canopy (potentially less overall influence of shadow due to vertical canopy structure); canopy gaps in tundra vegetation (especially at high latitudes) will inevitably lead to increased soil exposure. This is supported by Verrelst et al. [68], who reported significant reflectance and VI anisotropy in pine forest and meadow and concluded that the proportions of photosynthetic, non-photosynthetic and background compounds are stronger drivers than shadow. However, it is likely that shadows caused by tussocks, shrubs, hummocks at forward scatter view angles play a more significant role in directional reflectance than backgrounds and soils where there is 100% canopy cover due to their size relative to the individual canopy elements (i.e., leaves/needles and stems). The grade of the slope in which plots were situated/sampled may have contributed to the variability of the shadowing; however, steep slopes where not sampled. Combined results of the empirical analyses showed that canopy structure (i.e., PAI) is best retrieved at −55 • ; however, results were not substantially different for the other VZAs. This suggests that in low statured tundra canopies structural attributes (e.g., PAI and gap fractions) can be observed and extracted similarly from both forward and backscatter VZAs. Finally, since spectral measurements made at each VZA integrated a range of SZA and SAA values, future research (i.e., biophysical and biochemical retrievals using empirical modelling approaches) should examine the effect that fine-scale solar angular variability observed in both forward and backscatter domains has on model performance. Expressing sun-sensor geometries as phase angle may prove useful (e.g., [156]).

Physical Modelling: Assessment of PROSAIL and Comparison of Inversion Techniques
The PROSAIL model when inverted with the LUT approach achieved results that were comparable to the empirical models ( Table 8). The LM-based approach, however, did not achieve the same overall performance (Table 7). PAI was typically better modelled and estimated than LCC and CCC with the LUT approach, and like the field-based SR results, CCC retrievals were comparable to LCC when averaged across field sites and VZAs. It was expected that CCC results would be similar or better than LCC and PAI since it is derived from a combination of LCC and PAI and both are known to conjointly control canopy reflectance across the visible-near infrared spectrum [5,20,21,24,157]. It should be noted that CCC is correlated with plant density and community/assemblage as it uses both LCC and PAI in its formulation and is not depicting only PAI values. Where vegetation is denser (in a continuous or a discontinuous canopy) there will be greater quantities of chlorophyll aggregated across the canopy, given the same plant phenological conditions and assemblages/species between ground plots. We observed the same plants/communities existing at various densities-in these circumstances CCC is more a function of density (or PAI). Where different plant assemblages with different levels of chlorophyll exist in similar densities, CCC is more related to leaf chlorophyll, as this will be different across the canopy between the ground plots. PAI models were the poorest with the LM approach but reflectance error was only about ±0.02 on the red-edge and the left wing of the NIR, which are regions known to be sensitive to Arctic leaf and canopy structure [29,31]. The optimal LUT-based results differ from those studies that have shown difficulties in estimating LCC with PROSAIL in heterogeneous plant assemblages given prediction errors were within 15% [20,158,159], but align with studies that show LAI is easier to predict than LCC [20,160].
The overall poor results of the LM inversions may be explained by the inability of the one-dimension PROSAIL turbid medium model to accurately describe the directional radiation field of a complex heterogeneous canopy. However, the results of the LUT retrievals, the similarity of the simulated PROSAIL SR retrievals (Table 9) to the field-based SR retrievals (Table 4), and the good performance of all inversion models at the field site with the greatest PAI variability/separability (i.e., Banks Island) indicate that inversion errors may instead be due to the sampled variable's parameter space (i.e., the separation in PAI values) and model ill-posedness. Improvement in variable retrieval using either inversion approach may require consideration of canopy heterogeneity effects such as clumping, varying leaf and canopy structures (e.g., PAI), and pigment variability as it relates to species and/or functional groups present (e.g., [39,161]). Darvishzadeh et al. [20] showed that low-level, heterogeneous grassland canopies can be simulated with PROSAIL when stratification of plant heterogeneity (based on the number of species present) is integrated. A similar approach using the stratification of Arctic tundra vegetation based on species numbers would be a difficult undertaking due to the natural mixture of species at the subplot scale (1 m 2 ) and would not be practical to implement within an ecosystem monitoring program; however, it may be possible that measuring biochemical and biophysical plant traits from distinct functional groups/species assemblages may be a way of stratifying the sample measurement space for maximal parameter separation (see [162,163]). Nevertheless, the use of PROSAIL may be of value to monitoring monotypic shrub patches as their abundance increases across the Arctic [74,79] and as high spatial and spectral resolution sensors become more available.
Our explanations for better results at higher latitudes/lower biomass sites are as follows. (a) It is important to sample the full range of values (both spectral and biophysical/biochemical) at each site as this will allow for greater differentiation in the model predictions (i.e., separability in predicted values and/or classes). The Banks Island site contained greater PAI contrast between plots compared to other sites. (b) The physical model ill-posedness is related to the variability of the sampled parameters; i.e., where there is less spatial variability between spectral and biophysical/biochemical values, less differentiation will be seen between plots, thus leading to higher confusion in the inversion. It may be that the PROSAIL algorithms cannot distinguish the fine spectral differences between vegetation communities/parameters (e.g., where differences of LCC=34 µg/cm 2 vs. 35 µg/cm 2 occur) at a more uniform site than at a spatially heterogeneous site, as can be seen in Figure 5. We concluded that the one-dimensional PROSAIL turbid medium can be successfully applied in tundra environments where there is significant spatial variability between plots (when fine spatial scale variability has been accounted for/removed through proper field sampling approaches) and where broader classes of biophysical and biochemical variables are used that have a greater degree of separation. Indeed, when we grouped the LCC datasets into 5 µg/cm 2 classes and inverted the mean spectra from these classes, results were good for both LM and LUT inversion techniques (e.g., r 2 =0.89 to 0.92, respectively), though this should be interpreted with caution as these classes may not correspond to distinct vegetation types or assemblages that are geographically locatable. The implementation of a three-dimensional modelling approach (e.g., DART, [164,165]) may be more appropriate for high spatial resolution sensor applications in Arctic environments (e.g., [58]) but may have their own limitations (e.g., proper parameterization and inversion). PROSAIL is commonly applied and understood by many researchers so, in selecting an appropriate approach for this first such modelling study in Arctic environments, this study provides a baseline for which other biochemical and biophysical retrieval studies can be improved. From here, we expect others to try to improve on the results we found with refinements or different physical modelling approaches.
The poor results obtained by the LM inversion can typically be attributed to model ill-posedness; it is a well-known problem that multiple parameter combinations will lead to the same spectral reflectance [19,40]. LM was therefore considered extremely susceptible to ill-posedness as only one set of parameters was produced from each inversion. Attempts were made to account for ill-posedness and local minima by employing commonly used regularization strategies, including (a) setting the minimum and maximum values for each parameter based on field measurements and accepted ranges based on the literature, (b) setting algorithm start values based on average leaf/canopy conditions, and (c) fixing algorithm variables with known fixed parameters (e.g., VZAs).
The LUT approach was able to produce superior results to the LM approach when steps were taken to account for ill-posedness and should, therefore, be considered as the optimal approach for predicting biochemical and biochemical variables in tundra environments. The inversion of PROSAIL model against plot-scale ASD full spectrum and CHRIS M1 simulated hyperspectral data and a field-based validation dataset led to the following conclusions about cost functions, regularization options, and parameterization data. (a) The systematic evaluation of cost functions for inversion of biophysical and biochemical parameters in Arctic environments has not previously been conducted but have been evaluated in agricultural landscapes [59]. Divergence measures (i.e., RMSE and Neyman chi-square) were the best performing, with RMSE having the best results across all vegetation variables and Neyman chi-square performing best with PAI and CCC regardless of VZA. These results differ from that of Rivera et al. [59], who showed that LCC was best retrieved with LAE, while LAI was best retrieved with LSE. It is recommended that future inversion studies examine the use data normalization when implementing LUT inversion approaches. (b) The optimal inversion cost functions benefited from the addition of~9-12% Gaussian noise and the use of the mean of the top 10-12% of spectral matches. These results are similar to those achieved by Rivera et al. [59] who showed that the using 6% of multiple solutions and the addition of <20% noise yielded best results for retrieval of LCC and LAI in agricultural areas. The addition of noise and multiple best solutions provides further insight into why stratification is important to achieving optimal results in tundra environments. The degree of spectral separation between similar vegetation values (e.g., PAI values of 0.10 vs. 0.15) may be too small of a difference to overcome with current inversion techniques. (c) The use of prior information was shown to benefit the LUT-based inversion approach. For example, we compared the following: (i) the inversion of multiple view angles (i.e., all five angles) from the field spectroscopy data using one parametrized PROSAIL LUT that used one view angle (e.g., nadir); and (ii) the inversion of one field spectroscopy view angle (e.g., nadir) to multiple PROSAIL LUTs that were parameterized with all view zenith angles. In both cases, results were very similar to those achieved using the correctly parameterized angular parameters (VZA). This result further highlights the problem of ill-posedness that can be encountered within the realm of physical modelling. It also highlights that angular information based on view angle may not be as important as correctly parameterizing the other leaf/canopy variables in the PROSAIL model. The results of this simple test highlight the need to provide as many correctly parameterized values as possible in the inversion approach regardless of approach taken and regularization strategies used. It is recommended that inversion studies be accompanied by high-quality calibration/validation field data; however, where these are not available, the values provided here can serve as a baseline for tundra vegetation biochemical and biophysical parameters.
Spectral subsets using coarser spectral sampling (i.e., the M1 simulated CHRIS/PROBA, bandwidths of~10 nm) produced slightly less accurate predictions for all modelling approaches as compared to the full spectrum narrow bands (<2 nm). This outcome appears to support the notion that narrower bands provide better exploitation of subtle absorption and reflectance features of Arctic vegetation. However, the results of the predefined narrowband vegetation analysis of this study suggest that a bandwidth difference of~2 to~10 nm has no effect (as reported for the ROSAVI results). An explanation for better performance with the ASD dataset compared to the simulated CHRIS/PROBA M1 data may be that band center positions for the CHRIS sensor are not in the optimal locations for detecting Arctic tundra vegetation biochemical and biophysical properties. It is possible that full spectrum data that include short-wavelength infrared regions may provide additional benefits for retrieving and mapping biochemical and biophysical variables in Arctic environments and should be explored in future research (e.g., AVRISI-NG, EnMAP, PRISMA), especially with the physical modelling approach. In addition, the removal of spectral bands that are not well simulated by PROSAIL should be explored as this may increase inversion retrieval performance (e.g., [21]).

Conclusions
Arctic tundra vegetation properties of leaf and canopy chlorophyll content and plant area index were modelled using multi-angle field spectroscopy. Empirical parametric and non-parametric regression techniques were compared to the physically-based PROSAIL model. In empirical analysis, optimization of band selection for two-band vegetation indices was compared to pre-defined narrowband indices. Analyses were conducted using data acquired at four sites from the low to high Arctic, at two spatial scales (field and satellite remote sensing, and for two spectral datasets (full spectrum and simulated CHRIS M1 bands). Overall conclusions are as follows. (1) The simple ratio (SR) and revised optimized soil adjusted vegetation index (ROSAVI) were the optimal two-band and predefined narrowband multi-angle VIs, respectively. (2) The GPR algorithm produced favorable results when used in combination with a spectrally reduced dataset (i.e., CHRIS M1). (3) The optimal multi-angle spectral regions for variable retrievals were largely 450-550 nm and 650-750 nm. The −55 • view angle produced better empirical retrievals results. However, the overall conclusion drawn from the compiled empirical and physical modelling results, when examined across the field sites, was that a multi-angle approach (view angles within ±55 • ) provides comparable/viable characterizations of vegetated surfaces and are therefore equally suitable for biochemical and biophysical variable retrievals in the Arctic. (4) The PROSAIL model, when implemented with a field-based spectral data and inverted using a LUT approach (accompanied with appropriate cost functions and regularization approaches), was determined to be suitable for the prediction of LCC, PAI, and CCC in Arctic environments when vegetation variables are separated into distinct numerical classes (i.e., grouped ranges of values ordered in sequence).