Do recent NDVI trends demonstrate boreal forest decline in Alaska?

Remote sensing analyses of boreal forest regions have found widespread decreasing or increasing trends in normalized difference vegetation index (NDVI). Initially, these trends were attributed to climate change induced shifts in primary productivity. It is emerging, however, that fire disturbance and subsequent succession also strongly impact the optical properties of boreal forests. Here we use NDVI time series data from Landsat (1999–2018) paired with surveys of 102 forest stands with known recent fire history to investigate the relationship between NDVI and forest structure during succession. We found that NDVI varies systematically with stand age as a result of successional changes in forest structure and composition and that the proportion of deciduous (broad-leaved) trees in the upper canopy is a better predictor of NDVI than leaf area index. Recent fire disturbance led to strong NDVI decreases and early post-fire recovery of herbaceous and deciduous vegetation to strong NDVI increases. The mid-succession transition from deciduous to evergreen (needle-leaved) stands led to weak NDVI decreases, while mid-to-late succession thinning of evergreen canopies led to weak NDVI increases. Thus, both increasing and decreasing NDVI stands occur naturally across the landscape, and do not necessarily reflect a large-scale shift in boreal forest productivity.


Introduction
Boreal forests account for nearly one third of global forest area and store~30% of the total terrestrial ecosystem carbon pool, but climate change, increased fire, and human disturbance may be altering forest dynamics and productivity (Kasischke et al 1995, 2010, Pan et al 2011, Ma et al 2012, Bradshaw and Warkentin 2015. Remote sensing provides a valuable tool for studying the composition and productivity of this large and often inaccessible region, and many studies have utilized satellite-derived measurements of normalized difference vegetation index (NDVI) to infer landscape-scale productivity. Time series are used to calculate NDVI trends over time, with NDVI decreases interpreted as 'browning' and NDVI increases interpreted as 'greening' (Goetz et al 2005). Previous studies have reported various regional trends, with large spatial and intensity differences between datasets (e.g. GIMMS vs Landsat) and methodologies (Alcaraz-Segura et al 2010, de Jong et al 2011, Forkel et al 2013, Guay et al 2014. Previous studies comparing NDVI trends across datasets have found conflicting results, with some showing NDVI increases and others decreases for the same location (Alcaraz-Segura et al 2010). One source of discrepancy in NDVI analyses lies in the resolution of the dataset, which ranges from 8 km resolution for the Global Inventory Monitoring and Modeling Studies (GIMMS) Advanced Very High Resolution Radiometer (AVHRR) data to 30 m for Landsat data. Recent work has emphasized finer resolution datasets like Landsat to better diagnose the drivers of changing NDVI (Ju andMasek 2016, Sulla-Menashe et al 2018).
A second issue involves the interpretation of NDVI itself. Many studies have tied NDVI to field measurements of biophysical properties of vegetation and found positive correlations with vegetation cover, leaf area index (LAI), biomass, vegetation health, and chlorophyll concentration (Carlson and Ripley 1997, McMillan and Goulden 2008, Yang et al 2017. However, these biophysical properties are not fully independent, which challenges efforts to interpret recent NDVI trends in terms of changes in actual ecological or biophysical properties. Furthermore, relationships between biophysical properties and NDVI may not hold throughout the whole life history of a forest stand. For example, NDVI increases linearly with LAI for stands with low LAI before saturating at high LAI (Carlson et al 1990, Turner et al 1999. These ambiguities underscore the need for systematic field observations to better interpret observed changes in remotely sensed vegetation indices. We made field observations of vegetation structure and composition at over 100 forest stands of varying age in interior Alaska to investigate the relationship between recent NDVI trends and forest structure and composition. Locations were selected based on fire history and the 1999-2018 Landsat NDVI trend to sample four groups: (1) Stands that burned after 2005; (2) Stands that burned during 1980 and 2005; (3) Stands with no burn reported since 1980 and a decreasing NDVI; and (4) Stands with no burn reported since 1980 and an increasing NDVI. This collection of stands allowed us to explore four questions: (1) How are changes in forest structure and composition in recently burned stands associated with NDVI decreases? (2) How are changes in forest structure and composition in young stands associated with NDVI increases? (3) In the absence of recent fire, how are NDVI decreases in mid-tolate successional forests related to forest structure and composition? (4) In the absence of recent fire, how are NDVI increases in mid-to-late successional forests related to forest structure and composition?

Fire history
We used the 1940-2018 Alaskan fire perimeter database produced by the Alaska Fire Service (AFS) and acquired from the Alaska Interagency Coordination Center (https://fire.ak.blm.gov/incinfo/aklgfire.php) to identify sample locations. We subsequently determined the actual age of each stand with tree cores, so errors in the fire perimeter database and the lack of fire history information before 1940 would not affect our analysis.

NDVI extraction and temporal trend analysis
June-August Collection 1, Level 1 Landsat 5, 7, and 8 Surface Reflectance and Brightness Temperature images with less than 30% cloud cover for WRS2 Paths 66-72 and rows 14-16 were downloaded from the USGS (https://espa.cr.usgs.gov; 658 images total; 604 images after 1998) after reprojection at 0.0002695 • resolution (30 m latitude, variable longitude). The images were masked for snow, cloud and cloud shadow, and NDVI was calculated, homogenized across instruments and the summer median calculated (see details in Goulden and Bales 2019). Images were not masked for standing water. There was some residual effect of the Landsat 7 Scan Line Corrector (SLC) failure despite masking, presumably due to systematic alignment of gaps and resulting differences in scene availability. Subsequent analyses focused on 1999-2018 due to large gaps in the earlier record with the lack of a consistent downlink station in Alaska. Simple linear regressions were calculated from the resulting NDVI summer medians over the 20 years time period. Trends with a significant nonzero slope were calculated using a two-tailed t-test, at p ≤ 0.1 (slopes were considered significant for α ≤0.05 or ≥0.95).

Field sites
We surveyed 102 forest stands across interior Alaska, spanning a distance of over 425 km (figure 1) in August 2017 and 2018. The sites were selected before visiting the field to include locations with and without a fire since 1940. Recently burned sites were selected to span a range of years since fire, while sites without a recent fire were selected to include a range of Landsat NDVI trends.
Sites were further screened to avoid large topographic and land use gradients and to allow efficient access. Two perpendicular 60 m transects were established at each site, and stand characteristics were measured at 10 m intervals along each transect. Measurements included canopy cover, percent species composition of upper and lower canopies (measured as proportional cover of species over total area on interval), average height of upper and lower canopies, ground cover, standing and downed deciduous stems, and standing and downed evergreen stems. Any evidence of disturbance (fire, permafrost thaw, beetle damage, logging, or other human impacts) was noted. Fraction deciduous was subsequently calculated from the proportion of deciduous cover in the upper canopy over the total canopy coverage.
LAI (leaf area index) data were collected at each site using LAI-2000 plant canopy analyzer (LI-COR). Eight LAI-2000 measurements were taken under the canopy in each 30 m interval along the 60 m perpendicular transects and corresponding above canopy measurements were taken in the nearest canopy clearing. The measurements were then processed with the internal LI-COR LAI software to provide four estimations of LAI per site. The LAI-2000 is known to underestimate LAI in coniferous stands, so species specific corrections for black spruce (Picea mariana (Mill.) B.S.P.) were applied to black spruce site data (Chen 1996, Law et al 2008. Corrections for white spruce (Picea glauca (Moench) Voss) are not available, and these stands were frequently mixed with deciduous species; we consequently were unable to correct for clumping in the white spruce stands and note that LAI may be underestimated at these sites. The mean LAI was then calculated from the four estimations and used as the site level LAI.
To estimate tree age of each site, multiple tree cores were collected from representative trees of each species in the upper canopy (2-4 trees per species) with 4.3 mm outer diameter increment borers (Haglöf). Tree cores were mounted, sanded manually, and rings were counted using Lignovision software (Rinntech). The average tree age was subsequently determined by taking the mean ring count of cores analyzed. For recently burned and early fire recovery sites the site age was determined from the year of fire and confirmed by tree cores when possible.
Each successional stage had a distinct forest structure, and a characteristic NDVI trend (figure 2 and table S1 (available online at (stacks.iop.org/ERL/15/095007/mmedia))). Recently burned stands (age 0-13 years) had charred, dead trees with occasional live, evergreen trees in the upper canopy that survived the fire, and deciduous saplings in the lower canopy (figure 3). Herbaceous plants, moss, and charred material covered most of the ground. Early fire recovery stands (age 13-38 years) had upper canopies of evergreen trees that had survived fire, along with deciduous species such as paper birch (Betula papyrifera Marshall) and quaking aspen (Populus tremuloides Michx.), and lower canopies of birch and aspen, and mostly herbaceous ground cover (figure 3). Standing and fallen dead trees from the fire were still visible. Mid-to-late successional deciduous stands (age 30-180 years) had dense upper canopies of birch and aspen, and lower canopies of birch, aspen, and alder (Alnus sps.), with progressively more white spruce or black spruce in older stands (figure 3). The ground cover of deciduous stands consisted of bare ground and litter, herbaceous plants and small percentages of moss. Mid-to-late successional white spruce stands (age 50-200 years) often had multiple generations of white spruce in the upper canopy, with varying percentages of deciduous species in the upper canopy, and lower canopies of white spruce with dense thickets of alder (figure 3). Midto-late successional black spruce (age 30-200 years) had sparse black spruce in both the upper and lower canopies, and occasionally small percentages of birch and willow (Salix sp.) in the upper and lower canopies. Thick layers of moss, lichen, and herbaceous plants covered most of the ground in mid-to-late successional black spruce stands (figure 3).

NDVI trends across the landscape
There was a visible mosaic of increasing and decreasing NDVI trends across interior Alaska, which was Sample size (number of sites visited) is shown above each successional stage label, and p-value calculated using a two-tailed, one sample t-test on all trends for each group is shown above each boxplot. Null hypothesis is that the mean trend for each group is '0' . clearly linked with recent fire history (figures 4 and S1). Burn scars from recent decades show consistent behavior in NDVI trends, with trend perimeters typically matching fire perimeters. Timing of fire largely controlled NDVI trend direction and intensity, with recent burns  showing strong NDVI decreases and burns from  showing NDVI increases (figures 4, S1 and S2). Trend significance inside burn scar perimeters was variable, with some pixels showing significant trends (p-value < 0.1), and others insignificant within the same burn scar ( figure 4(B)).
There was heterogeneity in NDVI trends in the absence of recent fires, with some undisturbed patches showing increasing NDVI and others showing decreasing NDVI. These undisturbed patches often showed smaller magnitude trends than areas with reported fires, though 17.2% of these trends were significant (p-value < 0.1) (figure 4(B) and S2). Moreover, the unburned locations with significant trends were clustered in relatively homogenous patches (figures 4(B) and S1(b)), and were not intermixed with pixels of opposite slope as would be expected with a random distribution, implying that NDVI trends in the absence of recent fire may reflect ecological variation.

Successional effects on NDVI
Our data allowed us to create a 200 years long chronosequence of the boreal forest NDVI successional trajectory (figure 5), and to calculate the resulting 20 years chronosequence of NDVI trends related to recent fires (since 1940) (figure 6). Stands with recent fire displayed large ranges in NDVI due to the rapid changes in forest structure and composition following fire, but overall had the lowest average NDVIs ( figure 5(A)). Fire disturbance caused a sharp drop in NDVI due to the loss of vegetation, followed by rapid recovery of herbaceous and deciduous vegetation leading to strong NDVI increases (figures 3 and 5(B)). NDVI continued to increase until around 30 years since fire ( figure 5(A)), at which point many stands displayed maximum fraction deciduous in the upper canopy ( figure 3). The immediate effect of fire and subsequent recovery were large drivers of NDVI change in stands <30 years old (figures 5(B) and 6).
After stands reached 30 years of age, NDVI typically stabilized and displayed more gradual changes. Mid-to-late successional deciduous stands collectively had the highest summer median NDVI, followed by white spruce, then black spruce stands ( figure 5(A)). With increasing age, deciduous stands exhibited weak NDVI decreases and black spruce stands of all ages showed weak NDVI increases (figures 2 and 5(B)). White spruce stands displayed both decreasing and increasing trends, with younger white spruce stands typically showing decreases, and older white spruce stands showing increases (figures 2 and 5(B)).

Ecological effects on NDVI
The relationship between successional stage and NDVI can be explored further by examining the relationships between fraction deciduous in the upper canopy and NDVI, and LAI and NDVI. We found that NDVI was positively correlated with fraction deciduous in the upper canopy (r 2 = 0.63) and with LAI (r 2 = 0.15) (figure 7). Deciduous stands by definition had the highest fraction deciduous and highest NDVI. These stands showed marked variation in LAI, which was correlated with differences in NDVI, especially in sparsely vegetated, recent burns. Compared to deciduous stands, white spruce stands had a smaller fraction deciduous, a smaller range in LAI, and lower NDVI. Here, upper canopies typically consisted of a mixture of white spruce and deciduous species and an increasing fraction deciduous generally increased NDVI (figures 3 and 7(A)). Black spruce stands had the lowest NDVI, a relatively small range in LAI, and were less likely to have deciduous species in the upper canopy. Overall, fraction deciduous in the upper canopy was a stronger predictor of NDVI than LAI.

Discussion
The central concept addressed in this study is how changes in the biophysical properties of forests following fire and subsequent succession impact NDVI. Our results indicate that fire was the main driver of large changes in NDVI in our sites. Fires after 2005 were associated with large decreases in NDVI due to vegetation loss from fire, and fires between 1980 and 2005 with large increases in NDVI due to post fire recovery of deciduous vegetation (figures 2, 4 and 6) (see also Randerson et al 2006, Alcaraz-Segura et al 2010, Ju and Masek 2016. The NDVI increases from fire recovery last for roughly 30 years post-fire, at which point many stands displayed maximum fraction deciduous (figures 3, 5(A) and 7). This finding  is supported by (Beck and Goetz 2011) who also found that high NDVI coincided with high deciduous coverage around 20-40 years post-fire, and by (Amiro et al 2006) who found high albedo around 30 years post-fire that they linked to high deciduous coverage.
After this peak in NDVI and deciduous cover, mid-succession changes in forest structure and composition led to weak decreases in NDVI in deciduous stands (figures 2 and 5(B)). We believe the ongoing replacement of high-NDVI deciduous, broad-leaved trees with lower-NDVI evergreen, needle-leaved trees explains this trend, as evidenced by the changes seen in forest composition (figure 3) and the relationship between fraction deciduous and NDVI ( figure 7(A)). Our results confirm that NDVI saturates at high LAI (figure 7(B)) (Carlson et al 1990, Gong et al 1995, Turner et al 1999, and that fraction deciduous in the upper canopy is a better predictor of NDVI than LAI in boreal forests. Our data also indicate a mid-to-late succession NDVI increase, especially evident in black spruce stands and older white spruce stands (figures 2 and 5(B)). We believe this increasing trend in NDVI is due to the natural self-thinning and, or paludification of evergreen forests over time that reveals understories with higher NDVI, consisting of deciduous shrubs, grasses, and/or mosses (figure 3). These results are consistent with previous work showing that NDVI trends are related to cover type and that forests with sparse tree cover often show increasing NDVI trends Goetz 2006, Miles andEsau 2016). As such, our work further supports the conclusion that the understory of sparse forests plays a large role in remotely sensed NDVI and that increasing NDVI is a poor proxy of tree productivity (Caetano and Pereira 1996).
Many previous studies have documented widespread NDVI decreases across boreal forests, which were hypothesized to reflect climate change-induced drought stress (Goetz et al 2005, Bunn and Goetz 2006, Parent and Verbyla 2010 and later linked to changes in productivity via tree ring width analyses (Beck et al 2011, Bunn et al 2013. However, our study demonstrates that changes in NDVI also signal changes in forest species composition and structure that are not directly related to productivity and/or vegetation health. While the use of satellite-derived vegetation indices such as NDVI to analyze vegetation dynamics remotely has clear advantages, one evident drawback lies in the difficulty of directly ascribing changes in indices to the biophysical or ecological properties of forests and their understories on the ground.
One limitation of our study is that our geographic location did not cross any ecotones where climate change impacts are more likely to affect productivity (Goetz et al 2005, Andreu-Hayles et al 2011, Sulla-Menashe et al 2018. However, Landsat's 30 m resolution allowed us to better link NDVI to in situ forest properties. Our results indicate that NDVI varies systematically with stand age and successional stage, and that both fire and natural succession lead to composition and structure changes that cause changes in NDVI, and do not necessarily reflect climate change induced changes in tree productivity.

Conclusion
Our study linking surveys of forests with remotelysensed NDVI demonstrates how boreal forest fire and succession impact NDVI, and that fraction deciduous in the upper canopy is a better predictor of NDVI than LAI in boreal forests. We conclude that the immediate and legacy effects of fire are strong drivers of changes in boreal NDVI, and that forest patches with both increasing and decreasing NDVI occur naturally across the landscape, and do not unambiguously indicate large-scale changes in boreal tree productivity. While the response of forests to climate change remains unclear, our results underscore the need to deconvolve the legacy effect of disturbance and recovery when diagnosing recent trends in forest structure and function.