Future climate risks from stress, insects and fire across US forests.

Forests are currently a substantial carbon sink globally. Many climate change mitigation strategies leverage forest preservation and expansion, but rely on forests storing carbon for decades to centuries. Yet climate-driven disturbances pose critical risks to the long-term stability of forest carbon. We quantify the climate drivers that influence wildfire and climate stress-driven tree mortality, including a separate insect-driven tree mortality, for the contiguous United States for current (1984-2018) and project these future disturbance risks over the 21st century. We find that current risks are widespread and projected to increase across different emissions scenarios by a factor of >4 for fire and >1.3 for climate-stress mortality. These forest disturbance risks highlight pervasive climate-sensitive disturbance impacts on US forests and raise questions about the risk management approach taken by forest carbon offset policies. Our results provide US-wide risk maps of key climate-sensitive disturbances for improving carbon cycle modeling, conservation and climate policy.

INTRODUCTION

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 (Brodribb et al., 2020; Friedlingstein et al., 2014). 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 (Holland et al., 2019; Hurteau et al., 2009; Seidl et al., 2017). For example, an unprecedented and climate‐fuelled 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, Dymond, et al., 2008; Kurz, Stinson, et al., 2008). 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., 2010; Buotte et al., 2019; Hicke et al., 2012; 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 long‐term C storage (Anderegg et al., 2020; Clark et al., 2016).

Due to forests’ role as a C sink and important co‐benefits for biodiversity and ecosystem services, governments, corporations and non‐governmental organisations 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 (Cook‐Patton et al., 2020; Griscom et al., 2017; Roe et al., 2019). Yet significant scientific gaps remain that greatly limit the 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 prioritisation, and guiding climate policies and programmes relying on forest carbon uptake and storage. Continental‐scale risk assessment is currently lacking and urgently needed (Anderegg et al., 2020; Buotte et al., 2019; 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 programme 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, remote‐sensing 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.

METHODS

Overview and climate data

We constructed statistical models of climate risks from fire, (non‐fire) climate stress‐driven 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, pre‐processing, and climate downscaling are described in the Supplementary Information (SI) Methods.

Risk models

Fire

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, although updated with more recent data (through 2018 rather than 2010).

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 67 M ha with the smallest eight groups each representing <~2 M 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.76 M 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 km2 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 generalised linear models: a Binomial model with logit link function predicting zero versus 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 Figures S10 and S11).

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 (Figure 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 (Figure 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 visualisation purposes, predicted monthly burn areas from the model (fraction/month) were summed across months to estimate predicted burn area for each year (fraction/year).

Climate stress‐ and insects‐driven tree mortality

We constructed ‘climate stress’ and ‘insect’ tree mortality models using data from the US Forest Inventory and Analysis (FIA) dataset, which is a nationwide standardised 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 to 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 normalised 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 modelled ‘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 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 modelled 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 (Hember et al., 2017; Shaw et al., 2005; 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 generalised linear models: a Binomial model with logit link function predicting zero versus 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 non‐zero. 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., 2014; 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 self‐thinning 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 condition‐level 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 stabilised 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 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. deltaic << −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) for each forest type and mortality model, whereby one grid cell was held out from model training and a spatial buffer around that grid cell equal to the autocorrelation length was also removed from model training. The model was then fit on the remaining data and used to predict the hold‐out grid cell, and this was repeated 1000‐fold for each forest type and mortality model. Similar to the fire model, we examined model performance using cross‐validated area under the receiver operating curve (AUC) for the binary component of the hurdle model. We also considered the non‐zero‐value R for the beta‐regression part of the hurdle model (Figure S5), aggregating as above.

Finally, we imposed one further set of criteria on climate stress and insect models to incorporate climate‐dependence only where justified based on model performance. For all final model‐based analyses (i.e. Figures 1d, f, 3, 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 modelled 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.

FIGURE 1

Observed (left) and statistically modelled 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 modelled as burn area by wildfires, that is, fraction of a grid cell burned per year. Climate‐stress and insect‐driven tree mortality risk are modelled 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 estimates in climate stress and insect cases (c, e)

FIGURE 2

Projected 21st century risks for fire (a), non‐fire, climate stress‐driven tree mortality (b), and insect‐driven tree mortality (c) averaged across the US. Statistical risk models forced by simulations from each different climate model shown as transparent lines, coloured according to the three shared socioeconomic pathway (SSP) climate scenarios. The multi‐model mean for each SSP is shown opaque. Statistical model projection 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

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 visualise 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, coloured 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

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; Figure S2). For the insect mortality models, we compared our observed and modelled spatial patterns to an independent dataset of Aerial Detection Surveys (ADS) done by the US Forest Service to map bark beetle‐driven 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 modelled mortality and ADS observed mortality; Figure S3).

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) and the ensemble mean was taken across these baselines for each climate model and decade. All modelled and projection maps (e.g. Figures 1d/f, 3) were made on all conditions in FIA, regardless of number of censuses in the historical record, to cover all US forests, aggregated to 0.25° 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.

RESULTS

The fire risk model reliably predicted historical fires (cross‐validated AUC: 0.89), capturing interannual variability (cross‐validated R 2: 0.64), seasonal patterns (e.g. spring risk in the southeastern US and fall risk in the western US; cross‐validated R 2: 0.90), and spatial patterns (cross‐validated AUC: 0.78) (Figure 1a‐b; Figure S1). The model captured the spatial patterns of more prevalent fire across the western US, in particular in California and the northern Rocky Mountains.

Historical patterns of climate stress‐driven tree mortality, which is predominantly drought stress and includes biotic agents/insects, were highest across the western US and intermountain West, which was captured in the mortality model (cross‐validated spatial R 2: 0.18; Figure 1c–d). These patterns were consistent with the independent comparison dataset (Figure S2) and other recent studies (31). The inclusion of forest physiological metrics for drought tolerance, specifically community‐weighted plant hydraulic traits, substantially improved the predictive accuracy of the climate stress mortality models (ΔAkaike Information Criterion << −10), consistent with drought‐physiology studies (Anderegg et al., 2016). Observed historical permanence risks to US forests from insect‐driven mortality specifically were highest in the Rocky Mountains and modelled risks captured the key broad spatial patterns in risks (cross‐validated spatial R 2: 0.31; Figure 1e–f). These observed and modelled insect risks showed strong spatial agreement with an independent continent‐wide insect outbreak dataset (Figure S3).

Under all future shared socioeconomic pathway (SSP) climate scenarios, fire risks are projected to increase substantially throughout the 21st century (Figure 2a). Future risks increase similarly across scenarios through mid‐century, but diverge by 2050. By 2080–2099, the multi‐model mean projects 4‐fold (SSP2‐4.5), 9‐fold (SSP3‐7.0), and 14‐fold increases (SSP5‐8.5) in US‐averaged fire risk compared to historical (average 1990–2019) values. Projected climate stress risks increased substantially and varied by emissions scenario with average mortality increases by a factor of 1.3 in SSP2‐4.5, 1.5 in SSP3‐7.0, and 1.8 in SSP5‐8.5 by 2080–2099 (Figure 2b). Future US‐wide insect risk projections indicated increases of 1.2‐fold in SSP2‐4.5, 1.4‐fold in SSP3‐7.0, and 1.7‐fold in SSP5‐8.5 by 2080–2099 (Figure 2c). We note that the climate stress and insect mortality models are not independent and thus should not be considered additive. All three climate‐sensitive risks showed large differences across climate models, although the relative ranking of risk by SSP was consistent by the end of the century. The substantial differences between high and low emissions SSPs emphasises the critical importance of climate policy to mitigate climate risks to US forests.

We then conducted a risk assessment to quantify which regions and forests are likely to experience the highest risks in the 21st century (Figure 3, average of 2080–2099). By the end of the 21st century, high levels of fire risk, which were historically confined to pockets in California and the intermountain western US, are projected to expand across the entire western US (Figure 4). While these risks are substantially mitigated by emissions reductions (SSP2‐4.5, Figure 3a), risks are still projected to increase dramatically in regions like the Great Plains in the central US and southeastern US (Figure 4).

Future climate stress risks increased most across broad swaths of the intermountain and southwestern US, California, and western Texas, although parts of the eastern US and the upper midwestern US also exhibited increased climate stress mortality risk (Figures 3b, 4). Projected insect risk to forest permanence was highest across the Rocky Mountains in the intermountain western US, Sierra Nevada mountains in California, and parts of the northern Midwest (Figures 3c, 4). Climate stress and insect mortality model projections were only made for forest types 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.

DISCUSSION

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 (Bugmann et al., 2019; Fisher et al., 2018; Pugh 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) (Figure 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 (Figure 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 insect‐driven 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 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 visualised 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, large‐scale 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 (Abatzoglou et al., 2021; Barbero et al., 2015). 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‐characterised 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 (Figure 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 programmes. Further work could incorporate observed heterogeneity and future projections to better inform the construction of these climate programmes, such as by translating these risks into specific C loss estimates that could parameterise a better‐grounded buffer pool or other insurance programme. 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.

AUTHOR CONTRIBUTION

W.R.L.A., O.S.C., J.F., and J.J.H designed the project. W.R.L.A., O.S.C., G.B., A.T.T., J.F., and J.J.H conducted the analyses. W.R.L.A. and O.S.C. drafted the initial paper and all co‐authors provided input on writing and analyses.

PEER REVIEW

The peer review history for this article is available at https://publons.com/publon/10.1111/ele.14018.

Supplementary Material

Click here for additional data file.