2 Materials and methods

2.1 Geospatial classification

We used the Terrestrial Ecosystems of the World regional classification system hosted by the World Wildlife Fund (Olson et al., 2001; see supporting information). This dataset groups regions into a hierarchical system using eight realms at the coarsest level (e.g., the Indo-Malay Archipelago and associated lands in Southeast Asia), 14 biomes (e.g., temperate broadleaf and mixed forests), and, at the finest scale, 867 distinct ecoregions (e.g., Madagascar’s subhumid forest). These ecoregions represent distinct biota and potential habitats, but do not consider the real effects of human land-use practices such as agricultural clearing. We removed a total of 95 of the 867 original ecoregions from our analyses, including 76 ecoregions with no associated fire records (e.g., those existing in the Antarctic realm). In addition, we removed the mangroves biome type from all analyses, as this biome is comprised of only 19 small, discontinuous ecoregions with fire histories primarily driven by neighbouring biome types.

2.2 Identifying and assessing the ecoregion flammability threshold

We used estimated dead fine fuel moisture content (DFFMC: %) as drawn from the Canadian Fire Weather Index (FWI) system’s fine fuel moisture conde (Van Wagner, 1987) as a foundation for identifying ecoregion flammability thresholds (EFTs). We chose this estimation of DFFMC due to its established global applicability and ease of both calculation and interpretation (Wotton, 2008; Field, 2020). To calculate DFFMC using the FWI, we used the European Centre for Medium-Range Weather Forecast’s ERA5 atmospheric reanalysis data (Hersbach et al., 2020) due to its accessibility and worldwide coverage. It should be recognised that the current ERA5 reanalysis product holds regional biases in capturing observed moisture primarily within the tropics (Lavers et al., 2022). These data were used to calculate daily DFFMC representing noon local standard time between 1979 and 2019 at a gridded 0.25° spatial resolution. High-latitude overwintering periods were removed from the records to reduce the false identification of dry winter periods as highly flammable (McElhinny et al., 2020). See Figure 1 for a detailed workflow diagram.
To then identify each EFT controlling an ecoregion’s average ability to ignite and sustain a wildfire, we first calculated the cumulative proportion of burnt area from the Moderate Resolution Imaging Spectroradiometer (MODIS) MCD64CMQ product (2001 - 2021: Giglio et al., 2018) over the full range of potential EFTs by ecoregion. We constrained the upper bounds of potential EFTs to 70% DFFMC to allow for a safe degree of uncertainty above the known limits of the fibre saturation point (Fernandes et al., 2008). We then used a nonlinear least squares regression formula to fit a logistic curve to the relationship between the cumulative proportion of burnt area (BA) and the upper limit of the ecoregion’s driest DFFMC quartile using Eq. (1). This equation assumes the 25th percentile of an ecoregion’s monthly DFFMC is reflective of that ecoregion’s driest period – even where concrete seasonal climate patterns are more nebulous:
\(P(BA<DFFMC)=\frac{\phi_{1}}{1+e^{\frac{-\left(\phi_{2}\ -\ \text{DFFMC}\right)}{\phi_{3}}}}\), (1)
where ϕ1, ϕ2, and ϕ3 are the asymptote, curve inflection point, and scale parameter of the curve, respectively. We extracted ϕ2 from the model formulae as the inflection point where the greatest marginal increase in the cumulative proportion of burnt area occurred for a given reduction in DFFMC. The DFFMC value at ϕ2 is subsequently the EFT constraining fire potential for a given ecoregion. We also extracted the reciprocal of the model parameter ϕ3, which provides a proxy measure of the inflection point slope estimate (IPSE). The IPSE complements the EFT and provides insight into the strength of the relationship between DFFMC and fire activity in moisture-limited environments. A total of 74 ecoregions did not have associated fire data within the DFFMC constraints to support identifying an EFT or the associated IPSE. Instead of removing these ecoregions from our analyses, we chose to impute the missing values of those 74 ecoregions using the median EFT from the same hierarchical realm and biome classifications where available. For example, these 74 ecoregions include three boreal forests in either in the Nearctic or Palearctic realms. The two ecoregions in the Nearctic realm were imputed using the median EFT for Nearctic boreal forests (i.e., 14%), while the third from the Palearctic realm was imputed using the median EFT for Palearctic boreal forests (i.e., 16.3%). See supporting information all modelled ecoregional P(BA < DFFMC) curves as well as biome-level means.
We assessed the identified EFTs by first visualizing the EFT distributions by biome type, highlighting the median and quantile-based intervals. We statistically confirmed the apparent effects of median biome differences using a Kruskal-Wallis rank sum test (Kruskal and Wallis, 1952), and then identified the mean rank differences between specific biome types using a two-way Dunn’s test for pairwise multiple comparisons (Dunn, 1961). We extracted non-significant differences between biome types, highlighting the similarities in fuel-fire relationships dependent on the underlying biota and potential habitats.

2.3 Statistical analyses of the fuel moisture-fire relationships

To identify the eco-climatological factors driving the ecoregions’ EFTs, we used a combination of nonmetric multidimensional scaling (NMDS) and generalised additive modelling. This combination of variable ordination and inferential modelling provides a foundation to analyse biogeographic implications of EFTs across different ecoregions and vegetation types. First, we applied NMDS ordination to a suite of climatological and ecological data with known associations to wildfire ignition or spread and then plotted the final ordination against our identified EFTs. We used bioclimatic indicator means (1979 – 2018) for annual precipitation (BIO12: m s-1), annual temperature (BIO01: K), and precipitation seasonality (BIO15: %) from the ERA5 reanalysis dataset available via the Copernicus Climate Data Store (Hersbach et al., 2020; Wouters et al., 2021). We also used a combination of three additional MODIS products: Mean annual net primary productivity (NPP: t C ha-1 year-1) calculated between 2000 and 2015 (MOD17A3: Running et al., 2015), and two vegetation continuous fields representing the median percent of an ecoregion represented by both herbaceous and tree vegetation in 2020 (MOD44B: DiMiceli et al., 2015). This minimal set of variables were selected as they are commonly applied and have known relationships with fire activity. We tested both two- and three-dimensional nonmetric multidimensional scaling NMDS ordination based on exploratory stress scree and Shepard plotting, ultimately choosing the two-dimensional ordination with a stress index of 0.125 for simplicity. See supporting information for the exploratory NMDS diagnostics behind this rationale. The new, reduced NMDS dimensions were plotted using scatterplot variable ordination with the underlying distribution of EFTs to both explore the effects the different eco-climatological variables theoretically hold, as well as reduce the number of variables retained in the model.
Building off of the previously-identified association between intermediate productivity and fire activity (Pausas and Bradstock, 2007; Pausas and Ribeiro, 2013; Bowman et al., 2014; Ellis et al., 2022), we plotted global NPP as a function of mean fire activity underlain by both our identified EFTs and the associated IPSEs. This analysis uses 0.25° ERA5 grid cells for NPP and indexed fire activity as shown in Ellis et al. (2022), adding weight to the breadth and extent of both globally. Following those prior analyses, we calculated fire activity indices (FAI: Pausas and Ribeiro, 2013) using the mean annual (2002 – 2020) number of fire detections as recorded in the MODIS active fire database (MCD14DL: Giglio et al., 2003) while excluding permanent, anthropogenic heat sources. By showing both the identified EFT and its associated slope estimate, we highlight the global distribution of EFTs and its relationship with productivity and indexed fire activity, while the IPSE provides insight into the strength of the EFT as a constraint in a particular ecoregion.
Our use of generalised additive modelling employed a Bayesian framework to measure the effects relevant eco-climatological variables have in shaping EFTs. Our modelling focus is ultimately not guided by predictive power, but rather to verify the identified EFTs as a functioning, biological mechanism on ecoregion-level fire behaviour. Because our EFTs can naturally be interpreted as a percentage with a strong positive skew, we estimated the response likelihood function on a beta distribution. Note that while the Canadian FWI’s calculation of DFFMC can reach limits of up to 250%, our EFTs all fall under 100% due to the 70% constraint we applied during our identification of the EFTs. Additionally, it’s unlikely any potential EFT could extend above that maximum based on vegetation’s theoretical maximum fibre saturation point of 30% (Fernandes et al., 2008). Informed by our exploratory analysis of the EFTs and the NMDS ordination, we retained precipitation, temperature, precipitation seasonality, and herbaceous vegetation cover as fixed continuous effects in the model. Of the highly colinear variables, we chose annual precipitation over both NPP and percent tree cover due to the former’s reliance on precipitation, and the latter being inherently linked to herbaceous cover. As realm and biome type don’t reflect true vegetation, but rather a complex, overlapping network of theoretical habitats, we retained both realm and biome only as interactive random effects. See supporting information for additional model development details, including data preparation steps and specific model parameters.
Finally, to evaluate and explain our model within our inferential framework, we focused on a combination of variable importance using permuted RMSE dropout loss and Bayesian effects probability measures. First, we calculated the influence of individual variables in the model by using sampled (n = 1,000) change in the RMSE loss function over 100 permutations. We then evaluated the continuous fixed effects of our model using the sequential effect existence and significance testing framework. This framework succinctly describes the effects of model parameters, providing three interlinked probability measures of overall effect direction (i.e., existence), practical effect significance, and size (i.e., strength) while also being easy to interpret due to the rough equivalence to statistical significance (e.g., p < 0.05 or a probability greater than 95%). These statistical and visual tools work well in conjunction to describe the impacts of individual eco-climatological variables on effecting the identified EFTs while confirming the complexity of fuel moisture as a biological constraint on localised (i.e., ecoregion-level) fire behaviour.