Plant functional type aboveground biomass change within Alaska and northwest Canada mapped using a 35-year satellite time series from 1985 to 2020

Changes in vegetation distribution are underway in Arctic and boreal regions due to climate warming and associated fire disturbance. These changes have wide ranging downstream impacts—affecting wildlife habitat, nutrient cycling, climate feedbacks and fire regimes. It is thus critical to understand where these changes are occurring and what types of vegetation are affected, and to quantify the magnitude of the changes. In this study, we mapped live aboveground biomass for five common plant functional types (PFTs; deciduous shrubs, evergreen shrubs, forbs, graminoids and lichens) within Alaska and northwest Canada, every five years from 1985 to 2020. We employed a multi-scale approach, scaling from field harvest data and unmanned aerial vehicle-based biomass predictions to produce wall-to-wall maps based on climatological, topographic, phenological and Landsat spectral predictors. We found deciduous shrub and graminoid biomass were predicted best among PFTs. Our time-series analyses show increases in deciduous (37%) and evergreen shrub (7%) biomass, and decreases in graminoid (14%) and lichen (13%) biomass over a study area of approximately 500 000 km2. Fire was an important driver of recent changes in the study area, with the largest changes in biomass associated with historic fire perimeters. Decreases in lichen and graminoid biomass often corresponded with increasing shrub biomass. These findings illustrate the driving trends in vegetation change within the Arctic/boreal region. Understanding these changes and the impacts they in turn will have on Arctic and boreal ecosystems will be critical to understanding the trajectory of climate change in the region.

These advancements are important for several reasons. First, identifying which PFTs are expanding is important for understanding current and future vegetation and climate dynamics. For example, shrub expansion is not limited to deciduous shrubs, but also includes evergreen shrubs [48]. Discriminating between deciduous and evergreen shrub expansion is important because the implications of expansion are notably different-expansion of shorter stature evergreen shrubs is less likely to exert a positive feedback on climate warming [49]. Furthermore, increasing vegetation cover is not limited to shrubs, but can also include forbs and graminoids, especially grasses [4,5]. Conversely, increases in PFTs such as deciduous shrubs can cause decreases in lichens and other PFTs [4,[50][51][52][53]. Since lichen are critical winter forage, this has detrimental impacts on caribou, which are in decline across much of the ABZ [50,54]. Downstream effects of climate change such as changes in herbivory [55] and fire regimes [30,56] also affect PFTs differentially, which highlights the importance of precisely mapping where and how PFTs are changing in a warming ABZ. Moreover, changes in vegetation do not occur solely in a two-dimensional plane. Rather, increased shrub height and density are critical aspects of shrub expansion [4][5][6]57]. The AGB of individual PFTs has been mapped at local scales (mapped area 0.25 m 2 -12.5 km 2 ) using field-based biomass harvests combined with imagery from aircraft and/or unmanned aerial vehicles (UAVs). These efforts typically utilized red-green-blue and multispectral imagery to classify scenes and then relied on either Structure from Motion or lidar data to create canopy height models and predict biomass [58][59][60][61][62]. At the regional scale, Landsat satellite imagery has been used with airborne lidar to map boreal tree AGB in Canada [63], and with biomass harvest data to map arctic shrub AGB across the Alaska North Slope [46]. However, there are no regional biomass maps in Alaska and northwest Canada that parse biomass by multiple PFTs.
We mapped live AGB of non-tree PFTs (AGB PFT ) within Alaska and northwest Canada from 1985 to 2020 using biomass harvest, UAV, and Landsat satellite datasets. Individual PFTs included deciduous shrubs, evergreen shrubs, forbs, graminoids, and lichens. This effort focused on non-forested vegetation communities, but the study area included boreal regions. In areas where trees were present, our maps represent understory biomass only. We bridged the spatial gap between field data and satellite imagery using UAV based biomass maps [43,64,65]. Specifically, we: (a) Compared accuracy of AGB PFT modeled using satellite predictor data when trained using field vs. UAV based estimates of AGB PFT . (b) Assessed spatial variation in AGB PFT across the region. (c) Evaluated changes in AGB PFT from 1985 to 2020. (d) Examined how fires influenced changes in AGB PFT in recent decades.

Study area
We mapped AGB PFT across approximately 500 000 km 2 in Alaska and northwest Canada, corresponding to the extended ranges of the Porcupine and Fortymile caribou herds (figure 1). These caribou herds are culturally, ecologically and economically important across our study area. Since AGB PFT determines wildlife habitat and forage, mapping AGB PFT across this area provides valuable insight for wildlife management. Vegetation communities in this region range from boreal forest to dense shrublands to wet tundra and lichen barrens (table S1).

Field and UAV data collection and aggregation
During the summers of 2018 and 2019, we collected field data at 44 sites across the study area (table S1, figure 1). These sites capture the range in nonforest vegetation community types that occur across the region. At each site, we harvested live AGB at five 0.25 m 2 quadrats placed at 20 m intervals along a 100 m transect oriented parallel to the topographic contour. Harvested biomass was sorted into PFTs (deciduous shrub, evergreen shrub, forb, graminoid, lichen), oven-dried to a constant weight, and weighed with a precision of 0.0001 g. We also collected UAV imagery at each of the 44 sites. This UAV imagery was subsequently used to locally map AGB PFT . This mapping was done by modeling pixel-wise AGB PFT from pixel-wise volume, derived from Structure from Motion based canopy height models. For additional details on field data and UAV-based mapping see Orndahl et al [66]. We created two training/validation datasets, one aggregated directly from field biomass harvest data ('field-based data,' n = 266) and one aggregated from UAV-based AGB PFT maps ('UAV-based data,' n = 427), to assess how model performance varied based on differences in scale from training data to final map resolution. For both datasets, observations were 30 m pixels with biomass predictions derived from either field or UAV data (figures S1 and S2). Data aggregation details are provided in supplementary text 1.

Landsat spectral predictors and environmental predictors
We processed remotely sensed predictors in Google Earth Engine (GEE) [67]. To produce time series data, we analyzed all available Landsat 4/5/7/8 Collection 1 Tier 1 surface reflectance data from April 1st to September 30th for 1984-2020. Imagery was masked to exclude clouds, cloud shadows, snow, surface water, and gaps using the quality flags provided by CFMASK [68]. We also masked pixels with Normalized Differenced Snow Index (NDSI) >0.1 to remove some water bodies and subpixel snow. We normalized Landsat 8 reflectance and NDVI to Landsat 7 following [69], and Landsat 4 and 5 NDVI to Landsat 7 [37].
After quality-screening the Landsat imagery, we generated seasonal spectral predictors using the Continuous Change Detection and Classification (CCDC) algorithm [70]. First, we applied the CCDC algorithm to model the reflectance bands, thermal band, and NDVI (table S2). Using these model fits, we generated annual synthetic reflectance values for specific seasonal percentiles. Seasonal percentiles were defined as percentiles of the snow-free season. For example, summer was represented by the day of year corresponding to the 50th percentile of snow-free days (table S3). The other indices listed in table S2 were calculated from the synthetic reflectance estimates. We also calculated the daily rate of change for each season, and the overall mean, median and amplitude (maximum minus minimum) for each band/index. This resulted in 320 seasonal spectral predictors, listed in table S3. The CCDC modeling process filled gaps caused by clouds or other interference including Landsat 7 Scan Line Corrector gaps, and produced a continuous time series for each pixel, except for temporal gaps between time segment breaks. Temporal gaps were filled by iteratively assigning spectral data from the closest available year. We supplemented spectral predictors with a suite of environmental predictors related to climate, topography, permafrost, and vegetation (table S3). These procedures for Landsat data pre-processing and derivation of seasonal spectral predictors follow recent advances in vegetation mapping [45,70,71].

Modeling
We modeled AGB PFT as a function of the predictors in table S3. Because our predictor set was large, we employed a multi-stage approach to reduce the pool of available predictors before model fitting. First, we excluded predictors for which the standard deviation across our training data was less than 50% of the standard deviation across the full study area, as a model trained on this data would be subject to extrapolation when applied over the full study area. Then, we removed highly correlated predictors using hierarchical clustering across a range of distance thresholds. This produced a list of predictor subsets, Using the reduced predictor sets and AGB PFT as the response variable, we fit linear mixed effects models with LASSO regularization [72] using the glmmLasso package [73][74][75] in R v4.0.2 [76]. Site was included as a random effect, and the response variable was transformed to reduce heteroskedasticity and non-normality of the residuals. We tested both log and square root transformations and corrected predictions using Duan's smearing estimate [77]. During the model fitting process, the LASSO regularization scheme used by the glmmLasso package shrunk predictors based on their importance to the model. Some predictors were shrunk to zero and eliminated from the model.
We built a separate glmmLasso model for each PFT and used nested cross-validation to perform model selection and assess model performance. An outer leave-one-site-out cross-validation was used to assess model accuracy. An inner 10-fold crossvalidation, grouped by site, was used to select the model parameter λ (supplementary text 3), which controls how heavily predictors are penalized during the LASSO regularization (figure 2, green box). The metric used to assess model accuracy within the outer cross-validation was the root-mean-squarederror (RMSE) of predictions compared to a held-out test group.
We performed this modeling process for each predictor subset and evaluated model performance as it related to number of predictors. This allowed us to select, for each PFT, a predictor set that maximized the number of predictors while minimizing model fit and convergence issues. For each PFT, we selected the predictor subset and response variable transformation that produced the best fitting model and used these parameters to fit a final model on the full dataset. To select the final model, we used a composite metric that averaged the absolute value of normalized root-mean-squared-error (nRMSE), the absolute value of normalized mean-bias-error (nMBE) and one minus the correlation between the observed and predicted AGB PFT values.
During final model fitting, we propagated uncertainty using Monte Carlo simulations (n = 100). We randomly permuted field-based and UAV-based estimates of AGB PFT to capture uncertainty in sampling variability (field-based) and the model used to create UAV-based estimates (UAV-based). We also sampled each dataset with replacement to capture parameter uncertainty within the final model fit (figure 2, blue box) [78,79]. We aimed to follow best practice guidelines for biomass product validation [80]. Further details about error propagation are provided in supplementary text 4.
To assess predictor importance, we standardized the coefficients (supplementary text 5) and averaged them across all 100 model fits. Predictors with higher average standardized coefficient values were deemed more important.

Mapping and validation
We created maps of AGB PFT for eight years over a 35 year time series: 1985,1990,1995,2000,2005,2010, 2015 and 2020. Using the model coefficients from each of the 100 Monte Carlo simulations, we predicted AGB PFT across the study area for each time step. For each time step, the result was 100 AGB PFT estimates for every pixel. We considered the best estimate to be the 50th percentile (median) and also derived a 95% confidence interval (CI) based on the 2.5th and 97.5th percentiles (figure 2, orange box).
To reduce overestimation for log transformed models the input predictors were clamped to the minimum/maximum value found in the training data −/+ 10%. Additionally, AGB PFT was recorded as zero if: (a) the PFT was predicted as zero cover by Macander et al, 2022 [45], or (b) the land cover type was predicted as 'barren' by Wang et al [41].
To assess how training data source impacted modeling we compared the distribution of field-based and UAV-based data, produced final maps using both data sources, and compared them separately during validation.
We first validated our models using leave-onesite-out cross-validation. Then, final map products were compared to an existing map of shrub AGB [46] and independent harvest datasets from within [47] (n = 11) and outside [81,82] (n = 30-38, depending on the PFT) our study area ('external validation' , table S4). Based on results from model validation, we determined which training dataset (i.e. field or UAV) was best suited to predict each PFT. To make this determination, we compared actual vs predicted values and used a composite metric that averaged the absolute value of nRMSE, the absolute value of nMBE and one minus the R 2 (table 1). In cases where one training dataset was not definitively better than the other, we visually inspected the maps to make the final determination. We then produced a final regional map for each PFT using the best model.

Assessing AGB PFT distribution and change over time
To explore the spatial distribution of AGB PFT across the study area, we aggregated our AGB PFT predictions by land cover type [41] for the year 2010. This allowed us to summarize AGB PFT within commonly used categorical land cover classifications. To explore how AGB PFT changed over time, we used our 35 year time series, with eight time steps, to chart trends in AGB PFT across the study area. We also examined impacts of wildfires on AGB PFT by tracking changes in AGB PFT within burned areas. We delineated burned areas using fire polygons from the Alaska Interagency Coordination Center [83] and the Canadian National Fire Database [84] and excluded fires <200 hectares in size.

Model/map validation, impact of training data source, and feature importance
The distribution of training data was similar between field-based and UAV-based datasets (figure S4). Both datasets exhibited right skewed distributions for all PFTs. Because of the way the data were collected and aggregated, the field-based dataset had a lower number of observations (supplementary text 1) Based on cross-validation, external validation and visual inspection of the final map products, we found forb and lichen AGB PFT were predicted slightly better by field-based models, whereas deciduous shrub, evergreen shrub and graminoid AGB PFT were predicted slightly better by UAV-based models (figures S5 and S6, table 1). However, deciduous shrubs were the only PFT for which one model type (UAV-based) performed better than the other across all metrics (table 1). Final maps were produced using the best available model ( figure 3).
To compare mapping results between PFTs, we normalized RMSE (nRMSE) and MBE (nMBE) using the mean of the observed data and compared nRMSE, nMBE and R 2 values amongst PFTS. Considering all accuracy metrics, we found model predictions were best for deciduous shrubs and graminoids, and worst  S8). Predictor importance varied by PFT. For deciduous shrubs, cover was the most important predictor (β = 0.59). The best evergreen shrub predictors were autumn Normalized Burn Ratio (NBR, β = 21.16) and Enhanced Vegetation Index (EVI, β = 6.97), likely because they differentiate evergreen shrubs from plants that lose their leaves in autumn. For forbs, SWIR amplitude was most important (β = 7.87), possibly because forbs (especially horsetails) are bright green early in summer but are absent or inconspicuous in spring and senesce completely in autumn. Graminoid cover was the top predictor of graminoid biomass (β = 18.81). For lichens, predictors that captured invariance in seasonal lichen color (green CC change from spring to early summer, β = −26.34) or differences in seasonal lichen visibility (NBR change from late summer to autumn, β = 28.22), were most important. This is likely because, unlike vascular plants, lichens do not change color from spring to summer, and lichen visibility increases in autumn as leaf drop reduces the amount of overtopping by shrubs. A full predictor importance comparison is provided in figure S9.   (table 2). Most biomass was contained in Woodland (26%) and Tall Shrub (24%) land cover types, and biomass density was highest in Tall Shrubs (table S5). Deciduous shrubs accounted for 67% of regional non-tree/moss AGB (table 2) and constituted the majority of non-tree/moss AGB within most land cover types (figure S10). Deciduous shrub AGB density was highest in Tall Shrubs (701.3 g m −2 ); evergreen AGB density was highest in Bogs (122.3 g m −2 ); forb AGB density was highest in Littoral areas (7.3 g m −2 ); graminoid AGB density was highest in Herbaceous types (27.6 g m −2 ) and Tussock Tundra (25.2 g m −2 ); and lichen AGB density was highest in Open Shrubs (104.0 g m −2 ) (figure 5).
Fires strongly influenced trends in AGB PFT across the study area. Deciduous shrub, evergreen shrub and lichen AGB decreased, on average, 51%, 69% and 79% the year after fire, whereas graminoid and forb AGB increased, on average, 80% and 174% (figures 7 and S12). The largest changes in AGB PFT over the study period were often located in historic fire perimeters (figure S13). Deciduous and evergreen shrub AGB increased the most in . Independent validation results from best available models. For each plant functional type, the best model results are shown. Field-based models are used for forbs and lichens, UAV-based models are used for deciduous shrubs, evergreen shrubs and graminoids. Model predictions were compared to independent field biomass harvest data. Shapes denote the region the validation data comes from. Data from the Yukon North Slope (circles) fell within the study area. Data from the Alaska North Slope (squares) and Seward Peninsula (triangles) were outside the study area, but within similar ecosystems/vegetation community types. Biomass harvest data were compared to predictions from the biomass mapping year closest to field data collection (max difference = 2 years, mean difference = 1.1 years). Harvested area covered, at most, 0.06% of total site area and was therefore not expected to have a detectable influence on the Landsat spectral signal in the event that the harvest occurred before the mapping year. Solid lines depict best-fit regression models, while dotted lines depict 1:1 relationships. older (21)(22)(23)(24)(25)(26)(27)(28)(29)(30) year old) fire perimeters (175%, 36%), whereas forbs and graminoid AGB increased the most in newer (≤10 year old) fire perimeters (246%, 63%). Lichen AGB decreased for all fire ages, by 22% on average (figure S13).

Map accuracy and uncertainty
Field and UAV-based biomass data distributions were similar, lending confidence the UAV data accurately captured in situ conditions. Field-based data were measured directly through harvest, and thus more accurately represent on-the-ground conditions. Field data were also required to calibrate the UAVbased models. However, field data were only collected at five 0.25 m 2 quadrats per site and thus provide sparse spatial coverage relative to the spatial resolution of the Landsat sensors (900 m 2 ). This scale mismatch becomes more problematic as site heterogeneity increases [85]. Conversely, UAV-based data provide high spatial resolution biomass estimates across the full extent of each site, and can thus be directly aggregated to match Landsat's spatial resolution rather than extrapolated from a few small quadrats. However, UAV-based data are subject to additional sources of error, as uncertainty is introduced during the modeling process [66,80]. We found field and UAV-based data performed similarly, with fieldbased models performing slightly better for forbs and lichens, and UAV-based models performing slightly better for deciduous shrubs, evergreen shrubs and graminoids (see table 1 for accuracy assessment results). UAV-based modeling added a layer of complexity, thus it is important to discern whether the added complexity improves modeling results [64,80]. We found deciduous shrubs were the only PFT for which UAV-based modeling produced consistently better results (see table 1 for accuracy assessment results). Deciduous shrubs were likely modeled best by UAVbased data because they often occur in the canopy and have more vertical structure than other PFTs besides trees. UAV AGB PFT predictions were produced using Structure from Motion technology, which cannot penetrate vegetation canopies [66]. UAV-based biomass estimation of shrubs might therefore be improved with canopy-penetrating lidar technology that can produce more reliable mapping of the ground surface [64]. Independent evaluation of biomass maps is difficult due to a lack of biomass harvest datasets partitioned by PFT. To independently validate our results, we therefore relied on some data collected outside our study area [81,82]. We found graminoids were predicted best, comparing well during crossvalidation and external validation. Deciduous shrubs performed well during cross-validation, but less well during external validation. This is likely because our data set had limited observations where shrub AGB was high. For example, we had only one observation with deciduous shrub AGB > 4000 g m −2 . Our model therefore tended to underpredict deciduous shrub biomass as compared to external validation data.
Our AGB PFT estimates were likely less accurate in more heavily forested areas because (a) our training/validation site locations were concentrated in non-forested areas and (b) tree canopy obscures understory vegetation in satellite imagery. However, we took steps to facilitate prediction of understory AGB PFT in areas with trees. First, tree cover was included as a predictor in the AGB PFT models. This enabled prediction of understory AGB by allowing the model to differentiate AGB PFT spectral signals in areas with and without trees. Second, we used PFT cover maps [45] to mask our AGB PFT maps. Pixels with zero percent cover for a particular PFT were assigned zero AGB for that PFT. This masked out AGB PFT in some areas with dense tree cover. We acknowledge that despite these efforts, AGB PFT predictions in forested areas are subject to error. Notably, tree biomass might be erroneously associated with shrubs and/or AGB PFT might be underestimated where it is obscured by tree canopy.
Estimating uncertainty in map products is crucial to understanding their accuracy across time and space, and to facilitate comparison with other map products. We used a Monte Carlo approach to propagate major sources of uncertainty through our modeling process. Ideally, uncertainty estimation would incorporate error in the reference data including the modeling process used to produce UAV-based AGB PFT estimates, while also incorporating error in the predictor variables, sampling variability, and uncertainty in the final model building process [80]. We did not include error/uncertainty in the remotely sensed predictor data due to computational limitations. We also likely underestimated uncertainty in UAV modeling, which involved several stages that each added uncertainty [66]. Uncertainty in scaling from field data to UAV imagery is likely to be an important source of error, therefore the use of UAV data as an intermediate modeling step should be carefully considered in the context of mapping objectives. It remains challenging to account for the myriad sources of uncertainty that arise when modeling regional biomass, especially when incorporating field, UAV, satellite, and environmental datasets.

Land cover analysis
Deciduous shrubs made up the majority (67%) of total plant AGB (excluding trees and mosses) across the study area in 2020. On the Alaska North Slope, previous research estimated shrubs accounted for 43% of total AGB from 2007 to 2016 [46]. Our estimate is likely higher because shrub dominance tends to be higher in areas with warmer summers [46,86]. For example, across the Alaska North Slope, Berner et al 2018 found shrub dominance increased from about 30% to 50% between areas with the lowest June temperatures (1.4 • C) and the warmest June temperatures (11.3 • C). Conversely, graminoids made up only a small percentage (2%) of total AGB across the study area, despite being the second most common PFT by cover [45]. Deciduous shrub AGB constituted the majority of total AGB for most land cover types-even types named for dominance of a different PFT (e.g. graminoids in Tussock Tundra). This highlights not only the important role of deciduous shrubs in regional plant communities, but also the importance of quantifying PFT biomass in addition to cover [87].

AGB PFT change over time
Our analysis showed a 31% increase in total AGB (excluding trees and mosses) from 1985 to 2020. This is surprisingly consistent with an estimated 32% increase in total AGB from 1982 to 2010 across Northern Alaska derived from AVHRR NDVI [88], despite vastly different spatial resolutions (30 m vs 12.5 km). Furthermore, we found increasing AGB was primarily driven by deciduous shrubs (37%). Similarly, there have been large increases in deciduous shrub cover in recent decades across Alaska and the Yukon Territory [45]. These increases are consistent with trends of warming-induced deciduous shrub expansion across much of the Low Arctic [4,5,89].
Our study documents large regional increases in tundra AGB linked primarily to deciduous shrubs, while also revealing more subtle shifts in the AGB of other PFTs. We found evergreen shrub and forb AGB increased to a lesser extent than deciduous shrubs (7%, 14%), whereas graminoid and lichen AGB decreased (14%, 13%). Meta-analyses of warming experiments across the Arctic report a standardized mean difference (unitless) of +11 for deciduous shrubs, +24 for low shrubs, +36 for tall shrubs, +5 for evergreen shrubs, +5 for graminoids, +1 for forbs and −17 for lichens [4]. The direction and magnitude of changes reported by this meta-analysis were consistent with our results, with the exception of graminoids, for which we reported a decrease. Graminoid decreases we documented could be due to overtopping from expanding shrubs which could both obscure graminoid cover in satellite imagery and cause on-the-ground declines via light competition [45,53].
Heterogeneity is inherent in patterns of Arctic greening and vegetation change. For example, greening can be caused by shifts in relative abundance of different PFTs, or by more uniform increases in growth of existing vegetation [2,10,15]. The shrub AGB increases we report likely capture both of these greening pathways. Within PFTs, species are likely to respond differently to warming [5,44,90]. Our maps aggregate vegetation at the PFT level, but advances in spectral and spatial resolution of remotely sensed imagery might facilitate mapping of species level vegetation change [44,91]. Disturbances such as wildfire can cause abrupt and patchy vegetation change. Finally, landscape and environmental conditions, particularly soil moisture, influence where vegetation shifts occur [2,10,15]. We found the largest increases in deciduous shrub AGB typically occurred in wetter areas along riverbanks and in drainages, as well as within 11-30 year old fire perimeters, which is consistent with prior research in northern Alaska and Canada [10,92,93]. Evergreen shrub AGB similarly increased in intermediate age fire perimeters, and on the Alaska North Slope, which aligns with strong trends of greening recorded for the Alaska coastal plain [94]. Decreases in graminoid, forb and lichen AGB often occurred in areas with increasing shrub AGB, suggesting these declines might be driven by competition [53]. Graminoid and lichen AGB decrease were both widespread across the study area. However, lichen declines were often located within fire perimeters, whereas graminoid declines were less closely tied to fire and were prevalent along the Alaska North Slope. Changes in graminoid and forb AGB were pronounced in the Yukon Flats, although this should be interpreted with caution as we did not have training or validation data in this area.

Fire impacts on AGB PFT
Fire had a significant impact on AGB PFT distribution and trends. In general, deciduous shrub, evergreen shrub and lichen biomass decreased the year following fire whereas graminoid and forb biomass increased. Deciduous shrubs recovered relatively quickly after fire, generally returning to pre-fire AGB levels 5 years after fire and exceeding pre-fire conditions as succession continued, in line with overall trends of increasing shrub AGB. These trends were consistent with tundra fire research done in Alaska that suggests vegetation recovers to pre-fire levels approximately 3 years after fire [95,96]. However, other research suggests recovery periods for deciduous shrubs are generally longer (>10 years) [97,98]. A longer time-frame for shrub recovery is suggested by figure S13, where deciduous and evergreen shrub AGB increase was greatest for intermediate age fires. Within 11-30 year old fire perimeters, shrub AGB increase regularly exceeded 100%, demonstrating the ability of fires to facilitate deciduous shrub and tree growth [95,96,[99][100][101][102]. Boreal forest fires burn more severely and take longer to recover. However, they too seem to promote growth of shrubs and are often characterized by an intermediate stage of deciduous shrub growth and dominance [99][100][101][102].
Lichens experienced the largest relative declines in AGB post-fire. Lichen AGB decrease within fire perimeters regularly exceeded 25%, which was nearly double that experienced across the full study area. Lichens were also slowest to recover, in several instances not reaching pre-fire levels, even 30 years after fire. This is consistent with research documenting susceptibility of lichens to combustion, and slow recovery of lichens after fire [50,103,104]. Fire induced declines in lichen have serious implications for wildlife as lichen are an important food source for caribou. Caribou have been known to avoid burned areas up to 50 years post-fire [50,105].
We found graminoid and forb AGB increased after fire, then declined to pre-fire levels as succession continued. Sedges and grasses are able to re-sprout from surviving roots and rhizomes, and thus recover quickly following fire [95,97,98,106]. Graminoids then show increased productivity as they are released from competition with more slowly regenerating shrubs [106]. These factors might explain why graminoid AGB increased the year after fire. Forbs also respond positively to fire [102,107,108], and horsetails are known to survive fire due to deep rhizoids [109].
Within fire perimeters, forb and graminoid AGB increased in newer fire perimeters, but not those of intermediate age, likely because these PFTs are some of the first colonizers of post-fire landscapes [102,107,108]. Areas of forb and graminoid decrease within intermediate age fires were often correlated with areas of shrub increase, suggesting these PFTs were outcompeted as shrub growth accelerated.

Conclusion
The maps we present move beyond two-dimensions and capture the composition and structure of vegetation across Alaska and northwest Canada over the past four decades. Mapping biomass by PFT makes it possible to characterize spatial and temporal patterns of vegetation change in the Arctic. We report increases in deciduous and evergreen shrub biomass that are consistent with greening trends taking place across much of the Arctic. Conversely, we observed decreases in lichens and graminoid biomass, which were often coincident with shrub increase. Fires were influential in driving spatial and temporal patterns of AGB. Our maps capture patterns of post-fire succession consistent with previous research, notably early succession driven by an increase in graminoids and forbs, initial decreases in deciduous and evergreen shrubs, and long-lasting decreases in lichens. We present an example of a multi-scale modeling approach that could be applied to other regions, with appropriate training/validation data inputs. Mapping vegetation biomass remains challenging, due in large part to the difficulty of acquiring biomass harvest data. As biomass datasets continue to grow, a centralized database of harmonized AGB would facilitate future mapping efforts. Advancements in lidar technology (satellite, airborne and UAV-based) are also likely to further biomass mapping efforts as lidar is able to penetrate vegetation canopies and better reconstruct ground surfaces and canopy heights.

Data availability statement
The data that support the findings of this study are available upon reasonable request from the authors. DOI: 10.18739/A2Q52FF2B.