Evaluation of VEGETATION and PROBA-V Phenology Using PhenoCam and Eddy Covariance Data

High-quality retrieval of land surface phenology (LSP) is increasingly important for understanding the effects of climate change on ecosystem function and biosphere–atmosphere interactions. We analyzed four state-of-the-art phenology methods: threshold, logistic-function, moving-average and first derivative based approaches, and retrieved LSP in the North Hemisphere for the period 1999–2017 from Copernicus Global Land Service (CGLS) SPOT-VEGETATION and PROBA-V leaf area index (LAI) 1 km V2.0 time series. We validated the LSP estimates with near-surface PhenoCam and eddy covariance FLUXNET data over 80 sites of deciduous forests. Results showed a strong correlation (R2 > 0.7) between the satellite LSP and ground-based observations from both PhenoCam and FLUXNET for the timing of the start (SoS) and R2 > 0.5 for the end of season (EoS). The threshold-based method performed the best with a root mean square error of ~9 d with PhenoCam and ~7 d with FLUXNET for the timing of SoS (30th percentile of the annual amplitude), and ~12 d and ~10 d, respectively, for the timing of EoS (40th percentile).


Introduction
The study of vegetation phenology and its patterns on a global scale have become more important since the late twentieth century for analyzing the effects of climate change [1,2]. Remote sensing is a useful tool for characterizing land surface phenology (LSP) [3] and global changes of vegetation [4][5][6]. De Beurs et al. [7] analyzed a broad range of statistical methods designed to extract phenological metrics from satellite time series based on threshold percentiles [8][9][10], moving averages [11], first derivatives [12,13], smoothing functions [14] and fitted models [15].
Most literature on LSP has focused on the use of time series of vegetation indices mainly derived from MODIS data [6,16,17]. In previous studies, we showed the added value of using Copernicus Global Land Service (CGLS) leaf area index (LAI) time series derived from VEGETATION and PROBA-V data [10,18]. Bórnez et al. [18] found that the phenological metrics extracted from the CGLS LAI Version 2 (V2) time series agreed best with the available human-based ground observations of phenological transition dates for deciduous broadleaf forest in Europe (PEP727) and United States

Satellite Data: CGLS LAI
We used Copernicus land surface products (CGLS LAI V2) derived from SPOT-VEGETATION (1999 and PROBA-V (2014-2017) data. The spatial resolution is 1 km and the temporal frequency is 10 d [52].
The algorithm for LAI V2 product [53,54] capitalizes on the development and validation of already existing products (CYCLOPES version 3.1 and MODIS collection 5 products) and the use of neural networks [55]. The inputs of the neural networks are daily top of the canopy reflectances from VEGETATION and PROBA-V in the red, near-infrared and shortwave infrared spectral bands and the sun and view geometry. A multi-step procedure for filtering, temporal smoothing, gap-filling and compositing is then applied to the daily estimates to generate the final 10 d products [53].

PhenoCam Data
We used PhenoCam Dataset V1.0 [34,56]. It provides digital images every 30 min. In each image, a "region of interest" was defined manually based on the dominant vegetation type in the camera field of view [35] (Figure 2). The size of the ROI typically ranges from ~50 to ~500 m 2 [34,35]. The green chromatic coordinate (GCC) index [35] was calculated from the red (R), green (G) and blue (B) digital numbers (DN) as: (1) We used the 90th percentile GCC value and 1 d high quality composites to avoid noise from variations in meteorological, atmospheric or illumination conditions [57]. We manually filtered the poor-quality GCC observations and then gap-filled the missing data with a locally weighted scatterplot smoother (lowess)-based filter [35].
The algorithm for LAI V2 product [53,54] capitalizes on the development and validation of already existing products (CYCLOPES version 3.1 and MODIS collection 5 products) and the use of neural networks [55]. The inputs of the neural networks are daily top of the canopy reflectances from VEGETATION and PROBA-V in the red, near-infrared and shortwave infrared spectral bands and the sun and view geometry. A multi-step procedure for filtering, temporal smoothing, gap-filling and compositing is then applied to the daily estimates to generate the final 10 d products [53].

PhenoCam Data
We used PhenoCam Dataset V1.0 [34,56]. It provides digital images every 30 min. In each image, a "region of interest" was defined manually based on the dominant vegetation type in the camera field of view [35] (Figure 2). The size of the ROI typically ranges from~50 to~500 m 2 [34,35]. The green chromatic coordinate (GCC) index [35] was calculated from the red (R), green (G) and blue (B) digital numbers (DN) as: We used the 90th percentile G CC value and 1 d high quality composites to avoid noise from variations in meteorological, atmospheric or illumination conditions [57]. We manually filtered the poor-quality G CC observations and then gap-filled the missing data with a locally weighted scatter-plot smoother (lowess)-based filter [35].

FLUXNET Data
We used FLUXNET 2015 collection of gross primary production (GPP) flux data over 16 forest tower sites ( Figure 1) for the period 2003-2014 (110 site-years) [44]. We used the daily GPP (g C m-2 d-1) derived from half-hourly eddy covariance flux measurements using the night time based approach [58,59]. We smoothed the series of the daily GPP by using a Savitzky-Golay filter based on a second degree polynomial and a 30-day smoothing window [60][61][62].

Methods for Estimating Vegetation Phenology
We tested four methods for estimating phenology (Table 1 and Figure 3, [18]) from satellite (CGLS LAI time series) and ground-based data (PhenoCam GCC, FLUXNEX GPP) (e.g., Figure S2 in Supplementary Materials). The phenological metrics are the timing of the start of the growing season (SoS), the end of the growing season (EoS) and the length of the growing season (LoS). LoS was estimated as the length of time between the EoS and the SoS. The CGLS LAI 10 d time series were linearly interpolated at daily steps before phenological retrieval.

Method. Reference Principles and parameters Threshold based on percentiles
Verger et al. [10] SoS is defined as the first day of the year (DoY) when the vegetation variable exceeds a particular threshold. EoS is defined as the DoY when an index descends below a threshold. We established dynamic thresholds per pixel based on a percentile (10th, 25th, 30th, 40th and 50th) of the annual amplitude

FLUXNET Data
We used FLUXNET 2015 collection of gross primary production (GPP) flux data over 16 forest tower sites ( Figure 1) for the period 2003-2014 (110 site-years) [44]. We used the daily GPP (g C m −2 d −1 ) derived from half-hourly eddy covariance flux measurements using the night time based approach [58,59]. We smoothed the series of the daily GPP by using a Savitzky-Golay filter based on a second degree polynomial and a 30-day smoothing window [60][61][62].

Methods for Estimating Vegetation Phenology
We tested four methods for estimating phenology (Table 1 and Figure 3, [18]) from satellite (CGLS LAI time series) and ground-based data (PhenoCam GCC, FLUXNEX GPP) (e.g., Figure S2 in Supplementary Materials). The phenological metrics are the timing of the start of the growing season (SoS), the end of the growing season (EoS) and the length of the growing season (LoS). LoS was estimated as the length of time between the EoS and the SoS. The CGLS LAI 10 d time series were linearly interpolated at daily steps before phenological retrieval. Table 1. Description of the evaluated methods for the extraction of phenology metrics.

Method. Reference Principles and Parameters
Threshold based on percentiles Verger et al. [10] SoS is defined as the first day of the year (DoY) when the vegetation variable exceeds a particular threshold. EoS is defined as the DoY when an index descends below a threshold. We established dynamic thresholds per pixel based on a percentile (10th, 25th, 30th, 40th and 50th) of the annual amplitude

Logistic function
Zhang et al. [50] SoS is defined as the DoY with the first local maximum rate of change in the curvature of a logistic function fitted to the time series. EoS is defined as the DoY with the first local minimum rate of change in the curvature

First derivative
Tateishi and Ebata. [12] SoS is defined as the DoY of the maximum increase (maximum first derivative) in the curve. EoS is defined as the DoY of the maximum decrease in the curve

Autoregressive moving average
Reed et al. [11] A moving average is first computed at a given time lag (we tested 10-50 d and selected a 30 d time lag). SoS and EoS are then defined as the DoY when the moving-average curves cross the original time series of the vegetation index [12] maximum decrease in the curve

Autoregressive moving average
Reed et al. [11] A moving average is first computed at a given time lag (we tested 10-50 d and selected a 30 d time lag). SoS and EoS are then defined as the DoY when the moving-average curves cross the original time series of the vegetation index

Validation Approach
The LSP derived from VEGETATION and PROBA-V LAI V2 time series was compared with the LSP estimates using ground data from PhenoCam and FLUXNET when the same phenological extraction method was applied (Section 3.1). The statistical metrics used for assessing the performance are the root mean square error (RMSE), the mean error (bias); the coefficient of determination (R 2 ), slope and intercept of the Reduced Major Axis regression (RMA). Further the spatial patterns and latitudinal gradients of LSP estimates were assessed in Section 3.2. We used RStudio for the statistical analysis, Google Earth Engine (GEE, https://earthengine.google.org) for the

Validation Approach
The LSP derived from VEGETATION and PROBA-V LAI V2 time series was compared with the LSP estimates using ground data from PhenoCam and FLUXNET when the same phenological extraction method was applied (Section 3.1). The statistical metrics used for assessing the performance are the root mean square error (RMSE), the mean error (bias); the coefficient of determination (R 2 ), slope and intercept of the Reduced Major Axis regression (RMA). Further the spatial patterns and latitudinal gradients of LSP estimates were assessed in Section 3.2. We used RStudio for the statistical analysis, Google Earth Engine (GEE, https://earthengine.google.org) for the retrieval of LSP over the North Hemisphere, and ESRI ArcGIS 10.5 and gvSIG-desktop-2.3.1 for the graphs and maps.

Comparison of Satellite and Ground Phenologies
The coefficient of determination, R 2 , between the satellite-and ground-based estimates from PhenoCam and FLUXNET phenology ranges from 0.01 to 0.81 (p < 0.001) ( Table 2). The threshold-based method provided the best performances. The 30th percentile of annual amplitude was the best threshold for the SoS (RMSE < 9 d, bias < 2 d and R 2 = 0.74 with p < 0.001 for CGLS LAI V2 estimates compared to PhenoCam; and RMSE < 7 d, bias <4 d and R 2 = 0.81 with p < 0.001 when compared to FLUXNET) and the 40th percentile for the EoS (RMSE = 12 d, bias < 1 d and R 2 = 0.51 with p 0.001 compared to PhenoCam; RMSE < 10 d, bias < 5 d and R 2 = 0.53 with p < 0.001 compared to FLUXNET). Table 2. Statistics of the comparison between the SOS and EOS dates retrieved using the LAI, GCC, and GPP estimates for the four methods: thresholds, logistic function, derivative and moving average. * indicates significant correlations at p < 0.05; **, significant correlations at p < 0.001. The bold type highlights the best method. Evaluation over the 64 PhenoCam sites (356 samples (sites × years)) and 16 FLUXNET towers (110 samples (sites × years)) over deciduous forests in the North Hemisphere ( Figure 1).  5 show the scatter plots of the comparison of satellite and ground-based SoS and EoS retrievals for the four methods. The points were very close to the 1:1 line using the percentile and logistic-function methods while the derivative and moving-average methods produced worse results with more widely dispersed points, especially for the timing of EoS.

Metric
The satellite SoS (Figure 4, Figure 5 and Figure S1) retrieved with threshold and logistic-function methods showed RMSE < 11 d and bias < 2 d compared to PhenoCam, and RMSE < 8 d and bias < 4 d with FLUXNET (Table 2). Higher discrepancies for the SoS were found with the derivative and moving-average methods: RMSE of 19 d and 15 d, and bias of 2 d and < 1 d, respectively using PhenoCam estimates, and RMSE of 24 d and 16 d, and bias of -14 d and 2 d, respectively with FLUXNET estimates ( Table 2).
The EoS can be also robustly estimated using remote sensing observations (Figure 4, Figure 5 and Figure S1) although we observed a degradation of performances for all the methods for the estimation of the EoS as compared to the SoS: higher dispersion of points, higher RMSE and lower correlation for EoS ( Table 2). The EoS estimates from satellite time series of LAI agreed the best with GCC and GPP derived phenology metrics using the threshold method followed by the logistic function: RMSE of 12 d and 18 d, respectively, and biases < 1 d compared to PhenoCam, and RMSE of 10 d and bias < 5 d with FLUXNET (     Table 2.  Table 2. Figure S3 shows the spatial distribution of the GCLS LAI phenological estimates (SoS, EoS, and LoS) from 2000 to 2017 over the North Hemisphere using the threshold method. The length of the vegetation cycles regularly decreases from 220 days to 80 days when latitude increases from temperate to boreal regions. The SoS ranges widely from late march in the south to approximately mid-July in the north. The SoS is slightly earlier in central Europe than in North America for the same latitude. The EoS date ranges from early August to December.

Latitudinal Gradients of Satellite and Ground-Based Phenology
The latitudinal patterns of the timings of SoS and Eos derived from CGLS LAI V2 ( Figure S3) and PhenoCam GCC over deciduous forests from 30 • N to 53 • N in North America showed a very good agreement with a gradual decrease in the length of growing season of approximately five days per degree of latitude which resulted from symmetric variations of 2.5 days per degree of latitude in the start and end of season ( Figure 6). We found a correlation R 2 of 0.92 for the timing of SoS and 0.88 for the EoS when comparing the average satellite and PhenoCam phenology per latitude.
Remote Sens. 2020, 12, 3077 10 of 17 Figure S3 shows the spatial distribution of the GCLS LAI phenological estimates (SoS, EoS, and LoS) from 2000 to 2017 over the North Hemisphere using the threshold method. The length of the vegetation cycles regularly decreases from 220 days to 80 days when latitude increases from temperate to boreal regions. The SoS ranges widely from late march in the south to approximately mid-July in the north. The SoS is slightly earlier in central Europe than in North America for the same latitude. The EoS date ranges from early August to December.

Latitudinal Gradients of Satellite and Ground-Based Phenology
The latitudinal patterns of the timings of SoS and Eos derived from CGLS LAI V2 ( Figure S3) and PhenoCam GCC over deciduous forests from 30ºN to 53ºN in North America showed a very good agreement with a gradual decrease in the length of growing season of approximately five days per degree of latitude which resulted from symmetric variations of 2.5 days per degree of latitude in the start and end of season ( Figure 6). We found a correlation R 2 of 0.92 for the timing of SoS and 0.88 for the EoS when comparing the average satellite and PhenoCam phenology per latitude.  30-34ºN, 34-38ºN, 38-42ºN, 42-46ºN and 46-53°N.

Discussion
We assessed the agreement of the phenological metrics derived from satellite LAI (CGLS LAI V2 from VEGETATION and PROBA-V time series, 1999-2017) with those derived from PhenoCam (GCC) and FLUXNET flux towers (GPP) across 80 sites of deciduous forests mainly located in North America and Europe. The agreement between satellite and ground-based estimates depends on the method used to extract the transition dates. We compared four phenology methods: thresholds based on percentiles of the annual amplitude [10], first derivatives [12], autoregressive moving average [11] and a logistic function fitting approach [50]. Thresholds and logistic function resulted the most robust methods and the phenological metrics extracted from CGLS LAI V2 time series were strongly correlated with those derived from PhenoCam GCC and FLUXNET GPP. On the contrary the derivative and moving average methods showed higher discrepancies between satellite and ground estimates specifically for the timing of the EoS.
The threshold-based method performed the best in terms of accuracy of satellite estimates for the timing of the SoS and EoS: RMSE ~ 9 d and bias < 2d for the SoS, RMSE ~12 d and bias < 1d for

Discussion
We assessed the agreement of the phenological metrics derived from satellite LAI (CGLS LAI V2 from VEGETATION and PROBA-V time series, 1999-2017) with those derived from PhenoCam (GCC) and FLUXNET flux towers (GPP) across 80 sites of deciduous forests mainly located in North America and Europe. The agreement between satellite and ground-based estimates depends on the method used to extract the transition dates. We compared four phenology methods: thresholds based on percentiles of the annual amplitude [10], first derivatives [12], autoregressive moving average [11] and a logistic function fitting approach [50]. Thresholds and logistic function resulted the most robust methods and the phenological metrics extracted from CGLS LAI V2 time series were strongly correlated with those derived from PhenoCam GCC and FLUXNET GPP. On the contrary the derivative and moving average methods showed higher discrepancies between satellite and ground estimates specifically for the timing of the EoS.
The threshold-based method performed the best in terms of accuracy of satellite estimates for the timing of the SoS and EoS: RMSE~9 d and bias < 2d for the SoS, RMSE~12 d and bias < 1d for the EoS, and correlation of R 2~0 .7 compared to PhenoCam data; and RMSE < 7d and bias < 4d for the SoS, RMSE < 10 d and bias < 5d for the EoS, and correlation R 2~0 .8 compared to FLUXNET data. In both PhenoCam and FLUXNET comparison, the 30th percentile of the annual amplitude provided the best performances for the timing of the SoS and the 40th percentile for the EoS, confirming our previous findings [10,18]. These thresholds slightly outperformed 10th, 25th and 50th percentiles of the amplitude as proposed in PhenoCam Dataset V1.0 for the extraction of the phenological transition dates [35]. For the sites with concomitant measurements from the 3 sources of data: satellite LAI, PhenoCam GCC and FLUXNET GPP, we observed that the phenology derived from LAI V2 using percentiles 30 and 40 accurately reproduce the interannual variation of the SoS and EoS and usually provides an intermediate solution between PhenoCam and FLUXNET estimates with differences lower than 10 days (Figure 7). The latitudinal gradient in the northern hemisphere of the CGLS LAI V2 phenophases highly agree with PhenoCam observations with an advance (delay) of 2.5 days per degree of latitude from low to high latitudes in response to the South-North gradient of temperature and photoperiod [63,64]. These results are comparable to other studies [65]. The spatial variability in phenophases can be explained not only by the difference in climatic patterns but also by the elevation and soil conditions [66].
PhenoCam GCC and FLUXNET GPP, we observed that the phenology derived from LAI V2 using percentiles 30 and 40 accurately reproduce the interannual variation of the SoS and EoS and usually provides an intermediate solution between PhenoCam and FLUXNET estimates with differences lower than 10 days (Figure 7). The latitudinal gradient in the northern hemisphere of the CGLS LAI V2 phenophases highly agree with PhenoCam observations with an advance (delay) of 2.5 days per degree of latitude from low to high latitudes in response to the South-North gradient of temperature and photoperiod [63,64]. These results are comparable to other studies [65]. The spatial variability in phenophases can be explained not only by the difference in climatic patterns but also by the elevation and soil conditions [66].
Detecting phenology from carbon flux and PhenoCam data also faces some challenges [26,46]. The flux measurements are potentially ~20% biased due to the lack of energy balance closure, instrument response time, pathlength averaging and incomplete measurement of nocturnal CO2 exchange [46,67,68] which can lead uncertainties in phenology estimates. However, these errors are difficult to quantify and correct [26]. Furthermore, for some FLUXNET sites, there are substantial data gaps due to instrument malfunction or bad data quality [46]. In these cases, the gap-filling may lead to uncertainty in GPP time series and, consequently, in the phenological estimates [69].
The scale difference between ~1km VEGETATION and PROBA-V satellite pixels and the deca-/hectometric footprints of PhenoCam cameras and flux towers may introduce some difficulties for the comparison. This is partially minimized because our validation is limited to deciduous forests which tend to form large patches of the same vegetation type, reducing the influence of mixed or border pixels [26,70,71] The mixed signal due to multi-canopy layers may also introduce confounding effects since the understorey may have a different phenological cycle [70]. The emergence of forest understorey is interpreted in both ground-based and satellite observations as an increase in the greening signal. The agreement between the PhenoCam, FLUXNET and remotely sensed phenological metrics showed generally a higher accuracy for SoS than EoS, consistent with previous studies [18,29,64,[72][73][74][75]. Differences in the structure, ecophysiology and dynamics of the vegetation canopy at the start and end of the growing season [37] may partially explain this. The phenological dynamics for the Detecting phenology from carbon flux and PhenoCam data also faces some challenges [26,46]. The flux measurements are potentially~20% biased due to the lack of energy balance closure, instrument response time, pathlength averaging and incomplete measurement of nocturnal CO2 exchange [46,67,68] which can lead uncertainties in phenology estimates. However, these errors are difficult to quantify and correct [26]. Furthermore, for some FLUXNET sites, there are substantial data gaps due to instrument malfunction or bad data quality [46]. In these cases, the gap-filling may lead to uncertainty in GPP time series and, consequently, in the phenological estimates [69].
The scale difference between~1km VEGETATION and PROBA-V satellite pixels and the deca-/hectometric footprints of PhenoCam cameras and flux towers may introduce some difficulties for the comparison. This is partially minimized because our validation is limited to deciduous forests which tend to form large patches of the same vegetation type, reducing the influence of mixed or border pixels [26,70,71] The mixed signal due to multi-canopy layers may also introduce confounding effects since the understorey may have a different phenological cycle [70]. The emergence of forest understorey is interpreted in both ground-based and satellite observations as an increase in the greening signal.
Differences in the structure, ecophysiology and dynamics of the vegetation canopy at the start and end of the growing season [37] may partially explain this. The phenological dynamics for the timing of EoS tend to vary with species, age, dispersion and homogeneity, and can also differ across the same species, with differences of up to two weeks within the same ROI [26,43,71,76,77]. Further studies will use high resolution satellite data from Sentinel-2 to mitigate these issues and capture the spatial variability of EoS [78]. Note, however, that monitoring EoS is affected by the intrinsic uncertainties of satellite remote sensing at northern latitudes in autumn: atmospheric effects, snow and poor illumination conditions [79].

Conclusions
Phenological data from PhenoCam, FLUXNET and satellite remotely sensed data have become a broad resource for analyzing the relationships between global change and vegetation [29,35]. The network of PhenoCam webcams and eddy covariance towers cover only small areas around the camera or flux tower [40,80]. Satellite imagery has the advantage of providing continuous spatio-temporal coverage at the global scale. Near-surface digital cameras and flux towers have nevertheless become a good tool for characterizing local phenology and validate satellite estimates [26,37,57,81]. The high temporal frequency of PhenoCam and flux measurements provide continuous time series for applying the same phenology extraction methods to ground and satellite time series. This way, we avoid some of the issues identified in our previous research [18] related to the differences in the definition of satellite phenology metrics and ground phenophases when PEP725 and USA-NPN data were used for the validation (e.g.; the representativity and spatial distribution of the data as well as the gaps in the time series of ground measurements).
Results validate the land surface phenology estimated from CGLS LAI V2 time series, as well as the robustness of PhenoCam and FLUXNET data to analyze vegetation phenology. This study has put bounds on the uncertainty in satellite-derived phenological transitions, which should allow to analyze changes in the phenological distribution pattern and serve as a starting point for other studies that characterize anomalies and trends over vegetation phenology, as well as its possible relationship with changes in the climate pattern as a result of climate change.
Supplementary Materials: The following are available online at http://www.mdpi.com/2072-4292/12/18/3077/s1, Figure S1: "Boxplots of the bias errors of satellite-based minus the near-surface estimates of SoS (a) and EoS (b) over the 64 PhenoCam sites (a,b), and the 16 FLUXNET towers (c,d) for the four extraction methods: threshold method (the 30th percentile of annual amplitude for the SoS (a, c) and the 40th percentile for the EoS (b, d)), the logistic-function, first derivatives and moving-average. An elongated boxplot indicates a larger dispersion of the average bias". Figure S2: "Time series of CGLS LAI, PhenoCam GCC and FLUXNET GPP for the Harvard Forest site (42.5378N, -72.1715O) over the 2008-2012 period". Figure S3: "Maps of average SoS (a), EoS (b) and LoS (c) derived from CGLS LAI V2 time series (1999-2017) using the threshold method (30th percentile of annual LAI amplitude for SoS and 40th percentile for EoS). The maps show the estimated phenology in deciduous or mixed forest based on the annual C3S Global Land Cover for the year 2018 (http://maps.elie.ucl.ac.be/CCI/viewer/download.php). The continental areas in white are lakes, deserts, agricultural areas and evergreen forests. The phenology was not computed for pixels with very limited seasonality: when the annual amplitude ((max (LAI) -min (LAI)) was lower than the 30% of the median value in the time series (0.3 * LAI50th). For pixels with multiple growing seasons, we computed the phenological metrics for the growing season having the highest LAI amplitude". Table  S1: "Characteristics of PhenoCam [38] and FLUXNET fluxnet.fluxdata.org [44] sites. The Start date and End date indicates the period of available data. MAT is mean annual temperature and MAP is mean annual precipitation based on climate data are from WorldClim. Primary and secondary vegetation types are as follows: AG = agriculture; DB = deciduous broadleaf; DN = deciduous needleleaf; EB = evergreen broadleaf; EN = evergreen needleleaf; GR = grassland; MX = mixed vegetation (generally EN/DN, DB/EN, or DB/EB); SH = shrubs; TN = tundra (includes sedges, lichens, mosses, etc.); WL = wetland".