Changes in forest composition and associated biophysical climate feedbacks across boreal North America


 Deciduous tree cover is expected to increase in North American boreal forests with climate warming and wildfire occurrence. This shift in composition can generate biophysical cooling effects via increased land surface albedo. Here we use newly derived maps of continuous tree canopy and fractional deciduous cover to assess change over recent decades. We find on average a small net decrease in deciduous fraction cover from 2000 to 2015 across boreal North America, and from 1992 to 2015 across Canada, despite extensive fire disturbance that locally increased deciduous vegetation. We further find a near-neutral net biophysical change in radiative forcing across the domain due to relatively small net changes in albedo. Thus, while there have been widespread changes in forest composition over the past several decades across the domain, the net changes in composition and associated post-fire radiative forcing have not yet induced systematic negative feedbacks to climate warming.

about half of Eurasian boreal forests are deciduous (whether dominated by conifer or broadleaf species), post-re composition shifts in the NA boreal biome represent one of the largest potential regional energy budget feedbacks to the global climate system [19][20][21][22] .
Here we assess changes in deciduous cover over the past two decades across the NA boreal forest based on Landsat satellite observations coupled with an extensive eld sample plot network. Moreover, we quantify the magnitude of albedo climate forcings associated with the observed composition changes.
As part of this assessment, we developed novel maps of sub-pixel fractional deciduous cover at 30m resolution across the entire NA boreal biome as well as coincident maps of tree canopy cover. We used these spatially detailed, multi-temporal maps to assess the in uence of forest compositional changes on albedo-driven radiative forcing, both inside and outside of burn perimeters for res that occurred between 1950 and 2018. We nd res exerted a strong in uence on forest composition that persisted for decades, but there was coincident succession of deciduous to evergreen cover, thus relatively small net biomescale change in composition and albedo-driven radiative forcing when aggregated across the domain. Consequently, there have not yet been widespread systematic negative biophysical feedbacks to recent climate warming across boreal North America, despite highly dynamic composition changes locally and regionally.

Continental-scale fractional deciduous cover
Earth-observing satellites such as those from the Landsat series 24 provide global observational data useful for analyzing vegetation dynamics at ne spatial scale and relatively high temporal resolution 25 .
We used Landsat satellite data combined with in-situ forest inventory samples from multiple site networks (n = 27,494 plots) to develop continental-scale deciduous fractional cover and tree canopy cover maps for the NA boreal forest at multiple time-integrated epochs: 2000, 2005, 2010, and 2015. For the Canadian portion of the domain, which has greater Landsat image coverage, we were able to extent our mapping back to 1992. We used extensive forest inventory data that represents a robust sample of forest species and canopy conditions across broad bioclimatic gradients spanning the NA boreal domain ( Fig. 1 and Supplementary Table 1). We determined deciduous fractional cover for each forest inventory plot based on the ratio of basal area of the deciduous trees compared to all trees in the plot. We then used random forest regression models 26 , trained with multi-seasonal Landsat spectral bands and indices, to predict the per-pixel deciduous fractional cover and its uncertainty from plot locations to the entire NA boreal domain. We also developed models to map tree canopy cover for the same epochs. These canopy cover maps represent the fractional areal coverage of tree canopy within each 30 m Landsat pixel. For subsequent analysis, we considered pixels with tree canopy cover fraction (proportion) greater than 0.25 (25%) as tree-dominated based on visual assessment of high-resolution satellite imagery at 110 locations ( Supplementary Fig. 1). The combination of per-pixel deciduous fractional cover and tree canopy cover maps provide a comprehensive view of changing gradients of forest composition across the NA boreal domain over multiple time-integrated epochs.

Deciduous forest fraction cover change
We assessed changes in forest composition across the boreal forest in NA from 2000 to 2015 and Canada from 1992 to 2015 using the deciduous fractional cover and tree canopy cover layers (Fig. 2). It was not possible to assesses pre-2000 changes in Alaska due to insu cient availability of Landsat imagery. We found the average per-pixel deciduous fraction in the NA boreal domain decreased from 0.33 in 2000 to 0.31 in 2015. Deciduous fraction increased across 171.4 million hectares (Mha; 40%) and decreased across 246.9 Mha (58%) of the 427.6 Mha boreal forest domain (Table 1). Similarly, the average per-pixel deciduous fraction across the Canadian boreal forest decreased from 0.26 in 1992 to 0.24 in 2015 (Table 2). Our results from the 1992 -2015 (Canada only) and 2000 -2015 (boreal NA) time periods collectively indicate that, on average, the relative dominance of deciduous trees in the NA boreal forest slightly decreased during the past two to three decades. These ndings contrast with past reported increases in the extent of deciduous forest during recent decades across a smaller portion of the NA boreal forest domain based on land cover type classi cations 17 .
There were large regional differences in forest compositional change from 2000 to 2015. Average perpixel deciduous fraction decreased nearly twice as much in Alaska (-0.058) than Western (-0.019) or Eastern Canada (-0.023) ( Table 3). Examining changes in deciduous fraction across ecoregions 27 , we nd the Taiga Shield West was the only ecoregion with an increase in average per-pixel deciduous fraction (+0.009) whereas the largest average per-pixel decreases in deciduous fraction occurred in the Taiga Cordillera (-0.057) and Boreal Cordillera (-0.055) ecoregions (Table 3).
To assess the role of res in driving compositional change, we examined changes in deciduous fraction inside burned area perimeters of varying ages (Fig. 3). We obtained burned area perimeters from Alaskan and Canadian forest re agencies 27,28 and focused on the period from 1950 onward because of potential uncertainties in pre-1950 burned area perimeters. For all res that occurred after 1950, 59.4 Mha showed increases in deciduous fractional cover and 65.8 Mha showed decreases in deciduous cover during 2000 -2015. Over the same period, tree canopy cover increased in 51.7 Mha and decreased in 54.3 Mha within re perimeters (Fig. 3). Fires increased per-pixel boreal forest deciduous fraction by an average of 0.09, but this was highly dependent on when the re occurred (Table 1; Fig. 3).
We then considered three different time periods of re history: (i) recent res (1998 -2018), (ii) intermediate aged res (1978 -1998), and (iii) older res (1950 -1978). The majority (75%) of the area with post-re deciduous fraction increases occurred within burn perimeters of recent to intermediate ages, primarily due to changes from pre-re evergreen conifers to post-re early-successional deciduous trees and shrubs. During mid-succession (intermediate aged res), roughly equal increases and decreases in deciduous fraction were accompanied by an overall increase in tree canopy cover (Fig. 3a,b). Older res tended to show a continued decline in deciduous fraction as evergreen conifers established dominance ( Fig. 3a,b). Overall, our results suggest that "relay oristics" has been widespread after re 20 , with deciduous fraction initially increasing after re but then progressively declining as evergreen conifers establish dominance.
Effect of composition changes on seasonal albedo and radiative forcing Forest composition played a major role in changing stand albedo-induced radiative ux. Cooling effects associated with deciduous trees tend to be relatively small in summer but more prominent in spring and fall due to increased snow exposure under deciduous canopies 20 . We determined median seasonal albedo using daily "blue-sky" albedo derived from the MODIS satellites 28 and then trained random forest models to predict seasonal albedo based on deciduous fraction and tree canopy as predictors. We used the differences between predicted 2000 and 2015 seasonal albedo composite layers to derive seasonal shortwave radiative forcing for spring, summer, and fall. We did not include winter albedo in this study because winter forcings are smaller compared to other seasons at these latitudes and, crucially, the quality of albedo data from remote sensing tends to be unreliable due to short daylengths in the winter months. We note, however, that re-induced changes in winter albedo and radiative forcing tend to be in the same direction as changes in spring, summer, and fall 20 .
We assessed how changes in deciduous fraction and tree canopy cover between epochs impacted the seasonal albedo and associated radiative forcing. Speci cally, we used radiative forcing kernels for albedo 28 from the Community Atmosphere Model 5 (CAM5) 29 in the Community Earth System Model (CESM) to calculate seasonal forcing by multiplying differences in albedo values by mean spatiallydownscaled albedo radiative kernels. In addition to domain-wide albedo forcing for tree-dominated pixels, we also assessed forest re-induced composition alterations using burned area perimeters. Our random forest models explained signi cant variability in surface albedo with cross-validated R 2 values of 0.36. 0.44, and 0.29 for spring, summer, and fall seasons, respectively with p < 0.05 for all three seasons.
Spring and fall season models had higher explanatory power from tree canopy cover, whereas deciduous fraction explained most of the variability in summer ( Supplementary Fig. 2), which is consistent with the phenomenon of underlying snow dominating the albedo signal in non-summer months and tree species composition dominating the signal in summer months 21 .
Across boreal NA, changes in deciduous fraction and tree cover canopy from 2000 to 2015 generated a small net biophysical warming of 0.004 W m -2 in spring, 0.004 W m -2 in summer, and 0.001 W m -2 in fall (average of 0.003 W m -2 across non-winter months). These overall changes were the net result of regional heterogeneity in deciduous and evergreen forest changes due to the interplay between re disturbance history, climate, and other factors (Fig. 4). Although albedo forcings from res were locally strong within re perimeters, negative forcings from recent res and positive forcings from older res largely cancelled each other out, and thus contributed relatively little to net domain-wide outcomes (Fig. 3f, Fig. 4b-d). Most recent res (1998 -2018) generated large negative shortwave radiative forcings across all the non-winter seasons due to increases in surface albedo between 2000 and 2015, whereas re perimeters of intermediate age (1979 -1998) generated mostly positive seasonal forcings (Fig. 3e). In midsuccessional forests, the overall net change in the deciduous canopy was relatively low but increasing tree cover generated a biophysical warming effect. Similar patterns were observed for older res (1950 -1978), albeit with lower magnitude due to slower changes in forest composition in older stands. The overall radiative forcing inside mapped re perimeters averaged -0.058 W m -2 for spring, -0.012 W m -2 for summer, and -0.036 W m -2 for fall, leading to a relatively small cooling effect across non-winter months (Fig. 3f).

Discussion
Implications for forest -re -climate feedbacks in boreal North America Our results indicate that despite extensive shifts in forest composition over recent decades, there was little net change in the overall proportion of deciduous and evergreen cover across the North American boreal forest, thus nearly neutral non-winter net radiative forcings over the same time period. Similarly, res since 1950 contributed relatively little to net changes in forest composition, with increasing deciduous cover and biophysical cooling from recent res largely being offset by decreasing deciduous cover and biophysical warming from intermediate aged and older res.
All three major geographic zones showed little change in overall deciduous fraction, which was unexpected given our current understanding of how recent changes in climate and re regimes affect forest composition. A large body of literature indicates deciduous species are more competitive than evergreen conifers in severely burned areas, a phenomenon exacerbated by the warmer and drier climate Alaska and Canada have experienced over the previous 20 -40 years 1, 2, 7, 8 . Permafrost thaw also leads to more rapid mineralization and increases in plant-available nutrients, which are expected to favor fastgrowing deciduous species over conifers 15 . As a result, several studies have predicted increases in deciduous cover during the 21st century in Alaska 14,15 and associated decreases in landscape ammability 10,20 . While such large-scale increases in deciduous trees may still occur, our results suggest such a transition has not yet produced a systematic domain-wide shift.
Within burned areas, our results show slow but consistent increases in evergreen dominance in older re scars counteract increases in deciduous cover from recent res. This phenomenon is consistent with increasing evergreen cover outside mapped re scars, as our analysis only included re perimeters as old as 1950 and gradual successional changes can occur even 100 years after res 30 . Decades of re suppression around an expanding wildland urban interface in Alaska and Canada may have also played a role in increasing evergreen cover 31,32 . Other potential factors may include the impacts of drought and insect outbreaks on deciduous tree species [33][34][35] , an inability for deciduous trees to regenerate in thick organic soils without severe res, and the role of jack (Pinus banksiana) and lodgepole (Pinus contorta) pine in post-re succession throughout much of Canada 36-38 .
Future changes in re frequency, timing and severity could alter these patterns and successional cycles 4,6,9 . For example, with res in the boreal domain projected to increase with warming climate, the successional pathways of post-re deciduous species such as Populus tremuloides (quaking aspen) being replaced by evergreen conifers such as Picea mariana (black spruce), Picea glauca (white spruce) and Pinus banksiana (Jack pine) during intermediate and older re ages may occur less frequently, thereby prolonging the persistence of deciduous forest cover. More frequent re-burns also progressively decrease organic soil depth and expose greater areas of mineral soil, and deciduous trees are able to rapidly regenerate vegetatively, thereby increasing the proportion of deciduous trees across the landscape. Post-re regional cooling is also expected to increase in the future with an intensifying re regime 18 , although earlier snow melt will likely decrease the magnitude of cooling 22 . This effect is expected to be most pronounced in areas with high evergreen forest cover and closed canopies. The fact that we did not, on balance, observe such changes across the domain during the past 20-30 years highlights the potential for other mechanisms in conifer-deciduous dominance, such as some conifer species (e.g. black spruce) being replaced post-re by other evergreen conifers (e.g. jack and lodgepole pine) more tolerant of drought conditions on mineral soils.

Methods
Deciduous fractional layers: We developed deciduous fraction layers using random forest regression models in a supervised-learning approach (Supplementary Fig. 3). We used Landsat satellite data 24 from Landsat 5 TM, Landsat 7 ETM+, and Landsat 8 OLI sensors to derive deciduous fractional cover and tree canopy cover layers at 30-meter spatial resolution. While the Landsat datasets are publicly available to download through the United States Geological Survey (USGS), we used Google Earth Engine 39 (GEE) to pre-process these datasets via GEE python application programming interface (API). We used Landsat 5, 7, and 8 Collection 1 surface re ectance (SR) datasets for years 1987-1997, 1998-2002, 2003-2007, 2008-2012 These image composites were prepared for early spring, mid-summer, and fall seasons to identify key differences in deciduous and evergreen green-up amplitudes (Supplementary Table 2; Supplementary  Fig. 4). In addition to Landsat surface re ectance bands, we added several vegetation indices derived from Landsat bands as additional features for the random forest regression model: 1) normalized difference vegetation index (NDVI), 2) normalized difference wetness index (NDWI), 3) soil-adjusted vegetation index (SAVI), 4) variable atmospherically-resistant index (VARI), and 5) enhanced vegetation index (EVI). With the three seasonal composites and addition of slope, elevation, and aspect as topographical variables from GMTED2010 dataset 40 , the image composite layer-stack consisted of 36 bands or features ( Supplementary Fig. 5). The GMTED dataset was used to exclude any discontinuities in the topographical data due to latitude, or other factors.
We obtained ground samples for modeling deciduous fraction from 1) the Cooperative Alaska Forest Inventory (CAFI) 41 , 2) the Bonanza Creek Long Term Ecological Research (BNZ LTER) samples in Alaska 42 , 3) repeatedly-measured permanent sample plots (PSPs) 43 from the Canadian provinces of Yukon, Northwest Territories, British Columbia, Alberta, Saskatchewan, Manitoba, and Ontario, 4) ground sample plots from Northwest Territories 44 , and 5) Canadian National Forest Inventory 45 plot data for sites across Canada ( Fig. 1 and Supplementary Table 1). These ground samples include specie-speci c basal areas at each site. We converted the specie-speci c basal areas to deciduous fraction by taking ratio of deciduous tree basal area to total tree basal area to be used as the response variable in random forest 26 regression models.
We used the RandomForestRegressor module in scikit-learn python package as our random forest regression model. We divided the sample locations geographically into zones: 1) east and 2) west ( Supplementary Fig. 6) due to differences in spatial density of ground samples, signi cantly different model outputs, and large variations in importance of input features in the random forest regression models. In addition to geographically distinct random forest models, we trained two temporally distinct random forest models for each zone: 1) a 3-season model, using spring, summer, and fall seasonal composites and 2) a 1-season model, using summer seasonal image composite. The 1-season model used summer season image composite as input and was applied to pixels where all three seasons spring, summer, and fall were not available resulting in gaps in the 3-season model output. To reduce bias in the random forest models due to data imbalance because of the large number of evergreen sample plots, we resampled each zone's sample set by binned under-sampling to create a near uniform distribution of the response variable -deciduous fraction ( Supplementary Fig. 7). The random forest models were then parameterized using a multi-dimensional grid-search approach with number of trees (n_estimators), maximum number of features per tree (max_features), minimum number of samples per split node per tree (min_samp_split), and minimum number of samples at each leaf node in a tree (min_samp_leaf) as parameters (Supplementary Table 3). For each combination of the parameters, we parameterized 1,000 5fold cross-validated model runs up to a total of 30 million models to perform a grid search for parameter sets resulting in the highest mean cross-validated r-squared value, low root mean-squared error (RMSE) and low standard deviation in RMSE (Supplementary Table 3). The random forest regression models were then applied to tiled Landsat image composites. Landsat derived seasonal composites bands that are highly correlated (absolute value of correlation coe cient > 0.95; Supplementary Fig. 8) were not used in random forest models as inputs. We used variable importance plots to assess and validate bands used in the random forest models (Supplementary Fig. 9). We also used one standard deviation of all samples in output leaves across all the decision trees in the random forest regression model as a measure of uncertainty for deciduous fractional cover layers. Additionally, the samples not used in training the model were used to provide an estimate of variance explained by the RF model using adjusted r 2 value and model performance using root mean squared error (RMSE) value (Supplementary Table 3). Spatial autocorrelation effects in validation samples were reduced by resampling all the validation plot locations via gridding and using a minimum distance threshold of 1 km.

Tree canopy cover layers
The tree canopy cover layers were developed using a random forest modelling framework similar to deciduous fraction (Supplementary Fig. 3). The training sample were selected in a strati ed random manner, with a stratum for each tree canopy cover level and were extracted from the Global tree cover data for 2010 25 . We then used the tree canopy cover sample locations to extract band values from the 3season 36 band image layer-stack. Similar to deciduous fraction modelling approach, we used the RandomForestRegressor module in scikit-learn python package as our random forest regression model to derive tree canopy cover layers. We also used a 1-season model to ll in the gaps in output of the 3season model. However, we did not divide the study area into different zones for tree canopy cover. We used 30% of the ground samples that were selected to build the models exclusively for model validation for tree canopy cover layers. The parameterization of the model was performed similar to the deciduous fractional layer approach (Supplementary Table 3). To determine uncertainty in tree canopy cover value we again used one standard deviation of all the samples in output leaves across all the decision trees in the random forest regression model.
Assessing change and albedo-based radiative forcing: We used spring (March -May), summer (June -August), and fall (September -October) season image-composites of the blue-sky albedo product for 2000, 2005, 2010, and 2015 to assess changes in albedo along with deciduous fraction and tree canopy cover. In each seasonal mean-value composite, we used solar zenith angle > 70 deg in valid data values. Areas with low post-re tree cover were included in our analyses because they had higher pre-re tree cover values. To examine the strength of relationship between albedo, deciduous fractional cover, and tree canopy cover, we modelled albedo for spring (March -May), summer (June -August), and fall (September -October) using albedo data from 2000, 2005, and 2010 in a random forest supervisedlearning approach (Supplementary Fig. 3). We extracted strati ed random samples within each 0.1 albedo interval to train the random forest regression model from the blue-sky albedo product seasonal composites. We parameterized the random forest model by optimizing r-squared and RMSE values using a multi-dimensional grid-search approach to determine the best set of random forest regression parameters: n_estimators, max_features, min_samp_split, and min_samp_leaf (Supplementary Table 4). We multiplied the difference in our modeled albedo between the two epochs 2015 and 2000 by the "allsky" CESM-CAM5 albedo kernel 28 to compute radiative forcing. We extracted mean values of positive and negative change in deciduous fraction, tree canopy cover, blue-sky albedo composites, and radiative forcing to assess changes inside re perimeters and across the entire domain. We then used re boundaries obtained from Alaska Large Fire Database 46 and Canadian National Fire Database 47 for res that occurred between 1950 and 2018 to assess change. We used > 500,000 strati ed random samples to assess the validity of deciduous canopy layers by examining the relationships between deciduous fractional cover, tree canopy cover, and uncertainties ( Supplementary Fig. 10).

Caveats
Although our deciduous fractional cover and tree canopy layers capture the variation in tree-based forest composition across the boreal domain ( Supplementary Fig. 11), they do not speci cally consider the understory component of the boreal forest. In our analysis, we included trees that are above the height of 1.5 m. This excluded understory vegetation component may contribute to the unexplained variance in the deciduous fraction regression model. Another limitation of this analysis is the dependency on multispectral data to extrapolate deciduous fraction and tree canopy cover across the boreal domain. Since we trained the regression model on the spectral re ectance of tree-based deciduous fraction, there may be locations where the spectral re ectance of shrubs and other vegetation types resembles that of the treed components. While we used tree canopy cover as a masking layer to only include treed forest in our study, locations with high shrub density may be modelled as having higher tree canopy cover than the threshold, making them di cult to exclude. Further, per-pixel deciduous fraction does not categorically depict vegetation cover. While 30 m spatial resolution is ne enough to represent homogeneous vegetation, this pixel size may include mixed vegetation types such as deciduous shrubs and evergreen trees. Such mixing can occur for mid to lower tree canopy cover resulting in some pixels having a deciduous fraction value different from forest stands on the ground. Lastly, due to lack of remote sensing data for an entire re cycle, which is typically of the order of 80-150 years 19 , we analyzed variations in older re scars across a time-period of 15 years to assess changes in forest composition and forcing in early to midsuccessional post-re periods. We did not analyze res in late-successional period in this study due to lack of reliable re perimeter data from res that occurred 80-150 years ago. However, our assessment of changes in forest composition and forcing for the entire North American boreal domain provides an indication of variations in such late-successional re perimeters. We also were not able to integrate albedo and associated radiative forcing estimates for winter months. However, we note the direction of albedo changes and forcing in winter tends to follow that in spring and fall 20 , and, as such, we do not expect this to alter our conclusions. Finally, we were unable to account for potentially earlier snowmelt in spring due to warming during the analysis periods, which would otherwise lead to more biophysical warming 21 .

Software
We used Python 3.7.5 programming language for all work ows in this analysis. Data extraction from the raster layers was performed using custom packages and codes written in python programming language.
All statistical and machine learning analyses were performed using python packages numpy 1.15.0, scipy 1.2.0, and scikit-learn 0.23.2.

Data Availability
The data sets generated by this study will be made publicly available through the Oak Ridge National Laboratory Distributed Active Archive System (ORNL DAAC).

Code Availability
We developed custom python packages geosoup and geosoupML to process, extract, and manipulate geospatial data on multiple compute platforms and to use random forest and multiple linear regressor to extrapolate deciduous fraction and tree canopy cover from point samples over the entire boreal domain in our work ow and codes. The custom python packages are available at: https://github.com/masseyr/geosoup and https://github.com/masseyr/geosoupML, respectively. We also developed python package eehelper primarily to extract and manipulate geospatial data using GEE. This python package is available at: https://github.com/masseyr/eehelper. The python codes and bash scripts used in our work are available at: https://github.com/GEODE-Lab/DeciduousFraction.

Variables
Fire perimeters All re perimeters