Plant-water sensitivity regulates wildfire vulnerability

Extreme wildfires extensively impact human health and the environment. Increasing vapour pressure deficit (VPD) has led to a chronic increase in wildfire area in the western United States, yet some regions have been more affected than others. Here we show that for the same increase in VPD, burned area increases more in regions where vegetation moisture shows greater sensitivity to water limitation (plant-water sensitivity; R2 = 0.71). This has led to rapid increases in human exposure to wildfire risk, both because the population living in areas with high plant-water sensitivity grew 50% faster during 1990–2010 than in other wildland–urban interfaces and because VPD has risen most rapidly in these vulnerable areas. As plant-water sensitivity is strongly linked to wildfire vulnerability, accounting for ecophysiological controls should improve wildfire forecasts. If recent trends in VPD and demographic shifts continue, human wildfire risk will probably continue to increase.

W ildfire-burned area has increased rapidly in the western United States over the past four decades 1-3 , threatening humans and altering ecosystem function and vegetation growth patterns [4][5][6] . Several factors have contributed to the rise in burned area in the western US. Historic fire suppression has increased fuel loads in many areas 7 . Human populations have also expanded the wildland-urban interface (WUI)-a zone of transition between wildland and urban environments 8 -where human-caused ignitions are frequent and vegetation is abundant, posing high risk to human lives and structures [9][10][11] . Moreover, anthropogenic climate change and natural climate variability have combined to substantially increase the atmospheric aridity, which has contributed to a decline in fuel moisture and a resultant rise in burned area 2 . However, the sensitivity of burned area to atmospheric aridity can vary significantly across regions [12][13][14] . Understanding and predicting fire hazard as hot, dry conditions become more common [15][16][17] require accounting for the factors that modulate this sensitivity, including fuel availability and fuel moisture 18,19 .
As atmospheric aridity increases, live fuel moisture content (LFMC, measured as the mass of plant water per unit dry biomass) generally decreases 20,21 . However, the effect of atmospheric aridity on LFMC is regulated by a range of location-specific factors, including, but not limited to, topographic and soil controls on root-zone water availability, root water uptake and plant hydraulic traits that affect transpirational water loss. Specifically, plant hydraulic traits can cause up to 3-fold variation in LFMC, and thus affect fuel flammability 22 . For the same meteorological conditions (including both atmospheric aridity and precipitation), fuel moisture can vary widely, depending on plant species and hydraulic diversity [23][24][25] . However, a single fire can impact a wide range of species in a community, with a potentially diverse range of plant hydraulic traits within the area burned 26,27 . It is therefore unclear whether -and if so, to what extent -plant and soil hydraulic traits affect the spatial distribution of LFMC and fire hazard at large scales.
We investigate how the plant and soil features influencing LFMC's response to climate affect fire hazard. We consider the effects of climate on LFMC through a climate-derived moisture balance that considers both precipitation and vapour pressure deficit (VPD, calculated using the wetness of dead foliage and twigs in the litter; see Methods). We refer to the integrated sensitivity of LFMC to climate-derived moisture balance as plant-water sensitivity (PWS). Because of the significant memory in soil moisture 28 , LFMC depends on both current and previous meteorology. To account for this, we calculate PWS as the sum of the slopes of a multiple linear regression between LFMC anomalies and climate-derived moisture balance anomalies, with lags varying from 0 to 150 d in 15 d intervals.
We hypothesize that PWS regulates the effect of climate on burned area. Specifically, a unit decrease in climate-derived moisture balance will cause a larger decrease in LFMC in ecosystems where PWS is high. This decrease in LFMC may result in higher flammability 29 and, eventually, a larger burned area than in ecosystems with low PWS where the decrease in LFMC for the same decrease in climate-derived moisture balance is smaller (Extended Data Fig. 1). To test this hypothesis, we use LFMC maps derived from microwave remote sensing 30 . Furthermore, to assess whether current vegetation distributions will buffer or exacerbate future fire vulnerability in the western US, we test the interactions between PWS and other drivers of fire risk, including the rate of change in atmospheric aridity across the region, and the rate of growth of the WUI population.

Results
Link between PWS and wildfire vulnerability. The slope between timeseries of annual burned area and annual VPD ( is strongly linked to the PWS (R 2 = 0.71 and P < 0.0001; Fig. 1a). For different bins of PWS, spatially disparate locations with similar PWS are combined and the interannual variations of burned area and VPD is calculated across these locations (Fig. 1a).
is the slope between annual burned area and VPD for the years 2001-2020, and ranges from ~350 to 700 km 2 hPa −1 (Fig. 1c), depending on ecosystem PWS. We illustrate the simplified calculation of PWS for two pixels from high-and low-PWS ecosystems in Fig. 1b is robust in both shrublands (R 2 = 0.64, P = 0.005) and forests (R 2 = 0.69, P = 0.003), but uncertain in grasslands ( Supplementary Fig. 3). In grasslands, the slope between burned area and VPD is consistently low, resulting in a lower sensitivity to PWS (R 2 = 0.17, P = 0.24 increases sharply with PWS. Hereafter, we refer to this definition of 'hazard' when interpreting our results, but use the word 'risk' when the hazard coincides with human exposure and vulnerability. Note that the hazard zones as defined here only refer to the increase in burned area and are not related to the losses or benefits from the area burned. To verify that the effect of PWS on d(burned area) is causal, we test the correlation between PWS and several other potential biogeographical confounders that may correlate with PWS, but also affect d(burned area) d(VPD) . However, we do not find any strong correlation between PWS and mean or variance of VPD, mean normalized difference vegetation index (NDVI), dry-season NDVI, or dry-season length (see Supplementary Discussion 1 and Supplementary Fig. 4). Because each bin in Fig. 1a represents a collection of locations that are disparate in space ( Supplementary Fig. 2), the likelihood that other untested confounders exist is reduced.
Although we define PWS as the sensitivity of LFMC to climate-derived moisture balance, it is worth noting that this moisture balance will be influenced by VPD. Thus, VPD influences both the x axis and y axis of Fig. 1a, introducing the potential for circularity and artificially strong correlation. To evaluate this possibility, we compute a modified PWS that uses VPD as the predictor of LFMC rather than the climate-derived moisture balance. The correlation between this modified PWS and the original PWS is low (R 2 = 0.13; see Supplementary Discussion 2 and Fig. 5), suggesting that VPD only weakly affects PWS relative to other factors. Thus, this potential circularity has little effect on the strong correlation in Fig. 1a, which instead is probably driven by LFMC's effects on fire flammability and spread.
Antecedent precipitation (during the months of December-May, before the peak fire season) can spur vegetation growth, increasing fuel loads, and ultimately burned area, especially in arid ecosystems 12 . If precipitation-led growth in fuel availability is correlated with PWS, we may overestimate the role of PWS. However, we find no positive correlation between PWS and the sensitivity of fuel availability (represented by NDVI) to antecedent precipitation for any land cover ( Supplementary Fig. 6). Since the effect of antecedent precipitation on fire-season fuel availability can vary with mean precipitation, we also re-compute the previous calculation across the gradient of mean precipitation, but again find no significant correlations ( Supplementary Fig. 7). Overall, these sensitivity analyses suggest that the relationship in Fig. 1a is robust and probably driven by changes in LFMC dynamics across ecosystems with different PWS. prompts the question, what drives spatial variability in PWS?

Drivers of PWS. The strong link between PWS and
Using a random forest model (see Methods), we find that 58% of the variance in PWS is explained by 14 plant and soil hydraulic traits ( Supplementary Fig. 8). The most important drivers of PWS are average saturated soil hydraulic conductivity (K s , 20% importance), the shape of soil water retention curves (n, 13% importance) and root depth (9% importance) (Fig. 2). In total, soil and plant traits contribute to 55% and 45% of the explained importance, respectively. PWS is higher where VPD is increasing fastest. Given the role of PWS in regulating the sensitivity of burned area to VPD, we investigate how PWS has affected human exposure to wildfires. From 1980 to 2020, VPD increased across 91% of the western US, with a mean increase of 0.05 hPa yr −1 (Fig. 3a,d). We observe that PWS is higher in regions where VPD increased more quickly (Fig. 3a,b). The asymmetry in the joint distribution of VPD trends and PWS results in ~28% of the western US having both PWS and VPD trends that are greater than their respective median values (Fig. 3a). Results are similar when VPD trends are analysed relative to their long-term mean rather than as absolute values ( Supplementary Fig. 9). Since high PWS and high VPD growth both enhance fire hazard, their co-occurrence is likely to amplify increases in burned area. This effect of 'double-hazard' , where PWS is high (≥1.5) and VPD increased faster than average, is concentrated in regions such as the Sierra Nevada in California, eastern Oregon, the Great Basin in Nevada and the Mogollon Rim in Arizona (pink contour in Fig. 3c,d).
The WUI has expanded disproportionately in high PWS regions.
Beyond burned area, wildfire risk is also shaped by patterns of human exposure and vulnerability, both of which are strongly influenced by the distribution of the WUI. The WUI areas are proximal to large fuel loads, and the large number of humans and settlements in the WUI elevates both exposure and vulnerability. We thus investigate spatial patterns in WUI population growth as a function of PWS. Between 1990 and 2010, the population living in the WUI roughly doubled, growing from 10 million to 20.8 million (Fig. 4a), an increase of 108% (Fig. 4b). However, the WUI population rose most rapidly in regions with high PWS-driven wildfire hazard (160% rise). These are the same regions where burned area rose most rapidly relative to 2001 per unit rise in VPD. By contrast, in low-and medium-hazard regions, the WUI population grew by 107% and 95%, respectively.

Discussion
The relationship in Fig. 1a 18,36,37 . Although the relationship between interannual burned area and climate aridity is strong at ecoregion and landscape scales 2,38,39 , the role of vegetation in amplifying or dampening the effect of aridity on wildfires locally is still poorly understood. This is in part because LFMC is usually modelled primarily on the basis of meteorological conditions 40,41 . Given our empirical evidence that PWS -and thus, soil and plant hydraulic traits (Fig. 2) -regulate burned area variability (Fig. 1), a greater focus on the ecophysiological controls of burned area is needed to fully understand burned area drivers, both in the western US and elsewhere. Future studies that explicitly account for spatial variations in PWS or other ecophysiological controls may be better equipped for analysing and forecasting burned area. Global fire models parameterize vegetation controls on burned area using plant functional types 42,43 . However, the plant and soil hydraulic traits that control PWS (Fig. 2) are known to vary significantly within plant functional types 44,45 . This suggests that basing vegetation-fire relationships on functional types may lead to large errors 46 . The PWS metric (1) directly quantifies vegetation-climate sensitivity, (2) is related to burned area and (3) is scalable globally. The PWS metric thus offers a potential pathway to improve global vegetation-fire modelling.
Although PWS is not a widely recognized plant trait, we show that it is an indicator of whole-plant-water relations. It combines several plant and soil hydraulic traits that affect LFMC (Fig. 2). Notably, soil hydraulic traits explain a slightly larger fraction (55% vs 45%) of the variation in PWS than do plant hydraulic traits. This is to some extent expected; datasets of the explanatory variables have much coarser resolution than our PWS dataset, and because vegetation traits are more spatially heterogeneous than soil traits, this mismatch in scales adds noise and reduces the explanatory power of the plant hydraulic traits more than that of the soil hydraulic traits. Furthermore, because PWS is determined from LFMC, a property of the vegetation, even a co-variation with soil hydraulic traits merely reflects how plant water uptake is affected by soil hydraulics.

Fig. 2 | variable importance of plant and soil hydraulic traits to predict PWS.
Variable importance is estimated from average reduction in node impurity in random forests (see Methods). K s denotes saturated soil hydraulic conductivity, n denotes the shape parameter of soil water retention curves 74 , ψ 50 denotes xylem water potential at 50% loss in xylem conductivity, and g 1 denotes stomatal conductance slope parameter from ref. 78 , which is inversely proportional to the square root of water use efficiency. For description and data source of traits, see Supplementary Table 1.
Thus, we argue that our interpretation of PWS as a metric of the water sensitivity of plants is appropriate despite the influence of soil hydraulic traits (Supplementary Discussion 3). Vegetation drought acclimation may shift the PWS spatial pattern identified here, and thus its interaction with rising aridity and WUI population. Although drought-acclimated vegetation is likely to have lower PWS 47 , the pace of drought acclimation has probably been too slow to impact the spatial patterns of the double-hazard regions (Fig. 3) in the short term. For example, Trugman et al. 48 showed that ψ 50 of trees in the western US has reduced by ~0.001 MPa yr −1 between 2000 to 2019, but mean ψ 50 can vary from 1 to 5 MPa geographically. Thus, the spatial variation of ψ 50 is at least three orders of magnitude greater than its short-term temporal variation. In the medium-to-long term, however, major shifts in plant communities due to post-wildfire recruitment 7 could affect the overall spatial patterns in PWS.
Our results also show that in forests and shrublands, PWS is tightly linked to wildfire hazard. The shrublands in the arid ecosystems of the southwestern US are popularly conceptualized as fuel-limited, where wildfire hazard is expected to depend on fuel availability 21,37 . However, the moderately high PWS of shrublands, combined with the lack of any significant correlation between PWS and precipitation-driven fuel growth, suggests that arid shrublands in the southwest US may also be flammability-limited in some situations ( Fig. 3b and Supplementary Fig. 7). Our results thus provide a mechanistic explanation for empirical observations that previous droughts have led to a larger burned area even outside of forests in the southwestern US 49 . In grasslands, which are largely located in the Great Plains, the low slope between annual burned area and VPD suggests that other factors, such as availability of fuels, ignitions, strong winds and phenological stage, may have more dominant influences on burned area 50,51 . Since our data show that VPD played an insignificant role in governing burned area variability for grasslands ( Supplementary Fig. 3), and since the temporal resolution of the LFMC dataset of 15 d may be insufficient to capture rapid responses to water limitation in grasslands, we are unable has 29% unexplained variance (Fig. 1a). Studies quantifying the effect of dead fuel moisture response, vegetation demographic shifts and fire behaviour on d(burned area) d(VPD) may help to resolve this ambiguity.
The geographic co-occurrence of high PWS and high VPD suggests that the distribution of vegetation in the western US has amplified the effect of climate change on wildfire hazard (Fig. 3). Even though VPD has risen in most of the western US over the past few decades (Fig. 3d), PWS is higher in regions with larger VPD trends (Fig. 3b). Many of the resulting double-hazard regions with high PWS and high VPD trends are in the southwestern US. However, the PWS variations can be highly localized (Fig. 3c). For instance, while the southern Sierra Nevada has very high PWS and is among the double-hazard regions, the northern Sierra Nevada has very low PWS. Since future projections indicate continued VPD increases within the western US 16 , the spatial distribution of those VPD trends may interact with the spatial pattern of PWS to change the spatial distribution of overall fire hazard.
The expanding WUI has been viewed as a contributor to rising wildfire risk primarily due to increased human-caused ignitions 9-11,52,53 . However, as Fig. 4 demonstrates, the rise in wildfire risk to humans is also due to increasing population within high-hazard regions of the WUI. For instance, the increase in population residing in the high-hazard regions of the WUI between 1990 and 2010 (1.5 million people) is roughly equivalent to the combined current populations of San Francisco and Seattle, and represents an expansion of the share of total WUI population living in high-hazard regions from 9.7% in 1990 to 12.2% in 2010. Previous studies have recorded the growth of the WUI 8,54 , but we provide evidence that the WUI population has grown most in the most fire-vulnerable ecosystems in the western US. The disproportionate expansion of the WUI into the high PWS regions suggests that the increase in wildfire risk in the WUI is at least partly due to the higher vulnerability of vegetation to fires in particular areas of the WUI. This is corroborated by the evidence that the high-hazard zone witnessed the highest increase in relative WUI population and the highest relative d(burned area) d(VPD) (Fig. 4b). Further research is needed to understand the socioeconomic demographics of the populations occupying these high-hazard areas. However, at a minimum, disproportionate growth in WUI areas with high PWS and rapidly increasing VPD suggests that previous estimates of changing wildfire risk in the western US may have been conservative. While noting that PWS is not the only factor influencing wildfire risk, the PWS dataset could inform local and state-wide priorities related to wildland development, land use planning, vegetation management and home hardening solutions to curb wildfire risk.
Although wildfire risk in the western US has increased due to rising aridity and population in the WUI, the increase in risk has been significantly regulated by PWS to climate. This suggests that wildfire danger models, such as the National Fire Danger Rating System 40 , which rely on LFMC estimates derived from meteorology alone, misrepresent wildfire danger because they do not account for spatial variations in PWS. Although our analysis focuses on the western US, the effect of PWS on burned area identified here is probably also present elsewhere across the globe. Overall, the concurrence of high PWS, rising aridity and increasing WUI population has increased wildfire risk to people in many parts of the western US. In particular, the most sensitive ecosystems (PWS > 1.5) exhibit the most rapid increase in WUI population, along with widespread occurrence of above-average increase in VPD.
Accurate understanding and quantification of the processes shaping wildfire hazard, exposure and vulnerability are critical, given the growing losses from wildfire, including loss of habitat of vulnerable species 55 , loss of human life and structures 9,10 , direct and indirect economic costs 56,57 , and widespread impacts on public health beyond the area that burns 58 . Representations of wildfire risk that do not account for interactions between ecological, atmospheric and human drivers are thus susceptible to mischaracterizing wildfire risk. As we show here, examining the interplay between climate change, human population dynamics and the role of vegetation in regulating wildfire hazard can elucidate hidden interactions that lead to greater wildfire risk overall.

Methods
Calculating PWS. We defined PWS as the sum of slopes of a constrained linear regression between LFMC anomaly and lagged climate-derived moisture balance anomalies, with lags varying from 0 to 150 d in 15 d intervals, and a constraint that each slope be non-negative (equations 1 and 2). Because the slope between LFMC and climate-derived moisture balance is considered, PWS is an indicator of the degree to which the set of soil and plant traits in each location buffers the impact of recent and current climatic variations on LFMC. The unweighted sum of slopes was used to equally account for both concurrent and antecedent conditions. Unlike previous studies that use an unconstrained regression 59 , we constrained the regression coefficients to be non-negative. We thus did not allow an increase in climate-derived moisture balance anomaly to negatively affect an associated LFMC anomaly, which would be unphysical. Lags varied in 15 d intervals to match the temporal resolution of the LFMC data. Anomalies were computed by subtracting the pixel-specific seasonal cycle from the raw quantities to eliminate spatially varying seasonalities.
We used dead fuel moisture content (DFMC) to represent climate-derived moisture balance. The DFMC is an indicator of wetness and combines precipitation, temperature and humidity over lagged timescales into one index. It is sensitive to not just demand-side components of the water balance, such as atmospheric moisture demand, but also supply-side components, such as precipitation 60 . Here we used the 100 h DFMC, which is a meteorological estimate of wetness of twigs and branches 1-3 inches in diameter. It is calculated from the 24 h antecedent precipitation and the equilibrium moisture content corresponding to the drying-wetting potential of the atmosphere during the same period after adjusting for the duration of daylight 61 .  Fig. 1a). b, Percent change in burned area (relative to 2001) per unit rise in VPD versus percent change in WUI population from 1990-2010 for each hazard zone.
We calculated PWS as follows.
where, LFMC ′ s,t = Live fuel moisture content anomaly = LFMCs,t − LFMC s,doy(t) (3) DFMC′ is the dead fuel moisture content anomaly (computation same as LFMC′), t is the date between June to November for the years 2016-2020, doy is day of year, s is the pixel index, γ is the intercept, β is the slope and i is the lag between LFMC′ and DFMC′ in number of days. The i increases in steps of 15 d intervals (half-monthly). The LFMC s,doy(t) represents the mean LFMC at pixel s, on day of year t (across 2016-2020). According to equations (1) and (2), if LFMC purely reflects climate-derived moisture balance, the PWS is high, whereas if LFMC is buffered against climatic water limitation (for example, through regulatory mechanisms such as stomatal control, root water uptake and so on 62 ), PWS is low. We used data from 2016 to 2020 to estimate PWS. We limited our observations to June-November to focus on the period during which plant water dry-down occurs due to climatic water limitation. By limiting our observations to this period, the sensitivity of PWS to other processes, such as plant growth in the spring, is minimized. We commenced our analysis in 2016, the first year of LFMC data availability. The PWS was not computed for the 9.5% of pixels where more than 50% of the LFMC temporal record was missing due to cloud or snow cover, or where synthetic aperture radar data was absent. A maximum plant water 'memory' of 150 d was chosen due to the apparent hydrological memory timescales of ecosystems of the western US 63 . The linear regression between LFMC′ and DFMC′ has a mean R 2 of 0.3, with a wide distribution of fits (R 2 ranging from 0 to 0.85; Supplementary Fig. 10). A relatively low mean R 2 of 0.3 is consistent with expectations, as plant and soil hydraulic traits that contribute to plant water control are known to cause significant differences between LFMC and meteorology 62,64,65 . Since the regression coefficients were constrained to be non-negative, ~50% of the coefficients are zero ( Supplementary Fig. 11).
We used LFMC data from an artificial intelligence model trained and validated in ecosystems of the western US using microwave and optical remote sensing 30 . The LFMC dataset was trained on field-measured moisture content of live leaves, needles and thin branches from the National Fuel Moisture Database 66 . The dataset has an overall R 2 = 0.63, RMSE = 25% and bias = −1.9%, with consistent performance across all land cover types. We used DFMC data from GRIDMET 67 . We resampled the DFMC data to 15 d averages to match the start and end times of the LFMC maps. We rescaled the LFMC dataset from 250 m to 4 km using bilinear interpolation to match the spatial resolution of the DFMC dataset.
Linking PWS to the sensitivity of burned area to VPD. We separated the PWS map into 15 regions with equal-vegetated area on the basis of the histogram of PWS ( Supplementary Fig. 1). Vegetated areas included regions with trees, shrubs and grasslands obtained from the National Land Cover Database version 2016 68 . We then computed the slope between burned area and VPD in each region using linear regression. Combining different locations of similar PWS enabled calculating the burned area, which otherwise was not possible at the pixel-level. Our results were consistent when PWS was separated into 10 regions instead of 15 ( Supplementary Fig. 12). The linear fits between burned area and VPD have R 2 ranging from 0.48 to 0.58 and P values ranging from 0.0001 to 0.0010.
Although many aridity metrics have been used to study interannual changes in burned area 2,37 , we chose VPD to represent atmospheric aridity because of its parsimonious definition, its direct connection to atmospheric moisture demand, and its strong link to burned area compared with that of other meteorologically based aridity indicators 49,69 . Further, by not choosing the same meteorological indicator as the one chosen to assess the response of PWS to climate-derived moisture balance (DFMC), we minimized the introduction of artificial cross-correlations.
For each of the 15 PWS bins, we computed the slope between annual burned area and mean annual VPD, each calculated across April to March. We computed the slope using linear regression with an unconstrained intercept. We used burned area data for 2001-2020 from the Moderate Resolution Imaging Spectroradiometer burned area product MCD64A1 70 , and VPD data from the Parameter-elevation Regressions on Independent Slopes Model 71 . We aggregated burned area from 500 m to 4 km resolution to match the resolution of VPD. We then computed the slopes and their associated linear standard errors for each of the 15 PWS bins, followed by computation of a final linear regression between d(burned area) d (VPD) and PWS (Fig. 1a).

Spatial drivers of PWS.
To investigate what factors influence PWS, we regressed the PWS map against eight plant hydraulic traits and six soil hydraulic traits using an ensemble approach. For plant traits, we used canopy height from 2005 from ref. 72 and several hydraulic traits, including xylem capacitance, a stomatal conductance slope parameter that is inversely proportional to marginal water use efficiency (g 1 ), maximum xylem conductance, the ψ 50 , and a hydraulic functional type derived from K-means clustering of these plant hydraulic traits determined from ref. 41 . We also included isohydricity (a stomatal regulation trait) from ref. 45 , and maximum rooting depth from ref. 73 . For soil traits, we used average saturated hydraulic conductivity (K s ) the shape parameter describing the exponent of soil water retention curves (n), soil porosity from ref. 74 , and fraction of soil, silt and clay from ref. 75 (Supplementary Table 1). We rescaled all traits to 4 km resolution using nearest neighbour interpolation since all traits (except for canopy height and rooting depth) have a coarse resolution of 25 km.
We used the rescaled traits as explanatory variables in a random forest regression to assess the importance of each trait. We ignored all time-varying factors in our analysis and used static maps of all traits to explain the variance in the static map of PWS. We used 3-fold cross-validation to evaluate strength of fit (Supplementary Fig. 8). To compute variable importance, we measured the average decrease in node impurity due to splitting at each variable and normalized the output (Fig. 2). For the random forest ensembles, we chose a minimum of six samples per terminal node, a minimum of 5 × 10 −6 reduction in node impurity at each split, a minimum of two samples split per node, bootstrapped samples during node splits, and 50 trees within each forest. We used scikit-learn for the random forest regression 76 .
Estimating VPD trends. We estimated VPD trends using a pixel-specific linear regression of the annual-average VPD against year, for the period 1980-2020. For further investigation, we calculated relative VPD trends by dividing the absolute VPD trends by the pixel-specific 1980-2020 VPD mean ( Supplementary Fig. 9).
Estimating WUI expansion. We used maps of population density and WUI from Martinuzzi et al. 77 to track WUI expansion from 1990 to 2010 in regions with different wildfire hazard. Martinuzzi et al. identified WUI regions using pre-defined thresholds for density of human settlements in wildland areas, with different thresholds for intermix (≥6.18 houses per km 2 and ≥50% cover of wildland vegetation) and interface regions (≥6.18 houses per km 2 and <50% cover of vegetation, located <2.4 km of an area ≥5 km 2 in size that is ≥75% vegetated). We did not differentiate between WUI intermix and WUI interface regions since both represent regions with high wildfire risk. We rescaled the WUI maps from their original census block-level resolution to 4 km to match the resolution of PWS. We classified any 4 km pixel that overlaps with a block-level WUI polygon as WUI since it contains at least some area with neighbouring urban and vegetated areas ( Supplementary Fig. 13).
Reporting Summary. Further information on research design is available in the Nature Research Reporting Summary linked to this article.

code availability
The scripts required to reproduce our results are available at https://github.com/ kkraoj/wildfire_from_lfmc.
Corresponding author(s): Krishna Rao Last updated by author(s): Dec 11, 2021 Reporting Summary Nature Portfolio wishes to improve the reproducibility of the work that we publish. This form provides structure for consistency and transparency in reporting. For further information on Nature Portfolio policies, see our Editorial Policies and the Editorial Policy Checklist.

Statistics
For all statistical analyses, confirm that the following items are present in the figure legend, table legend, main text, or Methods section.

n/a Confirmed
The exact sample size (n) for each experimental group/condition, given as a discrete number and unit of measurement A statement on whether measurements were taken from distinct samples or whether the same sample was measured repeatedly The statistical test(s) used AND whether they are one-or two-sided Only common tests should be described solely by name; describe more complex techniques in the Methods section.
A description of all covariates tested A description of any assumptions or corrections, such as tests of normality and adjustment for multiple comparisons A full description of the statistical parameters including central tendency (e.g. means) or other basic estimates (e.g. regression coefficient) AND variation (e.g. standard deviation) or associated estimates of uncertainty (e.g. confidence intervals) For null hypothesis testing, the test statistic (e.g. F, t, r) with confidence intervals, effect sizes, degrees of freedom and P value noted

Software and code
Policy information about availability of computer code Data collection https://github.com/kkraoj/wildfire_from_lfmc/tree/master/misc Data analysis https://github.com/kkraoj/wildfire_from_lfmc/tree/master/analysis For manuscripts utilizing custom algorithms or software that are central to the research but not yet described in published literature, software must be made available to editors and reviewers. We strongly encourage code deposition in a community repository (e.g. GitHub). See the Nature Portfolio guidelines for submitting code & software for further information.

Data
Policy information about availability of data All manuscripts must include a data availability statement. This statement should provide the following information, where applicable: -Accession codes, unique identifiers, or web links for publicly available datasets -A description of any restrictions on data availability -For clinical datasets or third party data, please ensure that the statement adheres to our policy The derived plant-water sensitivity map along with the scripts required to reproduce the results presented in this manuscript can be found at https://github.com/ kkraoj/wildfire_from_lfmc. Live fuel moisture content maps are available from https://kkraoj.users.earthengine.app/view/live-fuel-moisture. Climate data from GRIDMET is available from http://www.climatologylab.org/gridmet.html. Wildland-urban interface maps are available from http://silvis.forest.wisc.edu/data/wuichange/.