Recent California tree mortality portends future increase in drought-driven forest die-off

Vegetation tolerance to drought depends on an array of site-specific environmental and plant physiological factors. This tolerance is poorly understood for many forest types despite its importance for predicting and managing vegetation stress. We analyzed the relationships between precipitation variability and forest die-off in California’s Sierra Nevada and introduce a new measure of drought tolerance that emphasizes plant access to subsurface moisture buffers. We applied this metric to California’s severe 2012–2015 drought, and show that it predicted the patterns of tree mortality. We then examined future climate scenarios, and found that the probability of droughts that lead to widespread die-off increases threefold by the end of the 21st century. Our analysis shows that tree mortality in the Sierra Nevada will likely accelerate in the coming decades and that forests in the Central and Northern Sierra Nevada that largely escaped mortality in 2012–2015 are vulnerable to die-off.


Introduction
Droughts are among the worst climate hazards society faces, creating economic losses of tens to hundreds of billion US dollars per year (Mishra andSingh 2010, European Commission 2012). Forests are especially vulnerable to drought in a warming world, as higher temperatures increase atmospheric moisture demand, evapotranspiration and soil drying and tree die-off (Allen et al 2010, Williams et al 2013, Allen et al 2015, Fettig et al 2019. Future droughts are projected to become longer, more severe and frequent in many regions (Wehner et al 2011, Seneviratne et al 2012, Madakumbura et al 2019, with consequent impacts on vegetation structure and function. An improved understanding of forest response to drought is needed to better predict the impact of climate change on forested ecosystems at large scales.
Such capabilities would help forest managers anticipate patterns of tree vulnerability to drought and proactively allocate resources and time at local and landscape scales (Heinimann 2010, Keenan 2015. Studies investigating the causes of tree die-off have advanced our knowledge of forest tolerance to drought, and the underlying interactions of biological, physiological and environmental factors. Robust statistical relationships between drought induced tree mortality and various metrics of drought intensity have often revealed ecological thresholds The patterns of die-off ultimately depend on both drought intensity and duration, and processes that occur over a range of timescales (Mcdowell et al 2008). The biological response to drought stress has been investigated intensively at the scale of individual plants based on theoretical (Mcdowell et al 2008) and experimental (Barbeta and Peñuelas 2016) analyses, but studies that focus on multidecadal records and at landscape scales and larger are scarce. Given potentially large spatial variations in drought tolerance and drought severity across landscapes, such larger scale studies are needed. For example, while large variation in vegetation stress and tree mortality were noted across the Sierra Nevada landscape during California's 2012-2015 drought it is unclear whether the variation arises from forest drought tolerance (e.g. local water availability or plant physiological factors) or drought severity.
Comparisons of drought duration and remotelysensed vegetation properties can provide an opportunity to investigate plant response and tolerance to drought at large spatial scales. The Normalized Difference Vegetation Index (NDVI) shows strong responses to droughts lasting 2-4 months in arid and humid biomes, and 8-10 months in semiarid and sub humid biomes in the world (Vicente-Serrano et al 2013). In a regional setting, a comparison between California forest water balance and canopy density, as measured by canopy water content (CWC) (Asner et al 2016), showed a comparatively long time scale response, with coniferous vegetation tolerating drought for several years before accelerated dieoff (Brodrick et al 2019). Despite the insights from such comparisons, the controlling factors (e.g. environmental and tree physiological properties) behind drought tolerance have not been fully explored. Efforts to forecast the long-term impact of drought are also complicated by the likelihood that climate, environmental conditions, and vegetation distributions may be different in the coming century (Kelly and Goulden 2008, Parks et al 2018, Holsinger et al 2019. Here we explore the timescale of drought response for coniferous forests in California's Sierra Nevada (Myers et al 2000, Bales et al 2011. This region experienced a severe drought in 2012-2015 that was followed by widespread die-off in 2015-2016. We used data from pre-2012 to identify the duration of drought that was best correlated with anomalies in canopy density time series as measured by the Normalized Difference Moisture Index (NDMI) (hereafter referred to as the drought sensitivity timescale or DST). We then investigate how the pre-2012 DST interacted with drought severity to produce the observed spatial and temporal patterns of die-off during the 2012-2015 drought. Finally, we use output from state-of-the-art global climate simulations to examine how changes in multiyear droughts may amplify future die-off episodes.

Vegetation indices
We used vegetation indices NDMI, NDVI, and CWC in this study. NDVI and NDMI were derived from Landsat 5, 7 and 8 surface reflectance and brightness temperature images. Data were obtained from USGS (https://espa.cr.usgs.gov) for the period 1984-2017 after being regridded to a resolution of 0.0002695 • (approximately 30 m). Snow-and cloud-affected pixels were removed using the Landsat collection 1 pixel-quality data layers. We used late growing season NDMI (August-October) in this study. Further details of the derivation of NDMI and NDVI can be found in Goulden and Bales (2019). Dry season (July-August) CWC data at 30 m resolution (derived in Brodrick et al 2019) was obtained from https://doi.pangaea.de/10.1594/PANGAEA.897276.
As a direct measurement of tree mortality, we used the number of dead trees from the Aerial Detection Survey data from the USFS (www.fs.usda. gov/detail/r5/forest-grasslandhealth/?cid=fsbdev3_ 046696). We re-projected and rasterized the dead trees per acre (DTPA) '1' product in this geodatabase.

Historical climate data
To derive drought indices, we used high resolution monthly precipitation (PR) data (Flint and Flint 2012), potential evapotranspiration (PET) and climatic water deficit (CWD) from the Basin Characterization Model (Flint et al 2013), a high resolution physically-based hydrologic model developed for California (https://ca.water. usgs.gov/projects/reg_hydro/basin-characterizationmodel.html#provisional). Data from 1980-2016 were used at 270 m resolution. Annual mean NDVI was used to calculate the total evapotranspiration (ET) of water year from an exponential relation between ET and NDVI. This relation is derived from in situ measurements of ET, as in Goulden and Bales (2019).

Future climate projections
To assess future drought conditions, we used monthly precipitation for the period 1850-2014 from historical simulations and from 2015-2100 for SSP2-4.5, SSP3-7.0 and SSP5-8.5 warming scenarios (O'Neill et al 2016) from ten of the state-of-the-art global climate models (GCMs) (supplementary table 1, available online at stacks.iop.org/ERL/15/124040/mmedia) participating in Coupled Model Intercomparison Project Phase 6 (CMIP6) (Eyring et al 2015). We selected available models with two or more ensemble members for all three scenarios at the time of analysis. In total, 90 ensemble members per scenario were used. Use of a large ensemble dataset allows the sampling of extreme conditions without using statistical resampling methods (Swain et al 2018).

Predictors of DTPA
We carried out a regression analysis to explain the spatial variations of dead trees per acre. Apart from climate and vegetation variables explained above, the following variables were extracted for the analyses: To represent resource competition (Young et al 2017), tree basal area data were obtained from the LEMMA group (https://lemma.forestry.oregonstate.edu/data/structu re-maps, variable: basal area of live trees ⩾2.5 cm diameter at breast height). Soil plant available water content (AWC) was obtained from the US General Soil Map database (http://websoilsurvey.nrcs.usda.gov/) by extracting the variable aws0150wta (units: cm). Subsurface drying can be directly linked to the amount of evapotranspiration exceeding the precipitation (Goulden and Bales 2019). To represent the vegetation-induced subsurface moisture use, we use the mean evapotranspiration as a predictor variable.

Conifer fraction and fire
We limited the analysis to conifer dominated forests in the Sierra Nevada. To mask out non-coniferdominated regions, we used the existing vegetation classification maps from the USFS (www.fs.usda. gov/detail/r5/landmanagement/resourcemanagement/ ?cid=stelprdb5347192). From the vegetation classification maps, we selected pixels with vegetation types: Sierran mixed conifer, Ponderosa pine, montane hardwood-conifer, lodgepole pine, red fir, white fir, subalpine conifer, Jeffrey pine, Douglas fir, Eastside pine and pinyon-Juniper. We removed fireimpacted pixels from 1980 through 2016 from the analysis by using the fire history data product from California's Fire and Resource Assessment Program (https://frap.fire.ca.gov/mapping/gis-data/).
All vector data used in this study were first rasterized and all data were reprojected to World Geodetic System 1984 using ArcGIS 10.7 (https://desktop.arcgis.com/en/) and the GDAL library (https://gdal.org/). All analyses (except for future projection from climate models) were done after bilinearly interpolating to the resolution of climate data (~270 m) using the Climate Data Operator library (Schulzweida 2017).

Calculation of drought indices
We consider four drought indices. Standardized Precipitation Index (SPI; Mckee et al 1993), cumulative precipitation minus evapotranspiration (PR-ET), standardized precipitation-evapotranspiration index (SPEI) and CWD. Using observed data, we calculated drought indices with integration periods from 1 to 6 years. The cutoff of 6 years was chosen considering the most significant multiyear droughts during the historical period in California (https://water.ca.gov/Water-Basics/Drought). All the years here are 'water years' , defined as starting from October of the previous year and ending in September of the corresponding year. For example, the 4 year SPI corresponding to 2015 would be calculated from the standardized PR from October-2011 to September-2015September- (i.e. water years 2012September- -2015. PR-ET and CWD were calculated in a similar manner to SPI and SPEI, but instead of standardizing, we calculated the cumulative PR-ET and CWD during the drought integration period. Such cumulative moisture deficits have been linked with physiological thresholds of tree mortality (Anderegg et al 2015, Goulden and Bales 2019).
For CMIP6 ensembles, historical and future time series of each ensemble member for each grid cell were first merged to obtain a time-series spanning 1850-2100. SPI for 4 year integration periods was calculated from these time series.

Calculation of DST
To determine DST, we calculated temporal correlations between vegetation and drought indices for each pixel (for the period 1983-2011), and retained the timescale with the maximum correlation: First, the vegetation anomaly was obtained for NDMI and CWC by removing the long term median (1984-2011for NDMI and 1990-2011. This can be considered to be the vegetation change (e.g. dNDMI and dCWC). For CWD, the long-term (1980-2011) median was also removed (but not for SPI, SPEI and PR-ET, since they are already in anomaly form). Considering the documented delayed response of NDMI and CWC to drought (Goulden and Bales 2019, supplementary figure 1), a lag of 1 year was imposed between drought and vegetation anomalies in the correlation calculations. For example, in the case of dNDMI and SPI, the 1984-2012 dNDMI timeseries was correlated with the 1983-2011 water year SPI timeseries. The drought anomaly time-series were obtained for various timescales from the four drought indices and the above steps were repeated, e.g. 4 year PR-ET in 2015 is the cumulative 2012-2015 PR-ET and the corresponding vegetation anomaly for that year is 2016 dNDMI and dCWC. Finally, the timescale giving the maximum positive correlation between vegetation anomalies and SPI, SPEI and PR-ET (or for CWD, the minimum negative correlation) was obtained as the drought timescale having maximum influence on vegetation, similar to the methodology used in Vicente-Serrano et al (2013). Only pixels where the correlation was significant at 80% were plotted in figures and used for analysis.

Random forest model
A random forest (RF) regression analysis was carried out to examine drivers of DTPA. RF is a nonparametric supervised machine learning algorithm (Breiman 2001). It has the advantage of handling nonlinear interactions between variables. RF is widely used in ecological studies to identify the key components of complex processes and their relative importance (

Statistical analysis and significance tests
The Pearson correlation coefficient was calculated in temporal correlation analyses. A two-tailed Student's t-test was applied to calculate p-values and statistical significance, including that of the difference in means of different scenarios of the multimodel GCM ensembles.
To obtain large geographic scale patterns of the DST from the pixel-wise analysis conducted using high resolution (~270 m) data, spatial smoothing was conducted. For spatial smoothing, K-nearest neighbor (KNN) regression was applied using the python package scikit-learn (https://scikit-learn.org/stable/). The KNN regression algorithm employs a userspecified distance measure and a threshold K to find the nearest K neighbors to each item. It then assigns the value based on the average. The K value (the distance between points during the spatial smoothing) was taken to be 250 pixels (an approximate distance of 70 km). This was chosen as a compromise between smoothing and retaining small scale spatial features (supplementary figure 4).

Spatial patterns and the interpretation of DST
The drought sensitivity timescale (DST) was calculated across the Sierra Nevada (see figure 2 for geographic domain) based on the relationship between SPI and NDMI ( figure 1(a)). Different forest types in the Sierra Nevada have different DST values, with mesic stands having DSTs around 3-4 years (see supplementary discussion). The longest drought time scales (>4 years) are seen in an elevation band around 1500-2500 m and a latitude range of 36-39 • N. Lower or higher elevation forests generally have shorter timescales (<2 years), depending on latitude. The timescale generally decreases with increasing elevation south of 38.5 • N. Alternative combinations of drought and canopy moisture indices revealed time scales that agreed with the spatial distribution obtained from the SPI-NDMI analysis (supplementary figures 5-7).
A simple interpretation of DST is that it represents the number of years required to empty a previously full, subsurface water storage under steady rate conditions. If this interpretation is correct, we should be able to multiply the mean drawdown rate by DST to predict subsurface water storage depth. We can quantify the subsurface storage depth as the maximum rootzone drying during the recent past ( . Carrying out the multiplication described above, we find that the simple interpretation of DST holds remarkably well for Sierran forests (supplementary figure 8). We further explore the drivers of DST and its interpretation in supplementary discussion. DST may be shaped by complicated environmental, climate and biological processes across temporal and spatial scales, but preliminary analysis shows that the factors thought to control plant water drawdown (Fellows and Goulden 2017) such as mean precipitation, evapotranspiration, potential evapotranspiration, plant available water content, can account for DST variation (from a random forest regression model with R 2 = 0.79, supplementary figure 9, see supplementary discussion). This confirms that DST reflects the plant accessible subsurface water buffer.   (2019). Error bars represent ±0.5 of the standard deviation. For panels (b)-(g), the elevation bin width is 50 m, and the x-axis is the mean elevation value of each elevation bin. In panels (b)-(g), only cold, mesic and dry forests (see supplementary discussion) were included for the analysis. The few grid cells where the mean PR-ET was less than zero were removed from data shown in panels (b)-(g). figure 1(a). In 2013, the year after the beginning of the drought, there were only modest losses of canopy moisture at all elevations (figure 1(e)). By 2014 canopy moisture loss increased for elevations where the normalized cumulative PR-ET was negative (figure 1(f)), and by 2015 extreme canopy loss was seen at these elevations (figure 1(g)).

Usefulness of DST as a vegetation stress predictor
In 2013, after 1 year of drought, some elevations showed moderate vegetation stress (figure 1(e)), but no systematic relationship between DST and vegetation stress. However, as the drought lengthened, vegetation stress became most apparent at those elevations where DST is lowest and drought severity is high, i.e. orange/reddish dots in figure 1(f) (Elevations above about 2700 m are characterized by a low DST, but they did not experience severe drought and hence remained unaffected). By 2015, after an even longer period of extreme precipitation deficit (at least 3 years), even elevations with the longest DST showed signs of extreme canopy moisture loss. Only the highest elevations avoided moisture loss (figures 1(b)-(d)). Canopy moisture loss precedes tree mortality and could ultimately trigger it (Brodrick and Asner 2017, Paz-Kagan et al 2017, Goulden and Bales 2019). The predictive power of DST for vegetation conditions is apparent in directly measured tree mortality, as shown in figures 1(e)-(g). Very low tree mortality is seen in 2013, similar to the background rate (Byer and Jin 2017). By 2014 values were slightly higher for those elevations with signs of vegetation stress, and relatively low DST values (orange/reddish dots in figure 1(f)). By 2015 a dramatic increase in tree mortality is seen everywhere except the highest elevations.
These results suggest that both DST and drought severity are useful predictors of the progression of vegetation stress during drought. To test this conclusion further, we examine the geographical variation in tree mortality at lower elevations of the Sierra Nevada ( figure 1(a)). Figure 2(b) shows the observed DTPA at the end of the drought. In the low elevations of the Southern Sierra, where DST is low ( figure 1(a), supplementary figure 6), the die-off is high. By contrast, the low elevations of the Northern Sierra, where DST is similarly low, little tree mortality occurred. Clearly the early signs of vegetation stress at low elevations (figure 1(f)), and the higher levels of eventual tree mortality ( figure 1(g)) came from the southern portion of the Sierra Nevada. This spatial pattern appears to reflect the greater precipitation shortfall in the south relative to the north ( figure 2(a)).
To test this hypothesis, we predicted the spatial distribution of tree die-off assuming the entire Sierra Nevada faced a precipitation deficit as large as that seen in the southern part (i.e. 4 year SPI = −2.5). The hypothetical distribution is generated with a model of DTPA based on random forest regression. As predictors, we use direct and indirect factors  contributing to drought severity (4 year SPI, mean potential evaporation during the drought, mean evapotranspiration, basal area) and drought tolerance to tree die-off (DST, plant available water capacity) (model R 2 = 0.64). The model predicts that mesic and dry forests in central and northwestern regions of the Sierra Nevada ( figure 1(a)) suffer marked tree die-off (figure 2(c)), consistent with their low DST values. Thus, differences in drought severity can account for the differences in tree mortality in those low elevation zones where DTS is comparably low.
To assess the potential of droughts similar to the 2012-2015 California drought in the future, we analyzed future projections of multiyear precipitation variability in California simulated by ten state-of the art GCMs (see supplementary table 1) participating in the CMIP6 (Eyring et al 2015). The simulated changes in 4 year SPI occurrence from historical  to future (2050-2100) under a 'no emissions reduction policy' warming scenario, SSP585 (Eyring et al 2015, O'Neill et al 2016, are shown in figure 3. Under climate change, the precipitation distribution shifts systematically, resulting in increased probability of severe multiyear droughts ( figure 3(a)). Examining scenarios associated with lower greenhouse gas emissions ( figure 3(b)), we see that the probability of severe droughts similar to 2012-2015 increases with growing emissions. But even for the lowest emissions scenarios, this increase in frequency is statistically significant at 99% compared to the historical simulations. We also note that the future may bring droughts even more severe than that in 2012-2015, as seen in the emergence of 4 year SPI anomalies more negative than −2.0 in the SSP585 scenario ( figure 3(a)). Such deep droughts would bring very dry conditions to large swaths of the Sierra Nevada, killing trees over areas with low DST, as in the hypothetical case of figure 2(c).

Discussion and conclusions
Using multiple drought and vegetation indices spanning many decades, we obtain a drought sensitivity timescale, DST, which can also be interpreted as the plant water buffer. We find that forests on the lowelevation western slopes and high-elevation eastern slopes of the Sierra Nevada have the shortest DST, and hence the least tolerance to drought. When DST is combined with drought severity, it can be used to map the progression of vegetation stress during the 2012-2015 drought. The low elevation slopes were the first to respond to the drought as it deepened, consistent with their low DST values. However, the largest ultimate response, and the greatest tree mortality, occurred in the southern low elevations. We demonstrate that this enhanced response can be explained by the spatial patterns of drought severity. The low elevation forests in the central and northern Sierra Nevada are also vulnerable to drought and would have likely experienced extensive dieback if the 2012-2015 drought had extended further north.
Previous work suggests interannual precipitation variability in California will increase in the future (Swain et al 2018). This work, combined with the results shown here demonstrating the connection between tree mortality and multiyear drought, raises the question of whether Sierra Nevada forests will experience greater tree mortality in the future. We found that multiyear droughts in California will increase with increasing greenhouse gasses based on state-of-the-art climate model simulations from the CMIP6 project (Eyring et al 2015). The distribution of multiyear precipitation anomalies shifts to a drier regime, and the likelihood of a 4 year drought as deep (or deeper) than the 2012-2015 event increases by up to a factor of three by the end of the century. These results imply a future increase in the likelihood of tree mortality in the Sierra Nevada, especially in areas with a low DST. The low-elevation central and northern Sierran forests did not exhibit die-off in the 2012-2015 event; but this does not mean they are not vulnerable to drought. Future droughts will almost certainly be distributed differently in space from the 2012-2015 event (as were the pre-2012 droughts which allowed us to diagnose the sensitivity of the low-elevation central and northern forests with low DST).
We note the warmer temperatures accompanying these projected droughts (Dai 2011(Dai , 2013, and that a limitation of our study is that our estimates of future drought ignore this warming effect. Warming will increase potential and actual evapotranspiration, which will further increase drought stress when precipitation is low; thus, our estimates represent a lower range of possible increases to drought intensity. Conversely, it is possible that plant physiological responses to elevated CO 2 may mitigate tree mortality impacts of drought, though it is unclear how much tolerance to mortality such responses confer (Swann et al 2016, Sperry et al 2019. Lastly, future climate projections are based on coarse resolution GCMs, and the spatial patterns of climate variables in this region of complex topography are inadequately represented. Future research can refine our analysis using downscaled climate model data, as well as state-ofthe-art land surface and vegetation models to simulate evapotranspiration changes and physiological responses to CO 2 .
Our results imply future changes in Sierra Nevada ecology. Previous work has shown that projected changes in mean climate would be associated with shifts in vegetation distribution (Holsinger et al 2019). Changes in mean conditions may lead to a decline in Sierra Nevada forests assuming current vegetation types migrate to new preferred climate zones (Parks et al 2018). These ecological transitions may be accelerated by changes in extreme events. Previous work has shown that tree mortality markedly impacts subsequent species composition and forest structure (e.g. Cobb et al 2017). Future increases in forest-die-off frequency and magnitude would likewise be associated with large impacts. The increase in tree mortality may increase fuel load and wildfire risk (Ruthrof et al 2016, Stephens et al 2018. Our results point to the low elevation western slopes of the Sierra Nevada as a hotspot of increasing die-off. Conversely, other sub-regions with high drought tolerance may be more tolerant to ecological change, and may even become refugia (Morelli et al 2016, Mclaughlin et al 2017. The drought timescale metric and results we show potentially provide information that may aid conservation planning.

Acknowledgments
We acknowledge the World Climate Research Programme's Working Group on Coupled Modelling, which is responsible for CMIP6, and we thank the climate modeling groups for producing and making available their model output. We also thank the Earth System Grid Federation (ESGF) for archiving the data and providing access, and various funding agencies who support CMIP6 and ESGF. This research was funded by the University of California Laboratory Fees Research Program. CDK is supported by the Department of Energy, Office of Science, Office of Biological and Environmental Research through the Early Career Research Program administered by the Regional and Global Model Analysis Program. CAN was supported by the Ridge to Reef Graduate Training Program funded by NSF-NRT award DGE-1735040. GDM would like to thank Neil Berg for discussions on early ideas of this paper.

Data availability statement
Climate model data were obtained from the CMIP6 archive, accessed via the ESGF at https://esgfnode.llnl.gov/projects/cmip6/. The data that support the findings of this study are available upon reasonable request from the authors.