Biophysical Determinants of Shifting Tundra Vegetation Productivity in the Beaufort Delta Region of Canada

Temperature increases across the circumpolar north have driven rapid increases in vegetation productivity, often described as ‘greening’. These changes have been widespread, but spatial variation in their pattern and magnitude suggests that biophysical factors also influence the response of tundra vegetation to climate warming. In this study, we used field sampling of soils and vegetation and random forests modeling to identify the determinants of trends in Landsat-derived Enhanced Vegetation Index, a surrogate for productivity, in the Beaufort Delta region of Canada between 1984 and 2016. This region has experienced notable change, with over 71% of the Tuktoyaktuk Coastlands and over 66% of the Yukon North Slope exhibiting statistically significant greening. Using both classification and regression random forests analyses, we show that increases in productivity have been more widespread and rapid at low-to-moderate elevations and in areas dominated by till blanket and glaciofluvial deposits, suggesting that nutrient and moisture availability mediate the impact of climate warming on tundra vegetation. Rapid greening in shrub-dominated vegetation types and observed increases in the cover of low and tall shrub cover (4.8% and 6.0%) also indicate that regional changes have been driven by shifts in the abundance of these functional groups. Our findings demonstrate the utility of random forests models for identifying regional drivers of tundra vegetation change. To obtain additional fine-grained insights on drivers of increased tundra productivity, we recommend future research combine spatially comprehensive time series satellite data (as used herein) with samples of high spatial resolution imagery and integrated field investigations.

ABSTRACT Temperature increases across the circumpolar north have driven rapid increases in vegetation productivity, often described as 'greening'. These changes have been widespread, but spatial variation in their pattern and magnitude suggests that biophysical factors also influence the response of tundra vegetation to climate warming. In this study, we used field sampling of soils and vegetation and random forests modeling to identify the determinants of trends in Landsat-derived Enhanced Vegetation Index, a surrogate for productivity, in the Beaufort Delta region of Canada between 1984 and 2016. This region has experienced notable change, with over 71% of the Tuktoyaktuk Coastlands and over 66% of the Yukon North Slope exhibiting statistically significant greening. Using both classification and regression random forests analyses, we show that increases in productivity have been more widespread and rapid at low-to-moderate elevations and in areas dominated by till blanket and glaciofluvial deposits, suggesting that nutrient and moisture availability mediate the impact of climate warming on tundra vegetation. Rapid greening in shrub-dominated vegetation types and observed increases in the cover of low and tall shrub cover (4.8% and 6.0%) also indicate that regional changes have been driven by shifts in the abundance of these functional groups. Our findings demonstrate the utility of random forests models for identifying regional drivers of tundra vegetation change. To obtain additional fine-grained insights on drivers of increased tundra productivity, we recommend future research combine spatially comprehensive time series satellite data (as used herein) with samples of high spatial resolution imagery and integrated field investigations.
Key words: Landsat; random forests; remote sensing; low arctic; arctic tundra; greening; EVI; climate change; land cover.  Lantz 2020). Vegetation productivity can also be measured at broad scales using multispectral satellite vegetation indices (Gao and others 2000) such as the Enhanced Vegetation Index (EVI) and the Normalized Difference Vegetation Index (NDVI). Changes in these indices have been observed across the Arctic with increasing productivity referred to as 'tundra greening' (Jia and others 2003;Bhatt and others 2010;Epstein and others 2012). Continental and pan-Arctic scale changes in vegetation productivity have generally been attributed to rapid temperature increases at high latitudes (Jia and others 2003;Bhatt and others 2010;Miller and Smith 2012;Fraser and others 2014a;Berner and others 2020). Plot-scale warming experiments and repeated observation also provide evidence that vegetation change has been caused by increasing temperature (Chapin and others 1995;Walker and others 2006;Hudson and Henry 2009;Elmendorf and others 2012).
Although widespread, observed increases in Arctic vegetation productivity have not been uniform, with some regions exhibiting stable or decreasing productivity (Jia and others 2006;Bhatt and others 2010;Epstein and others 2012;Ju and Masek 2016). Recent evidence suggests that variation in the response of tundra vegetation is related to both broad-scale (Wang and Friedl 2019; Berner and others 2020; Campbell and others 2021; Chen and others 2021) and fine-scale (Moffat and others 2016; Myers-Smith and others 2019; Bjorkman and others 2020) variation in biophysical factors, but few studies have explored linkages between these scales. In their conceptual model, Pearson and Dawson (2003) suggest that climate variables (long-term temperature and precipitation trends) have greater influence at a global and continental scales, whereas biophysical variables (such as soil moisture, surface topography and soil conditions, land cover, and land use) are likely to influence processes at landscape or local scales.
At the landscape-scale, research suggests that variability in soil moisture, land cover type, and landscape position are responsible for the heterogeneous response of Arctic vegetation productivity Tape and others 2012;Martin and others 2017;Bonney and others 2018;Campbell and others 2021). In this study, we explore the environmental factors driving heterogeneity in vegetation productivity trends across the Beaufort Delta region by combining plot-based fieldwork with an analysis of the Landsat satellite archive (Wulder and others 2019). We use a random forests (RF) ensemble decision tree algorithm to determine the environmental factors influencing spatial heterogeneity in tundra productivity trends in the Beaufort Delta region. Ultimately, we seek to understand how variation in environmental conditions influences the spatial patterns of vegetation productivity at a landscape-scale. An improved understanding of factors mediating tundra vegetation change will contribute to local and regional planning and inform earth system models that include feedbacks between vegetation growth and ecological processes such as permafrost dynamics, albedo, and evapotranspiration (Verseghy 1991;Verseghy and others 1993;Verseghy 2000;Bonan and others 2003;Quillet and others 2010).

STUDY AREA
This study focuses on the tundra ecosystems across the Yukon North Slope (0.54 MHa) and Tuktoyaktuk Coastlands (2.92 MHa). These Low Arctic ecosystems are located in the Beaufort Delta region of the western Canadian Arctic (Figure 1a). This coastal region borders the Beaufort Sea to the north and is bounded by the northern edge of the subarctic forest to the south (Timoney and others 1992). Both regions are located within the Inuvialuit Settlement Region and are significant to the communities of Tuktoyaktuk (population of 900), Inuvik (population of 3200), and Aklavik (population of 600), who use these lands for hunting, fishing, trapping, traditional harvesting of plants, and other cultural practices (Alunik and Morrison 2003;Murray and others 2005;Tyson and others 2016).
The Yukon North Slope extends from the foothills of the Richardson Mountains to the coast of the Beaufort Sea. Areas at high elevation drain to the Beaufort via the deeply incised canyons of the Babbage, Blow, and Big Fish Rivers among others (Rampton 1982;Yukon Ecoregions Working Group [YEWG] 2004). Rolling hills and hummocky terrain dominate the eastern portion of the Yukon North Slope that ends in a steep escarpment at the Mackenzie Delta (Rampton 1982). Average annual temperature between 1984 and 2016 at Shingle Point was -9.2°C with an average summer temperature (June through August) of 9.0°C (Environment and Climate Change Canada [ECCC] 2018b). On average, Shingle Point received 254 mm of precipitation annually (ECCC 2018b), about half of which fell as rain (Burn and Zhang 2009). The Tuktoyaktuk Coastlands extend from the eastern limit of the Mackenzie Delta to Cape Bathurst to the east. This gently rolling landscape is  (Fulton 1989) and c elevation in metres (Porter and others 2018) used in random forests analyses across the study area. scattered with lakes and ponds and is characterized by hummocky terrain in upland areas and polygonal terrain and wetlands at lower elevations (Rampton 1988;Ecosystem Classification Group [ECG] 2012). The climate at Tuktoyaktuk is similar to Shingle Point, with average annual and summer temperatures between 1984 and 2016 of -9.4°C and 8.9°C, respectively (ECCC 2018a). Average total precipitation during this period was 146 mm, about 40% of which fell as rain (ECCC 2018a). The Yukon North Slope and Tuktoyaktuk Coastlands are separated by the low-lying alluvial terrain in the Mackenzie Delta ecoregion. This dynamic ecosystem was excluded from our study area because much of it is forested and the tundra communities present are strongly influenced by hydrological dynamics (Gill 1972(Gill , 1973Pearce 1986;Burn and Kokelj 2009).
The Yukon North Slope hosts a diversity of terrain types including coastal beaches and estuaries, low-lying wetlands, upland tussock tundra, and shrub tundra (YEWG 2004; Wang and others 2019). The upland tundra in the foothills of the Richardson Mountains occurs on well-drained soils that support characteristic communities of tall and dwarf shrubs, lichens, graminoids and forbs (YEWG 2004). The Tuktoyaktuk Coastlands are largely covered by shrub and tussock tundra with wetter areas dominated by sedge and moss tundra (ECG 2012; Moffat and others 2016). In the southern portion of this region, scattered spruce woodlands are located along sheltered creeks and other lowlying areas (Lantz and others 2019).
Both the Yukon North Slope and Tuktoyaktuk Coastlands are underlain by continuous permafrost and are characterized by thermokarst features including polygonal terrain, earth hummocks, thaw slumps and pingos (Rampton 1982(Rampton , 1988YEWG 2004;ECG 2012). The Laurentide Ice Sheet covered most of this region during the Wisconsinan glaciation with the exception of land south and west of the Richardson Mountains in northern Yukon and the northern tip of the Tuktoyaktuk Peninsula and Cape Bathurst (Jessop 1971;Hughes and others 1981;Duk-Rodkin and Hughes 1995;YEWG 2004;ECG 2012).

METHODS
To investigate the drivers of vegetation change in the Beaufort Delta region, we combined RF modelling of regional EVI trends with multivariate analyses of plot-scale field data. We classified pixelbased trends in EVI (1984EVI ( -2016 as: (1) exhibiting significant increases in EVI or (2) un-trended. We used these binary classes (increasing EVI/un-tren-ded EVI) as response variables in a classification RF model to determine the factors facilitating and constraining increased productivity. To identify the factors influencing the magnitude of greening, we also performed a regression RF model that predicted the slope of significant positive EVI trends. To facilitate site selection for field sampling and multivariate analyses, pixels exhibiting significant increases in EVI were further classified as moderate or high magnitude greening.

EVI Trend Analysis
To document changes in the productivity of tundra vegetation across the study area, we tracked changes in EVI using the Landsat satellite archive (1984 -2016). EVI is a modified form of NDVI in which the blue band of the visible spectrum and satellite-specific correction terms (C 1 , C 2 , and L) account for soil (L) and atmospheric/aerosol (C 1 and C 2 ) influences (Gao and others 2000).
We used EVI in this analysis because it is sensitive to differences in vascular plant phytomass and vascular plant net primary productivity (Kushida and others 2015) and it performs well across a range of soil moisture conditions (Raynolds and Walker 2016). We obtained the EVI trend surface from Chen and others (2021) who generated the trend surface using annual composite Landsat images from 1984 to 2016. These images were obtained using the Composite2Change method for selecting best-available-pixels from the Landsat archive to produce a gap-free surface reflectance raster for each year (Hermosilla and others 2016). The imagery used to produce best-available-pixel composites was cross-calibrated among Landsat sensors (see: Markham and Helder 2012;Hermosilla and others 2016). Composite imagery was used to calculate EVI on a pixel basis (30-m spatial resolution) for each year and a pixel-based, nonparametric Theil-Sen regression was performed on the resulting time series (Theil 1950;Sen 1968) using the EcoGenetics package (Roser and others 2017) in R (R Core Team 2019). The significance of pixel-based slopes was assessed using a Mann-Kendall test for monotonic trends (Mann 1945;Kendall 1948) performed using the Kendall package in R (McLeod 2011). In our analysis, we used the resulting p-values to determine whether significant change (p < 0.05) had occurred in each pixel from 1984 to 2016.
We also used the raster surface derived from Theil-Sen regression of EVI trends to classify the study area into site types exhibiting: (1) high greening, (2) moderate greening, (3) no significant change (stable sites), and (4) browning (see histogram: Supplementary Fig. 1) for the purpose of our field investigation. Areas showing non-significant EVI trends (p > 0.05) were classified as stable. Pixels with significant increases in EVI were further classified based on the slope of the EVI trend. Pixels with slopes that were within one standard deviation (SD; 1.45 9 10 -3 per year) of the mean EVI trend (< 2.24 9 10 -3 + SD per year) were classified as moderate greening, and pixels with an increase greater than one standard deviation above the mean EVI trend (> 2.24 9 10 -3 + SD per year) were classified as high greening. Pixels with a significant negative slope (decline in EVI over time) were classified as browning. Since browning pixels only accounted for 0.63% of pixels in the study area, we did not consider this class as a category in this analysis.

Random Forests Analysis and Variable Importance
To identify the biophysical variables that best explain spatial variation in EVI trends, we used two RF models (Breiman 2001). RF models are decision-tree based, ensemble machine-learning methods involving the assembly of many regression or classification trees using a subset of available data to increase predictive capability (Breiman 2001;De'ath 2007). The RF method can also be used to determine variable importance by measuring the mean decrease in model accuracy upon removal of a given variable (Cutler and others 2007). In our first analysis, we used a classification RF to discriminate pixels showing a significant increase in EVI (p < 0.05, positive slope) from un-trended pixels (p > 0.05). In a second analysis, we used a regression RF to model the magnitude of the change in EVI (Theil-Sen regression slope) of pixels exhibiting significant increases in EVI (p < 0.05, positive slope) using the same suite of biophysical variables. We created both models using the ran-domForests package (Liaw and Wiener 2002) in R. Explanatory variables used in each model are presented in Table 1. We ran the classification RF model using a random selection of 40,000 pixels of each class to ensure balanced sampling. For the regression RF, we used a random sample of 1% (66,442 pixels) of all pixels in the study area. We ran each model using 1000 trees and calculated the variable importance using the importance function (Liaw and Wiener 2002). We assessed the influence of the four variables that had the largest impact on model accuracy using partial dependence plots showing the marginal effect of a given variable on the modelled parameter while keeping all other variables constant (Friedman 2001;Hastie and others 2009). For the classification RF, our partial dependence plots show the probability of significant greening. The partial dependence plots for a regression RF show the predicted magnitude of the slope in the EVI trend.

Explanatory Variables
Broad-scale biophysical data used in this study were obtained from a variety of sources and data processing and aggregation were completed using the R statistical software (Table 1; R Core Team 2019). We selected explanatory variables known to effect growth and productivity but were not able to directly assess the influence of microclimate because suitable data are not available to adequately capture climate variability at such fine scales. A digital elevation model (DEM), with 2-m spatial resolution (Figure 1c), was sourced from the Polar Geospatial Center's (PGC) ArcticDEM dataset (Porter and others 2018). We calculated slope, aspect, and terrain ruggedness index (TRI) from this DEM using the terrain function from the raster package (Hijmans 2020). TRI is the absolute difference between the elevation of a given cell and its surrounding neighbours (Riley and others 1999). The topographic wetness index (TWI) is a DEMdriven index of soil moisture in which potential soil moisture availability is calculated based on the slope and flow directions of the surrounding landscape (Kopecký and Č ížková 2010). Flat areas surrounded by upslope terrain will have higher TWI values (as moisture will tend to accumulate in these areas) compared to steep sloping cells that will have lower values (as moisture will tend to run off). TWI was calculated using the upslope.area function from the dynatopmodel package (Metcalfe and others 2018). Topographic position index (TPI) is another DEM-driven index that defines the relative elevation of a cell based on the mean elevation of surrounding cells (De Reu and others 2013).
Where TPI is positive, the cell has a higher elevation than the mean of the surrounding cells. A negative TPI indicates that the cell has a lower elevation than the mean of the surrounding cells. TPI was calculated using the tpi function in the spatialEco package (Evans 2020). We calculated solar insolation using the equation provided by Roberts and Cooper (1989): This index of solar insolation ranges from 0 to 1, with north-northeast facing aspects taking on values close to zero, and south-southwest aspects having values closer to 1. Data on land cover were obtained from the classification created by Wang and others (2019). This land cover classification includes 15 terrain types derived from a RF classification model for Landsat surface-reflectance at 30-m resolution using high-resolution field photos and imagery from the NASA ABoVE project to assign land cover classes for each year between 1984 and 2014 (Wang and others 2019). We used land cover data from 1984 as a predictor in the RF models to test the sensitivity of different vegetation classes to shifts in productivity. The average overall accuracy for these land cover data is 84.1 ± 4.1% (mean and 95% confidence interval across all years of the time series; Wang and others 2019). Surficial geology data were obtained from Fulton (1989), cropped to the extent of the study area and rasterized for further analysis (Figure 1b). We masked out any cells of the land cover and geology layers that were occupied by classes representing less than 1% of the study area (class descriptions can be found in Supplementary Tables 1 and 2). Spatial resolution of the data sources described above was matched to the EVI trend data (30 m). Where resolution of a dataset was finer than 30 m, cells were aggregated by taking the mean of sub-pixels. We resampled continuous data using a bilinear interpolation method while categorical data used the nearest neighbour method. These operations were performed with the resample function from the raster package (Hijmans 2020). Surficial geology and elevation data are shown in Figure 1b and c , and all other explanatory variables can be found in Supplementary Fig. 2.

Field Surveys
To assess the influence of biophysical variables on recent changes in tundra productivity, we measured biotic and abiotic variables at sites across the study area. Field sites were selected by randomly choosing 40 points in each of the productivity classes defined using the EVI trend (moderate greening, high greening, and stability). Randomization was constrained such that the pixel (30 m 2 ) containing the selected point was surrounded by the same greening class. We visited 21 stable sites, 32 moderate greening sites, and 26 high greening sites in July and August of 2019 (Figure 1a). At each site, we measured a suite of biotic and abiotic variables, as described below.
Vegetation surveys and soil sampling were conducted along two 30-m transects. Transects were orientated in north-south and east-west directions such that the 15-m mark on both transects was centered on the predetermined site coordinate. These dimensions were selected to create a plot corresponding to the 30-m resolution of one Landsat pixel. Thaw depth and soil moisture were measured at 5-m intervals along each transect as well as inside all vegetation quadrats at a site. We measured soil moisture at 3-5 cm below the surface using a handheld moisture probe (Delta-T Devices HH2 Moisture Meter with ML3 ThetaProbe Soil Moisture Sensor). Thaw depth was measured using a graduated metal probe inserted into the ground until the depth of refusal. We visually estimated the percent cover of plant species or species groups inside nested 5 m 2 and 1 m 2 quadrats at four locations at every site. We positioned quadrats using a random point located inside each of the four quadrants of the cross transect (northwest, northeast, southwest, and southeast). We estimated the cover of upright shrubs and trees using the larger quadrat and the cover of dwarf shrubs, graminoids, herbaceous species, lichens, and bryophytes using the smaller quadrat.
When vegetation cover estimates were completed, we collected a composite active layer soil sample from within each quadrat using a small shovel. The profile exposed during sample collection was also used to estimate the thickness of the moss layer and organic soil horizon. Soil samples were stored at -20°C before they were submitted for analysis of gravimetric soil moisture, macronutrients including nitrogen and phosphorus, and micronutrients including magnesium, sulphur and calcium. Chemical analyses were carried out using an inductively coupled plasma mass spectrometer (ICP-MS) at the Chemical Services Laboratory at the Pacific Forestry Centre of the Canadian Forest Service in Victoria, British Columbia.

Vegetation Community and Environmental Data Analysis
To explore differences in community composition among sites exhibiting different levels of landscape scale greening, we used non-metric multidimensional scaling (NMDS) ordination. We applied a log(x + 1) transformation to the raw vegetation percent cover data and ran the NMDS using a Bray-Curtis dissimilarity matrix of the transformed percent cover data. This analysis was set to repeat 100 times and select the best two-dimensional representation of the original data. The NMDS was performed using the metaMDS function from the vegan package (Oksanen and others 2019) in R. We used the analysis of similarities (ANOSIM; Clarke 1993) function (anosim) from vegan (Oksanen and others 2019) to test for statistically significant differences in vegetation community composition among the three sites types (stable, moderate greening, and high greening). To identify the species making the greatest contribution to differences among site types, a similarity percentage (SIMPER) analysis was conducted using the simper function in the vegan package (Oksanen and others 2019). We also compared the cover of vegetation classes in 1984 and 2014 using the multiyear land cover classification presented in Wang and others (2019). Specifically, we calculated the percent cover of land cover classes in 1984 and 2014 and compared the differences across the entire study area and between the Yukon North Slope and Tuktoyaktuk Coastlands.

Study Area Response Overview
Approximately 70% of the study area showed significant increases in EVI between 1984 and 2016 ( Figure 2, Table 2). Split across the study area, 71% of the Tuktoyaktuk Coastlands showed significant greening compared to 66% of the Yukon North Slope ( Table 2). The Tuktoyaktuk Coastlands also had a higher average EVI trend than the Yukon North Slope but was similar to the average across the entire study area (Table 2).

Classification Random Forests Analysis
The top four predictors in the classification RF were surficial geology, elevation, land cover in 1984, and TWI ( Figure 3). This model had user's accuracy of 72.05% and was capable of classifying pixels as stable and greening with user's accuracy of 73.5% and 70.6%, respectively. The area under the curve (AUC) of the receiver operating characteristic (ROC) is 0.796 (see Supplementary Fig. 3).
Partial dependence plots for surficial geology show that till blanket, glaciofluvial complex, glaciofluvial plain, and lacustrine sand were more likely to exhibit greening than colluvial and alluvial materials, which had a greater probability of being stable (Figure 4a). The partial dependence plot for elevation shows that the greatest probability of greening occurred at elevations between 0 to 60 m (Figure 4b). Elevations above approximately 60 m had a higher probability of being stable. This analysis also showed that significant greening was most likely in low shrub, tussock tundra, and sparsely vegetated land cover classes, while tall shrub, fen, and barren cover classes were more likely to be stable ( Figure 4c). With respect to topographic variables, the partial dependence plot for topographic wetness (TWI) indicates an increased probability of greening at moderate landscape wetness (Figure 4d). Values of TWI greater than 7.5 are largely associated with the margins of ponds and lakes as well as coastal tundra areas on the Tuktoyaktuk Peninsula.

Regression Random Forests Analysis
The four most important variables in the regression RF, which predicted the magnitude of greening, were surficial geology, elevation, and land cover, and TWI followed by a suite of variables related to topography ( Figure 5). This model explained 20.41% of the variance in EVI trend based on the out-of-bag data.
The partial dependence plots for the regression RF show that lacustrine sands, till blanket, and glaciofluvial complex cover are associated with more rapid greening (Figure 6a). Areas underlain by fine colluvium and till veneer had rates of EVI change lower than the average value (3.07 9 10 -3 per year) for significantly trended pixels (Figure 6a). Areas covered by till blanket comprised roughly 46% of the entire study area, of which, over 51% are significantly greening (Table 3). The partial dependence plot for elevation shows that greening was most rapid at lower elevations (Figure 6b). At elevations above 50 m, the calculated mean EVI trend was below the average EVI trend of significantly greening pixels (red dashed line; Figure 6b). Tussock tundra, low shrub, and sparsely   vegetated land cover classes exhibited the greatest predicted EVI trends (Figure 6c). The mean EVI trends across low shrub and tussock tundra classes were greater than the average EVI trend of significantly greening pixels (Table 4). The predicted EVI trend as a function of TWI decreased with higher index values, dropping below the average mean EVI trend of significantly greening pixels at approximately TWI of 8 (Figure 6d).
Analysis of data from Wang and others (2019) indicates that sparse and herbaceous cover decreased by 7.5% and 6.2% across the entire study area between 1985 and 2014. These data also show that low shrub and tall shrub cover increased by 4.8 and 6.0% over this period (Figure 8). The pattern of vegetation change was similar in both regions, but the magnitude of increases in shrub-dominant terrain were higher in the Tuktoyaktuk than the Yukon North Slope. The decline in herbaceous cover was also greater in the Tuktoyaktuk Coastlands compared to the Yukon North Slope (Figure 8).

DISCUSSION
Landscape-scale variation in surficial materials, topography, and vegetation structure were good predictors of changing tundra productivity because these factors influence access to both moisture and soil nutrients. Surficial materials govern the development of tundra soils and control nutrient and moisture availability. Parent materials impact soil formation through differential weathering and the mineral and chemical composition of the substrate (Brady and Weil 1996;Walker 2000). Poorly sorted materials such as till veneer are often associated with exposed bedrock and are found at higher elevations in the foothills of the Richardson Mountains and on the Anderson Plain (Fulton 1989). These thin tills contain larger gravel and boulders and result in shallow, poorly developed soils that are unable to retain moisture and nutrients (Brady and Weil 1996). Till veneer was associated with stable vegetation in our study area, while the deeper soils that typically develop in areas of till blanket or glaciofluvial complexes were associated with greening ( Figure 4a) and had higher EVI trend values than the average across the study area (Table 3). Finer-grained soils at lower elevation have a larger surface area, are less susceptible to leaching, and are able retain more nutrients for plant uptake (Walker and Everett 1991;Brady and Weil 1996).
The importance of elevation in our RF models indicates that microclimatic variation associated with landscape position also impacts tundra vege-  tation productivity. Our observation that lower rates of greening and a reduced probability of greening were associated with higher elevations suggests that cold temperatures, dry soils, and reduced snowpack on hilltops limit the effects of increasing regional temperatures on tundra productivity. Soils at higher elevation also tend to have more unconsolidated sediments and exposed bed-rock, which likely limit productivity because of reduced moisture and nutrients.
Differences in relative elevation drive variation in soil moisture that likely influence vegetation community development and productivity. Previous research shows that landscape scale variation in temperature, soil moisture and snow accumulation can impact the establishment and growth of tundra plants (Sturm and others 2001; Myers-Smith and others 2011; Niittynen and others 2020a; Niittynen and others 2020b). Our NMDS analysis also showed an association between elevated soil moisture and plant community composition at moderate and high greening sites (Figure 7). A number of recent studies have demonstrated that soil moisture influences tundra vegetation growth (Myers-Smith and others 2015; Cameron and Lantz 2016; Ackerman and others 2017; Bjorkman and others 2018) and indices of vegetation productivity (Campbell and others 2021; Chen and others 2021). Topographic variation in microclimate and soil moisture have also been shown to be better predictors of tundra greening at fine scales than broader-scale climaterelated factors such as the total length of growing season (Gamon and others 2013). Soil moisture levels influenced by microclimate can drive increased nutrient mineralization and are likely an important mechanism of increased productivity (Chapin and others 1988; Deslippe and Simard 2011; Deslippe and others 2012; Mekonnen and others 2021). More rapid greening associated with moderate levels of topographic wetness indicate that tundra vegetation on mesic to moist soils is most sensitive to regional warming. This is likely because moderate soil moisture limits the negative effects of temperature-induced moisture stress on tundra growth and productivity (Johnson and Caldwell 1975;Dagg and Lafleur 2011;Myers-Smith and others 2015;Ackerman and others 2017). Snowmelt associated with the onset of spring can also increase soil moisture in areas of large drifts, contributing to an influx of moisture to tundra soils. Gamon and others (2013) found that earlier snowmelt was associated with drier soils and lower mid-season NDVI which highlights the complex interaction between moisture availability and growing season. Winter conditions (notably snow cover) are also important to consider since winter is a crucial period in the development of vegetation communities and their functional diversity (Niittynen and others 2020a; Niittynen and others 2020b).
Vegetation type was a good predictor of EVI trends because some functional groups, such as Figure 6. Partial dependence plots for the top four variables used in the regression random forest predicting EVI trend (per year): a Quaternary surficial geology (Fulton 1989), b elevation in metres (Porter and others 2018), c land cover in 1984 (Wang and others 2019) and d topographic wetness index. Note differing yaxis ranges among plots A-D. The dashed line on the yaxis in plots B and D highlights the mean EVI trend value for significantly greening pixels (3.07 9 10 -3 per year). Plots B and D are presented with decile rug marks of training data.
shrubs and graminoids, are more responsive to changes in temperature and soil conditions. Specifically, many species of deciduous shrub are adapted to respond to changes in nutrients availability, the length of the growing season, and moisture availability at the onset of spring (Chapin and others 1995;Hobbie and Chapin 1998;Bret-Harte and others 2001;Myers-Smith and others 2011;Kelsey and others 2020). Increased productivity in shrub tundra communities was likely also driven by the ability of deciduous shrubs to rapidly allocate resources to secondary growth and asexual reproduction (Bret-Harte and others 2001; Wiedmer and Senn-Irlet 2006; Ropars and Boudreau Figure 7. Non-metric multidimensional scaling (NMDS) ordination of vegetation community composition across the a Yukon North Slope (stress = 0.15) and b Tuktoyaktuk Coastlands (stress = 0.13). Sites split by classification with stable, moderate greening, and high greening shown as blue, green, and red points and polygons, respectively with abiotic vectors (p < 0.1) overlaid. The ordinations also plot associations between species and functional type (lichen) to cumulative contribution of 65% and NMDS score. Species abbreviations are defined in Table 5. 2012). Taller, shrub dominated vegetation may also increase soil moisture and nutrient availability via increased snow capture (Sturm and others 2001;Wipf and Rixen 2010;Leffler and others 2016;Niittynen and others 2020a). Observed increases in the productivity of tussock tundra may be related to variability of substrate properties including grain size and soil chemistry that are predicted to favour development of moist-acidic tundra with climate warming (Walker and others 1998;  The table shows the species and species groups accounting for 80% of the total dissimilarity and is ranked by the contribution to dissimilarity. The site type showing higher cover for a given species is shown in bold. The NMDS Code column lists the abbreviations used in Figure 7. This analysis was not completed for the Yukon North Slope because there were no differences in community composition among sites. others 2004a; Epstein and others 2012). The tussock growth form of species such as Eriophorum vaginatum is also highly conducive to growth in nutrient poor environments through nutrient cycling within tussocks (Cholewa and Griffith 2004) and growth of deeper roots (Chapin and others 1988). Our results are also consistent with a number of studies linking tundra greening with productive cover types (Jia and others 2006;McManus and others 2012;Frost and others 2014;Campbell and others 2021;Chen and others 2021) and shrub proliferation in the low Arctic in particular (Tape and others 2006;Ropars and Boudreau 2012;Lantz and others 2013;Frost and others 2014;Moffat and others 2016). Our comparisons of land cover over time also show that the proportion of shrub and tussock tundra have increased since 1984 in our study regions. Higher cover of shrubby vegetation types in the Tuktoyaktuk Coastlands likely drove differences in the extent and distribution of greening between regions. Past research in this region has also shown an association between gains in shrub and herbaceous cover and increases in NDVI . The Yukon North Slope has experienced less of a transition to shrub tundra than the Tuktoyaktuk Coastlands ( Figure 8) with greater overall similarity in community composition across the region (Figure 7) and we attribute the difference in the extent of observed greening between these two regions ( Table 2) to differences in the intensity of shrub proliferation. Average annual and summer temperatures across the Beaufort Delta region have increased by 3. 5°C and 1.9°C, respectively, between 1926(Travers-Smith and Lantz 2020 with largely homogenous warming across the region (Vincent and others 2015). Widespread increases in tundra productivity are likely a response to the direct and indirect effects of this warming. Productivity responses to warming are well documented, but environmental limitations also drive differences in the response of vegetation across the Arctic (Walker and others 2006;Hudson and Henry 2009;Elmendorf and others 2012;Myers-Smith and others 2015). NDVI trends across the Arctic Coastal Plain of Alaska also show complex responses to changes in temperature and precipitation across terrain types, suggesting that climatic drivers are mediated by regional environmental factors (Lara and others 2018). Underlying climate drivers, landscape and regional variability in soils and topography are key determinants of spatial heterogeneity in vegetation responses (Raynolds and others 2008;McManus and others 2012;Lara and others 2018). This conclusion is consistent with recent studies on spatial patterns of productivity trends across northwestern North America (Chen and others 2021) and on Banks Island (Campbell and others 2021), and complements previous research highlighting terrain variability driving these regional patterns (Jia and others 2006;Tape and others 2006;Walker and others 2009; McManus and others 2012; Tape and others 2012; Berner and others 2020; Niittynen and others 2020a). Our results using EVI trends to map tundra greening are comparable to those of studies using NDVI responses; suggesting similar spatial patterns of greening (Myers-Smith and others 2020) and mechanisms driving vegetation change (Jia and others 2006;Raynolds and others 2008;McManus and others 2012). Further, EVI uses additional spectral and angular information not used by NDVI, implemented to address known issues with NDVI related to solar incidence and atmospheric conditions present (Liu and Huete 1995).
Although some researchers recommend accounting for temporal autocorrelation in time series through pre-whitening methods (Guay and others 2014; Berner and others 2020), we do not believe that this is an issue in our analyses given the temporal revisit rate of the data used and the rates of phenological development present in tundra ecosystems. Weak evidence of temporal autocorrelation of NDVI has also been documented in areas dominated by deciduous vegetation compared to evergreen vegetation due to seasonal foliage replacement and other differences in reliance on previous-year nutrient storage (Berner and others 2011). Additionally, our methods are consistent with other studies using vegetation index time series that do not implement pre-whitening procedures (Fraser and others 2014b;Nitze and Grosse 2016;Raynolds and Walker 2016;Lara and others 2018).
The Tuktoyaktuk Coastlands and Yukon North Slope provide habitat to a diversity of mammals including caribou, muskox, bears, wolves, Dall's sheep, red fox, and wolverine (Russell and others 1993;YEWG 2004;Rickbeil and others 2018). Our observation that tussock and dwarf shrub tundra at lower elevations are most prone to increased vegetation productivity suggest that caribou who utilize this habitat type while avoiding upright shrub tundra (Russell and others 1993;Johnstone and others 2002;Rickbeil and others 2018) will be significantly impacted by ongoing vegetation change. The impacts of vegetation change on habitat use in this region should be assessed by combining landscape-scale data on vegetation change with satellite telemetry data for caribou and other important species (see: Rickbeil and others 2018). Range expansion of moose (Tape and others 2016) and beavers (Jung and others 2016) may also be related to vegetation change and should be explored using systematic surveys and satellite telemetry.
Several recent analyses suggest that coarse resolution imagery can mute the complexities of finerscale ecological patterns resulting in disagreement among coarse-scale remote sensing platforms (Guay and others 2014; Berner and others 2020; Myers-Smith and others 2020). Our findings show that landscape scale variation in biophysical characteristics strongly influences vegetation dynamics evident in moderate-resolution (30 m) Landsat imagery. Because tundra ecosystems exhibit heterogeneity at scales of 1-2 m (Epstein and others 2004b; Lantz and others 2010; Assmann and others 2020; Myers-Smith and others 2020), it is possible that change detection using higher-resolution sensors (such as 0.5 -3 m resolution imagery from WorldView, QuickBird, and PlanetScope) could account for some of the unexplained variation in our models. We suggest that future research utilize high-resolution sensors calibrated to temporally contemporaneous Landsat (Markham and Helder 2012;Belward and Skøien 2015) or Sentinel-2 (Drusch and others 2012) imagery. This integration of measurements from a greater variety of remote sensing platforms and the use of Sentinel-2 imagery cross-calibrated with Landsat offers a compelling advance in remote sensing capabilities (Wulder and others 2015; Claverie and others 2018) to avail upon samples of fine scale earth observations and improve models or offer insights on the nature of local vegetation heterogeneity.
Recent work using remotely piloted aircraft systems (RPAS; or drones) also emphasizes the influence of fine-scale variability in topography on the composition and productivity of tundra plant communities (Assmann and others 2020; Cunliffe and others 2020; Myers-Smith and others 2020). Increased affordability of RPAS systems has facilitated the collection of higher spatial resolution data representing increasingly large areas and future research can take advantage of these tools to monitor changes at fine scales (Fraser and others 2016; Assmann and others 2020; Cunliffe and others 2020). We encourage the use of RPAS surveying under standardized flight and recording protocols, such as through the High Latitude Drone Ecology Network (HiLDEN; https://arcticdrones.org), to address this gap in spatial data. These data can help further explain the complex interactions be-tween vegetation and environmental and climatic influences across spatial scales while covering greater area than is possible with ground surveys. Vegetation indices derived from satellite imagery (such as EVI or NDVI) may be influenced by factors altering surface reflectance like standing water or vegetative stress (Roy 1989;Ollinger 2011). Higher resolution imagery will make it possible to assess the influence of soil moisture, physiological stress, and disease on vegetation indices.

CONCLUSIONS
Global climate change is driving rapid increases in vegetation productivity across northern latitudes. The spatial heterogeneity in productivity trends highlighted in this study is the result of finer-scale, landscape processes that mediate the effects of regional climate warming. Surficial geology and topography are among the best predictors of spatial pattern in tundra greening because they influence soil conditions and moisture. Vegetation type also strongly influences changes in productivity because deciduous shrubs can respond rapidly to changes in moisture, nutrients, and temperature. Tundra vegetation change will impact wildlife habitat (Rickbeil and others 2018), surface energy balance, and carbon storage (McGuire and others 2006; Schaefer and others 2014) and understanding the drivers of landscape scale variation in vegetation shifts is critical to accurately characterize these relationships. We encourage continued research using random forests modelling to identify the role of regional processes as they relate to broader environmental change across spatial and temporal scales.
and WestGrid (www.westgrid.ca). We would also like to thank the Aklavik Hunters and Trapper Committee, Michelle Gruben, Dennis Arey, Kiyo Campbell, Tracey Proverbs, Angel Chen, Nicola Shipman, Zander Chila, and Hana Travers-Smith for their assistance in the field and for their valuable insight and discussions on these topics. We would also like to thank two anonymous reviewers whose insightful comments helped improve this manuscript.

De clarations
Conflict of interest The authors declare that they have no conflicts of interest.

OP E N A C C ES S
This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit h ttp://creativecommons.org/licenses/by/4.0/.