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.