• Developing land surface models for climate-fire interactions requires estimating and overcoming uncertainty in fuel accumulation processes • Models that simulate fuel accumulation differ in how they parameterize and represent fuel decomposition; assumptions are often hard coded • Sensitivity to parameter and model structure uncertainty increases with climate warming and decreases with increasing precipitation Supporting Information: Supporting Information may be found in the online version of this article.

litter and fine woody debris <7.6 cm in diameter; Table S1 in Supporting Information S1) and how it changes through time (i.e., fine fuel accumulation). Because fine surface fuel loading is a key driver of fire behavior, including its intensity, spread, and flame length (Kreye et al., 2013;McCaw et al., 2008;Rothermel, 1972;Sullivan et al., 2018;Thaxton & Platt, 2006;Weise et al., 2005), accurately projecting fine fuel accumulation is crucial for forecasting future fire regimes (Ren et al., 2021), optimizing fuel treatments, including their cost and longevity (Calkin & Gebert, 2006;Tinkham et al., 2016;Vaillant et al., 2015), and determining how wildfire will affect carbon (C) stocks across spatial and temporal scales (Campbell & Ager, 2013). Yet, projections of fuel accumulation and its contribution to wildfire remain highly uncertain (Dupuy et al., 2020;Varner & Keyes, 2009).

Fine fuel accumulation results from a balance between fuel inputs (i.e., productivity then phenology/mortality) and losses (i.e., combustion and decomposition), both of which are influenced by climate change (Figures 1 and 2; Agee & Huff, 1987). Land surface models, including those that developed from the CENTURY lineage (e.g., CLM, RHESSys, BiomeBGC (Lawrence et al., 2019;Tague & Band, 2004;Thornton et al., 2005), are capable of representing many of these processes and, in recent years, some land surface models have been adapted to incorporate fuel-wildfire feedbacks (e.g., LandClim; Gaillard et al., 2014;FireBGC;Keane et al., 2011;and RHESSys-WMFire;Kennedy et al., 2017). While such models provide a powerful tool for projecting climate effects on future fire regimes, adapting these models beyond their original domain to accommodate complex and nonstationary feedbacks among climate, vegetation, fuels, and wildfire is not a trivial task.

Figure 1. The parameters, processes, and state variables driving fire across spatial and temporal scales. This is an adaptation and extension of the conceptual figure developed by Moritz et al. (2005), which expanded the fire triangle concept to incorporate the feedbacks among fire drivers and processes at multiple scales, ranging from flames to fire regimes. Dominant drivers at each scale are identified along the sides of each triangle. Here, we illustrate the processes and feedbacks that are directly relevant to fine fuel accumulation, which controls fuel dynamics represented by the small green triangles at each scale. We use the term O-horizon to refer to litter (Oi horizon) and duff (Oe and Oa horizons). PFT stands for plant functional type and SA/V stands for the surface area to volume ratio.

Fine fuel accumulation is the balance between the input and removal of fuels (Figure 2; Table S1 in Supporting Information S1). The classic Olson (1963) model assumes that fuel accumulation is a function of the balance between the rate of fuel deposition and the rate at which it decays, represented as a simple curve of fuel density over time. However, in reality, fuel accumulation is constrained by a complex balance among several deposition and loss processes. For example, deposition is a function of net primary productivity (NPP), leaf and stem turnover, and vegetation mortality, including background mortality and pulses of mortality due to disturbances. Fuel losses can occur through decomposition, combustion, erosion, and herbivory. To further complicate this balance, climate change and wildfire can alter both deposition and losses at multiple spatial and temporal scales.

Many studies have focused on how decomposition rates influence the net exchange of C between ecosystems and the atmosphere (net ecosystem exchange; NEE; e.g., Kramer et al., 2017;Melillo et al., 1982;Schlesinger & Andrews, 2000), or on how decomposition influences nutrient cycling and NPP (Lal, 2004), but not on how it may modify fine surface fuel loading, fire spread, and associated feedbacks with greenhouse gas fluxes. Harris et al. (2016) reviewed many of the input processes controlling fuel deposition and their effects on fire regimes. Here, we focus on the role of decomposition in driving fine fuel accumulation and the fundamental equations used to represent it. Decomposition is expected to accelerate under future warming (Hopkins et al., 2012), but its response to increasing temperature and drought remains highly uncertain.

There are often inconsistencies between fuel layers defined in land surface models (e.g., litter and woody debris) and those defined in fire behavior models. To clarify (Table S1 in Supporting Information S1), we refer to the litter layer as comprised of fallen dead leaves and needles. Woody fuels (WD) include fine woody debris size classes (<7.6 cm, e.g., 1-, 10-, and 100-hr fuels; FWD) and coarse woody debris (>7.6 cm in diameter, e.g., 1,000-hr and greater; CWD). Fine surface fuels then comprise litter and fine woody fuels.

Although coarse and fine surface fuels categories are represented at varying levels of detail in fire regime models (Figure 3), most of our current theoretical understanding and models of decomposition are based on work involving litter and soil organic matter (SOM; Keane et al., 2011;Tague & Band, 2004). In most land surface models, this theoretical understanding is applied to simulate the decomposition of dead plant leaf litter. Woody fuel loss, on the other hand, is modeled using either a constant scalar (e.g., He et al., 1999;Rebain et al., 2015) or by including both fragmentation to smaller size classes and the estimated litter decomposition rates (e.g., Keane et al., 2011). Given that the parameters that define litter decomposition are also used in WD decomposition (Figure 3), we focus on how models of decomposition contribute to uncertainty in modeling litter mass loss, and how this uncertainty propagates to projections of climate-related changes to future fire regimes. Any model uncertainty that we identify for litter decomposition would also be present for WD decomposition. This analysis does not consider whether it is appropriate to use the same model structures for litter and WD decomposition, which is an important subject for future research.

Decomposition is controlled by three overarching factors: (a) environmental conditions, particularly temperature and moisture, (b) the amount and quality of substrate available for decomposers, and (c) microbial community structure and function (Chapin et al., 2011;Melillo et al., 1982). C cycling models then incorporate these factors into an exponential decay-like model. Environmental constraints on decomposition are represented by scalars associated with moisture (r w ) and temperature (r T ). A rate constant (k) that is empirically estimated and primarily associated with litter quality is typically used, either as a single value and/or as multiple sequential pools (l1, l2, l3, l4) that are increasingly recalcitrant (k varies by pool). This implicitly accounts for the role of microbes using a single, static parameter. Overall this representation uses first-order kinetics, meaning microbial decomposition processes are modeled using a single, first-order equation that is controlled by the size of each C pool (Parton et al., 1998;Running & Coughlan, 1988;Tague & Band, 2004;Thornton et al., 2005). The first-order model can then be written as

Figure 3. Representation of fine surface fuels in three example models (RHESSys, LandClim, FireBGCv2). Litter may be represented as a single pool or multiple pools (l1, l2, l3, l4). Woody debris (WD) may also be in a single pool, partitioned into fine (FWD) and coarse size classes (CWD), or the FWD may be further portioned into 1-, 10-, and 100-hr fuels, with CWD representing ≥1,000-hr + fuels. Decay rates for each of these pools may depend on environmental factors (temperature r T ; moisture r w ), litter quality (kl), or size class (k frag , s w ). Colored text in equations shows where common elements are used across fuel size classes. Equation details are given in Section 3.

Next, we describe the current understanding of processes underlying environmental constraints (e.g., r w and r T ), litter quality (e.g., k), and microbial dynamics.

Physical environmental conditions in an ecosystem or landscape influence decomposition in large part through their effects on temperature and moisture. Temperature regulates decomposition directly through its effects on soil microbial activity and indirectly through its effects on litter and soil moisture. Increasing temperatures increase microbial respiration rates exponentially across biomes. For example, in warm tropical forests, litter pools are small despite high rates of NPP, whereas in temperate coniferous forests, litter pools can be large even though NPP is much slower (Chapin et al., 2011;Lieth, 1975). Because temperature affects NPP and decomposition at different rates (Kirschbaum, 1995), it is crucial to understand mechanistic relationships between warming and litter decay to accurately predict fine fuel accumulation.

Traditionally, C cycling models have used empirically fitted temperature sensitivity functions (i.e., Q10) to describe how decomposition rates increase with warming (e.g., Davidson et al., 2006;Luo et al., 2001;Reichstein et al., 2003). Q10 is a measure of the extent to which 10°C rise in the temperature increases the rate of a chemical reaction. However, fitting Q10 functions to soil respiration data has yielded highly variable temperature sensitivities (Davidson et al., 2006). For example, Q10 can vary with season (Janssens & Pilegaard, 2003), soil organic matter content and quality (Reichstein et al., 2005), soil moisture (Meyer et al., 2018), land cover (Yuste et al., 2004), elevation (Wang et al., 2013), and latitude (Zhou et al., 2009). Modeling the effects of temperature on decomposition is extremely difficult because these environmental constraints can obscure the intrinsic temperature sensitivities of various substrates, and these constraints may themselves be sensitive to climate (Davidson & Janssens, 2006).

One of the biggest constraints on decomposition is moisture availability. Similar to plants, decomposers are most productive in warm moist environments where they are neither oxygen nor diffusion-limited. However, soil microbes are also less sensitive than plants are to drought (Hanan et al., 2017;Jackson et al., 1988;Parker & Schimel, 2011). Thus, in relatively xeric locations, warming and drying may decrease NPP and fine surface fuel inputs while increasing decomposition, thereby reducing fuel loadings and fire hazard (Batllori et al., 2013;Kennedy, Bart, et al., 2021;Kennedy, Johnson, et al., 2021). Furthermore, drying-rewetting cycles may become more frequent with climate change and can stimulate decomposition of labile substrates while slowing rates for recalcitrant ones (Haynes, 1986).

Temperature and moisture respond to both top-down climate drivers and bottom-up environmental drivers-such as topography, soil properties, and vegetation cover-and they influence decomposition both directly and indirectly. Therefore, wildfire modeling requires predicting future temperature and moisture regimes not only for their direct effect on wildfire behavior and spread but also how they will interact to drive fine fuel accumulation (Figure 1). Despite the clear need to account for temperature and moisture variability in C cycling models, there are several uncertainties that still must be resolved for future projections to be reliable. For example, the extent of future drought remains highly uncertain (Cook et al., 2020). While it is clear that temperatures and evapotranspiration (ET) will continue to increase, future precipitation is less predictable and thus for ecosystems that exist near the threshold of flammability to fuel-limitation, improved projections of future aridity will be extremely valuable for predicting fire hazard (Hanan et al., 2021).

Another limitation to modeling the effects of future aridity on decomposition comes from uncertainty in a model structure. Models that represent moisture controls on decomposition tend to focus more on soil moisture than litter moisture. For example, in RHESSys-WMFire and FATES-SPITFIRE, the moisture controls influencing fine fuel decomposition are based on soil water content and soil matric potential, respectively (Andren & Paustian, 1987;Tague & Band, 2004;Thonicke et al., 2010). Similarly, the moisture controls influencing decomposition in LANDCLIM are a function of evapotranspiration (ET; Gaillard et al., 2014). These variables do not always operate on the same timescales as fine fuel moisture (Hatton et al., 1988). However, modeling fine fuel moisture is challenging because there are limited empirical studies elucidating the mechanisms driving water adsorption by plant litter. In one study, Talhelm and Smith (2018) observed relationships between water adsorption and the structure and chemistry of leaf litter. Notably, it was shown that litter with high concentrations of heat content and lignin exhibited lower water adsorption (Talhelm & Smith, 2018). Further research is needed to improve models of fine fuel moisture and its effects on decomposition.

Finally, temperature and moisture can interact in complex ways, and these interactions may not be multiplicative, which can lead to possible equifinality when attempting to estimate their individual contributions through lab experiments (Tang & Riley, 2020). This is evident when comparing historical and future projections for different C cycling models. In many cases, C cycling models can have convergent projections over the historical period and highly divergent projections in the future (Z. Luo et al., 2015). We know this is problematic for slow cycling soil C stores, but it has not been tested extensively for litter/fine surface fuels.

At a given temperature and moisture regime, decomposition rates can vary by several orders of magnitude due to differences in litter quality (Silver & Miya, 2001). Litter quality refers to the relative proportions of labile metabolic compounds in litter stores, such as sugars, amino acids, moderately labile compounds such as cellulose and hemicellulose, and recalcitrant compounds such as lignin (Chapin et al., 2011). Two common indices for litter quality are its C:N ratio and its lignin:N ratio (B. R. Taylor, 1989). Litter with relatively high N tends to be composed of more labile C compounds and less structural material and will therefore decompose more quickly (Hobbie, 2000;Melillo et al., 1982). Litter quality also decreases rapidly with age because labile materials decompose quickly. Belowground resource availability is a key factor influencing litter quality-vegetation in high resource sites produces litter that decomposes quickly because the physiological traits that lead to high NPP, such as high surface to volume ratio and low C:N, also tend to favor rapid decomposition.

In multipool, first-order models, each litter pool has a single k value that is static through time and a single set of temperature and moisture reducing functions (Georgiou et al., 2017, Equation 1). An implicit assumption in first-order models is that the response functions do not change with the composition or size of the microbial community (Schimel, 2001). Thus, first-order models do not account for how such changes can modify decomposition rates. Research over recent decades, however, has shown that first-order model assumptions can be problematic, particularly for slow cycling soil C pools, which can experience accelerated decomposition when inoculated with heterotrophic microbes (Z. Luo et al., 2015) or which can experience decelerated decomposition when sorptive mineral surfaces in the soil interact with substrates, enzymes, and microbes (Tang & Riley, 2015). As a result, first-order models are potentially inadequate for representing processes such as priming, where the decomposition of organic C can be enhanced through plant root exudates or elevated CO 2 concentrations that stimulate the heterotrophic microbial community (Hungate et al., 1997).

More recently, models have attempted to capture the role of soil microbes in mediating decomposition and/or sorption of enzymes and organic matter to mineral surfaces (e.g., Hararuk et al., 2015;Kaiser et al., 2014;Wieder et al., 2013). Such models may explicitly represent enzymatic degradation of soil and litter C (i.e., through Michaelis & Menton, [1913] kinetics). In these models, decomposition rates depend on the sizes of both C and microbial pools. While such models may be needed to simulate the decomposition of recalcitrant soil organic matter pools at regional scales (Todd-Brown et al., 2013), they have not been tested in the context of fine surface fuels and wildfire. Furthermore, wildfire can dramatically reduce microbial biomass (e.g., Goodridge et al., 2018;Hanan, D'Antonio, et al., 2016;Hanan, Schimel, et al., 2016;Knicker, 2007) and alter the microbial function and enzyme activity over decadal timescales (Pellegrini et al., 2020). These feedbacks are also poorly represented in biogeochemical models.

Fire models range in complexity from simple empirical models that can be used to classify large-scale fire regimes (e.g., Littell et al., 2018) to fully physical models that have the potential to predict individual wildfire spread with precision (e.g., Mell et al., 2007). However, many of the existing fire models-at all levels of complexity-inadequately represent the full system of interactions and feedbacks that influence fine fuel loadings (Figure 1). For example, empirical models have shown how climate change is increasing regional-scale fire occurrence and area burned (Abatzoglou & Williams, 2016;Guyette et al., 2012;Littell et al., 2018;McKenzie & Littell, 2017), but these models do not account for dynamic vegetation and fuel-fire feedbacks, such as wildfire self-limitation. More complex fire behavior models often rely on generalized fuel schemes, such as the Scott and Burgan (2005) 40 stylized fuel models, the Australian Bushfire Fuel Classification (Cruz et al., 2018), or the Canadian Forest Fire Behavior Prediction System (Forestry Canada, 1992). However, many of these models do not accommodate dynamic fuel loadings and therefore ignore novel fuel beds that can arise from climate change, land management, and feedbacks with past fires. Because the relationships among climate, fuels, and wildfire that we observe today may not persist in the future, there are important limitations in applying empirical and physical fire models to projecting future fire regimes (Newman et al., 2019).

Land surface models, on the other hand, have the potential to simulate dynamic fuel loadings through bidirectional couplings among climate, vegetation, and fuels. Such models do not typically include the level of physical detail used for modeling individual fire spread but instead represent processes and feedbacks driving long-term fire regimes. In this paper, we are concerned with how the fundamental decomposition mechanisms outlined in the previous section are incorporated into land surface models that have been adapted to include fully coupled feedbacks with wildfire. These fire regime models can be used to investigate how climate change will alter future wildfire regimes, which is one of the most pressing questions in global change and forest management. We focus on three exemplar landscape-scale North American models (FireBGCv2, RHESSys-WMFire, LandClim) to demonstrate uncertainties that pertain to models across a wide range of scales. Below, we provide general model descriptions, then we detail the decomposition routines of these three models. We chose these three models to explore the consequences of underlying model assumptions that are commonly found in fire regime models (Figure 3). RHESSys-WMFire and FireBGCv2 have very similar decomposition model structures, with a few key differences. LandClim takes an entirely different approach, relying on relatively simple empirical relationships that deviate from the more common framing of Equation 1.

FireBGCv2, RHESSys-WMFire, and LandClim simulate how interacting ecosystem processes pertaining to climate, vegetation, soils, hydrology, and disturbance influence C fluxes (Gaillard et al., 2014;Keane et al., 2011;Kennedy et al., 2017). However, they differ in the set of processes they emphasize and in the scales that they represent.

FireBGCv2 is adapted from BIOME-BGC to represent individual-tree-based succession and wildfire (Keane et al., 2011). It merges biogeochemical processes from Biome-BGC (Running & Coughlan, 1988) with the FIRE-SUM gap model (Keane et al., 1989). FireBGCv2 operates at five distinct spatial scales, ranging from individual trees to entire landscapes and operates on a daily time step. Physiological processes such as photosynthesis, respiration, and decomposition are calculated at the finest scales, whereas fire is implemented stochastically at a landscape scale.

RHESSys-WMFire is unique in that it fully couples a biogeochemical model with a spatially explicit hydrologic model to simulate processes such as streamflow, evapotranspiration, NPP, respiration, mineralization, nitrification, and C and N export to streams (Tague & Band, 2004). Most processes are modeled at a patch scale, which typically varies between 30-m and 270-m resolution. Subsurface and surface water are routed laterally between patches within subbasins to produce streamflow. The largest spatial unit is the basin, which aggregates subbasins and is a closed drainage area encompassing a single stream network. Like, FireBGCv2, RHESSys-WMFire also operates at a daily time step.

LandClim is a spatially explicit, stochastic landscape model that developed from LANDIS to incorporate largescale disturbances such as fire and feedbacks with climate change (Gaillard et al., 2014;He et al., 1999). Land-Clim represents stand-scale (i.e., 25-m) vegetation as the number and biomass of trees in cohorts. Processes such as growth and mortality are simulated at an annual time step, and landscape-scale processes, such as fire, wind, and seed dispersal are simulated at a decadal time step (Gaillard et al., 2014).

These models also differ in the degree of complexity they use to represent fire (Figure 3). Both FireBGCv2 and LandClim simulate ignition and spread based on moisture, wind, and topography, given a threshold value for fuel presence. FireBGCv2 scales the probability of spread by a user-specified fire return interval, which is a surrogate for fuel accumulation that does not respond to changing climate and vegetation conditions. Fire behavior in FireBGCv2 is based on either Rothermel (1972) or Albini (1976) equations, which depend on intrinsic fuel properties and on fuel loading of different size classes. Fire effects are calculated using the FOFEM model (Reinhardt et al., 2001).

LandClim simulates random ignitions and spread is based on woody fuel presence and a drought index. It calculates fire intensity as a function of fuel load and moisture (Schumacher et al., 2006). Fire size in both LandClim and FireBGCv2 is limited by a user-specified maximum, a common simplification in many fire models to ensure that wildfire spread (e.g., by a specified maximum duration of spread or a maximum fire size) is limited. In such representations, the effects of fine fuel accumulation on wildfire area burned and feedbacks with wildfire activity may not be emergent from model projections, even though accounting for fire self-limitation can moderate projections of future area burned (Abatzoglou et al., 2021;Hurteau et al., 2019).

RHESSys-WMFire produces fire spread maps over randomized ignitions and stochastic spread, providing probability distributions of fire activity over time. In addition to topography, wind, and climate (as in LandClim and FireBGCv2), fire spread and effects also respond to dynamic changes in fuel loading (Bart et al., 2020;Kennedy et al., 2017); RHESSys-WMFire is therefore robust to climate nonstationarity and the positive and negative feedbacks that influence fuel dynamics fire regimes over time (Hanan et al., 2021).

The basic equations governing daily decomposition and mass loss of litter and woody fuels are similar between RHESSys-WMFire and FireBGCv2, but with key differences in how they represent C pools and the temperature and decomposition scalars. Fuels in both models are represented by litter and coarse woody debris pools. Litter is partitioned into multiple pools, such as labile C, cellulose, shielded cellulose, and lignin/remaining mass. In RHESSys-WMFire, woody fuels are aggregated into a CWD pool. FireBGCv2 further partitions the CWD pool into size classes (1-, 10-, 100-, and 1,000-hr fuels; Figure 3).

Both RHESSys-WMFire and FireBGCv2 use scalars to account for temperature (r T ) and moisture (r w ) limitations on decomposition. Woody debris mass loss includes those scalars and a fragmentation constant (k frag ). For example, in RHESSys-WMFire, the woody debris pool is fragmented into the litter pool at the following rate (r wd ):

In FireBGCv2, fragmentation differs by size class and is modified by a size class scalar (s w ). C in each size class fragments to the next lower size class. For example, fragmentation from 1000-hr to 100-hr fuels occurs at a rate (r wd.1000 ):

𝑟𝑟wd.1000 = 𝑘𝑘frag𝑟𝑟𝑚𝑚𝑟𝑟𝑇𝑇 𝑠𝑠𝑤𝑤𝑤1000

(3)

Litter in each pool (l1, l2, l3, and l4) then decomposes at the following rate (r l.i ):

where kl differs by litter pool. RHESSys-WMFire also includes a third scalar to represent nitrogen limitation by calculating the fraction of potential nitrogen immobilization (f; Tague & Band, 2004), so that the final decomposition rate is

Note that, within each model system, the temperature (r T ) and moisture (r w ) scalars are the same in the series of equations governing woody fuel and litter loss (Figure 3). Both FireBGCv2 and RHESSys-WMFire share the same temperature scalar, which depends nonlinearly on the soil temperature (T soil ):

10.1029/2021MS002818 9 of 22 𝑟𝑟𝑇𝑇 = 𝑒𝑒 308.56 * ( 1 71.02 -1 𝑇𝑇 soil +273.15-227.13 ) (6)

They do, however, have a key difference in their moisture scalars. The moisture scalar in FireBGCv2 (r w.B ) is calculated as in Biome-BGC (Thornton, 1998) and depends on the soil water potential ( 𝐴𝐴 𝐴𝐴 ) relative to the range of possible soil water potentials (min, max):

The RHESSys-WMFire moisture scalar (r w.R ) follows a modified NGAS model, where 𝐴𝐴 𝐴𝐴 is soil water content (Parton et al., 1996).

LANDCLIM follows a model lineage distinct from FireBGCv2 and RHESSys-WMFire. It also represents fuels as litter and woody debris but partitions its woody debris into 2 size classes (Figure 3): fine (<7.6 cm) and coarse (>7.6 cm). This is in contrast to the single pool in RHESSys-WMFire and the four pools in FireBGCv2. The decomposition rate of woody debris in LandClim is based on the empirical model of Mackensen et al. (2003), who fit a curve between coarse woody decomposition rate and temperature, based on a cross section of different studies and locations. LandClim then applies a scalar to this equation so that fine wood is lost 5x faster of course wood (Schumacher et al., 2006). The general relationship between fine wood decomposition rate (r w ) and air temperature (T a ) is

The scalar m is the increase in the decomposition rate for fine wood (m = 5); d 1 and d 2 are the Mackensen et al. (2003) estimated coefficients (0.0166 and 0.093, respectively).

Calculation of the litter decomposition rate in LandClim is achieved using an empirical regression equation estimated by Meentemeyer (1978) using data from multiple sources to estimate general relationships between annual foliage litter decomposition rate (r l ), annual actual evapotranspiration (AET), and percent lignin. The best fit synthesis model for foliage litter decomposition rate explained 70% of the variability and included AET as a main effect and an interaction between AET and lignin (represented by the ratio AET/lignin):

The models described above and larger-scale global fire models (e.g., FATES-SPITFIRE; Thonicke et al., 2010) are in part adaptations of existing models that were not originally developed to simulate wildfire (Hantson et al., 2016;Rabin et al., 2017). For example, decomposition in RHESSys-WMFire and the ELM or CLM component of FATES-SPITFIRE grew out of algorithms developed for the biogeochemical model CENTURY, which was originally designed to simulate biogeochemical fluxes in natural and managed ecosystems (Holm et al., 2020;Koven et al., 2020;Parton, 1996;Tague & Band, 2004). However, there has not been a detailed assessment or validation of their prediction of surface dead biomass, which can play an important role in projected wildfire activity.

Next, we demonstrate how the potential parameter and model structure uncertainty in decomposition rate calculations can propagate to uncertainty in fine surface fuel loading. We start with an exploration of sensitivity to parameter estimation uncertainty using the empirical fine woody fuel decomposition model used in LandClim.

We then compare predicted litter decomposition rates among the three model structures described above. We chose three models as examples of current state-of-the-art fire regime models, not to imply that these models are particularly problematic in this regard, but rather to illustrate potential uncertainties that apply to many similar models of this type. Methods for the sensitivity analysis are detailed in Supplementary Text (Sections S1 and S2 in Supporting Information S1).

To explore potential uncertainty in predicted decomposition rates due to uncertainty in parameter (coefficient) estimates, we use the LandClim equation for fine woody fuel decomposition. As empirical (regression) estimates, each coefficient in Equation 9 has an associated standard error. The curve that was estimated explained 34% of the variability in the decomposition rate, and there was increasing variability in the decomposition rate as the temperature increased (Mackensen et al., 2003). This might be of particular concern in climate scenarios with the increasing temperature. At the maximum temperature of 25°C, observed decomposition rates for coarse wood varied from ∼0 to ∼0.6. Given the fine wood multiplier (m) of 5 in LandClim, this would propagate to fine wood decomposition rates of around 0 to 3.0.

To explore the consequences of uncertainty in coefficient estimates on decomposition rates and fuel loadings, we conducted a simple sensitivity analysis (SA) by systematically varying coefficient values in the underlying equation (m, d 1 , d 2 ; Equation 9), decoupled from other model processes. Unfortunately, standard errors were not given in the source material, making it difficult to determine plausible bounds of uncertainty. Therefore, we evaluated ranges of coefficients +/-33% for the empirical estimates, which yielded a range of decomposition rates at 25 C that were similar to those in the data used to estimate the curve, thereby characterizing the expected variability in the decomposition rate at that temperature. We recorded both decomposition rates (Figure 4, top) and percent of initial fuel loading remaining assuming no fuel inputs (Figure 4, bottom). The latter is akin to a common litter bag experiment, where a portion of litter is left in a given environment and its mass tracked over time.

Given that the relationship between temperature and the rate of fine wood fuel decomposition is exponential, the sensitivity of that relationship to the exponent must also be nonlinear (Figure 4, top; Text S1 in Supporting Information S1), as is the effect on future woody fuel loading (Figure 4, bottom). The sensitivity of decomposition rates to model coefficients increases with increasing temperature, with the widest uncertainty bounds at the highest temperatures. Thus, as the climate warms, we can expect increased decomposition rate parameter uncertainty.

We used RHESSys-WMFire as a base model to investigate the sensitivity of litter decomposition and projected wildfire to model structure uncertainty. We programed the litter decomposition models from FireBGCv2 and LandClim in RHESSys-WMFire and then conducted three sets of simulations, one with each decomposition model. These simulations were identical except for the decomposition model. This analysis provides an example of how decomposition model structure can modify future projections for a single model system (RHESSys-WM-Fire). We also considered two climate change scenarios using GCM data from the Coupled Model Inter-comparison Project 5 (CMIP 5; K. E. Taylor et al., 2012) that have been statistically downscaled across using the Multivariate Adaptive Constructed Analogs (MACA, Abatzoglou, 2013;Abatzoglou & Brown, 2012; Text S1 in Supporting Information S1). These included a hotter and drier future (ProDrought) and a hotter and wetter future (ProVeg).

For each climate change scenario, we simulated associated baseline conditions (Text S1 in Supporting Information S1). This enabled us to consider short-term (2040s) and long-term (2070s) futures against a baseline for each time period. The baseline scenarios experienced wildfire and decomposition for the same amount of time as the climate change scenarios and thus differences from baseline could be attributed to climate change. We also conducted a simulated litter bag experiment, as if a litter bag was left in the watershed over each of the time periods and for each of the scenarios. From that experiment, we calculated litter mass loss over 15 years for each climate scenario.

We found large differences in decomposition rates among the three model structures (Figure 5). RHESSys-WM-Fire consistently predicted higher decomposition rates than both FireBGCv2 and LandClim. LandClim predicted the lowest decomposition rates. Both RHESSys-WMFire and FireBGCv2 predicted increasing decomposition over time in the climate change scenarios (Figure 5), largely associated with increasing temperature. FireBGCv2 had the most variability in projected decomposition rates, responding to interannual variability in moisture associated with precipitation. RHESSys-WMFire was less sensitive to changing precipitation than FireBGCv2.

Decomposition rate differences associated with each model depended on precipitation and temperature (Figure 6). The largest differences between FireBGCv2 and RHESSys-WMFire occurred at the lowest precipitation and highest temperature, likely due to the greater sensitivity to moisture limitation shown by FireBGCv2. The largest differences between LandClim and the other two models occurred at the highest temperatures, because decomposition was not sensitive to temperature in LandClim. LandClim had the smallest decomposition rates at the highest temperatures.

Differences in decomposition rate propagated to litter mass loss across the three models and climate change scenarios (Figures 7 and 8; Figure S1 in Supporting Information S1). We found that model structure uncertainty was inconsistent between historical and future scenarios and also differed between dry vs. wet scenarios. In general, litter mass loss was fastest for RHESSys-WMFire, leading to fuel limitation and the lowest area burned. LandClim on the other hand had the slowest litter mass loss and greatest area burned. RHESSys-WMFire and FireBGCv2 both projected increases in litter mass loss under the climate change scenarios (Figure 7).

Different model sensitivities of decomposition rates to changes in temperature and precipitation also propagated to sometimes large differences in projected wildfire area burned (Figure 8). Relative to baseline for all model 9. The middle (solid) curve is the hardcoded value in the model. The "multiplication factor" for this curve is 0.0166 from Equation 3 multiplied by 5. The upper and lower curves illustrate how projected decomposition rates might vary if the components of the multiplication coefficient each increased or decreased by 33%. The green line indicates the point where the fraction lost is 1, or in other words where inputs are equal to outputs. The relationship between the decomposition rate exponent and decomposition rates are shown in separate panels for 4 different mean annual temperatures (that is, 5-, 15-, and 25-°C. Model sensitivity to parameter uncertainty increases with increasing temperature. Bottom: Sensitivity of LandClim percent of fine woody fuel remaining at year 10 to the coefficient values in Equation 10. The middle line is the hard-coded value in the model. The upper and lower lines illustrate how projected decomposition rates might vary if the components of the multiplication coefficient each increased or decreased by 33%. Model sensitivity to parameter uncertainty increases with increasing temperature.

structures, this was associated with increases in annual area burned in the near future (i.e., the 2040s) for the wet future scenario (ProVeg) and similar annual area burned in the dry scenario (ProDrought; Figure 8). In the distant future (2070s), patterns interacted with the decomposition model structure. For FireBGCv2 and RHESSys-WM-Fire, the burned area remained similar to baseline for the ProVeg scenario and decreased from baseline in the ProDrought Scenario. For LandClim, annual area burned increased from baseline much more substantially for all climate change scenarios than it did in the other two models, and the increase was largest in the 2070s and under drought scenarios (Figure 8). RHESSys-WMFire generally predicted lower area burned than FireBGCv2, although the differences were smaller than those from LandClim.

As Earth's climate continues to change, we need insights from both experiments and models to understand how fine surface fuel loading and its properties will vary over space and time, and how they will affect fire behavior and fire regimes. Land surface models are valuable tools that, when fully coupled with models of fire behavior and spread, provide projections of fire regimes that incorporate feedbacks among climate, vegetation, fuels, and wildfire. It is crucial to understand uncertainty in projections of wildfire regimes associated with the coupling of land surface and fire behavior models.

There are a vast number of fire models in existence, including empirical, mechanistic, stochastic, and various combinations of the three (Rabin et al., 2017;Reinhardt et al., 2001;Sullivan, 2009aSullivan, , 2009b;;Sullivan & Sullivan, 2009). These models are designed to target different spatial and temporal scales of fire forecasting, ranging from the physics of individual flames to fire regimes (Figure 1; Harris et al., 2016;Keane et al., 2004). For example, empirical models can provide valuable insight into climate-wildfire relationships at regional scales (e.g., Abatzoglou & Williams, 2016;Guyette et al., 2012;Littell et al., 2018;McKenzie and Littell, 2017). However, because these models do not explicitly represent fuel dynamics, projecting climate-wildfire relationships into the future would implicitly assume that vegetation and fuels will remain stationary. Empirical resource gradient models have addressed this constraint by using climate indices to identify the relative roles of fuel moisture vs. fuel loading on wildfire activity (Krawchuk & Moritz, 2011;Mann et al., 2016;Parisien et al., 2016). However, 0 25 50 75 100 0 2 4 6 8 10 12 Years Percent of litter pool remaining Baseline A. 0 25 50 75 100 0 2 4 6 8 10 12 Years Climate change B. proDrought 0 25 50 75 100 0 2 4 6 8 10 12 Years Percent of litter pool remaining Baseline C. 0 25 50 75 100 0 2 4 6 8 10 12 Years Climate change D.

RHESSys-WMFire LandClim proVeg because empirical models rely on pattern-matching, they still have limited utility in projecting future wildfire under novel climate and fuel bed conditions (McKenzie & Perera, 2015).

While fire spread models can provide more mechanistic representations of the relationships between fuel and fire behavior, they also have limitations under climate change. Such models typically rely on stylized fuel classifications that require calibration to obtain expected fire behavior for a given fuel and vegetation type (Cruz & Fernandes, 2008). These classifications can sometimes be dynamic (e.g., depending on predicted stand conditions as in FFE-FVS; Rebain et al., 2015); however, they typically do not represent novel fuel beds that arise from climate change-driven shifts in decomposition or fuel treatments (Johnson et al., 2011;Kennedy, Bart, et al., 2021;Kennedy, Johnson, et al., 2021;Varner & Keyes, 2009). Data assimilation models, on the other hand, have the potential to account for changing properties of novel systems (Y. Luo & Schuur, 2020). Such models can enable researchers to disentangle specific drivers of wildfire activity (e.g., Kelley et al., 2021); however, observational data on fuel properties must be available.

Alternatively, land surface models are designed to simulate climate change effects on vegetation productivity and decomposition and many of these models also include algorithms for prescribing fire effects (e.g., CENTURY/ DAYCENT; Parton, 1996;Parton et al., 1998), BIOME-BGC; Thornton et al., 2005, and RHESSys;Tague and Band, 2004). However, in these model systems, fire may be applied as an exogenous driver and is not always represented as an emergent property of the fuel landscape. Although such model systems can still provide a powerful framework for mechanistically simulating climate-vegetation feedbacks following fire (Gathany & Burke, 2012;Hanan et al., 2017Hanan et al., , 2018;;Hudiburg et al., 2017), they do not include fuel-fire feedbacks that are needed to simulate decadal-scale fire regimes.

There are many models that do incorporate bidirectional couplings to represent climate-fuel-fire relationships, many of which are reviewed and classified by Keane et al. (2004). In these models, climate, vegetation, and dynamic fuels inform wildfire spread, behavior, and effects using varying degrees of abstraction for the system of feedbacks represented in Figure 1. These models include mechanistic representations of fuel moisture and fuel loading, which support applications under climate change scenarios. However, there are still large uncertainties in how these models represent fuels.

Here, we examined two types of model uncertainty (parameter and model structure) in three state-of-the-art fire regime models (LandClim, FireBGCv2, and RHESSys). We found that the sensitivity of projected decomposition to both types of uncertainty increased with climate warming and decreased with increasing precipitation (Figures 456). The sensitivity of decomposition to the model structure was highest under hot, dry conditions and lowest under cool, wet conditions (Figure 6). Additionally, in LandClim, temperature is not included as a direct driver of decomposition and therefore differences in decomposition projections between LandClim and other models increased with warming (Figure 6).

Previous studies focused on SOM pools have found that the temperature and moisture sensitivities of decomposition can vary over space and time, interact in complex ways, and these interactions may not be multiplicative (Dijkstra et al., 2011;Steinweg et al., 2008). We found that in FireBGCv2 and RHESSys-WMFire, the sensitivity of decomposition to temperature and precipitation can also interact. For example, over time in the warming scenario, RHESSys-WMFire and FireBGCv2 decomposition rates were closer on average with increasing temperatures (Figures 5 and 6). But there was substantial interannual variability in FireBGCv2 decomposition rates that tracked with precipitation (Figure 5). This led to periodically large differences in decomposition rates between the two models. Such interactions can lead to equifinality (i.e., that a given end state can be reached by multiple paths) when developing the model structure and parameterizations from lab experiments (Tang & Riley, 2020), which can be problematic when projecting future fire regimes under novel climates.

Another source of uncertainty comes from the representation of fuels themselves. For example, some fire models that managers use for forest planning (e.g., FARSITE and BEHAVE; Andrews, 2007;Finney, 1998) only include the woody fuels. A prevailing challenge is that woody fuel decomposition and the interactions with fire are not well studied (Hyde et al., 2011(Hyde et al., , 2012)), in part because measuring mass loss of coarse woody fuels (due to fire or decomposition) can be challenging (Fry et al., 2018). When woody fuel decomposition is incorporated in models, it is often based on a constant value (e.g., FFE-FVS; Rebain et al., 2015), or a value adapted from litter models (e.g., FireBGCv2; Keane et al., 2011). Thus, in many models, uncertainty in the decomposition rate propagates to uncertainty in the more "management-relevant" fuel layers. Other challenges that arise with modeling climatefuel-fire feedbacks include the incorporation of processes such as snag-fall rates (Kennedy, Bart, et al., 2021;Kennedy, Johnson, et al., 2021), snag decomposition (Stenzel et al., 2019), and delayed litterfall from scorched trees that otherwise survive fires (Espinosa et al., 2018;Keane, 2008). To improve fire management in the future, we need to not only improve our models of litter decomposition, we also must develop better theories and models for the controls on fine woody fuel deposition.

Process-based fire regime models can account for feedbacks among climate, fuels, and wildfire (Figure 1), which enables us to evaluate how fire regimes and fire effects will be transformed in response to climate change and management actions. However, to appropriately account for such feedbacks, we need to evaluate and improve our understanding of the fundamental processes and parameters we use to simulate fine fuel accumulation. We described several uncertainties in the model structure and parameters used to represent decomposition, which may lead to large uncertainties in projecting future fire under climate change.

For process models to be reliable, they must be continually confronted with observations and empirical data, including data for parameterization, validation, evaluating uncertainty, and improving the way we represent various mechanisms. Empirical studies can help improve our representation of litter turnover but there are disconnects between our empirical understanding and ability to model processes over fire-relevant scales. These disconnects arise because empirical studies typically focus on individual scales and rarely account for feedbacks that occur across scales-such as the effects of climate change on the microbial processes regulating fine fuel decomposition, its subsequent effects on fire, and feedbacks to soil biogeochemical processes (Figure 1). Understanding these complex climate-fuel-fire feedbacks is critical for land surface models that forecast future fire regimes.

To refine our modeling approaches, future research should (a) leverage existing and continue to implement new long-term monitoring studies of fine fuel accumulation and compare model predictions to observed, (b) quantify and understand fuel accumulation-related parameter and model structural uncertainty, and (c) consider fuel dynamics and feedbacks when assessing climate-wildfire relationships.

Even though decomposition is a key component of landscape, regional, and global C budgets, litter decomposition in land surface and Earth system models has not been thoroughly evaluated and most studies have focused on soil organic C stores rather than fine surface fuel loading. Long-term empirical studies can assist with model evaluation. For example, the long-term intersite decomposition experiment (LIDET; Harmon, 2016) provided a 10-year study of litter decomposition at multiple locations across North and Central America. These data have been used to evaluate simulated litter mass loss for different models. For example, Bonan et al., 2013 used LIDET data to constrain temperature and moisture effects on decomposition in the community land model version 4 (CLM4; Lawrence et al., 2012) and found that simulated C loss was more rapid than the observations across all sites. The large discrepancies between the laboratory microcosm studies used to parameterize the CLM4 litter decomposition and the LIDET field study likely resulted from poorly constrained temperature, moisture, and nitrogen controls (Bonan et al., 2013).

While the LIDET study provided valuable in situ benchmarks for improving our process representation in models, it does not necessarily account for feedbacks between fire and fuel decomposition dynamics. Penman and York (2010) used a 22-year data set to examine the relative influence of climate and fire history on rates of litterfall, decomposition, and fuel loading, in a coastal Eucalypt forest in Southeastern Australia and found that litterfall and decomposition were both influenced by temperature, recent rainfall, and fire history. However, such feedbacks are not currently well understood or represented in models. While long-term studies are extremely valuable for evaluating and improving models, they are relatively rare-we need many more long-term decomposition studies across climates and fire regimes to better evaluate and improve our mechanistic representation of fine fuel accumulation in biogeochemical models-this must include studies of both litter and fine woody fuel decomposition. In many respects, long-term decomposition studies could follow the "body farm" design (Bass et al., 2004), where examples of woody debris and litter from different species commonly found in a given fire regime are tracked over the long term with associated factors such as fire intensity, microclimate variabilities, aspect, etc (e.g., Cornelissen et al., 2017;Trettin et al., 2021). Ideally, these sites should be adjacent to sites where long-term data relevant to fires and ecosystems are also being collected, such as National Ecological Observatory Network, Critical Zone Observatory, Long Term Ecological Research Network, or the Smithsonian Forest Global Earth Observatory (ForestGEO) locations.

Future research should also consider fuel dynamics and feedbacks when assessing climate-wildfire relationships. Decomposition and fire have typically been studied separately, even though they can strongly interact (Cornelissen et al., 2017;Hyde et al., 2011). For example, repeated, low-intensity fires can reduce microbial CO 2 respiration rates and extracellular enzyme activity in coniferous forests, which may promote mineral soil C storage (Pellegrini et al., 2021). Additionally, decomposition is highly sensitive to nutrient availability and prescribed burning can deplete N and P litter stoichiometry, further slowing litter decay (Butler et al., 2019). However, such feedbacks are not well represented in biogeochemical and land surface models, which may cause us to overestimate decomposition in areas that experience increasing fire frequency or severity.

In this paper, we focused on fine fuel accumulation, but, as recently highlighted by the sixth Intergovernmental Panel on Climate Change Assessment, an additional limitation of using fire-enabled Earth system models to assess the carbon cycles is the limited information on dynamic changes to vegetation strata following fires (IPCC, 2021). Extreme wildfire behavior, effects, and changing fire regimes are also driven by understory and canopy fuel dynamics including shrub cover, canopy base height, and canopy bulk density (Parsons et al., 2016;Peterson et al., 2005;Van Wagner, 1977). Similar analyses of model uncertainty and improvement are also required for the representation of these dynamic fuel properties in fire regime models.

Uncertainties associated with existing model structure and parameters must be thoroughly documented. Given that many of the governing decomposition equations are based on individual case studies from a single location, and because key parameters are often hard coded, a great deal of model structural uncertainty is currently ignored and difficult to characterize. To understand future climate-fuel-fire feedbacks, it is essential to be transparent about what model choices are being made, the reasons for those choices, and the associated uncertainty. This is particularly critical as the domain of biogeochemical and land surface models is expanded to include the evaluation of future wildfire regimes, wildfire effects, and how we can mitigate the effects of climate change on wildfire through management.