Earth's forests play a fundamental role in the global carbon (C) cycle and currently are a substantial carbon sink, sequestering up to 25% of human carbon dioxide emissions annually (Bonan 2008;Pan et al. 2011). Yet the future of forests in a rapidly changing climate is highly uncertain (Friedlingstein et al. 2014;Brodribb et al. 2020). In particular, increasing climate stresses and disturbance could compromise forest C storage, yielding manifold impacts on biodiversity, ecosystem services, and carbon cycle feedbacks and undermining the potential of forests as a climate solution (Hurteau et al. 2009;Seidl et al. 2017;Holland et al. 2019). For example, an unprecedented and climate-fueled bark beetle outbreak in Canada drove immense swaths of tree mortality and reversed an entire region of boreal forest from a C sink to a C source over a decade with large implications for climate policy (Kurz et al. 2008(Kurz et al. , 2009)). In addition to insect outbreaks, wildfires and climate stress have been widely documented as prominent risks because they strongly regulate forest C stocks and are likely to increase in future climates (Bentz et al. 2009;Hicke et al. 2012;Buotte et al. 2019;Wang et al. 2021). Thus, it is essential to rigorously quantify and understand drivers of historical risks and use this understanding to project future climate-driven risks for forest ecosystem functions and services, including longterm C storage (Clark et al. 2017;Anderegg et al. 2020).
Due to forests' role as a C sink and important co-benefits for biodiversity and ecosystem services, governments, corporations, and non-governmental organizations have shown widespread and growing interest in leveraging forests as "nature-based climate solutions" to sequester and store C as part of meeting climate policy goals (Griscom et al. 2017;Roe et al. effective use of forest-based climate solutions in an evidence-based climate policy framework.
Crucially, to be used for climate mitigation, forests must achieve some level of "permanence" whereby a management or policy action leads to more ecosystem C storage, averaged over time, compared to a rigorous baseline (Hurteau et al. 2009;Ruseva et al. 2017). Although fossil C emissions persist in the atmosphere for hundreds to thousands of years (Archer et al. 2009), many public and private carbon markets only require increased forest C storage to last at most for up to 100 years (Ruseva et al. 2017) -and sometimes only up to a few decades.
Rigorous forest climate risk assessment is crucial for understanding climate impacts on ecosystems and biodiversity, informing conservation and management prioritization, and guiding climate policies and programs relying on forest carbon uptake and storage. Continental-scale risk assessment is currently lacking and urgently needed (Buotte et al. 2019;Anderegg et al. 2020;Lecina-Diaz et al. 2021). Spatial quantification of risks can inform forest protocols in climate policy by ensuring that climate risks are adequately considered in program design -for example, through the construction of "buffer pools" and other insurance mechanisms -and can inform forest project development and conservation (Hurteau et al. 2013). However, current forest offset protocols tend to include fixed, spatially invariant risks that do not incorporate future climate impacts and likely underestimate the integrated risks to forests over long time scales (e.g. the 100 year horizon used by the Climate Action Reserve (Anderegg et al. 2020)).
In this paper, we combine forest inventory data across United States (US) forests, remotesensing data of wildfires, high resolution climate data, and downscaled climate model projections to assess climate-sensitive risks for forest C stocks in the US. We first quantify how forest structure and climate anomalies mediate major climate-related risks to US forests from wildfire and non-fire, climate stress-mediated tree mortality (defined here as tree mortality that is sensitive to climate, excluding fire-driven mortality). We then model the spatial patterns and magnitudes of these risks over the historical record. Finally, we use downscaled future climate data to project how these risks might evolve in the future due to climate change, revealing where forests are likely to be the most vulnerable in the 21st century.
We constructed statistical models of climate risks from fire, (non-fire) climate stressdriven tree mortality, and insect-driven tree mortality using high-resolution historical climate data, satellite data for fire burn area, and forest inventory plot data for tree mortality. We performed extensive cross-validation and comparisons against independent datasets over the historical period. We then developed a high-resolution downscaled climate dataset from six climate models to project these climate risks across the US for three future climate scenarios.
Statistical risk models and validation are described below, and full input dataset details, preprocessing, and climate downscaling are described in the Supplementary Information (SI) Methods.
We developed a statistical model to create gridded (4-km spatial resolution), monthly predictions of burn area as a function of climatic variables. This model built on previous fire risk estimation efforts (Barbero et al. 2014). Many of the methods are similar, though updated with more recent data (through 2018 rather than 2010).
Forest permanence risks -Manuscript -6
The model was fit to historical fire data from the Monitoring Trends in Burn Severity (MTBS) database (SI Methods), which comprises 30-m annual rasters of burn severity as well as burn area boundary polygons for individual fires (Eidenshink et al. 2007). The dataset covers fires from 1984-2018 and includes fires larger than 202 ha (404 ha in the Western US) for the continental US. As predictors, we considered both temporally-varying climatic variables as well as time-invariant vegetation variables. Our primary climatic variables were monthly temperature, precipitation, and climatic water deficit (CWD) derived from the TerraClimate dataset (SI Methods) for the historical period. For vegetation, we used the National Forest Type Dataset.
The forest groups ranged in area between 51k to 67M ha with the smallest eight groups each representing less than ~2M ha. To limit the number of variables and prevent overfitting, spatially sparse forest groups were aggregated into supersets by combining the smallest forest groups with the most spatially similar, larger forest group (we chose an area threshold of 1.76M ha). This consolidation decreased the number of forest groups from 25 to 17 and had little effect on model behavior. We converted these forest group maps into 17 binary, gridded maps, each of which became a predictor in the model. Every pixel was assigned to one forest group.
We fit a "hurdle" regression model to predict burn area as a function of climate and vegetation variables. This model jointly predicts the probability of a non-zero value and, if a non-zero value is present, its continuous value (Cragg 1971). Intuitively, this model can be thought of as combining a classifier ("was there fire?") and a regression ("if there was fire, how large was the burn area?"). For computational reasons, all datasets were aggregated to a 16 km 2 grid for fitting. The model was then applied to create predictions on the 4-km grid. We formally represented the hurdle model using a sequence of two generalized linear models: a Binomial model with logit link function predicting zero vs non-zero values, and a linear Gaussian model with normal link function predicting burn area in the locations where it was non-zero. We implemented the hurdle model in Python using scikit-learn by combining the LogisticRegression and LinearRegression methods (Pedregosa et al. 2012).
In addition to the variables described above, we included two additional predictors to better capture inter-annual trends. To create these predictors we first calculated two timeseries representing a conterminous US-average monthly temperature and precipitation and then calculated a 12-month rolling maximum of each of the two timeseries. Conceptually, these two predictors provide a measure of longer-term drought stress when conterminous US-wide high temperatures and sharp precipitation regimes occur simultaneously. In practice, including these extra predictors improved overall model performance only slightly, but allowed the model to better reproduce both monthly trends, interannual variability, and the observed increase in burn area over the observation period (see Fig S10-11).
For the full model, we assessed accuracy using area-under-the-ROC-curve (AUC) from the output of the logistic regression portion of the hurdle model. We consider this AUC the primary metric of interest given the sparse and nearly binary nature of the training data. We report these AUC values obtained using split-halves cross-validation, where the held out set was constructed by sampling years independently (Fig S1). We also assessed performance through the model's ability to reproduce three patterns: annual, seasonal, and spatial. For annual and seasonal trends, we computed an R 2 between the value computed directly from the data and the model's prediction allowing for a constant offset difference (Fig S1). For the spatial trends we computed an AUC just as we did for the full model, except for first averaging over time, which we consider the appropriate metric given the sparse nature of the data. While these are not the metrics on which the model is trained, they provide an indication of how well the model captures important and observable patterns in the data. For visualization purposes, predicted monthly burn areas from the model (fraction / month) were summed across months to estimate predicted burn area for each year (fraction / year).
We constructed 'climate stress' and 'insect' tree mortality models using data from the US Forest Inventory and Analysis (FIA) dataset, which is a nationwide standardized network of >100,000 long-term forest monitoring plots that track growth, mortality, and overall health of US forests. We used FIA data from 2000-2018. We aggregated FIA forest plot data on live and dead basal area from a tree-level to a 'condition' level, grouping together conditions representing repeated inventories of the same location. To construct climate stress and insect risk models, we screened for plots that had at least 2 or more inventory measurements, which enables the estimation of a true mortality rate. We next screened out plots that had a "fire" or "human" disturbance code or a "cutting" treatment code to remove major confounding disturbances.
We estimated the fraction of mortality based on the concept of a census interval, which we define as a pair of measurements in two measurement years (t₀, t₁). The fraction of mortality is defined as the ratio of new dead basal area in t₁ to the total live basal area in t₀, which was then normalized by the census length to give annual mortality rates. We computed this ratio separately for each condition. Given that many FIA plots only had one repeated measure (only one census interval), we used the first census interval for all conditions. We modeled 'climate stress-driven' mortality as the mortality that occurred during this census interval (with other confounding mortality drivers excluded, see above) and 'insect-driven' mortality using the "agent code" (AGENTCD) tree-level data, where codes of 10-19 indicate insects as the primary Forest permanence risks -Manuscript -9 causal agent of death. We note that the climate stress mortality models include mortality from insects, which was a deliberate decision because insects and climate stress such as drought often co-occur and interact to kill trees in many forests across the US and thus cannot be clearly separated (Anderegg et al. 2015), although we performed a sensitivity analysis of the climate stress model when excluding mortality of trees with insect agent codes and observed very similar modeled mortality patterns (R 2 =0.60, p<0.0001). Drought and other climate stress-driven mortality does not have a clear or widely used agent code in the FIA database; instead, climate/drought-driven mortality is often attributed to a wide array of more proximate agents (including insects, disease, weather, and other/unknown); see Anderegg et al. (2015) for a detailed discussion. Thus, our attribution is that this mortality is likely driven by "climate stress" broadly defined, as we have aimed to remove other major drivers of mortality, notably fire, human disturbance, and account for stand and self-thinning dynamics in model construction (see below). This mortality attribution approach has uncertainties but is generally reasonable and is the standard approach that has been widely used in numerous climate-related FIA studies (Shaw et al. 2005;Hember et al. 2017;Stanke et al. 2021).
We fit a statistical model predicting mortality as a function of climatic and stand variables. We formally represented the hurdle model using a sequence of two generalized linear models: a Binomial model with logit link function predicting zero vs non-zero values, and a linear model with beta-distributed link function, which is used for modeling proportions where values are between zero and one, for predicting mortality in the conditions where it was nonzero. The beta-distributed link function for the linear regression was chosen based on inspecting the behavior of the raw data distributions. We implemented the hurdle model in R using the glm function in the default 'stats' package and the 'betareg' package (Zeileis et al. 2010).
For each condition, we extracted the mean, minimum, and maximum over the census interval for six annual climate variables that were selected based on their importance in the drought and insect mortality literature: precipitation, temperature, Palmer Drought Severity Index (PDSI), potential evapotranspiration (PET), climatic water deficit (CWD), and vapor pressure deficit (VPD) (Bentz et al. 2010;Creeden et al. 2013;Williams et al. 2013). We also extracted the stand age for each condition from FIA and the community-weighted mean and range of the functional trait of the water potential at 50% loss of hydraulic conductivity (P50) from maps published in a recent study (Trugman et al. 2020) scaled to 0.25 degree. P50 has been widely linked to drought-driven mortality risk in site-level (Nardini et al. 2013) studies and meta-analyses (Anderegg et al. 2016). We also included in the mortality models two stand variables, age or age-squared, to account for background ecological dynamics such as selfthinning and background mortality, following Hember et al. (2017). All predictor variables were z-scored across the full dataset for that variable to ensure that variable ranges did not drive model outputs.
Climate stress and insect mortality models were fit independently to each FIA "forest type," which was chosen as an intermediate compromise of capturing the diversity of responses across US forests but aggregating above a species-level to enable adequate estimation of mortality levels. To ensure that each forest type had 50 or more condition measurements, we aggregated some sparse forest types into more common ones (59 were so aggregated out of the initial set of 171), leading to 112 initial forest types in our dataset. We aggregated conditionlevel mortality rates, age, climate data, and functional traits to a 0.25 degree grid for each forest type. This grid size was chosen through sensitivity analyses to determine the optimal aggregation where the coefficient of variation of mortality rate stabilized but large-scale climate variation was preserved. All climate stress and insect mortality models were fit using this 0.25 degree gridded data for each forest type.
We considered collinearity among predictor variables by examining variance inflation factors. We found that variance inflation factors were too high for comparing mean/min/max of the same variable (e.g. mean vs min vs max annual temperature), but were generally within acceptable levels (<5) across the six predictor climate variables, stand age, and P50 hydraulic trait. Thus, we conducted a stratified model selection analysis where we fit all possible model combinations with one of each predictor variable (i.e. varying all possible combinations of mean vs min vs max of each climate variable, age vs age-squared, P50 mean vs P50 range) and selected the most parsimonious model via Akaike Information Criterion (Burnham & Anderson 2004). For all analyses in this paper, we fit the same predictor variables across all forest types to reduce complexity. Thus, individual forest types were not allowed to have separate predictor variables. Model selection analyses were done separately on climate stress and insect mortality dependent variables.
We examined optimal model complexity by comparing the AIC and R 2 of nested sets of models. We compared climate stress and insect mortality models as a function of: i) a null model of forest type-only (i.e. each forest type would receive only its mean mortality), ii) a null model of mortality as a function of forest type and age only (i.e. no climate predictors), iii) mortality as a function of forest type, age, and climate predictors, and iv) mortality as a function of forest type, age, climate, and functional traits. We observed that climate variables significantly improved (i.e. deltaAIC << -3) model performance beyond both null models for both climate stress and insects and that the range of P50 significantly improved climate stress mortality models, but not insect models.
We assessed model performance with cross-validation and used two primary metrics that reflect the performance of different parts of the hurdle model. We first tested for spatial autocorrelation using Moran's I for each forest type and each of the climate stress and insect mortality models. For forest types and mortality models where significant autocorrelation was detected, we used a comparison of Moran's I by distance bin, using the 'correlog' function in the pgirmess package in R (Giraudoux et al. 2018), to determine the autocorrelation length. We set the spatial autocorrelation length as the midpoint between the last significant bin and the first non-significant bin. We then conducted spatial hold-out cross-validation (Ploton et al. 2020) 4), we identified forest types where cross-validated AUC was greater than 0.6 and the forest type had >20 grid cells with mortality observed in the historical record, based on a recent analysis of stability and information criteria in regression models (Jenkins & Quintana-Ascencio 2020). This led to risks being modeled with climate variables and projected for 30 forest types in the insect models and 61 forest types in the climate stress models out of 112 possible forest types. For all forest types that did not meet these criteria, we projected mortality simply as the mean of historical mortality for that forest type, and thus set all future mortality to that value. We note that this is a very conservative decision and is likely to underestimate future risks.
We further performed two evaluations of our mortality models against independent metrics or datasets. For the climate stress mortality models, we compared our observed mortality rates by species against the recent "Forest Stability Index" for eight major western US forest species (Stanke et al. 2020) and observed a strong relationship (R 2 =0.72; Fig S2). For the insect mortality models, we compared our observed and modeled spatial patterns to an independent dataset of Aerial Detection Surveys (ADS) done by the US Forest Service to map bark beetledriven mortality across the US (Williams et al. 2016). Despite large differences in types of dataset (e.g. aerial versus plot; "bark beetle-driven" mortality versus "insect-driven" plot agent codes) and spatial scales, we found strong agreement between our models and that independent dataset (AUC=0.79 and R 2 =0.29 comparing our modeled mortality and ADS observed mortality;
For future CMIP6 projection-driven risk models, we used the same climate variables as chosen by the final model selection analysis and projected climate stress and insect risk (i.e. % basal area killed per year) over each decade from 1950-2100 in different climate models and scenarios. Future decadal climate variables were z-scored against three decadal baselines (1990-1999; 2000-2009; 2010-2019) all conditions in FIA, regardless of number of censuses in the historical record, to cover all US forests, aggregated to 0.25 degrees by forest type, and then aggregated across forest types as described above. For future projections, we used a constant stand age and P50 functional trait based on the 2000-2018 historical values due to uncertainties about future forest dynamics and composition. This is an assumption and uncertainty, but a full exploration of stand age dynamics, species distribution and composition shifts, and demography is beyond the scope of this current analysis.
The fire risk model reliably predicted historical fires (cross-validated AUC: 0. where models showed skillful cross-validated performance (i.e. AUC>0.6) and thus lower risk in some regions (e.g. southern pine beetle risk in the southeastern US (Weed et al. 2013)) may reflect data and model limitations rather than inherently lower risks.
Our results provide a synthesis of fire-, climate stress-, and insect-driven climate risks to forests in an open-source dataset available at continental scales. Climate-sensitive risks to US forests have major impacts on forest C cycling and climate change feedbacks, and thus quantifying forest permanence risks is important for conservation and climate policy efforts.
Tree mortality and disturbance are large uncertainties in current land surface and vegetation models (Fisher et al. 2018;Pugh et al. 2019;Bugmann et al. 2019) and better large-scale historical datasets are needed for benchmarking and improving these models. Thus, the disturbance risk and mortality maps and their climate sensitivities derived here can help advance C cycle models. Our results reveal that US forests are very likely to experience increasing risks from climate change that undermine their C sequestration potential, an important factor that should be considered in climate change mitigation policy.
The spatial patterns in our risk models -both historical and future risks -broadly agree with other similar efforts of individual disturbances in the literature, such as the burn area patterns of large fires (e.g. Barbero et al. 2014) or drought risks (Buotte et al. 2019). The spatial patterns of insect model projections are consistent with previous projections for several major insect species (Bentz et al. 2010) and overall magnitude is similar to coarse-level ecoregion projections in parts of the western US (McNellis et al. 2021) (Fig 4). Further, the climate sensitivities of insect mortality for several western US pine species with the highest historical insect-driven mortality were consistent with estimates in the literature (Fig. S4) (Bentz et al. 2010;Creeden et al. 2014).
Our statistical modeling approach with static vegetation for estimating future climate risks to forests due to fire, climate stress, and insects is an informative 'first-step' analysis and complementary approach to current process-based methods for several reasons. First, accurately capturing mortality due to fire, climate stress and insects is quite challenging and often not included in state-of-the-art land surface models used in CMIP6 projections. Mortality of any type is a major uncertainty in process-based models, and efforts to predict drought mortality from first principles are nascent and still need substantial work (Bugmann et al. 2019). Prognostic insectdriven mortality is completely absent in CMIP6 mechanistic models currently (Fisher et al. 2018). Prognostic fire is present in some CMIP6 models, but none are able to capture the extent of current extreme fire events (Fisher et al. 2018). The statistical approach presented here is rigorously validated against historical observations and likely provides an upper bound of the extent of future disturbance given the lack of vegetation-disturbance feedbacks and dampening factors, further discussed among the several important caveats and limitations below. These models and approaches could be applied in other regions or countries by leveraging global fire Forest permanence risks -Manuscript -18 data from MODIS and forest loss/disturbance data from Landsat, bringing in ground plot networks where possible and accounting for direct human land-use change.
These climate-sensitive risk maps and projections provide spatial quantification and uncertainty assessment across climate models, climate scenarios, and risk models that can inform risk management and conservation decisions. To support these aims, all data and code underlying these models are publicly available, and can be easily accessed and visualized via a web portal (https://carbonplan.org/research/forest-risks). As with any analysis, these projections are subject to several uncertainties and caveats. In addition to uncertainties in underlying Earth system models and statistical climate downscaling approaches (SI Methods), these projections use empirical models based on static forest composition and structure over the 1984-2018 period.
Thus, these projections do not account for shifts in forest composition or distribution, interactions among risks, and carbon dioxide effects on plant drought stress. In particular, largescale impacts of fires, drought, or insects could substantially reduce biomass, and thus risk, although this is not likely to exert a material influence before the 2050s (Barbero et al. 2015;Abatzoglou et al. 2021). The risk projections also do not include impacts of land-use management, considered to be a strong potential lever in fire risk (Smith et al. 2016) and to a lesser degree climate stress and insect risks. We note, however, that there are also many reasons that these risk projections may be conservative or underestimates for insects and climate stress mortality, including projections made only for strong historical models, frequent non-linear impacts of drought and insects that may not be well-characterized in inventory data, and novel pests and pathogens. A comparison of our risk projections with mechanistic land surface models with prognostic fire from CMIP6 results revealed strong and consistent spatial correlations (Fig Forest permanence risks -Manuscript -19 S7), providing additional confidence in the patterns of future risks and their impact on forest carbon sequestration.
Our results clearly show both spatial heterogeneity and future increases in risk across broad swaths of the continental US. While some current forest offset protocols incorporate risk, for example through the construction of "buffer pools" or related insurance-like mechanisms, current risk estimates do not incorporate either form of variability (Anderegg et al. 2020). Thus, our findings raise serious questions about the integrity of these programs. Further work could incorporate observed heterogeneity and future projections to better inform the construction of these climate programs, such as by translating these risks into specific C loss estimates that could parameterize a better-grounded buffer pool or other insurance program. Our results provide a critical starting point in quantifying risks over space and time and can inform management, conservation, and policy actions. Taken in sum, our results increase the urgency and magnitude of response needed for reducing greenhouse gas emissions to mitigate climate change given the increasing risks of climate change to nature-based climate solutions.
Figure 1: Observed (left) and statistically modeled historical (right) risk maps for fire (A&B), non-fire, climate stress-driven tree mortality (C&D), and insect-driven tree mortality (E&F) reveal widespread and spatially varying risks. Fire risk is modeled as burn area by wildfires, i.e. fraction of a grid cell burned per year. Climate-stress and insect-driven tree mortality risk are modeled as basal area mortality per year. For each impact risk model (B, D, F), anywhere shaded is considered forested. The forest mask for fire differs slightly from those used for climate stress and insects due to different input data. Data gaps in forest inventory in WY and OK preclude observed risk calculations in climate stress and insect models (C, E). Figure 3: Risks for fire, climate stress-driven tree mortality, and insect-driven tree mortality (rows) averaged over the 2080-2099 period, separated by shared socioeconomic pathway (SSP) climate scenario (columns). Note that color-bars are substantially expanded relative to those in Figure 1 in order to visualize future projections that exceed historical risks.
Figure 4. Regionally-averaged time series of forest risks from three impacts: fire (top), non-fire, climate stress-driven tree mortality (middle), and insect-driven tree mortality (bottom). Regions of interest shown as gray boxes in the maps at the top of each column. As in Figure 2, statistical risk model simulations from each climate model simulation are shown as transparent lines, colored according to the three shared socioeconomic pathway (SSP) climate scenarios. The multi-model mean for each SSP is shown opaque. Statistical model projections for each region driven by historical meteorological data (rather than meteorological data derived from a climate model) are shown in black. Fire risks are calculated with a 10-year centered moving average, while climate stress-related and insect-related are presented as decadal averages. steps described in detail below. For all climate data processing, we used the years 1970-2000 (inclusive) as our climate normal period.
Quantile mapping is widely-used in the climate downscaling literature with a variety of implementations depending on the application (Wood et al. 2004;Li et al. 2010;Sobie et al. 2017;Maraun et al. 2013). For this study, we performed a simple form of single-variable trendpreserving quantile mapping (Wood et al. 2004;Stoner et al. 2013) to correct for persistent biases in the raw CMIP6 climate model output, matching the statistical distribution of the TerraClimate dataset on a gridcell by gridcell basis. Using the Scikit-downscale package (Hamman and Kent, 2020), the quantile mapping procedure was implemented for each variable and grid cell as follows. The process begins by extracting the three timeseries: X: historical (i.e. CMIP) y: observed (i.e. TerraClimate) X ': modeled (i.e. CMIP or ScenarioMIP) and by removing their trends. Empirical CDFs are then calculated for all three timeseries using Cunnane plotting positions. The CDF from X' is then interpolated onto the CDF of X using onedimensional piecewise linear interpolation, thus correcting for statistical differences between the X and X'. The mapping from X to y is performed by again using one-dimensional piecewise linear interpolation. Finally, the trend from X' is reimposed to the resulting vector of quantile mapped data.
Though this method is performed on monthly data, we still take care to handle extreme values in a reasonable way. As implemented in our study, we take the common approach of applying the proportional differences between X' and X for all values outside the bounds of the training data.
We used multiple hydroclimatic variables from TerraClimate in the development of the individual risk models (described below). Here we reimplemented the core components of the simple hydrologic model from Abatzoglou et al. (2018), the code for which is available here [https://github.com/carbonplan/cmip6-downscaling/blob/main/cmip6_downscaling/disagg/terraclimate.py]. We used the PDSI implementation of the climate_indices software package (Adams 2017). We used climatological wind from TerraClimate for derivation of hydrologic variables. Surface elevation and soil water holding capacity data was taken directly from the TerraClimate dataset.
As an additional form of cross-validation, we considered the ability of fire model to extrapolate beyond the period of which it was fit. For each of several years Y, we fit the model on years between 1984 and Y, and predicted on years between Y + 1 and 2018, and computed an R2 as above. Note that, because our calculation of R2 allows for an arbitrary offset, it can be high even when the absolute magnitude of the prediction is off (i.e. when there is a bias), which is especially likely when extrapolating model fits. For that reason, when computing these statistics, we separately compute and report the bias as the fractional difference between the average risk between prediction and observation (Fig S1 ).
We performed a sensitivity analysis to ensure that our climate-stress and insect-driven mortality models were robust to consideration of design weights that are assessed for FIA plots.
FIA plots are often assigned weights based on strata and estimation units for given inventories.
Following previous studies (61), we used the rFIA package (Stanke et al. 2020) to extract weights for all our conditions and constructed both climate stress and insect mortality models with and without weighting. We found very strong agreement between the two approaches in both modeled mortality patterns (R 2 =0.95, R 2 =0.98, respectively; p<<0.00001 for both) and model coefficients (R 2 =0.88, R 2 =0.90, respectively; p<<0.00001 for both), indicating that our approach was robust to weighting decisions (Fig S9). We further checked to ensure that mortality driven by invasive species was not driving our mortality model patterns by verifying that our dominant and high-sensitivity forest types were not those known to be widely impacted by invasive species.
We note that future projections of hydrological variables do not account for CO2-related water savings effects from declines in stomatal conductance, although our results are consistent with studies that do include these effects (Sperry et al. 2019;Qin et al. 2020). We verified that the mortality models predict reasonable mortality patterns and magnitudes when driven by
In addition, we further compared our risk model projections to results from process-based vegetation models that include fire from nine Earth system models included in the CMIP6 archive. We downloaded the surface upward mass flux of carbon dioxide expressed as carbon due to emission from fires excluding anthropogenic land use change ("fFire"; kgC m -2 sec -1 ) from the nine ESMs that have prognostic fire models and archived this variable for both We observed a strong and non-linear correlation (quadratic regression, r=0.68, p<0.00001) whereby larger increases in projected burn area were associated with dramatically smaller or even negative changes in tree carbon (Fig S7). While these two analyses are a simple, first-order comparison, they provide additional support that our statistical risk models 1) reasonably capture spatial patterns in fire risks from process-based vegetation models and 2) are tightly linked with projected changes in carbon in these mechanistic vegetation models. We further provide a figure comparing projected SSP3-7.0 risks from our risk models with current carbon stored in US forests (Fig S8).
We performed all analysis using open source software and publicly available data. Our team used the Pangeo cloud environment (Robinson et al. 2020)
The fitted model predicts patterns of historical fire with a cross-validated AUC of 0.89.
Comparing spatial patterns of observed burn area and model predictions for the historical period 1984 through 2018 reveals moderate agreement (Figure 1a) with a cross-validated spatial R 2 of 0.057. To support visualization of higher burn area regions, the same data is also shown with a -validation performance and shuffle statistics for the fire model. Each panel shows a single performance metric with two different forms of cross-validation (black dot, splithalves; gray dot, extrapolation) compared to a null distribution generated by shuffling the relationship between training and testing (thick gray line, median; thin gray lines, 10th and 90th percentiles). The shuffle distribution provides an indication of the performance that might be expected due to fitting noise. ROC/AUC scores of 0.5 are roughly random chance and 1 are perfect prediction (values >0.65 are typically considered strong performance). For the annual and seasonal R2 we allow for an arbitrary offset, so we separately computed the bias when crossvalidating these predictions (median 2.4% bias for split-halves and -52% for extrapolation). The substantial negative bias for extrapolation suggests that our projected estimates of fire risk are, if anything, conservative.
Figure S2: Comparison of forest type-level average mortality rates (% basal area killed / year) for major western US forest types from our data to a recently published "forest stability index" (FSI) 84 at a species level. Tree species are Pinus contorta (PICO; lodgepole pine forest type in this analysis), Picea engelmanii (PIEN; Engelmann spruce forest type), Populus tremuloides (POTR; aspen forest type), Abies lasiocarpa (ABLA; subalpine fir forest type), Pinus edulis (PIED; pinon-juniper woodland forest type), Juniperus osteosperma (JUOS; juniper woodland forest type), Pinus ponderosa (PIPO; ponderosa pine forest type), and Pseudotsuga menziesii (PSME; douglas fir forest type). Gray line shows a linear regression.
Figure S3: Comparison of FIA-derived observed (A) and modeled (B) insect-driven tree mortality (% basal area killed / year) to an independent bark-beetle tree mortality dataset (C) from aerial detection surveys (% of grid cell area affected / year) 71 .
Figure S4: Insect-driven tree mortality climate sensitivities are consistent with previous studies finding connections between temperature and drought stress and insect-driven tree mortality 36,85 . Model coefficients for the binomial (black) and beta (gray) portions of the hurdle model are shown for two example forest types, limber pine (A) and lodgepole pine (B).
Figure S5: Drought and insect model performance. Cross-validated area under receiver operator curve (AUC) for insect (A) and climate stress (B) models. ROC/AUC scores of 0.5 are roughly random chance and 1 are perfect prediction (values >0.65 are typically considered strong performance). Cross-validated R 2 values for grid cells with non-zero mortality for insect (C) and climate stress (D) models. We separately performed an aggregated spatial cross-validation at the continental scale, which yielded R 2 values of 0.18 (90% CI: 0.16 to 0.20) for drought models and 0.31 (90% CI: 0.23 to 0.34) for insect models. Figure S9: Mortality models are highly consistent across Forest Inventory and Analysis (FIA) stratum weighting methods. (A, B) Continental climate stress and insect coefficients from regression models with or without including FIA weights (R 2 =0.90, R 2 =0.88, p<<0.00001 for both). (C, D) Continental climate stress and insect mortality predictions from regression models with or without including FIA weights (R 2 =0.98, R 2 =0.95, p<<0.00001 for both).
Figure S10. Seasonal patterns of burn area from model predictions, averaged over the 1984-2018 historical period. The model captures the observed evolution of spatial burn area patterns throughout the year. For example, modeled burn area in the Southeastern US peaks in the Spring, followed by increases in the intermountain west in the summer. Modeled burn area in Autumn is predominantly within California.