Projecting future fire regimes in semiarid systems of inland northwestern U.S.: interactions among climate change, vegetation productivity, and fuel dynamics

Fire regimes are influenced by both exogenous drivers (e.g., increases in atmospheric CO2; and climate change) and endogenous drivers (e.g., vegetation and soil/litter moisture), which constrain fuel loads and fuel aridity. Herein, we identified how exogenous and endogenous drivers can interact to affect fuels and fire regimes in a semiarid watershed in the inland northwestern U.S. throughout the 21st century. We used a coupled ecohydrologic and fire regime model to examine how climate change and CO2 scenarios influence fire regimes over space and time. In this semiarid watershed we found that, in the mid-21st century (2040s), the CO2 fertilization effect on vegetation productivity outstripped the effects of climate change-induced fuel decreases, resulting in greater fuel loading and, thus, a net increase in fire size and burn probability; however, by the late-21st century (2070s), climatic warming dominated over CO2 fertilization, thus reducing fuel loading and fire activity. We also found that, under future climate change scenarios, fire regimes will shift progressively from being flammability to fuel-limited, and we identified a metric to quantify this shift: the ratio of the change in fuel loading to the change in its aridity. The threshold value for which this metric indicates a flammability versus fuel-limited regime differed between grasses and woody species but remained stationary over time. Our results suggest that identifying these thresholds in other systems requires narrowing uncertainty in exogenous drivers, such as future precipitation patterns and CO2 effects on vegetation.

 For a given vegetation type, there are temporally stable thresholds-specifically, the ratio of changes in fuel loading to changes in fuel aridity-that determine whether a location is fuel   Figure 2). The mean annual precipitation in the area is around 980 mm, of 139 which 60% is snow (Frenzel, 1989). Trail Creek is characterized by cold, wet winters and warm, Creek are mainly coarse, permeable alluvium (Smith, 1960). No  154 Table 1. Fire regime groups and corresponding characteristics ).  and the 10th percentile. Based on the ranking shown in Figure S1, we selected four models to     each year), and number of fire starts (fires that surpass 30 burned patches in 30 years' assessment 332 period). We then compared distributions of these variables among the model scenarios.

333
We also considered patch-level summaries of burn probability (Pburn), fuel aridity, and 334 fuel loading. For each patch, the burn probability (Pburn) was calculated as: 336 Fuel aridity is calculated by WMFire as relative deficit and was calculated monthly as: We calculated patch-level mean burn probability, fuel loading, and fuel aridity to understand the and CO2 fertilization effects on fire activity, while "climate change effect" represents the climate 344 change effects without CO2 fertilization. increased decomposition rates, which further reduced fuel loading ( Figure S6 and Figure S9).

384
Faster decomposition of litter and reduced vegetation productivity are the dominate driver than   scenario see Figure S10.

423
In the 2070s, changes in Pburn were more consistent across storylines and fuel conditions. converged with Region A (grass-dominated), where the watershed as a whole became intensely fuel-limited (caused by climate change effect on vegetation growth and litter decomposition).

427
The CO2 fertilization effect tempered the magnitude of decreasing fuel loads in the 2070s but

432
In all cases, the CO2 fertilization effect had an upper limit (i.e., around 610 ppm, which        indicates an increase (e.g., increases in productivity cause increases in fuel loading; "-" 598 indicates a decrease (e.g., increases in litter decomposition lead to decreases in fuel loading).

599
"Mesic" and "arid" indicates that the specified affect is location dependent and occurs in mesic

662
Therefore, even though we observed a non-monotonic trend in fire response to climate change 663 from the 2040s to 2070s (Figure 4 and Figure 5), here we found that the ∆fuel loading:∆ fuel meaning that increases in fuel aridity dominated over the decreases in fuel loading (Figure 11 a).

669
In grass ecosystems, fuel loading was the only factor that separated increases or decreases in  allocation of limited resources to reduce fire risks and societal vulnerability will need to take 731 these changes into account. Our modeling approach demonstrates that fire regimes will likely change differentially across watersheds, but there is still substantial uncertainty in the models. As scientists work to reduce key uncertainties-including future precipitation patterns, future decomposition rates, and species-specific responses to CO2 fertilization-future projections will 735 continue to be refined, providing better support to wildfire management decisions.

741
Code and data availability

742
The coupled RHESSys-WMFire model code is available online at:  Text S1. Model descriptions S1.1 Fire spread model WMFire is a stochastic fire-spread model designed to be coupled with RHESSys (Kennedy et al., 2017). It takes output variables from RHESSys and uses them to predict fire spread. Because it is an intermediate-complexity stochastic model, WMFire is not designed to predict the perimeters and timing of individual fires. Instead, the model can be used to predict aggregate spatial and temporal characteristics of fire spread across basins over time (i.e., the fire regime). A successful ignition occurs when there is an ignition source on the landscape, and it successfully starts a wildfire. The successful start of fire (Pi(l,d), Eqn 1) is calculated from the probabilities associate with litter load and relative deficit (Pi(l)) and Pi(d), respectively), where deficit is calculated as 1-ET/PET (here, we use relative deficit as a surrogate for fuel aridity). The probability of a successful spread (Ps(l,d,S,w)) is calculated based on the probabilities associate with litter load (l), relative deficit (d), topographic slope (S) and the orientation of spread relative to wind direction (W), giving Ps(l), Ps(d), Ps(S), and Ps(w), respectively (Eqn 2). The probability of spread (Ps) increases with increasing fuel load and relative deficit, is highest in the direction of wind, increases in the uphill direction, and decreases in the downhill direction. After a successful ignition, WMFire tests the orthogonal neighbors of that patch against the probability of spread to determine if there is a successful spread. WMFire model has demonstrated accuracy in the Santa Fe (New Mexico) and HJ Andrews (Oregon) watersheds (Kennedy et al., 2017).
Pi (l,d)=Pi (l)×Pi (d) Equation (S1) Ps (l,d,S,w)=Ps (l)×Ps (d)×Ps (w)×Ps (S) Equation (S2) S1.2 Fire effects model The fire effects model is built to match the complexity of the coupled RHESSys-WMFire model (Bart et al., 2020). The fire effects model uses the probability of spread (Ps) as an index of fire intensity but also accounts for canopy structure, which links fire intensity and spread with fire severity. Fire consumes C in the litter and coarse woody debris (CWD) pools based on the CONSUME model (an empirical model developed using statistical relationships derived from measured woody fuel consumption data; . Passive crown fire is spread from the understory to the overstory. For the understory, Ps is a proxy for fire intensity, and the firecaused mortality is a function of intensity (currently it is a 1:1 linear relationship, but can be changed to account for different landscapes). For the overstory, mortality and consumption of fuel are based on how much litter and understory fuel are consumed. There is also a set of parameters to account for the understory height threshold, overstory height threshold, and functional form representing the relationship between mortality and consumption. The fuel that is not consumed moves into the litter and CWD pools.

Text S2. Model parameterization
We used a Monte Carlo approach to calibrate six groundwater-related parameters: saturated hydraulic conductivity (Ksat), decay of Ksat with depth (m), air-entry pressure (φae), pore size index (b), bypass flow to deeper groundwater stores (gw1) and groundwater drainage rates to the stream (gw2). We selected the best parameter set for Trail Creek by comparing observed and modeled streamflow using the Nash-Sutcliffe efficiency metric (NSE), R2 for the correlation between daily observed and modeled flow, and percent error in annual flow estimates as well MODIS ET and NPP . The Monthly NSE is 0.94 with a percent error of 2.6 for calibration period. The model estimates ET and NPP in a reasonable range. A detailed description of hydrologic calibration is described by Ren et al., (2021). The fire spread model WMFire has previously been shown to replicate expected spatial patterns of fire spread for a wildfire in the Pacific Northwest (Kennedy et al., 2017) and fire regime characteristics for two different watersheds with contrasting historical fire regimes (LANDFIRE, ). The model is robust to applications in both historically low severity frequent-fire regimes in the southwest and mixed to high severity infrequent fire regimes in the Pacific Northwest. We selected three criteria (spatial distribution of fire spread, fire seasonality, fire return interval (FRI) to validate the fire spread model against LANDFIRE estimates for Trail Creek. Spatial distribution and seasonality of fire were not sensitive to ignition rates and agreed with LANDFIRE estimates. FRI, on the other hand, was sensitive to ignition rates. We adjusted ignition rates according to the size of Trail Creek, which enabled RHESSys-WMFire to simulate spatial variation in FRIs that agreed with LANDFIRE estimates. A detailed description of WMFire model calibration is described in (Hanan et al., 2021). Figure S1. The ranking of fire-related climate change variables among GCM models for Trail Creek according to the future greenhouse gas emission RCP8.5 scenarios. These variables include the 1971 -2000 vs. 2040 -2069 changes in the monthly mean of a) actual evaporation (∆AET) and b) water deficit (∆DEF), the standard deviation of monthly c) AET (σAET) and d) DEF (σDEF); annual number of days per year where 100-hour dead fuel moisture is below the historical (1971 -2000) value for the e) 3rd percentile and f) 10th percentile. ProDrought had the largest increase in the standard deviation of DEF from the historical (1971 -2000) period; this represents a storyline where drought would be promoted. ProVeg had the smallest increase in the mean of DEF, which would increase productivity and limit fire; and ProFire had significant variability in fire-related metrics; MultiMean is the GCM model that was closest to the multimodel mean. Figure S2. Basin scale fuel load (i.e., litter carbon) response to climate change and atmospheric CO2 fertilization effect under four different GCMs scenarios for 2040s. Baseline scenario used historical climate data for the future scenario spin-up. These scenarios are without fire model on.