*2.4. Multi-Model Very-Large Fire Predictions*

The structural uncertainties arising from training PETs to observed meteorological data are compounded by structural uncertainties arising from the application of these models to long-term climate forecasts. To improve the quality of the probability estimates and VLF occurrence forecasts, multi-model averages of fire event probabilities are used to integrate both sources of structural uncertainty: the choice of the PET and selection of the climate model. The final probability estimates are assumed to be an average of predictions from all combination pairs of the 100 PETs within each forest and 13 modeled weather datasets; a total of 1300 individual predictions for each region, month, and probability of interest. The modeled weather data used to make PET predictions come from the second version of regional Multivariate Adaptive Constructed Analogs (MACA) [47] dataset that was trained on gridMET, and downscaled with 13 GCMs: bcc-csm1-1-m, BNU-ESM, CanESM2, CCSM4, CNRM-CM5, CSIRO-Mk3-6-0, GFDL-ESM2M, HadGEM2-ES365, inmcm4, IPSL-CM5A-MR, MIROC5, MRI-CGCM3, NorESM1-M.

The multi-model averages are calculated using both unweighted and weighted approaches. The unweighted approach assigns equal weight to each individual prediction and assumes that each PET and climate model is equally credible. The weighted approach assigns unequal weights to predictions from each pairing of PETs and climate models, to try to bias-correct the multi-model averages and optimize predictive performance. The final weight applied to each individual prediction is the product of two independent components: a climate model weight and a PET weight. For a specified region and month, we can write the estimated probabilities as,

$$\begin{aligned} \vec{p}\_{LF} &= \sum\_{i=1}^{100} \sum\_{j=1}^{13} u\_{LF,i} \upsilon\_j p\_{LF,i,j}, \\\\ \vec{p}\_{VLF} &= \sum\_{i=1}^{100} \sum\_{j=1}^{13} u\_{VLF,i} \upsilon\_j p\_{VLF,i,j}. \end{aligned}$$

Here, *u*∗,*<sup>i</sup>* represents the weight applied to the predictions from PET *i*, *vj* represents the weight applied to predictions utilizing climate model *j*, and *p*∗,*i*,*<sup>j</sup>* is the prediction obtained from PET *i* utilizing climate model *j*. We estimate the weight components using a fully Bayesian approach that incorporates fire occurrence and modeled climate forcings from 1984–2005, as well as probabilistic representations of possible parameter estimates. Via Bayes rule, we know that the posterior of the model weight components, *θ* = −→*u LF*, −→*u VLF*, −→*v* , is proportional to the product of the likelihood and prior probability distributions. The likelihood component represents the probability of observing the fire occurrence time series, *<sup>X</sup>* <sup>=</sup> −→*<sup>X</sup> LF*, −→*<sup>X</sup> VLF*, assuming that they were generated from a Bernoulli process parameterized with our weighted multi-model averages:

$$\begin{array}{c} \mathcal{X}\_{LF,t} \sim \text{Bernoulli}(\vec{p}\_{LF,t}), \\\\ \mathcal{X}\_{VLF,t} \sim \text{Bernoulli}(\vec{p}\_{VLF,t}\mathcal{X}\_{LF,t}). \end{array}$$

The prior component, *p*(*θ*), which is a probability density function representing our a priori belief regarding the parameter values, is defined using independent Dirichlet priors with uninformative concentration parameters:

$$
\overrightarrow{u} \underset{LF\prime}{\rightrightarrows} \underset{VLF\prime}{\rightrightarrows} \sim \text{Dirichlet}(\overrightarrow{1}^\*).
$$

The posterior was approximated using Just Another Gibbs Sampler (JAGS) software (Plummer 2003 [41]) and the runjags package in R (R Development Core Team 2008 [53]). An initial set of 30,000 samples were generated from three parallel Markov Chain Monte Carlo chains using a burn-in interval of 10,000 steps, adaptive phase of 10,000 steps, and thinning interval of 100. Calculations were performed on a MacBook Pro (Quanta Computer, Inc, Shanghai, China) with a 2.7 GHz Intel Core i7 processor (Hillsboro, Oregon, USA). Convergence was monitored visually, and also via the calculation of the potential scale reduction factor, using the range of the central 90th percentile of the marginal posteriors as a test statistic (Brooks and Gelman 1998 [54]). We assume that the second half of the chain has approximately converged if the maximum potential scale reduction factor fell below 1.01. If the chain had not converged, then it was continued in batches of 1000 iterations until the maximum potential scale reduction factor was less than 1.01 (Gelman and Shirley 2011 [55]). To provide guidance to future analysts looking to perform similar analyses, an informal computational comparison of JAGS (Version 4.3.0) and Stan software (Version 2.17.3) was completed, which is described in the supplementary materials. The final VLF probabilities can then be calculated by averaging both probabilities of interest—either with point estimate averages of the model weights or using the full posterior in the weighted approach—and applying Bayes rule (Figure 2). The resulting distribution of VLF probability time series is then used to estimate the changes in VLF frequency in the future by finding the difference between the expected number of VLFs in the historical climate (1956–2005), and the expected number of VLFs in the future (2050–2099) under moderate and severe warming scenarios, RCP 4.5 and RCP 8.5, respectively.

**Figure 2.** Workflow of very-large fire probability calculations. In the prediction phase, probability estimates are calculated for every combination pair of climate model and probability estimation tree. The model-averaging phase combines these predictions into probability estimates using either a weighted or unweighted averages. The final phase uses Bayes rule to calculate the very-large fire probability as the product of both components from the model-averaging phase.
