Vegetation Indices Do Not Capture Forest Cover Variation in Upland Siberian Larch Forests

Michael M. Loranty 1,* , Sergey P. Davydov 2, Heather Kropp 1, Heather D. Alexander 3, Michelle C. Mack 4, Susan M. Natali 5 and Nikita S. Zimov 2 1 Department of Geography, Colgate University, Hamilton, NY 13346, USA; hkropp@colgate.edu 2 Northeast Science Station, Pacific Institute for Geography, Far East Branch, Russian Academy of Sciences, Cherskiy 678830, Russia; davydoffs@mail.ru (S.P.D.); nzimov@mail.ru (N.S.Z.) 3 Department of Forestry, Forest and Wildlife Research Center, Mississippi State University, Starkville, MS 39759, USA; heather.alexander@msstate.edu 4 Center for Ecosystem Science and Society, Northern Arizona University, Flagstaff, AZ 86011, USA; Michelle.Mack@nau.edu 5 Woods Hole Research Center, Falmouth, MA 02540, USA; snatali@whrc.org * Correspondence: mloranty@colgate.edu; Tel.: +1-315-228-6057


Introduction
Boreal forest responses to warming temperatures will feed back to climate in a variety of important ways.Boreal ecosystems store globally significant amounts of terrestrial carbon in vegetation and soils.The balance between vegetation uptake and soil carbon losses is particularly important yet also represents a large source of uncertainty in high-latitude carbon cycle climate feedback, particularly in the context of permafrost thaw [1,2].In addition to carbon cycle effects, vegetation change may also feed back to climate via altered ecosystem water and energy dynamics.Increases in evapotranspiration cool the land surface and increase atmospheric water vapor [3,4].Increasing canopy cover reduces land surface albedo during the spring when vegetation masks snow, leading to surface warming [5].When evergreen trees are replaced by deciduous species, by contrast, reductions in snow masking lead to increased albedo.Changes in light attenuation with canopy composition and density may also alter soil temperatures [6], which have important implications for soil carbon cycling.The magnitude of climate feedback associated with each of these processes depends on the rate and spatial extent of vegetation change.Consequently, there has been a great deal of effort to understand how boreal forest vegetation is responding to climate [7][8][9][10][11][12].
Optical remote sensing has emerged as one of the key tools used to understand ecosystem responses to climate.The existence of long-term satellite records such as the Advanced Very High Resolution Radiometer Global Inventory Modeling and Mapping Studies (AVHRR GIMMS) have enabled spatially and temporally continuous analyses of vegetation productivity that were among the first indications of biome-wide vegetation responses to climate [13].Numerous Vegetation Indices (VIs) have been developed to use spectral information as a proxy for vegetation productivity, and the Normalized Difference Vegetation Index (NDVI) is the most commonly used VI.Long-term trends in NDVI are interpreted as increases or decreases in boreal forest productivity, and have been linked to tree growth and changes in forest composition [14].However, comparison of multiple data sources that include products from new moderate resolution satellite sensors have revealed inconsistencies in NDVI trends [15,16].Interpreting the ecological drivers of changes in NDVI is a key challenge for understanding high-latitude ecosystem responses to climate.
Recent increases in spatial and temporal resolution of satellite data, along with increases in computational efficiency are enabling investigation of spatial and temporal variation in high-latitude vegetation at finer scales.Landsat data with 30 m spatial resolution is increasingly being used to link the magnitude and direction of NDVI trends with different vegetation types [17,18] and to corroborate AVHRR trends [19].Very high-resolution multispectral satellite data (i.e., <~4 m) from commercial platforms, such as IKONOS, QuickBird, and WorldView, provide even finer spatial detail and can be combined with historical aerial photographs to confirm hypothesized vegetation changes such as shrub and tree expansion into tundra [20], which in turn can be linked with NDVI trends [21].Though very high-resolution satellite data has relatively poor temporal resolution, it has proven useful for mapping fine-scale variation in vegetation or ecosystem types and vegetation indices, primarily in tundra ecosystems [22][23][24][25].Spatial relationships between vegetation properties and NDVI in Arctic tundra ecosystems are often used to support interpretations of temporal trends in NDVI [26].To date there have been relatively few efforts to utilize very high-resolution vegetation indices and ecosystem properties as a tool to aid ecological interpretation of NDVI trends in boreal forests.
In this study we examine the effects of variability in forest cover on spatial and temporal variability in vegetation indices derived from very high-resolution multispectral satellite data with high temporal resolution.By comparing satellite data with field observations from boreal forests in Siberia that range from sparse to dense canopy cover we seek to answer the following questions: (1) Do values of commonly used vegetation indices increase with forest cover?and (2) Are there differences in overstory and understory phenology that correspond to differences in vegetation indices?

Study Site
The study was conducted along the Kolyma River in northeastern Siberia near the Northeast Science Station (NESS) located in the town of Cherskiy, Sakha Republic, Russia (68.74 • N, 161.40 • E).Mean annual, January, and July temperatures are −10 • C, −32 • C, and 13 • C, respectively, for the 1986-2015 period.Boreal forests in the area are composed exclusively of Larix cajanderi, a deciduous conifer hereafter referred to as larch.The understory is dominated by deciduous shrubs (Betula nana exilis, B. divaricata, Salix spp., and Vaccinium uliginosum), small evergreen shrubs (Vaccinium vitis-idaea, Empetrum nigrum, and Ledum decumbens (Ait.)),mosses (Aulacomnium turgidum, Dicranum spp., and Polytrichum spp.), and lichens (Cetraria spp.and Cladonia spp.).Forests of varying density are dominant in uplands to the east of the Kolyma River.Lowland floodplain ecosystems in the area contain forests but also large extensive shrublands dominated by Betula nana exilis, Salix spp., and Alnus fruticosa (Duschekia fruticosa), and also graminoid-dominated areas dominated by Eriophorum polystachion, Calamagrostis purpurea, and Carex spp.[27].For this study, we examined a ~25 km 2 area located to the east of NESS and the town of Cherskiy, as shown in Figure 1.

Field Data
We characterized tree density and canopy cover for 46 stands throughout the study area during the 2010-2017 field seasons [27,28].There were no major disturbances or changes in land cover across the study area during this period.At each stand we established three 20 m long belt-transects spaced approximately 10 m apart that ranged from 1 to 4 m in width depending on tree density.Along each transect we recorded the number of stems and the diameter at breast height (DBH) for each stem.Existing allometric equations were used to calculate aboveground biomass from DBH, and multiplied by the proportion of carbon in live biomass (0.48) to determine aboveground larch biomass carbon content for each tree [27,28].Larch canopy cover was measured using a hemispherical densitometer by taking four measurements in each cardinal direction at the beginning, middle, and end of each transect.Stand-level estimates for each variable were calculated as the mean of the three transects.
In June 2016 we installed in situ sensors to measure NDVI (SRS, Decagon Devices, Pullman WA, USA) at one high-and one low-density stand, as shown in Figure 2. At each stand, four downward facing sensors measured spectral radiance in the red (650 nm) and near infrared (810 nm) wavelengths.Of the four sensors at each stand, two were mounted at 7 m height on a tower mast above the larch canopy, and two were mounted below the larch canopy at 1 m above the ground surface.These downward facing sensors had a 36° field of view (fov) and were mounted at a 45°

Field Data
We characterized tree density and canopy cover for 46 stands throughout the study area during the 2010-2017 field seasons [27,28].There were no major disturbances or changes in land cover across the study area during this period.At each stand we established three 20 m long belt-transects spaced approximately 10 m apart that ranged from 1 to 4 m in width depending on tree density.Along each transect we recorded the number of stems and the diameter at breast height (DBH) for each stem.Existing allometric equations were used to calculate aboveground biomass from DBH, and multiplied by the proportion of carbon in live biomass (0.48) to determine aboveground larch biomass carbon content for each tree [27,28].Larch canopy cover was measured using a hemispherical densitometer by taking four measurements in each cardinal direction at the beginning, middle, and end of each transect.Stand-level estimates for each variable were calculated as the mean of the three transects.
In June 2016 we installed in situ sensors to measure NDVI (SRS, Decagon Devices, Pullman WA, USA) at one high-and one low-density stand, as shown in Figure 2. At each stand, four downward facing sensors measured spectral radiance in the red (650 nm) and near infrared (810 nm) wavelengths.Of the four sensors at each stand, two were mounted at 7 m height on a tower mast above the larch canopy, and two were mounted below the larch canopy at 1 m above the ground surface.These downward facing sensors had a 36 • field of view (fov) and were mounted at a 45 • angle facing eastward and westward.In addition, three upward facing sensors with a 180 • fov measured solar irradiance in the same red and near infrared wavelengths; one of these was mounted at 7 m height at the high-density site and the remaining two were mounted at 1 m height at the high-and low-density sites.Red and near infrared reflectance (R and NIR, respectively) were calculated by dividing radiance measured with the downward facing sensors by irradiance measured by collocated upward facing sensors.We assumed that the single upward facing sensor at 7 m was representative for both stands since they were located less than 2 km apart.We calculated NDVI for each of the eight downward facing sensors from R and NIR as follows: Remote Sens. 2018, 10, x FOR PEER REVIEW 4 of 15 angle facing eastward and westward.In addition, three upward facing sensors with a 180° fov measured solar irradiance in the same red and near infrared wavelengths; one of these was mounted at 7 m height at the high-density site and the remaining two were mounted at 1 m height at the highand low-density sites.Red and near infrared reflectance (R and NIR, respectively) were calculated by dividing radiance measured with the downward facing sensors by irradiance measured by collocated upward facing sensors.We assumed that the single upward facing sensor at 7 m was representative for both stands since they were located less than 2 km apart.We calculated NDVI for each of the eight downward facing sensors from R and NIR as follows: (A) (B) Measurements made at 60-s intervals were averaged every 30 min and recorded to a Decagon EM50 data logger.We used precipitation data from the Cherskiy airport to exclude data from days with rain.Mean NDVI values from 12:00 to 14:00 local time were used for analysis.
We observed phenology for overstory larch trees and understory shrubs (Betula and Salix) for the 2000-2017 time period across the landscape that included all of the study sites.We relied on visual observations made in forests near NESS because the phenophases described below occur rapidly and in concert across the study area.Start of the growing season (SOS) was defined as the timing of leaf emergence beginning with budburst and ending with leaf unfolding.This typically occurred rapidly over several days at the very end of May or the first half of June, and was followed by continued leaf growth through the end of June.End of the growing season (EOS) was defined as the visible onset of fall senescence as indicated by changes in leaf coloration (i.e., from green to yellow or red).This phenophase typically occurred rapidly over a period of several days in mid to late August.We defined SOS and EOS as the mid-point of the several day period over which leaf emergence and yellowing occurred, respectively.Growing season length (GSL) was defined as the difference between EOS and SOS.

Satellite Data
We utilized imagery from several satellite sensors to examine relationships between forest cover and VIs across the study area.To characterize forest cover across the study area we used a panchromatic WorldView2 image acquired on 8 March 2015.The 12-bit image had ~0.75 m spatial resolution and was converted to top-of-atmosphere (TOA) reflectance.This image was provided by the Polar Geospatial Center at the University of Minnesota (https://www.pgc.umn.edu/).To examine NDVI and EVI across the study area at high spatial and temporal resolution, we used data from the PlanetScope constellation of satellite sensors [29].The PlanetScope constellation consists of Measurements made at 60-s intervals were averaged every 30 min and recorded to a Decagon EM50 data logger.We used precipitation data from the Cherskiy airport to exclude data from days with rain.Mean NDVI values from 12:00 to 14:00 local time were used for analysis.
We observed phenology for overstory larch trees and understory shrubs (Betula and Salix) for the 2000-2017 time period across the landscape that included all of the study sites.We relied on visual observations made in forests near NESS because the phenophases described below occur rapidly and in concert across the study area.Start of the growing season (SOS) was defined as the timing of leaf emergence beginning with budburst and ending with leaf unfolding.This typically occurred rapidly over several days at the very end of May or the first half of June, and was followed by continued leaf growth through the end of June.End of the growing season (EOS) was defined as the visible onset of fall senescence as indicated by changes in leaf coloration (i.e., from green to yellow or red).This phenophase typically occurred rapidly over a period of several days in mid to late August.We defined SOS and EOS as the mid-point of the several day period over which leaf emergence and yellowing occurred, respectively.Growing season length (GSL) was defined as the difference between EOS and SOS.

Satellite Data
We utilized imagery from several satellite sensors to examine relationships between forest cover and VIs across the study area.To characterize forest cover across the study area we used a panchromatic WorldView2 image acquired on 8 March 2015.The 12-bit image had ~0.75 m spatial resolution and was converted to top-of-atmosphere (TOA) reflectance.This image was provided by the Polar Geospatial Center at the University of Minnesota (https://www.pgc.umn.edu/).To examine NDVI and EVI across the study area at high spatial and temporal resolution, we used data from the PlanetScope constellation of satellite sensors [29].The PlanetScope constellation consists of approximately 120 CubeSats operated by Planet, a commercial satellite company.The PlanetScope satellites acquire imagery in four spectral bands (blue, green, red, and near infrared) with a ground sample distance (GSD) of 3-4 m.We used the Planet Surface Reflectance (SR) product collected using satellites deployed to a Sun-synchronous orbit, which is produced at 3 m spatial resolution and derived by first processing the Planet Analytic Radiance Product to TOA reflectance and then atmospherically correcting this to bottom-of-atmosphere or surface reflectance.Atmospheric corrections are preformed using the 6SV2.1 radiative transfer model [30] informed with atmospheric optical and water vapor information acquired from concurrent MODIS scenes.The Planet SR product does not account for haze or low cirrus clouds, so we manually screened images for clarity.In total, we utilized images acquired on 5 clear days between 15 June and 31 August 2017.We also used a single Landsat 8 OLI image for the study site acquired on 28 July 2017 to compare with the Planet SR data and to determine if relationships between VI and canopy cover were consistent between the two products.We acquired the Landsat Surface Reflectance Level-2 data product from the USGS EarthExplorer interface.We used the included quality band to mask clouds and cloud shadows from the Landsat image.Examination of the resulting NDVI image revealed that not all cloud shadows were excluded by the mask, so we used an NDVI threshold of 0.5 to further mask cloud shadows and anthropogenic surfaces.In addition to NDVI, we used the Enhanced Vegetation Index (EVI).EVI incorporates blue reflectance (B), in addition to the NIR and R, that helps to reduce impacts of the atmosphere and soil background [31].We calculated EVI as follows:

Data Analyses
We processed field data and satellite images, and performed all data analyses with R version 3.4.3 [32].Manipulation and processing of geospatial data were performed using the Raster [33] and sp [34] packages in R. To compare field observations with satellite data, we used the geolocation of the center of the middle transect from each site, and then extracted data from all pixels within 15 m of the center point, except for Landsat data, where we used a single 30 m pixel.
To create maps of canopy cover, we developed a relationship between observed canopy cover and the 12-bit digital number representing TOA reflectance (DN) in the panchromatic image.DN values for each stand were determined by taking the mean of all pixels within 15 m of stand center.We standardized DN (DN S ) values between 450 and 1000 to a range of 1-2 to avoid convergence issues when fitting an exponential model using the nls function in R. We fit an exponential model in the form of a*exp(-b*DN S ) using data from all 46 stands and used the parameters with a map of DN S to create a map of canopy cover for the study area.We tested for significant differences between tree and shrub phenology metrics using a student's t-test in R. To examine relationships between canopy cover and VIs we performed linear regression analyses using the lm function in R. All code used for analyses can be found at the following url: https://github.com/mloranty/ness_phenology/releases/tag/v2.0 [35].

Landscape Spatial Variability in Vegetation Dynamics
Canopy cover ranged from 0% to 93% for the 46 stands sampled across the study area.Above ground larch biomass carbon ranged from 0 to 2750 g C m −2 , and was closely linearly related with canopy cover, as shown in Figure 3, (r 2 = 0.81, slope = 27.1, intercept = 37.5, p < 0.001).Canopy cover was inversely related to the standardized digital number for each stand (DN S ) from the panchromatic satellite image, as shown in Figure 4 (a*exp(−b*DN S ), a = 2560, b = 3.4, p < 0.001, RMSE = 14.1).We used this relationship to predict canopy cover across the study area, as shown in Figure 5. Mean predicted canopy cover for raster cells within a 15 m of stand center was closely correlated with observed values of canopy cover, as shown in Figure 6A (r 2 = 0.69, slope = 0.71, intercept = 9.2, p < 0.001), and RMSE of predicted canopy cover for all stands was 14.1%.Residual canopy cover was inversely related to observed canopy cover, as shown in Figure 6B (r 2 = 0.26, slope = −0.29,intercept = 9.2, p < 0.001), meaning that there was a tendency towards over-predictions at low canopy cover (<~20%) and under-predictions at high canopy cover (>~50%).
Remote Sens. 2018, 10, x FOR PEER REVIEW 6 of 15 < 0.001), and RMSE of predicted canopy cover for all stands was 14.1%.Residual canopy cover was inversely related to observed canopy cover, as shown in Figure 6B (r 2 = 0.26, slope = −0.29,intercept = 9.2, p < 0.001), meaning that there was a tendency towards over-predictions at low canopy cover (<~20%) and under-predictions at high canopy cover (>~50%).We observed differences in tree and shrub phenology across the study area, as shown in Table 1.Start of season (SOS) in larch occurred in the first week of June, on average three days before shrubs, but these differences were not significant, as shown in Table 1 (p = 0.1).End of season (EOS), indicated by onset of leaf senescence, occurred significantly later in larch (p < 0.01) with larch needles remaining green for an average of six days later in August than shrub leaves.On average, GSL was significantly longer by 8 days (p < 0.01) for larch trees than for the dominant understory shrubs.) r2 = 0.81 Remote Sens. 2018, 10, x FOR PEER REVIEW 6 of 15 < 0.001), and RMSE of predicted canopy cover for all stands was 14.1%.Residual canopy cover was inversely related to observed canopy cover, as shown in Figure 6B (r 2 = 0.26, slope = −0.29,intercept = 9.2, p < 0.001), meaning that there was a tendency towards over-predictions at low canopy cover (<~20%) and under-predictions at high canopy cover (>~50%).We observed differences in tree and shrub phenology across the study area, as shown in Table 1.Start of season (SOS) in larch occurred in the first week of June, on average three days before shrubs, but these differences were not significant, as shown in Table 1 (p = 0.1).End of season (EOS), indicated by onset of leaf senescence, occurred significantly later in larch (p < 0.01) with larch needles remaining green for an average of six days later in August than shrub leaves.On average, GSL was significantly longer by 8 days (p < 0.01) for larch trees than for the dominant understory shrubs.We observed differences in tree and shrub phenology across the study area, as shown in Table 1.Start of season (SOS) in larch occurred in the first week of June, on average three days before shrubs, but these differences were not significant, as shown in Table 1 (p = 0.1).End of season (EOS), indicated by onset of leaf senescence, occurred significantly later in larch (p < 0.01) with larch needles remaining green for an average of six days later in August than shrub leaves.On average, GSL was significantly longer by 8 days (p < 0.01) for larch trees than for the dominant understory shrubs.

Variability in NDVI
Field observations of NDVI displayed expected seasonal patterns, as shown in Figure 7. Maximum overstory NDVI ranged from 0.5 to 0.7 in 2016, as shown in Figure 7A, and 0.45 to 0.65 in 2017, as shown in Figure 7C, with differences of approximately 0.15 between east and west facing sensors at the low-density stand during the middle of the growing season, while within stand differences at the high-density stand were notably smaller, as shown in Figure 7A,C.At the low-density site understory NDVI values, as shown in Figure 7B,D, were similar to overstory values at both sites.Understory NDVI at the high-density site was notably lower than the low-density understory and displayed the least seasonal variability, as shown in Figure 7B,D.

Variability in NDVI
Field observations of NDVI displayed expected seasonal patterns, as shown in Figure 7. Maximum overstory NDVI ranged from 0.5 to 0.7 in 2016, as shown in Figure 7A, and 0.45 to 0.65 in 2017, as shown in Figure 7C, with differences of approximately 0.15 between east and west facing sensors at the low-density stand during the middle of the growing season, while within stand differences at the high-density stand were notably smaller, as shown in Figure 7A,C.At the low-density site understory NDVI values, as shown in Figure 7B,D, were similar to overstory values at both sites.Understory NDVI at the high-density site was notably lower than the low-density understory and displayed the least seasonal variability, as shown in Figure 7B,D.Mean PlanetScope NDVI or EVI for each stand from five dates throughout the growing season, as shown in Figure 8A,B, was not positively related to observed canopy cover.The only statistically significant relationships between these vegetation indices and observed canopy cover were negative, and occurred late in the growing season.Mean growing season NDVI and EVI were negatively correlated with observed canopy cover, as shown in Figure 8C.Pixel-wise comparison of predicted tree cover with NDVI and EVI across the entire landscape yielded similar results, as shown in Figure 8D-F; the only significant relationships were negative associations between canopy cover and EVI late in the growing season.Landsat and PlanetScope exhibited similar geographic patterns in NDVI across the study area, as shown in Figure 9. Landsat NDVI was consistently higher than PlanetScope NDVI.Landsat NDVI was not significantly related to canopy cover, but Landsat EVI exhibited a significant negative correlation with canopy cover, as shown in Figure 10.Mean PlanetScope NDVI or EVI for each stand from five dates throughout the growing season, as shown in Figure 8A,B, was not positively related to observed canopy cover.The only statistically significant relationships between these vegetation indices and observed canopy cover were negative, and occurred late in the growing season.Mean growing season NDVI and EVI were negatively correlated with observed canopy cover, as shown in Figure 8C.Pixel-wise comparison of predicted tree cover with NDVI and EVI across the entire landscape yielded similar results, as shown in Figure 8D-F; the only significant relationships were negative associations between canopy cover and EVI late in the growing season.Landsat and PlanetScope exhibited similar geographic patterns in NDVI across the study area, as shown in Figure 9. Landsat NDVI was consistently higher than PlanetScope NDVI.Landsat NDVI was not significantly related to canopy cover, but Landsat EVI exhibited a significant negative correlation with canopy cover, as shown in Figure 10.

What Drives Variability in Vegetation Indices?
NDVI and EVI did not increase appreciably with forest cover across a landscape with stands ranging from 0% to upwards of 90% canopy cover.The only significant relationships that we found were negative correlations between observed canopy cover and VIs late in the growing season, and when VIs were averaged across the growing season.Using modeled canopy cover across the entire landscape yielded the same results.The consistency of these results bolsters interpretation of analyses with modeled results, despite the inherent error associated with the relationship between DN and observed canopy cover.Moreover, given that the lowest VI values occurred in pixels with predicted canopy cover around 60%, and that the model was biased towards under-prediction at

What Drives Variability in Vegetation Indices?
NDVI and EVI did not increase appreciably with forest cover across a landscape with stands ranging from 0% to upwards of 90% canopy cover.The only significant relationships that we found were negative correlations between observed canopy cover and VIs late in the growing season, and when VIs were averaged across the growing season.Using modeled canopy cover across the entire landscape yielded the same results.The consistency of these results bolsters interpretation of analyses with modeled results, despite the inherent error associated with the relationship between DN and observed canopy cover.Moreover, given that the lowest VI values occurred in pixels with predicted canopy cover around 60%, and that the model was biased towards under-prediction at

What Drives Variability in Vegetation Indices?
NDVI and EVI did not increase appreciably with forest cover across a landscape with stands ranging from 0% to upwards of 90% canopy cover.The only significant relationships that we found were negative correlations between observed canopy cover and VIs late in the growing season, and when VIs were averaged across the growing season.Using modeled canopy cover across the entire landscape yielded the same results.The consistency of these results bolsters interpretation of analyses with modeled results, despite the inherent error associated with the relationship between DN and observed canopy cover.Moreover, given that the lowest VI values occurred in pixels with predicted canopy cover around 60%, and that the model was biased towards under-prediction at high canopy cover, these low VI values may actually represent areas with higher canopy cover.While PlanetScope VI values were consistently lower than Landsat, potentially due to wider spectral bands relative to Landsat [36] or differences in imaging geometry [37], geographic patterns and relationships between canopy cover and VIs were consistent across the two platforms.
Phenological differences between larch and the dominant understory shrubs were not clearly reflected in field or satellite VIs.If the phenological differences we observed contributed strongly to variation in VIs then areas with greater canopy cover should have had higher VI values during the start and end of the growing season because overstory trees green-up earlier and senesce later than dominant understory shrubs.However, mean growing season VIs were still lower in denser forests despite the fact that canopies there would have been green for longer than in low-density forests where understory plants contribute more to stand-wide measures of VIs.However, mean growing season VIs were inversely related to canopy cover.Despite the relatively high temporal resolution of PlanetScope data, we were not able to examine phenology differences during green-up or senescence, and the spatial distribution of our field measurements of NDVI were limited.Previous analyses have shown that seasonal NDVI dynamics are dominated by a combination of snow-melt as well as overand under-story canopy phenology effects [38], and this is apparent in our 2017 field observations of NDVI.This may also be related to the fact that our phenology observations document changes in leaf coloration, but not senescence.Additionally, even 3 m pixels in all but the lowest and highest density stands likely contain a combination of overstory and understory reflectance.Further, Betula and Salix shrubs are only one component of understory vegetation communities in these ecosystems [27].The presence of evergreen and non-vascular plants in the understory may dampen the phenological signal from trees and deciduous shrubs.These results highlight the difficulties of interpreting phenological metrics in Siberian larch forests with heterogeneous mixtures of needleleaf, broadleaf, and non-vascular plants with both deciduous and evergreen leaf habits.
Our field observations show that understory NDVI in low-density stands may exceed overstory NDVI, as shown in Figure 7, though these results should be interpreted with caution due to potential effects of different combinations of illumination and viewing geometry [37].While the effect of understory vegetation has been recognized as an important factor influencing boreal forest reflectance [39], the relative invariance in VI across a large gradient in forest cover is unexpected.It may be the case that higher stem to leaf ratios in larch trees relative to understory vegetation components result in tree stems and branches having a strong negative effect on the VIs in high density forests.Lowland shrublands directly south across the river from NESS that were not included in our analyses had higher NDVI values than upland forested areas, as shown in Figure 9. Similarly, the dark dendritic features within the study area (i.e., high NDVI), as shown in Figure 9, are likely water tracks with dense tall shrub canopies.Broadleaf shrub canopies may mask the underlying woody plant components better than larch needles.In addition, other observations have shown a high degree of variability in understory reflectance [40] suggesting that variability in understory vegetation communities may contribute to spatial patterns in NDVI independently of overstory composition.Moreover, moss NDVI varies over short timescales with moisture conditions [41], which may help to explain daily and seasonal variability in understory NDVI observations, as shown in Figure 7B,D.

Implications for Interpreting NDVI Trends
Northeastern Siberia exhibits consistent greening trends across analyses of long-term NDVI data sets [15,42].The magnitude of observed NDVI trends is relatively small in comparison to the range of NDVI values that we observed across a fairly large range of forest canopy cover and biomass values.Our observations are representative of forest cover across a much larger area [43] where 75% of the study regions had less than 1900 g m −2 in larch biomass.Seasonally integrated GIMMS NDVI is positively related to canopy cover across the 100,000 km 2 Kolyma River watershed [9], however, this likely reflects geographic patterns in growing season length.On the other hand, tree ring indicators of larch growth were positively related to July NDVI [9], consistent with other analyses of larch growth and NDVI [44,45].This may indicate that larch and understory vegetation are responding similarly to interannual climate variability, and that NDVI captures changes in ecosystem productivity in the absence of changes in forest cover.
While relationships between larch radial growth and NDVI trends indicate links between these two variables [44,45], our results illustrate the potential for relatively large increases in forest cover with no discernable NDVI response.Because larch density across this region is low [46], it is likely that NDVI trends are the result of combined understory and overstory responses to climate, or even dominated by understory responses.This highlights the need to separate overstory and understory spectral dynamics [47][48][49] in boreal larch forests.It is even possible that browning trends could result from increases in forest cover in these ecosystems.Indeed, the highest biomass stands in our study area are all approximately 50 years old, having originated from a stand-replacing fire that occurred approximately 80 years ago.The NDVI signal associated with this high biomass regrowth was most likely influenced by the vegetation communities preceding forest growth after fire.Our sites may also serve as useful treeline analogs because many of the understory species at our sites are common to tundra ecosystems as well.As such, our results indicate that ecotonal changes such as treeline advance may not be detected by NDVI trends.

Conclusions
We did not observe variation in NDVI or EVI across an upland boreal forest landscape in northeastern Siberia with forest canopy cover ranging from 0% to over 90%.Temporal differences in overstory and understory canopy phenology were not captured by satellite data, likely due to a lack of temporal resolution.Our results highlight the importance of understory vegetation in boreal VI dynamics, particularly in relation to interpretation of long-term trends.In these ecosystems, which dominate a large portion of boreal Siberia, invariant or declining NDVI and EVI with increasing forest cover highlights the potential for large changes in forest cover with no change in VIs.This is not unique to this region [50], but nonetheless poses a challenge for interpreting the potential climate feedback consequences associated with vegetation productivity trends inferred from NDVI in Siberian larch forests.Thus, long-term trends in satellite derived VIs must be interpreted with caution.Our results highlight the utility of high-resolution satellite data for understanding how VIs vary with ecosystem composition, and the continued need to investigate the contribution of understory vegetation to ecosystem VI dynamics.
Remote Sens. 2018, 10, x FOR PEER REVIEW 3 of 15 spp., and Polytrichum spp.), and lichens (Cetraria spp.and Cladonia spp.).Forests of varying density are dominant in uplands to the east of the Kolyma River.Lowland floodplain ecosystems in the area contain forests but also large extensive shrublands dominated by Betula nana exilis, Salix spp., and Alnus fruticosa (Duschekia fruticosa), and also graminoid-dominated areas dominated by Eriophorum polystachion, Calamagrostis purpurea, and Carex spp.[27].For this study, we examined a ~25 km 2 area located to the east of NESS and the town of Cherskiy, as shown in Figure 1.

Figure 1 .
Figure 1.Map of the study area in northeastern Siberia.The blue polygon denotes the upland larch forested area used for analysis in this study.The red star denotes the location of the Northeast Science Station.Colored points indicate forest stand sample locations and associated canopy cover values.The background image is a panchromatic WorldView2 satellite image from 8 March 2015.

Figure 1 .
Figure 1.Map of the study area in northeastern Siberia.The blue polygon denotes the upland larch forested area used for analysis in this study.The red star denotes the location of the Northeast Science Station.Colored points indicate forest stand sample locations and associated canopy cover values.The background image is a panchromatic WorldView2 satellite image from 8 March 2015.

Figure 2 .
Figure 2. Aerial photographs of the high-density (A) and low-density (B) stands instrumented for in situ Normalized Difference Vegetation Index (NDVI).Photographs were both taken during leaf-out on 1 June 2017.Tower masts are located in the upper right corner of the photographs.Photographs: M. Loranty.

Figure 2 .
Figure 2. Aerial photographs of the high-density (A) and low-density (B) stands instrumented for in situ Normalized Difference Vegetation Index (NDVI).Photographs were both taken during leaf-out on 1 June 2017.Tower masts are located in the upper right corner of the photographs.Photographs: M. Loranty.

Figure 3 .
Figure 3. Linear relationship between canopy cover and aboveground larch biomass carbon for 46 stands sampled across an upland landscape.Whiskers indicate standard error of the mean (SE).

Figure 4 .
Figure 4. Inverse non-linear relationship between stand-level canopy cover and digital number (DN) from a March 2015 panchromatic WorldView2 image.Note the model was fit using standardized digital numbers (DNS) and data were back-transformed for plotting.Vertical whiskers indicate standard error (SE) of mean canopy cover, and horizontal whiskers are one standard deviation (SD) of mean stand DN.

Figure 3 .
Figure 3. Linear relationship between canopy cover and aboveground larch biomass carbon for 46 stands sampled across an upland landscape.Whiskers indicate standard error of the mean (SE).

Figure 3 .
Figure 3. Linear relationship between canopy cover and aboveground larch biomass carbon for 46 stands sampled across an upland landscape.Whiskers indicate standard error of the mean (SE).

Figure 4 .
Figure 4. Inverse non-linear relationship between stand-level canopy cover and digital number (DN) from a March 2015 panchromatic WorldView2 image.Note the model was fit using standardized digital numbers (DNS) and data were back-transformed for plotting.Vertical whiskers indicate standard error (SE) of mean canopy cover, and horizontal whiskers are one standard deviation (SD) of mean stand DN.

Figure 4 .
Figure 4. Inverse non-linear relationship between stand-level canopy cover and digital number (DN) from a March 2015 panchromatic WorldView2 image.Note the model was fit using standardized digital numbers (DN S ) and data were back-transformed for plotting.Vertical whiskers indicate standard error (SE) of mean canopy cover, and horizontal whiskers are one standard deviation (SD) of mean stand DN.

Figure 5 .
Figure 5. Map of canopy cover across the study area predicted using the relationship between DN and field observations shown in Figure 4.

Figure 5 . 15 Figure 5 .
Figure 5. Map of canopy cover across the study area predicted using the relationship between DN and field observations shown in Figure 4.

Figure 6 .
Figure 6.Relationship between observed and predicted canopy cover (A), and between observed and residual canopy cover (B).Dashed lines show regression fits, thin solid lines indicate 1:1 lines (A), and abscissa (B).

Figure 6 .
Figure 6.Relationship between observed and predicted canopy cover (A), and between observed and residual canopy cover (B).Dashed lines show regression fits, thin solid lines indicate 1:1 lines (A), and abscissa (B).

Figure 7 .
Figure 7. Field observations of NDVI from overstory (A) and understory (B) sensors for 2016, and from overstory (C) and understory (D) sensors in 2017.Blue and black symbols represent measurements from the low-and high-density stands, respectively.Open and closed circles indicate replicate (east-and west-facing) sensors at each stand.Points represent daily observations and lines are a 5-day running mean.

Figure 7 .
Figure 7. Field observations of NDVI from overstory (A) and understory (B) sensors for 2016, and from overstory (C) and understory (D) sensors in 2017.Blue and black symbols represent measurements from the low-and high-density stands, respectively.Open and closed circles indicate replicate (east-and west-facing) sensors at each stand.Points represent daily observations and lines are a 5-day running mean.

Figure 8 .
Figure 8. Relationships between observed canopy cover and PlanetScope NDVI (A) and EVI (B) across all 46 forest stands for five days during the growing season, and also for both indices averaged across the growing season (C).Plot whiskers show the standard error of mean observed canopy cover, and the standard deviation of NDVI/EVI for all pixels with centers located within 15 m of the stand center.Plots (D-F) show the same, but with modeled canopy cover aggregated into bins of 10%.Shaded polygons represent the standard deviation of NDVI/EVI within each canopy cover bin.In all plots, trend lines show significant (p < 0.05) linear relationships between canopy cover and NDVI/EVI.Absence of trendlines indicates no significant relationship.

Figure 8 .
Figure 8. Relationships between observed canopy cover and PlanetScope NDVI (A) and EVI (B) across all 46 forest stands for five days during the growing season, and also for both indices averaged across the growing season (C).Plot whiskers show the standard error of mean observed canopy cover, and the standard deviation of NDVI/EVI for all pixels with centers located within 15 m of the stand center.Plots (D-F) show the same, but with modeled canopy cover aggregated into bins of 10%.Shaded polygons represent the standard deviation of NDVI/EVI within each canopy cover bin.In all plots, trend lines show significant (p < 0.05) linear relationships between canopy cover and NDVI/EVI.Absence of trendlines indicates no significant relationship.

Figure 9 .
Figure 9. NDVI maps across the study area from PlanetScope on 26 July 2017 (A) and Landsat on 28 July 2017 (B).Note the maps use the same scale, indicating that Landsat has higher NDVI values than PlanetScope.

Figure 10 .
Figure 10.Relationships between Landsat Vegetation Indices (VI) from 28 July 2017 image and observed canopy cover.Dashed lines show statistically significant relationship (p < 0.05).Horizontal whiskers indicate standard error (SE) of mean canopy cover.

Figure 9 .
Figure 9. NDVI maps across the study area from PlanetScope on 26 July 2017 (A) and Landsat on 28 July 2017 (B).Note the maps use the same scale, indicating that Landsat has higher NDVI values than PlanetScope.

Figure 9 .
Figure 9. NDVI maps across the study area from PlanetScope on 26 July 2017 (A) and Landsat on 28 July 2017 (B).Note the maps use the same scale, indicating that Landsat has higher NDVI values than PlanetScope.

Figure 10 .
Figure 10.Relationships between Landsat Vegetation Indices (VI) from 28 July 2017 image and observed canopy cover.Dashed lines show statistically significant relationship (p < 0.05).Horizontal whiskers indicate standard error (SE) of mean canopy cover.

Table 1 .
Phenology differences between overstory larch trees and understory Betula and Salix shrubs for the 2000-2017 period.Values in parentheses represent one standard deviation.
a significantly different at p < 0.01; b significantly different at p < 0.01.SOS: Start of season; EOS: End of season; GSL: Growing season length.

Table 1 .
Phenology differences between overstory larch trees and understory Betula and Salix shrubs for the 2000-2017 period.Values in parentheses represent one standard deviation.
a significantly different at p < 0.01; b significantly different at p < 0.01.SOS: Start of season; EOS: End of season; GSL: Growing season length.
Relationships between Landsat Vegetation Indices (VI) from 28 July 2017 image and observed canopy cover.Dashed lines show statistically significant relationship (p < 0.05).Horizontal whiskers indicate standard error (SE) of mean canopy cover.