**1. Introduction**

Although representing only a small fraction of the total number of fires, very-large fires (VLFs) are events often associated with dramatic economic, human health, and environmental risks that are unlike most other wildfires. The most salient and immediate economic impacts are suppression costs and property losses (Barrett 2018, [1]), which are often relatively large in VLFs compared to other smaller events (González-Cabán 1983 [2], Stephens et al., 2014 [3]). In addition to these direct costs, there is a suite of indirect economic impacts—such as damages from post-fire hazards, rehabilitation costs, lost tax and business revenue from community evacuations (Dale 2009 [4])—that are increasingly probable and costly in larger wildfires (Neary et al., 2003 [5], Peppin et al., 2011 [6], Beverly and Bothwell 2011 [7], Beverly et al., 2011 [8]). VLFs have the potential to burn large areas of vegetation and emit tremendous quantities of smoke within a short duration of time, which can adversely impact air quality for months at a time (Stephens et al., 2014 [3]), even at long distances from any active burning (Forster et al., 2001 [9], Val Martin et al., 2013 [10]). The large areas of active burning and sudden increase in air pollutants pose numerous risks to public health (Reid et al., 2016 [11]) and safety (Achtemeier 2009 [12], Stephens et al., 2014 [3]); and hospital admissions and treatment costs are expected to increase during VLFs (Moeltner et al., 2013 [13]). Although there are some ecological benefits of fire, VLFs have also been associated with significant, deleterious, and sometimes irreversible environmental changes. These include the production of environmental conditions conducive to the establishment of invasive species (Crawford et al., 2001 [14]), loss of ecosystem services (Rocca et al., 2014 [15]), and long-term modifications to forest structure (Haffey et al., 2018 [16]). Given the disproportionate costs VLFs pose to economic, social and environmental values, there is a broad need across disciplines, to better understand patterns and trends of their occurrence to mitigate future hazards. There is even greater urgency for this given that multiple lines of evidence suggest that VLF frequency has increased (Williams 2013 [17], Dennision et al., 2014 [18], Barbero et al., 2014 [19]) and will continue to increase into the future (Stavros et al., 2014 [20], Barbero et al., 2015 [21]). Still, there are several challenges to obtaining reliable predictions of future wildfire activity, many of which are related to various kinds of scientific unknowns or uncertainties. These uncertainties can come from many sources (Chen et al., 2018 [22]) and can be associated with a particular quantity of interest—what will the future frequency of very-large fire events be?—or associated with model structures—how are environmental conditions related to very-large fire frequency? The latter is referred to as model or structural uncertainties, which can have significant effects on the conclusions one draws from an analysis (Morgan et al., 1992 [23], Syphard et al., 2018 [24]). In the context of forecasting future wildfire activities, these structural uncertainties can arise from the selection of vegetation models (Sitch et al., 2008 [25], Syphard et al., 2018 [24]), assumed anthropogenic effects (Westerling et al., 2011 [26]), as well as greenhouse gas emmissions and their effects on the environment.

Structural uncertainties associated with the characteristics of the future environment are commonly accounted for in long-term climate impact studies through the use of multiple General Circulation Models (GCMs). Each GCM predicts climatological responses based on unique assumptions regarding chemical and physical interactions between a suite of factors including land, water, atmosphere, and the cryosphere. These models can be used to forecast future climate by forcing the models to observe historical atmospheric conditions and running the models forward using representative concentration pathways (RCPs) as plausible carbon emission scenarios. There are four future RCP scenarios, which are labeled RCP 8.5, RCP 6, RCP 4.5, and RCP 2.6; each assumes various levels of fossil fuel use and economic activity. The suffix labels correspond to the approximate 2100 radiative forcing levels. For instance, the high-emission (RCP 8.5) scenario corresponds to an approximate radiative forcing increase of 8.5 W/m<sup>2</sup> by 2100 compared to pre-industrial conditions. Since the choice of GCM and RCP can be thought of as competing plausible models of the future environment, and the precise climatological response and quantities of greenhouse gases that will be emitted and sequestered prior to 2100 are unknown, it makes sense to interpret them as structural uncertainties (Taylor et al., 2012 [27]).

In addition to the structural uncertainty arising from relating greenhouse gas emission scenarios to climatological impacts, structural uncertainty further arises when relating climatological variables to an impact of interest. In the context of fire, these relationships could include weather patterns such as temperature, precipitation, atmospheric moisture, winds, and clouds. Identifying and describing the relationship between these variables and VLFs is an immensely complex and subjective process, as there can be many competing hypotheses. For instance, temperature controls landscape flammability, is associated with thunderstorm activity and by extension ignition frequency (Flannigan et al., 2009 [28]), mediates tree mortality through drought (Allen et al., 2010 [29]) and insect pests (Bentz et al., 2010 [30]), and influences the length of the snow-free season (Westerling 2016 [31]). Additionally, the timing and amount of precipitation can also influence wildfire behavior in parallel with temperature by controlling the availability of fine fuels (Meyn et al., 2007 [32]), fuel moisture

(Flannigan et al., 2016 [33]), and distribution of flammable species (Bradley et al., 2016 [34]). Multiple weather variables may adequately measure a common phenomenon associated with wildfire risk, such as drought (Zargar et al., 2011 [35]), resulting in highly correlated covariates that when utilized in wildfire risk prediction, produce models with near-identical goodness-of-fit. Hence, although statistical models can be a useful tool for representing and identifying the relative importance of relationships between the environment and fire, they are still only approximations to reality. Confounding of unmeasured variables like vegetation management (Holsinger et al., 2016 [36]) may influence the predicted importance of some weather variables, which can make model selection challenging. In some cases, the same suite of covariates can be used to make predictions with multiple mathematical representations (Mallick and Gelfand 1994 [37]), presenting an additional level of uncertainty that is easily overlooked. Given the frequency with which we face structural uncertainties when modeling highly complex phenomena like wildfire, it is extremely risky to select any single model as an approximation of this phenomena, and a more robust approach should explore results from multiple models (Morgan et al., 1992 [23], Littel et al., 2011 [38]).

Bayesian model averaging is flexible and a commonly used method of accounting for structural uncertainties like these, lessening many of the risks of traditional model selection techniques and improving performance across a variety of metrics (Raftery and Zheng 2003 [39]). Within this framework, model weights, which are assumed to be uncertain, are used to combine covariates or predictions from multiple sources into a single probability distribution. The uncertainties in the model weights are represented using the posterior, which is a probability distribution representing the belief in the model parameters conditional on the observed data. While posteriors provide a natural framework for interpreting uncertainties, in practice, closed form expressions of the posterior are non-trivial and direct calculation is often impossible. Hence, simulation methods like Markov Chain Monte Carlo are typically used to generate samples from the distribution, which are in turn used to approximate the quantities that are of interest to the analyst (Fragoso et al., 2018 [40]). Computational barriers to Markov Chain Monte Carlo techniques have diminished greatly since they were first introduced, and a range of recently developed software options such as JAGS (Plummer 2003 [41]), Stan (Carpenter et al., 2017 [42]), and Integrated Nested Laplace Approximations (Rue et al., 2009 [43]) have facilitated the application of these methods in novel and previously infeasible contexts (Monnahan et al., 2017 [44]).

Hence, in this paper, we account for both kinds of structural uncertainty—uncertainty from the climate models and uncertainty from the choice of VLF models—using Bayesian model averaging to generate predictions of event frequency in the last half of the 21st century in the Continental United States. In Section 2, we present the methods used to produce robust predictions of future wildfire activity using GCMs and multiple fire occurrence models. In Section 3, the results of this analysis are available and demonstrate that increases in VLF activity should be expected in many regions in the Continental United States at the end of the century. In Section 4, we close the paper with a discussion of the implications of this analysis to decision-makers and researchers.

### **2. Methods**

#### *2.1. Fire Occurrence Data*

Data from the Monitoring Trends in Burn Severity (MTBS) project [45] , which describes individual fire size and severity based on changes in satellite imagery, are used to measure monthly fire occurrence. The original dataset included all detected fire events within the continental United States for the years 1984–2015 and is further filtered to remove all events 404 hectares or smaller, or that were non-wildfire. The filtered data are grouped into 18 regions with broadly similar climate and vegetation characteristics using a geospatial dataset of ecosystem divisions (Figure 1, Bailey 2016 [46]). For each region, two binary time series were constructed: one representing large fire (LF) occurrence, and another representing very-large fire (VLF) occurrence. The LF occurrence time series, *XLF*,*t*, reports "1" if at least 1 event of at least 404 hectares is recorded during that month; otherwise, it reports

"0". The VLF occurrence time series, *XVLF*,*t*, records "1" if at least 1 VLF—one that exceeds the 95th percentile of the region's filtered MTBS burn area records (Table 1)—is recorded during that month and region; otherwise, it reports "0". The Marine Division had one fire event and the Subtropical Regime Mountain had no VLF events during 1984–2005, and both were dropped from further consideration, yielding a total of 16 independent ecoregional analyses. All of the time series are split into a tuning dataset (1984–2005) and a training dataset (2006–2015).

**Figure 1.** Map of multi-model mean temperature changes between 1955–2005 and 2050–2099 over the relevant Bailey's divisions under the representative concentration pathway (RCP) 4.5 (top) and RCP 8.5 emission scenarios (bottom). Regions with insufficient fire occurrence data for analysis are colored white.


**Table 1.** Very-large fire cutoffs, the number of fires, large fire months, and very-large fire months for two time periods: 1984–2005 and 2006–2015.
