Time-series maps reveal widespread change in plant functional type cover across Arctic and boreal Alaska and Yukon

Widespread changes in the distribution and abundance of plant functional types (PFTs) are occurring in Arctic and boreal ecosystems due to the intensification of disturbances, such as fire, and climate-driven vegetation dynamics, such as tundra shrub expansion. To understand how these changes affect boreal and tundra ecosystems, we need to first quantify change for multiple PFTs across recent years. While landscape patches are generally composed of a mixture of PFTs, most previous moderate resolution (30 m) remote sensing analyses have mapped vegetation distribution and change within land cover categories that are based on the dominant PFT; or else the continuous distribution of one or a few PFTs, but for a single point in time. Here we map a 35 year time-series (1985–2020) of top cover (TC) for seven PFTs across a 1.77 × 106 km2 study area in northern and central Alaska and northwestern Canada. We improve on previous methods of detecting vegetation change by modeling TC, a continuous measure of plant abundance. The PFTs collectively include all vascular plants within the study area as well as light macrolichens, a nonvascular class of high importance to caribou management. We identified net increases in deciduous shrubs (66 × 103 km2), evergreen shrubs (20 × 103 km2), broadleaf trees (17 × 103 km2), and conifer trees (16 × 103 km2), and net decreases in graminoids (−40 × 103 km2) and light macrolichens (−13 × 103 km2) over the full map area, with similar patterns across Arctic, oroarctic, and boreal bioclimatic zones. Model performance was assessed using spatially blocked, nested five-fold cross-validation with overall root mean square errors ranging from 8.3% to 19.0%. Most net change occurred as succession or plant expansion within areas undisturbed by recent fire, though PFT TC change also clearly resulted from fire disturbance. These maps have important applications for assessment of surface energy budgets, permafrost changes, nutrient cycling, and wildlife management and movement analysis.


Introduction
Large and accelerating changes in vegetation composition and structure have been documented in the Arctic, oroarctic, and boreal bioclimatic zones over the past 50 years (Esper and Schweingruber 2004, Tape et   landscape to regional scales, and trends in northern vegetation have also been examined at circumpolar scales using generalized spectral metrics such as the normalized difference vegetation index (e.g. Berner et al 2020). However, robust, multi-temporal quantitative modeling of the cover of specific plant functional types (PFTs) has not been undertaken over broad spatial extents, despite the utility of such products for application in physical science, ecosystem services, and wildlife management. Recently a moderate resolution (30 m) remote sensing image analysis was applied to map the distribution and change in land cover categories related to the dominant PFT (Wang et al 2020). Other efforts have applied similar data to map the continuous abundance of selected PFTs or other units (Beck et al 2011, Macander et al 2017, Nawrocki et al 2020, usually for a single time step. We improve on previous methods by modeling a continuous measure of plant abundance-top cover (TC, Karl et al 2017), the proportion of area for which a vegetation class forms the top of the canopy (0%-100%)-for seven ecologically informative PFTs through time. Our goal for these continuous measures is to provide the spatiotemporal resolution necessary to detect changes in the absolute magnitude of specific plant groups. PFTs observed or modeled in observational and experimental studies can then be more directly connected to any observed changes in abundance through time. Additionally, land managers can better manage current resources and plan for future conditions with these types of maps through improved understanding of changes in forage availability and habitat changes for animals, fuel for fires, or shifts towards or away from PFTs that include species important for human subsistence. Our quantitative maps for individual PFTs improve on traditional categorical vegetation maps, which suffer limitations for multi-temporal change detection because many vegetation shifts are likely to occur within the range of variability of a given map class rather than with a transition from one class to another.
We selected seven mutually exclusive PFT categories that collectively include all vascular plants within the study area as well as light macrolichens, a nonvascular class of high importance to caribou management (table 1). The PFTs, which largely follow Nelson et al (2022), are separated by growth form and leaf habit to optimize detectability in optical remote sensing and to characterize ecologically important distinctions related to vegetation dynamics and wildlife habitats. Dynamics and ecological feedback mechanisms for these PFTs have been extensively studied in recent decades. Conifer trees compose the primary fuel source for boreal forest fire and also are expanding in spatial distribution with increasing temperatures along the forest-tundra ecotone (e.g. Lantz et al 2019, Terskaia et al 2020, Maher et al 2021. Broadleaf trees form a smaller spatial footprint but are important in regeneration after fire. There is growing evidence that the abundance of broadleaf trees is increasing in Alaska and western Canada with intensified wildfire (e.g. Baltzer et al 2021, Mack et al 2021. Deciduous shrubs are expanding in the Arctic and oroarctic (e.g. Tape et al 2006, Myers-Smith et al 2011, which can lead to changes in snowpack, nutrient cycling, hydrology, and the abundance of lichens and other understory vegetation (e.g. Sturm et al 2001, Myers-Smith et al 2011, Mekonnen et al 2021. Notably, some deciduous shrub taxa, such as feltleaf willow, are key wildlife forage species (scientific names for all common plant names are provided in table S5 available online at stacks.iop.org/ERL/17/054042/mmedia). They often co-occur with trees in the boreal forest and can be dominant during mid-term post-fire succession. Tall, low, and dwarf deciduous shrubs are included because all can be canopy-forming and many Arctic taxa occur across growth forms depending on environmental factors. Evergreen shrubs are generally darker, occur in low or dwarf growth forms, and are green in the shoulder seasons adjacent to snow cover, allowing remote sensing time-series analyses to distinguish them from deciduous vegetation. Models have predicted that understory PFTs such as graminoids and light macrolichens could decline due to light competition from shrubs (Mekonnen et al 2018).

Study area
Our study area followed the southern boundary for North American Beringia from Nawrocki et al (2021), excluding temperate portions of Alaska (figure 1). We extended the map area eastward over most of the Yukon and small portions of British Columbia and Northwest Territories to include the full ranges of several Canadian caribou herds. The maps were gridded at 30 m resolution in the NASA Arctic Boreal Vulnerability Experiment (ABoVE) coordinate system (EPSG:102001; Loboda et al 2017).

Covariates
We compiled a suite of environmental and spectral covariates for the training and prediction of the timeseries PFT models. The environmental covariates represent topographic, climatic, permafrost, hydrographic, and phenological gradients across the study area; and were constant for all models. The spectral covariates are based on Landsat thematic mapper (TM), enhanced thematic mapper plus (ETM+), and operational land imager (OLI) data collected over 1984-2020 and follow the approach used to derive the ABoVE land cover product (Wang et al 2019). We used the continuous change detection and classification algorithm (CCDC; Zhu and Woodcock 2014) to derive a suite of annual seasonal reflectance and spectral metric covariates. Because PFT TC can change gradually over time, we used annual spectral covariates rather than aggregating over multiyear temporal segments, as Wang et al (2020) had done. We also extended the CCDC model to cover additional years (1984-2020, instead of 1984-2014) and incorporated Landsat 4 and Landsat 8 data. Detailed information on the definitions and development of the environmental and spectral covariates is provided in the supplemental materials.

Vegetation cover reference data
We compiled and normalized field-based vegetation cover data collected during 1994-2019 from a variety of sources in Alaska and Yukon, including resource management agencies, academic researchers, industry, and consultants. Data included both ocular and quantitative vegetation cover estimates collected from the ground or from aerial platforms (typically helicopters). Fractional cover derived from classified unmanned aerial systems (UASs) imagery were also included. Details are provided in the supplemental materials.
A subset of plots reported only total cover values, which includes overtopped understory vegetation. Total cover values are higher than TC values when vegetation layers overlap. TC can be directly observed by airborne or spaceborne optical sensors, while overtopped understory vegetation generally cannot be observed directly. To avoid bias from mixing different cover metrics in our modeling, we selected TC as a property directly observable by optical sensors, and normalized the total cover values to TC using an empirical model described in the supplemental materials.

Spatial autocorrelation assessment and blocking
Failure to address spatial autocorrelation in model calibration and assessment can lead to poorly specified models and overly optimistic assessments of model performance (Ploton et al 2020). However, prior Arctic and boreal continuous cover mapping efforts have generally not addressed spatial autocorrelation. To mitigate against these issues, we assessed the spatial autocorrelation of response variables before modeling and applied a spatially blocked cross-validation approach (Roberts et al 2017) to mitigate bias from spatial autocorrelation in model optimization and assessment. We calculated empirical semi-variograms to assess levels of autocorrelation for different PFTs at different lags (Ploton et al 2020). Through permutation of observations, a semi-variogram envelope was calculated to compare against observed semi-variance. We then assigned each observation to a spatial block at a scale larger than the autocorrelation sill. To assign spatial blocks, we calculated the X and Y coordinates for each observation (in meters, EPSG:102001), divided by the number of meters corresponding to the selected scale, and took the floor integer value. Each unique combination of X and Y integers was assigned to a spatial block.

Statistical models
We applied two stochastic gradient-boosting models, implemented in LightGBM (Ke et al 2017), to map PFT distributions based on the training data and spatial predictors. Stochastic gradient boosting is a form of machine learning where sequential weak learners, in our case short decision trees, fit the residuals remaining from prior learners; it is robust to outliers and collinear variables and incorporates interactions and nonlinear responses (Friedman 2002). A binary probability model was applied to map PFT distribution and a regression model was applied to map PFT abundance. The two models were combined for a final prediction of PFT cover. The statistical modeling approach follows Nawrocki et al (2020), Nawrocki et al (2021) except with five inner and outer crossvalidation folds that were spatially blocked. We predicted PFT cover at eight time steps spaced at five year intervals covering the period 1985-2020. Predictions for 1985,1990 and 2020 are extrapolated based on the models developed using reference data and covariates for 1994-2019.
2.6. Accuracy assessment 2.6.1. Pixel scale Model performance was assessed from the predictions in the independent test folds. We calculated mean absolute error (MAE), root mean square error (RMSE), bias, and R 2 from the observed and predicted cover values with reference data from all years pooled. These metrics were all assessed relative to the 1:1 line, rather than a regression line fit through observations and predictions. We also assessed accuracy and area under the receiver operating characteristic curve (AUC) for the presence/absence predictions. We assumed that model performance was similar for the extrapolated years, since the extrapolation range was short; we also calculated performance metrics for different epochs of reference data.

Uncertainty mapping
We assessed how RMSE varied with the magnitude of predicted TC and used this information to characterize the spatial distribution of uncertainty. We summarized the RMSE for each 1% bin of predicted cover for each PFT, based on the predictions in the test folds. We used the value for the absence bin (0% cover) directly and performed locally weighted smoothing through the range of values from 1% to the highest predicted cover. Map predictions above the highest predicted cover from the test folds were very rare but could occur, so we filled the RMSE estimate from the highest predicted value from the cross-validation folds to cover the full range to 100% TC. After linking RMSE to the magnitude of predicted cover, we were able to generate spatial maps of error.

Maps and analyses of PFT distribution and changes
We prepared TC and RMSE maps for each PFT at each time step. Since categorical land cover maps are a common and widely used tool to characterize the distribution and change of PFTs on the landscape (e.g. Wang et al 2020), we compared the information contained in the ABoVE land cover map from Wang et al (2019) to our PFT maps. We summarized the areal extent of PFTs within each land cover class for the year 2010.
We summarized and plotted the areal extent of all PFTs for each five year time step as the percent cover multiplied by the pixel area. Summary units included the full study area and three bioclimatic zones (Arctic, oroarctic, and boreal). The Arctic zone was derived from the Circumpolar Arctic Vegetation Map (Walker et al 2005); the oroarctic zone consists of northern, largely treeless highlands south of the Arctic where altitude (rather than latitude) controls climate and vegetation patterns and was derived from Virtanen et al (2016).
We assessed change over the 1985-2020 period using the model RMSE to stratify percentage changes into categories based on the direction and statistical significance of change. We assessed change significance using the unequal variances t-test (Welch 1947), where: We binned significance levels to three categories (p > 0.2, ⩽0.2 and >0.05, ⩽0.05) and summarized change for each PFT by the significance level, direction, and zone using weighted areas. To visualize changes, we summarized the net areas changed by PFT for 960 m superpixels (32 × 32 map pixels).

Fire disturbance effects
We stratified results by fire history to distinguish patterns of PFT abundance and change in areas that were unburned and that burned in different epochs before and during the time frame covered by our maps. We separately summarized PFT trajectories for the 20 largest fire years recorded in the study area based on the annual area burned to assess the impacts of those years and also to qualitatively assess the responsiveness of the five year analysis interval to annual disturbances (results in supplemental materials). For fire data, we used fire perimeters through the 2020 fire year derived from the Alaska Fire History Perimeter Polygons (https://fire.ak.blm.gov/predsvcs/maps.php) and Canadian National Fire Database (CNFDB, version dated 07-07-2021).

Spatial autocorrelation assessment and blocking
Empirical semivariograms for PFT TC (figure 2) indicate that spatial autocorrelation is significant out to a range of approximately 90 km or less for most PFTs. The shaded ribbons depict the variogram envelope, which is the range of semivariance based on randomly permuting the data values across the spatial locations 99 times. Where the semivariance levels out, and ideally crosses the lower envelope, spatial autocorrelation is insignificant. Semivariance for conifer trees intersects the envelope at a lag of 5-10 km but then decreases and remains close to but below the envelope until about 300 km. Evergreen shrubs remain consistently below the envelope across the 500 km range we assessed, which is approximately half the diameter of the study region. All other PFTs have clearly insignificant spatial autocorrelation within a range of 100 km. Based on this analysis, we assigned the vegetation plot reference data to 100 km blocks that were maintained in the five-fold cross-validation.

Accuracy assessment
Model performance, assessed using spatially blocked, nested fve-fold cross-validation, showed overall RMSE from 8.3%-10.7% for most PFTs (14.6% for graminoids, and 19.0% for deciduous shrubs; figure 3, table 2). R 2 was highest for trees, intermediate for shrubs, graminoids, and lichen, and low for forbs (table 2). The presence/absence models had the same pattern with accuracy and AUC highest for trees and lowest for forbs. There was little spatial structure in the model residuals (figure 2), evidence that the spatial blocking cross-validation approach was successful. Model performance metrics for different epochs of reference data are presented in supplemental table S6. The analysis of model error by predicted cover (figure 4) showed that RMSE was lowest for most PFTs at or near a predicted cover of zero. However, deciduous shrubs and graminoids had high RMSE at a predicted cover of zero and lower RMSE at predicted cover values of about 1%-25%. RMSE increased for all PFTs as cover increased from 1% to 40%-60%, and then declined at higher predicted TC values. In general, TC was underpredicted for high values across all PFTs. High values (at or near 100%) were rare in the reference data, especially for non-tree PFTs.

Maps and summaries of PFTs
All PFT top cover and RMSE maps are available at https://doi.org/10.3334/ORNLDAAC/2032 (Macander and Nelson 2022). Mapped light macrolichen cover for 2020 (figure 5) showed extensive areas with little or no lichen TC, large areas with moderate lichen TC, and areas of high lichen TC scattered throughout the study area but concentrated in southwest Alaska and interior uplands. Mapped deciduous shrub cover (figure 6) was nearly ubiquitous outside of high mountains and water bodies and was highest in south-central Alaska. Deciduous shrubs were the dominant PFT in the study area (table 3), followed by graminoids, evergreen shrubs, and conifer trees.
Conifer tree and broadleaf tree TC (figures S3 and S4) was mostly restricted to the boreal zone and was patchy, related to fire history. Evergreen shrub TC (figure S5) was very widespread, was highest in uplands of intermediate elevations, and was low in major river valleys. Graminoid TC (figure S6) was highest in the Arctic zone and was low at higher elevations and most of the boreal zone. Forb TC (figure S7) was low in most areas and was high on some river deltas and mountain slopes.

Comparison of continuous PFT maps with land cover categories
The ABoVE land cover classes (Wang et al 2019) are intended to capture the dominant PFT in each pixel. A cross-tabulation of the PFT cover maps for 2010 and the land cover map from the same year provides additional detail about the composition of the land cover classes in aggregate, in particular the sub-or co-dominant PFTs (figure 7). Land cover classes included large components of sub-and codominant PFT cover, showing that the land cover classes represent mixtures of PFTs. In addition to being the most abundant PFT overall, deciduous shrubs were the dominant PFT within several land cover classes associated with both forest and tundra ecosystems. While this analysis focuses on the weighted area distribution, the frequency distribution of PFT cover values within land cover classes provides additional information (figure S8).

Maps and analyses of PFT changes
The largest changes in PFT TC were associated with deciduous shrubs, which increased by 65 682 km 2 from 1985 to 2020 from 17.4% to 21.1% of the mapped area (table 3). Graminoids showed the second biggest change and the largest loss, decreasing by 39 850 km 2 from 11.9% to 9.7% of the mapped area. Conifer trees, broadleaf trees, and evergreen shrubs increased by over 15 000 km 2 each while lichens declined by 12 892 km 2 overall. The rates of change were mostly constant over the study period (figure 8), except between 2000 and 2005 in the oroarctic and boreal zones. That interval included the 2 largest fire years on record for the mapped area.
We summarized PFT weighted area by year within the six most extensive land cover strata from Wang et al (2020) (figure 9). These strata are based on aggregated land cover and include categories for both stable and changed land cover. The greatest relative change in PFT TC occurred in the changed land cover strata (Evergreen Forest Loss and Shrub Gain) and show the expected patterns. Deciduous shrub TC increase, and graminoid TC decrease, is observed in the four stable strata.
Combining the net change over 1985-2020 and the estimated RMSE at both timesteps, we estimated the overall magnitude and 30 m pixel scale statistical significance of the observed changes by bioclimatic zone (figure S9) and land cover strata (figures S10 and S11). For example, within stable evergreen forest, the PFT cover of conifer trees and deciduous shrubs showed large net increases and corresponding net decreases for evergreen shrubs, graminoids, and lichen. For the regional summaries, net change overall is the most useful summary, while the pixel-scale significance is relevant for detailed assessment of    to recent firescars. Increases in shrubs are widespread in all bioclimatic zones; increases of trees are widespread in the boreal and oroarctic. Graminoid TC decrease is very widespread, though some areas, such as the northernmost Arctic, show an increase. Lichen TC is mostly stable or declining. More conservative maps of change were prepared by masking out 30 m pixels where the change magnitude was less than the  overall model RMSE for the PFT, before aggregation (figure S14); the overall patterns are similar, but the spatial distribution of change is more limited, especially for shrubs, forbs, and graminoids.

Fire disturbance effects
PFT cover trajectories were summarized by fire history (figure 11) using 1 unburned category, 2 intervals that burned in 1985 or earlier (with no prefire PFT mapping), and five intervals that burned over the 1985-2020 period of the PFT mapping. The unburned category represents most of the area and shows linear increases in TC for tree and shrub PFTs and declines in lichen and especially graminoid TC. For burns that occurred during our mapping period (1986 or later), there are abrupt decreases in most PFTs (except forbs and graminoids) immediately following fire. Forb and graminoid TC increased after fire, then began to decline after 5-10 years. Shrub TC initially decreased, but rapidly increased after 5-10 years. Beginning about ten years post-fire, tree TC tends to increase as graminoids and forbs decrease. The two categories that burned in 1985 or earlier show different stages of post-fire succession at the start of our mapping period with the 1971-1985 firescars starting with higher graminoid TC and  to 1998, so much of the reflectance model data was interpolated between 1986 and 1999 with few reflectance observations. Therefore, our models did not always capture when fire occurred, and several years of post-fire recovery may have already occurred by the time suitable data became available. There was usually insufficient time to reinitialize the CCDC model after fires that occurred in 2020. Therefore, the cover estimates for 2020 were in most cases filled from prior maps (e.g. 2015). We note that a quality assurance raster associated with each map allows endusers to identify pixels where TC was filled from a different year.

PFT dynamics
Deciduous shrubs increased across the study domain in general and after fire in particular. Warminginduced increases in shrub recruitment and growth reported elsewhere connect with our finding of increasing deciduous shrub TC (Tape et al 2006, Myers-Smith et al 2011. We examined how the location of changes in deciduous shrub TC overlapped with changes in other PFTs ( figure 12). For all PFTs except graminoids, our models predicted extensive areas where deciduous shrub TC increased while TC of the other PFT remained stable (light green); these are mostly in the Arctic and oroarctic. Extensive areas of simultaneous decrease in graminoid TC and increase in deciduous shrub TC occurred (light orange), especially in the Arctic and in Southwest Alaska, where Arctic, oroarctic, and boreal bioclimatic zones converge. Simultaneous increase of both deciduous shrub TC and TC of the other PFT (dark green) was widespread in the boreal bioclimatic zone for woody PFTs (trees and shrubs), reflecting both post-fire recovery and woody plant expansion. Evergreen shrubs also increased across the study domain but to a far lesser extent and magnitude than deciduous shrubs. Deciduous shrub increases in the context of increases in other woody PFTs and decreases in herbaceous and nonvascular PFTs provide a clear indication of an ongoing long-term shift towards larger statured vegetation and away from tundra and graminoids in much of northwest North America.
The widespread loss of graminoid TC across the domain was an important finding and does not appear to be strongly associated with fire. Models projecting PFT changes due to warming (Mekonnen et al 2018) have connected graminoid decline to light competition with shrubs, which may explain the overall trend we observed. While some overtopped graminoids can persist in the understory, the strength of the decline in graminoid TC suggests a decline overall, not just at the top of the canopy. The simultaneous decline of graminoid TC and increase in shrub TC appears to overlap regions that have already lost or are losing underlying permafrost (SNAP 2022). Thus, declines in graminoid TC and its replacement by deciduous shrubs may have both causes and feedback cycles related to below-ground processes, such as permafrost degradation and nutrient cycling, in addition to surface properties, such as albedo and snow cover (Liston et al 2002, Sturm et al 2005. Lichens decreased across the study domain after fire, but most of the area with lichen TC loss corresponded to places with increased deciduous shrub Environ. Res. Lett. 17 (2022)  TC (figure 12). Competition by large statured plants, such as shrubs, has been shown in experiments to decrease lichen cover through time, and declines in lichen TC throughout the study area agree with observed and predicted outcomes of lichen loss due to increased shrub cover (e.g. Fraser et al 2014). Increased shrub TC in these areas indicate some loss of lichen may be related to competition through shading and smothering by leaf litter. Our results suggest that the magnitude and extent of lichen loss due to competition with shrubs may exceed the loss of lichen after fire. Future work will examine trends in lichen and other forage in relation to caribou population dynamics.
The landscape-scale to regional trends that our models predicted corresponded with some site-scale studies of multi-decadal vegetation change in Alaska (e.g. Racine et al 2004) and conflicted with others (e.g. lichens in Villarreal et al 2012). Similar to the findings of Roland et al (2019), discrepancies between our results and site-scale studies suggest that local vegetation changes do not necessarily indicate the broader landscape-scale or regional patterns of change. An important application of our maps is to provide landscape-scale and regional context for future analyses of local vegetation change. Comparison of our proportional PFT maps to the Wang et al (2019) land cover classes demonstrate that the proportional PFT maps represent greater complexity in vegetation change patterns, which can improve relevance to plot-driven metrics of local vegetation change.

Uncertainty
Model fit varied substantially between PFTs, leading to differences in the interpretability of resulting maps and in some cases limiting our ability to demonstrate significant change at the 30 m pixel scale. Bias was generally low, however, with all PFTs having a bias of ⩽1.9% except deciduous shrub (3.6%). This suggests that while pixel-level predictions have relatively large errors, predictions aggregated over larger areas (i.e. superpixels) are more reliable with little systematic offset between observed and predicted values.
Conifer and broadleaf tree TC were the best models and also the largest statured PFTs. Observed tree TC was also at or near 100% at some plots, which provided an endmember with unmixed spectral properties of that PFT that were generally lacking for lower-statured PFTs. In addition, the tree PFTs each consisted of a low number of species. Low species diversity likely both reduced the niche breadth relevant to abiotic gradients represented by our covariates and reduced spectral variability; we speculate that both factors contributed to higher model performance than for other PFTs.
The light macrolichen and graminoid models had intermediate relative performance, with deciduous and evergreen shrubs performing slightly worse. There is likely some confusion between dense deciduous shrubs and broadleaf trees, which can appear similar when broadleaf trees are young; future work should attempt to quantify that confusion. Lichen model error metrics MAE and RMSE reported here were similar to other studies in Alaska and western Canada (Kennedy et al 2020 while the R 2 reported here was somewhat lower. Forbs had the poorest fitting models, which was not surprising given that this PFT includes a large number of species that fill numerous and varied niches, are small, and generally occur at low TC values intermixed with other vegetation. Our models showed some difficulty distinguishing absences from very low cover values (e.g. 1%-3%). Thus, given that forbs often occur at low cover values, the reduced performance of the forb PFT may result from the lack of variation to model. Future PFT mapping efforts could consider lumping graminoids and forbs into a single herbaceous class or, alternatively, breaking forbs out into smaller units aligned to similarities in ecological niche and mapping distribution only.
The reference data included ground plots with both quantitative cover estimates and ocular estimates, aerial ocular estimates, and aggregated UAS classifications. Most of the ground plots were smaller than the 30 m pixels, and there are errors due to scale mismatch and observer variability. Estimates from the classified UAS data are based on full 30 m pixels and should have minimal scaling errors compared to small ground plots.
We reviewed the best practices recommended for validation of aboveground woody biomass product validation (Duncanson et al 2021) to assess approaches to characterize uncertainty in the current maps, and to consider approaches to more robustly characterize uncertainty in future PFT mapping refinements. Because our overall reference data were not from a probability sample, design-based or hybrid methods to estimate error were not possible, although we note some of the individual datasets included do have a probability sample design. There is also considerable error and observer variability in ocular vegetation cover estimates. The optimal approach would be a Monte Carlo simulation with error propagation that included error distributions in the reference data observations, the TC normalization model, covariates, and predictive model residuals. 30 m pixels where change was less than the magnitude of model RMSE were set to zero before aggregation, in order to emphasize areas with higher confidence of change.
This would be a major computational effort with the current models, which require several hours of compute time for each model training and annual map generation step.
For the current effort, we modeled RMSE error as a function of prediction magnitude. We observed the general pattern that error increased as predicted cover increased and then plateaued or slightly declined at the highest cover predictions (figure 4). For deciduous shrubs and graminoids, there were larger errors at zero TC than at low TC (e.g. 1%-10%). This was related to false negatives (i.e. predictions of absence when the PFT was present at low to moderate TC). These are represented in the binned scatterplots (figure 3) as a horizontal row of non-zero observed cover below the x-axis in each panel, which represents the zero predicted cover bin. Exploratory modeling using major shrub genera (e.g. willow, alder) did not show a similar effect, so less generalized PFTs may have narrower ecological niches such that their distribution is more readily modeled. Taller shrubs (⩾ ∼ 1.5 m) and trees can also be resolved using airborne lidar, though this covers very limited extents and does not provide taxonomic precision. More detailed modeling of deciduous shrub cover based on reference data that incorporates both lidar-derived structural information and taxonomic information from other sources will be explored in future work.

Conclusions
Maps of subpixel abundance and change for multiple PFTs for a large Arctic and boreal study area, including 6 vascular PFTs and light macrolichens, represent an important improvement from previous mapping efforts that focused on classifying the dominant PFT class, or mapping one or a few PFTs for a single snapshot in time. Despite extensive wildfire over the 1985-2020 study period, our results indicated that both broadleaf and coniferous tree cover increased throughout the domain. Rapid increases in deciduous shrub cover after fire and throughout the domain support both post-fire pulses in shrub abundance and long-term shrub increase and canopy dominance in the absence of fire disturbance. Graminoid TC decreased throughout the domain, likely driven by competition from deciduous and evergreen shrubs, evidence over a large area that supports predictions from process models (Mekonnen et al 2018). Slow, steady loss of lichen is related to both increases in shrubs and associated losses from fire.
Considerable challenges remain in resolving continuous cover of PFTs through time. Errors in our models stemmed from a wide range of sources. High resolution PFT maps derived from high resolution remote sensing imagery and hyperspectral imagery could provide extensive additional training data to further refine our maps in the future. Leveraging lidar to directly incorporate vegetation structure is also promising, especially to partition the deciduous shrub PFT into tall and low-stature groups. The continued refinement of quantitative PFT models will support a wide range of science applications and will provide valuable context for interpreting general signals of 'greening' and 'browning' evident in other long-term satellite records (e.g. Berner et al 2020, Potter and Alexander 2020). Our maps will enable PFTs observed or modeled in observational and experimental studies to be more directly connected to observed changes in abundance through time. Additionally, land managers can better assess current resources and plan for future conditions through understanding changes in forage availability and habitat changes for animals, fuel for fires, or shifts towards or away from PFTs that include species important for human subsistence. Our quantitative continuous cover maps for individual PFTs improve on traditional categorical vegetation maps, which suffer limitations for multi-temporal change detection because many vegetation shifts are likely to occur within the range of variability of a given map class, rather than with a transition from one class to another.

Data availability statement
The data that support the findings of this study are openly available at the following URL/DOI: https://doi.org/10.3334/ORNLDAAC/2032.