The role of land cover change in Arctic-Boreal greening and browning trends

Many studies have used time series of satellite-derived vegetation indices to identify so-called greening and browning trends across the northern high-latitudes and to suggest that the productivity of Arctic-Boreal ecosystems is changing in response to climate forcing at local and continental scales. However, disturbances that alter land cover are prevalent in Arctic-Boreal ecosystems, and changes in Arctic-Boreal land cover, which complicate interpretation of trends in vegetation indices, have mostly been ignored in previous studies. Here we use a new land cover change dataset derived from Landsat imagery to explore the extent to which land cover and land cover change influence trends in the normalized difference vegetation index (NDVI) over a large (3.76 M km2) area of NASA’s Arctic Boreal Vulnerability Experiment, which spans much of northwestern Canada and Alaska. Between 1984 and 2012, 21.2% of the study domain experienced land cover change and 42.7% had significant NDVI trends. Land cover change occurred in 27.6% of locations with significant NDVI trends during this period and resulted in greening and browning rates 48%–128% higher than in areas of stable land cover. While the majority of land cover change areas experienced significant NDVI trends, more than half of areas with stable land cover did not. Further, the extent and magnitude of browning and greening trends varied substantially as a function of land cover class and land cover change type. Forest disturbance from fire and timber harvest drove over one third of statistically significant NDVI trends and created complex mosaics of recent forest loss (as browning) and post-disturbance recovery (as greening) at both landscape and continental scale. Our results demonstrate the importance of land cover changes in highly disturbed high-latitude ecosystems for interpreting trends of NDVI and productivity across multiple spatial scales.


Introduction
Climate change is warming Boreal and Arctic ecosystems twice as rapidly as the global mean warming (Pithan and Mauritsen 2014) and causing changes in disturbance regimes and ecosystem function in northern high-latitudes (Keeling et al 1996, Kasischke and Turetsky 2006, Kasischke et al 2010, Graven et al 2013. Arctic-Boreal 'greening' and 'browning' trends, commonly inferred via time series of satellite-derived normalized difference vegetation index (NDVI) data, have been widely interpreted to indicate changes in northern high-latitude ecosystem function. Specifically, greening (positive NDVI trends) has been broadly interpreted as reflecting increased ecosystem productivity or vegetation cover, while browning (negative NDVI trends) has been interpreted as reduced productivity or vegetation cover due to insect infestations, drought, and other sources of stress (Zhang et al 2008, Parent and Verbyla 2010, Beck and Goetz 2011, Verbyla 2011, Rogers et al 2018.
However, recent studies have cast doubt on such interpretations of NDVI trends (Ju and Masek 2016, Sulla-Menashe et al 2018). Specifically, NDVI trends can reflect numerous other changes in land surface properties, and NDVI has low sensitivity to changes in land surface properties at high values (e.g. NDVI≈0.8) (Myneni and Williams 1994, Carlson and Ripley 1997, Buermann et al 2002, Pastick et al 2019. Further, NDVI trends are sensitive to fire disturbance (Sulla-Menashe et al 2018) and changes in surface moisture (Raynolds and Walker 2016), and NDVI imagery from sensors with coarse (i.e. >500 m) spatial resolution, such as the commonly-used Advanced Very High Resolution Radiometer (AVHRR), do not resolve landscape-scale disturbances (Guay et al 2014, Ju andMasek 2016). Challenges with AVHRR data also include low radiometric quality, orbital drift, and sensor cross-calibration issues ( Despite the prevalence of disturbance and land cover change in high-latitude ecosystems (White et al 2017, Wang et al 2019b, long-term land cover datasets with sufficient spatial resolution are lacking, and most studies of high-latitude NDVI trends have therefore not considered how land cover changes influence observed trends (e.g. Park et al 2016). More generally, land cover is a fundamental ecosystem attribute that is widely used to parameterize vegetation properties in data products and ecosystem models (Myneni et al 1997, Running et al 2004, Zhang et al 2008, Luus and Lin 2015, Zhu et al 2017. Therefore, improved understanding of land cover dynamics is required to accurately interpret the satellite record and characterize high-latitude ecosystem change.
In this paper, we quantify the extent to which land cover dynamics are associated with high-latitude greening and browning trends derived from time series of 30 m Landsat imagery, which, unlike imagery from coarse spatial resolution sensors, can resolve landscape-scale patterns in disturbance and regrowth.
To do this, we analyze joint variation in multidecadal land cover change and NDVI trends from Landsatbased datasets spanning the Core Study Domain of NASA's Arctic-Boreal Vulnerability Experiment (ABoVE), which encompasses most of Alaska and northwestern Canada. Specifically, we address the following questions: (1) To what extent are NDVI trends related to land cover change in the ABoVE Core Study Domain?
(2) How do different land cover classes and land cover change types control the area, sign, and magnitude of greening and browning trends?
( Land cover and land cover change maps The land cover dataset generated by Wang et al (2019a) provides maps of annual land cover created using the continuous change detection and classification (CCDC) algorithm (Zhu and Woodcock 2014), which identifies land cover change as statistically significant breaks in time series of Landsat surface reflectance data at each pixel. Land cover class labels were assigned to time segments between breaks using Random Forest classification (Breiman 2001 The land cover dataset maps ten land cover classes. For this work, we focused on six classes that represent dominant vegetation functional types: Evergreen Forests, Deciduous Forests, Shrubs, Herbaceous, Sparse Vegetation, and Barren land (table 1). The remaining classes were aggregated into an 'Other' class, which we do not analyze. For land cover change, we defined six major types based on the land cover change types analyzed in (Wang et al 2019b): Evergreen Forest Loss, Shrub Loss, Deciduous Forest Gain, Evergreen Forest Gain, Herbaceous Gain (from Sparse Vegetation or Barren), and Shrub Gain (from Herbaceous, Sparse Vegetation, or Barren). Pixels exhibiting cyclic land cover changes (e.g. from Evergreen Forest to Shrub and back to Evergreen Forest) occurred at negligible rates but are still considered to have changed. To avoid double counting areas with conflicting transitions, the change types are explicitly designed not to allow overlapping change types (e.g. Evergreen Forest becoming Shrub is Evergreen Forest Loss and not Shrub Gain). Each pixel is thus assigned a single stable land cover class or land cover change type, and approximate areas of loss do not sum to areas of gain because land cover is not in steady state in the domain (Wang et al 2019b). Annual land cover class maps had an overall accuracy of 84.1%, specific land cover change types had an overall accuracy of 76.1%, and land cover status (i.e. changing versus stable) had an accuracy of 86.6%. The land cover change dataset is described in detail in Wang et al (2019b).

Greening and browning data
To evaluate spatial patterns in greening and browning, we used the NDVI trend dataset created by Ju and Masek (2018), which was derived from time series of Landsat 5 and 7 surface reflectance imagery for the period 1984-2012 and spans most of Canada and Alaska. This dataset provides estimates of NDVI trends at each pixel by fitting an ordinary least squares linear model to peak summer (July and August) NDVI as a function of time and includes the magnitude (slope) and statistical significance (p-value) of the estimated model. We consider statistically significant trends at α=0.05.
Our analysis excludes areas that are mapped as water and small sub-areas of the ABoVE Core Study Domain that were not included in the Ju and Masek dataset. The resulting study area covered 3.86×10 6 km 2 . The Ju and Masek (2018) and Wang et al (2019a) datasets cover slightly different time periods (1984-2014 and 1984-2012, respectively). However, only 1776 km 2 (0.046% of the domain) experienced land cover change in 2013 and 2014, and we therefore assumed this difference in time periods had negligible impact on our analyses.

Disturbance maps
To characterize relationships between disturbance and NDVI trends, we used existing datasets that catalogue the location and timing of fires and timber harvest throughout the study domain. For fires, we combined the Alaskan Large Fire Database (Kasischke et al 2002) and the Canadian Fire Database (Stocks et al 2002) to generate a database of fire perimeters for the ABoVE Core Study Domain. For timber harvest, we used Landsat-derived datasets for the timing and location of harvest in Canada that span 1985-2010 (White et al 2017). Timber harvest did not occur at significant rates in the Alaskan portion of the study domain (note that the ABoVE Core Study Domain excludes the southern portion of Alaska) (Smith 2002).

Analysis
Our analysis includes four elements. First, we intersected the land cover dataset with the NDVI trend dataset to quantify the extent to which greening and browning is co-located with land cover change, stable land cover classes and specific land cover change types. Second, we sampled pixels with significant trends from across the study domain, equivalent to 1.5% of the study domain, and analyzed how the magnitudes of NDVI trends were distributed within and across each stable land cover class and land cover change type. The final sample included 2383 873 pixels (2145 km 2 ; table 2). Third, to characterize the spatial heterogeneity of land cover change and NDVI trends at landscape scale, we selected two sub-regions with complex landscape composition caused by wildfire or timber harvest. Specifically, we selected an area with extensive harvest in British Columbia and an area with extensive fires in the Northwest Territories, and analyzed spatial relationships between NDVI trends and land cover change in each. Finally, to better understand relationships between NDVI and land cover, we sampled (n=1750 000) from a publicly available Landsat-based dataset of peak growing season NDVI (July and August) that spans the ABoVE Core Study Domain from 1984 to 2014 (Melaas et al 2019). We sampled from stable land cover pixels and characterized the distribution of summer 2010 NDVI values associated with each land cover class. Table 1. Land cover classes and land cover change types considered for this study, adapted from Wang et al (2019b). Land cover change types are indicated by italics, in contrast to the stable land cover classes. Extents listed are estimates for the particular sub-domain (i.e. the intersection of the ABoVE domain and the Ju and Masek (2018) NDVI trend map), rather than for the entire ABoVE domain as described in Wang et al (2019b).

Land cover/change type Description
Extent ( 1(a)). Concurrently, just over one fifth of the domain (816 148 km 2 ; 21.2%) experienced land cover change between 1984 and 2014. A large proportion of significant NDVI trends were co-located with land cover change. Among areas with significant NDVI trends (p<0.05), over a quarter (454 630 km 2 ; 27.6% of significant trends) also experienced land cover change. The majority (55.7%) of areas of land cover change experienced significant NDVI trends, while only a minority (38.9%) of areas with stable land cover showed significant NDVI trends ( figure 1(b)). Greening and browning occurred both in areas of stable land cover and in areas of land cover change.
The magnitude of greening and browning varied significantly by land cover change status. Areas of stable land cover experienced modest NDVI trends (+0.0023±0.0015 year −1 for greening and −0.0021±0.0019 year −1 for browning; media-n±interquartile range). In contrast, areas of land cover change showed larger magnitude trends (+0.0034±0.0028 year −1 for greening and −0.0048±0.004 year −1 for browning) (figure 1(c)). On average, the rate of greening was 47.8% higher and the rate of browning was 128.6% higher in areas of land cover change compared to stable areas.
Regionally, both significant NDVI trends and land cover change were spatially extensive, with the Arctic biome experiencing predominantly greening and the Boreal biome experiencing greening, browning, and disturbance (figure 2). NDVI trends associated with land cover change had higher magnitude and were concentrated in localized patches, particularly in areas of Evergreen Forest Loss and fire disturbance (figures 2 and 3). In contrast, NDVI trends associated with stable land cover change were widespread, but of relatively lower magnitude (figure 3).

Greening and browning by land cover
We observed substantial variability in the area and magnitude of NDVI trends as a function of land cover class and land cover change type (figures 4 and 5). Except for Sparse Vegetation, the majority of pixels in each stable land cover class showed no significant NDVI trend (figure 4(c)). Significant NDVI trends in stable land cover tended to be positive, with notable ( i.e. more than 5%) browning only occurring in Deciduous Forests (11.7%), Evergreen Forests (7.5%), and Sparse Vegetation (6.9%) (figure 4(c)). In contrast, the majority of pixels in all but one land cover change type, Shrub Loss, showed significant trends ( figure 4(d)). Browning was only a large component in Evergreen Forest Loss, in which 43.7% of the areas showed browning and 20.6% showed greening. In other land cover change types, greening was considerably more prevalent than browning.
Trend magnitudes varied substantially by land cover class and land cover change type. Among stable land cover classes, Shrub, Herbaceous, and Barren land cover types showed the largest NDVI trends

Landscape-scale complexity in disturbed regions
The timing of disturbances strongly influences the sign and magnitude of NDVI trends, which can vary substantially at landscape scale. Timber harvest creates complex patterns, with localized patches of NDVI trends on the order of 100's of meters that vary in both sign and magnitude (figure 6). Browning is common where harvest occurred later in the time series (i.e. after 1995), and is generally mapped as Evergreen Forest Loss (figures 6(a)-(c)). Greening is common where harvest occurred early in the time series, which provides time to observe regrowth. In contrast, NDVI    trends in stable areas with no harvest or land cover change were dominated by low magnitude greening trends ( figure 6(d)).
Fires impact both land cover and NDVI trends throughout the Boreal part of the study domain ( figure 2). At landscape scale, fires create complex  mosaics of NDVI trends and land cover change, with both sign and magnitude of trends influenced by the timing of fire ( figure 7). Browning tends to occur where fires happened later in the time series (i.e. 1990 and later), while earlier fires were associated with greening, representing forest recovery (figures 7(a)-(c)). Burned areas with stable land cover or Evergreen Forest Gain showed substantial and high-magnitude greening, reflecting ongoing recovery from fires that occurred before the start of the Landsat time series ( figure 7(d)).
Together, fire and timber harvest occurred in over one third (34.3%) of significant NDVI trends. Fire events occurred in 5.5% of greening pixels and in 24.4% of browning pixels, while harvest occurred in 1.0% of greening pixels and 3.4% of browning pixels. The largest magnitude NDVI trends were generally found in these disturbed pixels. Pixels with land cover change not attributable to either fire or harvest account for 19% of positive trends and 24% of negative trends.

NDVI trends and land cover change
Our results suggest that spatial patterns in NDVI trends are strongly, but not exclusively, influenced by land cover change and disturbance processes at both regional and landscape scale (figures 2, 3, 6 and 7). Land cover change has a disproportionate effect on NDVI trends-over one quarter of pixels with significant NDVI trends also experienced land cover change, even though land cover change occurred in only one fifth of the domain. Further, the highest magnitude NDVI trends tended to co-occur with land cover change, with areas of change experiencing median rates of greening and browning that were 47%-128% larger than corresponding rates in areas of stable land cover. NDVI trends in areas of stable land cover were predominantly positive and were widely distributed throughout the study domain. In contrast, areas of land cover change exhibited both greening and browning, and because they were associated with disturbance, tended to be more localized. These results imply that the largest changes in NDVI are caused by changes in land cover rather than ecological changes within land cover classes (e.g. increased leaf area), which has important ramifications for the interpretation of greening and browning trends. For example, Evergreen Forest Gain and Deciduous Forest Gain exhibit similar trend magnitudes (figure 5) but have different rates of carbon cycling and sensitivities to climate, resulting in different ecological outcomes despite similar trends (Zimov et al 1999, Welp et al 2007, Augusto et al 2015. Furthermore, ecosystem models are often parameterized using land cover in combination with vegetation indices that serve as proxies for photosynthetic activity (Running et al 2004, Mahadevan et al 2008. Hence, since many NDVI changes at high-latitudes are caused by land cover change, studies that use greening and browning trends to characterize changes in productivity must also account for land cover change. Boreal browning and arctic greening trends Browning was mostly prevalent in the Boreal biome, which is often interpreted as reflecting stress-induced forest decline. However, Boreal browning may actually reflect a combination several different ecological processes. For example, Evergreen Forest Loss was associated with a quarter (25.7%) of pixels showing browning trends and was largely caused by disturbances like fire and harvest, rather than stress-induced forest decline. This is comparable in scale to the contribution of stable forests to browning. Roughly one third of browning trends were in stable forests, including both Evergreen Forest (25.1% of browning trends) and Deciduous Forest (8.3% of browning trends). Browning trends in stable forests likely reflect subtle modes of forest decline arising from drought However, browning in some areas may reflect successional shifts in plant functional types. Counterintuitively, browning occurred in 7.1% of Evergreen Forest Gain areas, which is not generally compatible with stress or decreasing productivity ( figure 4(d)). An alternative explanation is that browning in these forests reflects successional shifts from Deciduous Forests, which typically have high NDVI, towards later successional Evergreen Forests with lower NDVI (figure 8). Similar processes may be driving browning in mature Deciduous Forests that are transitioning to Evergreen Forests. The ecological interpretation (i.e. succession versus climate driven forest decline) is very different between these scenarios, and information on land cover change is essential for interpreting NDVI trends.
In contrast, NDVI trends in Arctic areas above the tree line showed widespread greening with negligible browning and disturbance (figure 2). Greening in the Arctic was often associated with areas of Shrub Gain and Herbaceous Gain, as well as stable Herbaceous and Shrub cover, potentially reflecting increases in vegetation cover and leaf area that are consistent with field studies (Tape et al 2006, Myers-Smith et al 2011. Challenges using NDVI to characterize ecosystem changes We expect that most land cover changes are accompanied by changes in NDVI. However, while most land cover change pixels showed statistically significant NDVI trends, a substantial proportion (44.3%, figure 1(b)) did not. Non-significant NDVI trends were found in all land cover change types, with the highest proportion in areas of Shrub Loss ( figure 4(d)). The magnitude of NDVI trends is strongly influenced by the timing of disturbance, where present (Sulla-Menashe et al 2018). NDVI trends in areas that experienced disturbance-driven reductions in NDVI in the middle of the time series are generally smaller because NDVI tends to recover after disturbance events, and changes in NDVI that occur in the middle of time series have low leverage over estimated trends relative to those near the beginning or end of the time series. Hence, non-significant NDVI trends can be a misleading indicator that surface properties have been stable over the time series, when in fact the opposite is true.
Detecting ecosystem change using NDVI can also be challenging because NDVI saturates over areas with high leaf area (Myneni and Williams 1994, Carlson and Ripley 1997, Buermann et al 2002. Hence, changes in leaf area or biomass are difficult to detect in locations with high NDVI (i.e. >0.8), such as Deciduous Forests (figure 8). Indeed, greening in Deciduous Forests was relatively limited (18.5% of Deciduous Forests, figure 4(c)) and was of generally lower magnitude compared to other land cover classes (figure 5). The low sensitivity of NDVI to Deciduous Forest growth is particularly significant because this land cover type has expanded in area by as much as 15% since 1984 (Wang et al 2019b). Hence, the growth response to climate change in a growing area of Boreal forests may be masked by this saturation effect.
Similarly, detecting land cover transitions between Evergreen Forests and Shrubs can be challenging since they have comparable growing season NDVI values (figure 8). This may explain why significant NDVI trends are absent in nearly one third of areas with Evergreen Forest Loss, which typically transition to Shrub or Herbaceous following disturbance ( figure 4(d)). Further, treeline expansion into tundra, which has been previously documented (Lescop-Sinclair and Payette 1995, Harsch et al 2009), may not lead to detectable NDVI trends since transitions from Shrub to Evergreen Forest are subtle. Hence, relying solely on NDVI can lead to the incorrect conclusion that relatively little change has occurred in some areas. Figures 6 and 7 illustrate the spatial heterogeneity in both NDVI trends and disturbances, like fires and harvest, in Boreal ecosystems. Coarse spatial resolution (i.e.500 m) instruments like AVHRR or MODIS provide valuable global imagery, but the scale and complexity of high-latitude ecosystem limit the utility of such data for studying these systems (Ju and Masek 2016). Furthermore, these are slow growing systems that sensors with shorter data records, such as MODIS, may not accurately characterize. Moderate spatial resolution instruments, such as those onboard Landsat (Woodcock et al 2008) or Sentinel 2 (Drusch et al 2012), resolve landscape-scale changes, and are therefore required to interpret NDVI trends in heterogeneous and dynamic high-latitude ecosystems. In fact, using moderate resolution imagery, we find that fires and harvest drive 27.8% of browning trends and 6.5% of greening trends. Advances in computational resources and availability (e.g. Gorelick et al 2017) are facilitating analysis of moderate resolution imagery across large high-latitude regions, although lingering issues remain regarding cross-sensor calibration (Guay et al 2014, Ju and Masek 2016), lower temporal frequency, and the relative paucity of Landsat data in some regions, such as Siberia.

Conclusions
Studies of high-latitude greening and browning suggest that large-scale changes in Arctic-Boreal ecosystem productivity have occurred in recent decades. However, reliance solely on NDVI trends as representing shifts in ecosystem productivity substantially oversimplifies these changes. Over a quarter of significant NDVI trends in the ABoVE Core Study Domain were associated with land cover change, and land cover change drove the largest trends, primarily due to forest disturbance and regrowth. Hence, understanding land cover composition and history is crucial to understanding the changes in ecosystem productivity and climate feedbacks that NDVI time series are indirectly capturing.
We suggest that studies characterizing greening and browning trends in high-latitude ecosystems must also consider disturbance and land cover change at both regional and landscape scales. The relatively long record and moderate spatial resolution of Landsat imagery provide previously unexplored opportunities to better characterize and understand high-latitude land change processes. By combining time series of Landsat imagery with growing sources of ground data from experiments such as ABoVE, increasing availability of computational resources, and new sources of imagery such as from Sentinel-2, more nuanced characterizations of recent high-latitude ecosystem change are both necessary and increasingly possible.

Acknowledgments
This study was part of the Arctic-Boreal Vulnerability Experiment. Resources supporting this work were provided by the NASA High-End Computing Program through the NASA Center for Climate Simulation at Goddard Space Flight Center. This work was supported by a National Science Foundation Graduate Research Fellowship under Grant No. DGE-1247312, by NASA ABoVE grant No. NNX15AU63A, and by NASA Grant No. 80NSSC18K0994. We gratefully acknowledge Lucy Hutyra, Curtis Woodcock, Oliver Sonnentag, and two anonymous reviewers for providing valuable feedback that improved the manuscript.

Data availability statement
The data used in this study are publicly available at the Oak Ridge National Laboratory Distributed Active Archive Center (ORNL DAAC). The annual land cover data from Wang et al (2019a) are available at https:// doi.org/10.3334/ORNLDAAC/1691, the NDVI trend data from Ju and Masek (2018) are available at https:// doi.org/10.3334/ORNLDAAC/1576, and the peak growing season NDVI data from Melaas et al (2019) are available at https://doi.org/10.3334/ ORNLDAAC/1698. The land cover change dataset is available upon request.