Global Patterns and Dynamics of Burned Area and Burn Severity

: It is a widespread assumption that burned area and severity are increasing worldwide due to climate change. This issue has motivated former analysis based on satellite imagery, revealing a decreasing trend in global burned areas. However, few studies have addressed burn severity trends, rarely relating them to climate variables, and none of them at the global scale. Within this context, we characterized the spatiotemporal patterns of burned area and severity by biomes and continents and we analyzed their relationships with climate over 17 years. African ﬂooded and non-ﬂooded grasslands and savannas were the most ﬁre-prone biomes on Earth, whereas taiga and tundra exhibited the highest burn severity. Our temporal analysis updated the evidence of a decreasing trend in the global burned area ( − 1.50% year − 1 ; p < 0.01) and revealed increases in the fraction of burned area affected by high severity (0.95% year − 1 ; p < 0.05). Likewise, the regions with signiﬁcant increases in mean burn severity, and burned areas at high severity outnumbered those with signiﬁcant decreases. Among them, increases in severely burned areas in the temperate broadleaf and mixed forests of South America and tropical moist broadleaf forests of Australia were particularly intense. Although the spatial patterns of burned area and severity are clearly driven by climate, we did not ﬁnd climate warming to increase burned area and burn severity over time, suggesting other factors as the primary drivers of current shifts in ﬁre regimes at the planetary scale.


Introduction
Fire activity plays a key role in shaping ecological communities, biogeochemical cycles, climate, and human lives and assets [1,2].More than half of the land surface on Earth is prone to fire [3], with about 3% burning annually [4].As a result, landscape fires generate annual emissions estimated in 2.2 PgC, an equivalent of 23% of global CO 2 from fossil fuels [5], influence global albedo [6], and cause premature deaths from poor air quality, dozens of direct fatalities, and annual economic losses estimated at above US$2500 million [7].Despite these global figures, burned area (BA) and its directly related fire regime variable, the fire frequency [8], are largely heterogeneous across space and time because of differences in their main determining factors.Among the BA determining factors there are fuel load, which is linked to primary productivity and herbivory [1,9,10]; fuel connectivity, which is enabled by non-fragmented landscapes and non-sparse vegetation [11,12]; flammability of fuels, closely linked to weather and vegetation properties [2,13,14]; and ignition sources, which might be natural or anthropogenic, the last exceeding 90% of ignition events in most terrestrial biomes [8].In the same way, the ecological consequences of fire vary depending on the intrinsic ecosystem traits, post-fire environmental conditions, and fire characteristics, mainly burn severity (BS).BS is closely linked to fire intensity and here is defined as the immediate degree of overall environmental change caused by a fire in ecosystems, including biomass loss, vegetation mortality, and biochemical and physical impacts on soil [15].Moreover, BS determines the post-fire processes, playing a key role in the future ecosystem pathways [15][16][17][18], and its characterization is essential for the refinement of carbon emission models [5,15,19].
Nowadays, the exploitation of some types of satellite imagery provides globally consistent BA products based on the detection of active fires and the spectral changes in surface reflectance [4,20].Remote sensing analysis revealed grasslands and savannas as the most fire-prone biomes on Earth, mainly those located in Africa and Northern Australia.In these regions many areas burn annually or biennially, exhibiting a large quantitative difference from the rest of the biomes, where generally less than 2% of the area burns every year [4,5].Satellite imagery also allows the characterization of spatial patterns of BS [19,21,22].BS has been customarily assessed by spectral indices based on the near-infrared and shortwave infrared reflectance, such as the difference in the normalized burn ratio (dNBR) [16,22,23].The dNBR accurately matches field measurements of the overall environmental change caused by a fire in multiple ecosystems (mainly calculated through composite indices combining several ecosystem metrics) and has been adopted as a standard BS metric by the EFFIS and MTBS programs of the European Union and the United States of America, respectively [19].Despite the lack of comprehensive global spatiotemporal BS characterizations [7], current knowledge of maximum fire intensities [20], biomass burned fractions [5] and BS drivers (mainly the fuel load available to burn and fuel continuity) [24][25][26] anticipate a large spatial variability in BS across the globe.
The evolution of BA and BS arouses great interest in the media and in the scientific community, which frequently warned about the increase or worsening of fire activity during recent years [7], often attributing that assumed trend to climate change (see examples in [17,[27][28][29]).In this sense, climate warming is expected to cause increases in fire weather danger in many regions [14], which is a driver of BA in a large proportion of Earth's land surface [13] and influences BS [24,25].However, former empirical evidence of BA and BS increases are often based on constrained observations, in terms of timescales or spatial coverage [7].In fact, global quantitative BA analysis has shown significant decreases in the global BA between 2003 and 2015, mainly concentrated in the tropical savannas of Africa, South America, and Asian steppes, albeit BA increases have been detected in many regions such as the Siberian Arctic [30] and others dominated by closed canopy forests [4], including Australia where BA increases have been fostered by climate change [31].Likewise, little is known about BS trends at the global scale as most of the literature focuses on USA, Australia, and Mediterranean Europe, which might lead to a global "Western" biased perception [7].For instance, it is well documented a shift from low to high severity fire regimes in southwestern US forests, caused by the implementation of fire suppression policies after the European settlement [1], and extensive spatiotemporal studies have revealed a generalized increase in high-severity fires in some parts of Western US between 1984-2015 [21].In Australia, an increase in the proportion of BA at high severity has been detected between 2013 and 2017 [17], and in Europe, an increased prevalence of extreme wildfire events attributed to climate change and human alterations of landscapes has been reported [29].
Here, we characterize the spatial patterns and temporal trends of BA, BS, and BA by BS levels at the global scale and by biomes and continents for the period 2003-2019 (17 years).In addition, we studied the potential relationships between spatiotemporal patterns of BA, BS, and BA by severity levels with climate variables (mean annual temperature, annual temperature range, annual precipitation, and annual precipitation range) to identify the former role of climate change on fire activity.To achieve these objectives, we used NASA-MODIS-derived products at 500 m spatial resolution and ERA5 ECMWF re-analysis climatic data.

Data Sources
The BA and BS data were obtained from the MOSEV database [19].MOSEV has been developed based on MODIS information from the MCD64A1 collection 6 [32], which is the standard NASA global BA product and probably the most used by the scientific community [19]; and from the Terra MOD09A1 and Aqua MYD09A1 version 6 products, which provide the highest quality surface reflectance observations in eight-day periods.MOSEV offers, among other information, monthly global data at 500 m spatial resolution from 2000 to 2019 of date of burning and BS measured by the dNBR spectral index ranging from −2000 to 2000 [19].Burn severity in the MOSEV database refers to short-term BS, also called initial assessment of BS as the closest pre-and post-fire MOD09A1 and MYD09A1 scenes are used for calculations.The alternative to initial BS assessments, are the extended assessments that are often made one year after burning and thus are influenced by ecosystem recovery [15,16], BS in the MOSEV database is directly related to short-term BS calculated with Landsat imagery (R = 0.74 for dNBR; R = 0.42 for RdNBR) [19].Here, we used the dNBR instead of the relativized index RdNBR, also available in MOSEV, according to (I) its better performance when using MODIS; (II) to its similar or better fit to field composite metrics of BS such as the Composite Burn Index in different types of ecosystems; and (III) to be consistent with our definition of BS, which is understood as the overall degree of environmental change caused by fire [19].For the present study, the period 2003-2019 was selected for consistency because Aqua products were not available for entire years before 2003.Mean temperature and precipitation data were extracted for each month of the time series from the ERA5 ECMWF reanalysis.ERA5 was downloaded from the Copernicus Climate Data Store (CDS), using the CDS Python API.

Data Preparation
All data were extracted globally and by the regions conformed by the intersection of terrestrial biomes (i.e., global vegetation units) [33] with the continent boundaries [34].The extraction of the global and regional data from both MOSEV and ERA5 was designed in the R programming language [35].All the routines were implemented in the supercomputing facilities of the Spanish Research Council (CSIC).
For each year, we computed globally and by regions the BA, the percentage of the land burned, the mean BS, the percentage of land burned at different BS levels, and the percentage of the BA burned at different BS levels.The number of BS levels and threshold values depend on the user's objectives [16,19,36].The differentiation of three categories (low, moderate, and high BS) is probably the most common approach [19,37,38] and many studies have followed the pioneering proposal by Key and Benson [16] to differentiate levels.In the present study, we differentiated low, moderate, and high BS levels by the standard <270, 270-440, and >440 dNBR ranges, respectively, following the thresholds proposed in [16].We have used the dNBR value of >440 to differentiate the high severity rather than other common higher values (e.g., 660) as those have been found to be excessive by multiple studies [37,38], and because values above 660 are scarce in the MOSEV database [19].Moreover, we highlight that we have applied the same thresholds to all the regions for consistency; thus, same categories reflect similar overall spectral change with respect to the pre-burn situation regardless of the region.
Likewise, for each year we computed the mean annual temperature, the annual temperature range as the difference between the hottest and coldest months, the annual precipitation, and the annual precipitation range as the difference between the wettest and driest months.This raw database with the annual values of all study variables was checked for coherence before performing the statistics based on regional data.Thus, biomes corresponding to lakes, water bodies, rock, and ice were removed, as well as those regions with BA registered less than 10 years of the study period (<60% of the time).
To perform the temporal trend analyses (2003-2019; n = 17), we normalized the global and regional annual data of fire variables to the corresponding 2003 values (with values of 2003 set to 100%), to facilitate interpretation and intercomparison.Afterward, we calculated the Sen robust regression estimator with 95% confidence intervals using the zyp.sen and confint.zypfunctions of the zyp package [39], and the Mann-Kendall significance using the cor.test function.The Sen robust regression estimator is a non-parametric statistical analysis computed from the median of the slopes of all lines through pairs of points, and thus it is insensitive to outliers.The non-parametric Mann-Kendall significance indicates if a variable consistently decreases or increases over time.
The influence of climate variables (independent) on each fire variable (dependent) was studied with the data by regions (n = 64; n = 63 for BA at high severity analyses because of the absence of this BS level in one region) using two complementary statistical approaches: first, the relationship of each single climate variable with each fire variable was shown by fitting the local polynomial regression (loess), and the Spearman's rank correlation coefficients and significances were calculated using the cor.test function; secondly, conditional inference regression trees were implemented to account for complex combined effects of climate variables.All mean 2003-2019 BA data were square root transformed for these statistical analyses and to facilitate outputs visualization.Regression trees were built with the ctree function of the partykit package [40], including as predictors the four climate variables together.Regression trees include the most significant partitions up to p < 0.10 calculated with the approximated finite-sample distribution of Monte Carlo for each node.The variance explained by the regression tree predictions against our data was also shown.The analyses were run to analyze both, the synchronous relationships between fire and climate by using the 2003-2019 means; and the influence of temporal trends, using the Sen slopes as input data.
To analyze the relationships between pairs of fire variables, as well as between pairs of climate variables, correlation matrices with the loess fit and Spearman's rank correlation coefficients and significances were computed for the synchronous spatial data (mean 2003-2019 values) and trend data (Sen slopes) by regions (n = 64) using the pairs.panelsfunction of the psych package [41].

Spatial Patterns of Burned Area and Burn Severity
We observed that 4.78 ± 0.12 Mkm 2 (mean ± standard error), and around 3.26 ± 0.08% of Earth's land surface burned annually in the period 2003-2019 (Table 1), similar to NASA-MODIS MCD64A1 BA data reported for 2003-2015 [4].The spatial patterns of absolute BA were highly heterogeneous among regions (biomes and continents) (Table A1).The largest contributors to global BA were the tropical and subtropical grasslands, savannas, and shrublands of Africa (2.79 ± 0.07 Mkm 2 year −1 ), Australia (0.36 ± 0.03 Mkm 2 year −1 ) and South America (0.26 ± 0.02 Mkm 2 year −1 ), a consequence of both their proneness to fire and their vast extent.
Relativizing the BA to the total extent of each region (Figure 1; Table A1), we found that the most fire-prone regions were the flooded grasslands and savannas of Africa (26.97 ± 1.01% year −1 ), as well as the tropical and subtropical grasslands, savannas and shrublands of Africa (20.07 ± 0.49% year −1 ), and Australia (17.01 ± 1.32% year −1 ).There, the seasonally dry climate enables positive feedback interactions between primary pro-ductivity, fuel connectivity, and fire [1,12,42] which are essential for maintaining the open structure of these biomes [9,43].Among grasslands and savannas, the highest potential to burn in those seasonally flooded has been attributed to greater fuel accumulation rates as consequence of a higher productivity and the low herbivory of their less palatable vegetation [10].We found other biomes to be particularly fire-prone with more than 5% of their area burned annually (Figure 1; Table A1), including the tropical and subtropical dry broadleaf forests of Africa and montane grasslands and shrublands of Africa and Australia.Likewise, a surprisingly high BA fraction (4.38% to 3.42% year −1 ) was detected in tropical and subtropical moist broadleaf forests of Africa and Australia as fire is not an intrinsic ecological process in rainforests [1,31].The BA fraction was below 1.20% year −1 in all temperate forests and boreal biomes, except for the Australian temperate broadleaf and mixed forests (1.72% year −1 ).Although Mediterranean biomes are considered among the most fire-prone on Earth [3], their BA ranged from 0.45% year −1 in Europe to 1.56% year −1 in North America.Our results revealed that most biomes burned more in Africa and Australia than on the other continents.In Africa, extensive burning is facilitated by the low population density and decreases in grazing by bulk-feeding herbivores [12], whereas the proneness of Australian biomes to fire is fostered by the extreme climate seasonality, the high number of dry lightning events, aboriginal fire management practices, the physiognomy of eucalypts vegetation [31,42,44], and even by the colonization of exotic grasses [45].
Table 1.Global mean values and normalized trends (value in 2003 = 100%) in burned area (BA), burn severity (BS), and BA by BS levels at the global scale for the period 2003-2019.Both, BA, and BA by BS levels are presented in absolute area values and in percentage with respect to the total extent of the global land surface, BA by BS levels are also presented in percentage with respect to the global BA.Note that BA and BA by BS levels, whether calculated as absolute area or as percentage with respect to the total extent of the region results in the same trend values.The entire land area of the world was considered in calculations, including regions with scarce or null fire occurrence.Significant trends (p < 0.05) are denoted by boldface and the number of asterisks denotes the level of statistical significance (p < 0.05, p < 0.01, and p < 0.001).SE: standard error; CL: confidence limit.

Fire Variable
Global BS spatial patterns were closely linked to biomes, being more consistent among continents than BA (Figure 1, Table A1).Generally, the BS patterns detected in the present study are in line with those of fire radiative power emitted by fires detected in previous work [20], revealing at the global scale the assumed relationship between BS and fire intensity [15].In general, we found the highest mean severities in the taiga, followed by tundra, temperate coniferous forests, and Mediterranean biomes.BS was particularly high in North American taiga, which is characterized by high-intensity stand-replacing crown fire regimes [25], whereas lower mortality (42% lower) and combustion completeness (36% lower) characterize the Eurasian boreal forests, due to different traits of the dominant species [2,46].The high severity of tundra fires is a consequence of the burning completeness of vascular vegetation and moss cover, as well as of the upper soil organic strata accumulated for decades, which contributes to permafrost degradation [47].It is important to note that biomes at high latitudes have the largest topsoil carbon stocks on Earth (0-30 cm) [48], and seasonally become available to severely burn.In Mediterranean biomes, high severity results from a marked crown fire regime [1,3].The lowest severities were found in tropical coniferous forests, mangroves, and in deserts and xeric shrublands; in the first, combustion completeness is flammability-limited, and the last BS is constrained by fuel-limitations.Tropical and subtropical forests in Africa exhibited lower BS than in Asia and North America, which can be a consequence of burning at smaller fire patches [46] as fire size covariates with BS [17], and potentially because of the predominance of burning in lower tree cover ranges (25-75%) [46].Low severity values were also found in non-flooded grassy biomes, and intermediate severities generally corresponded to temperate forests.These global patterns confirmed the tradeoffs between BA fractions and BS (ρ = −0.25;p < 0.05) (Figure A1) previously detected for fire intensity [20], at the time that supports the fuel load available to burn as the primary driver of BS [25,26,49].
The analysis of BA by different severity levels, particularly at high severity, is useful to acknowledge the territorial magnitude of fire impacts [21].Globally, we revealed that 0.21 ± 0.00 MKm 2 (0.14 ± 0.00%) of land surface burned annually at high severity (Table 1).In absolute terms, the extent of BA affected by high severity was significantly linked to the BA extent (ρ = 0.87; p < 0.01) (Figure A1), the largest contributors to the global high severity BA are the tropical and subtropical grasslands, savannas, and shrublands of Africa, South America, and the Asian taiga (Table A1).However, relativized to the total extent of each region, we found that the regions most affected by high-severity burning were the montane grasslands and shrublands of Australia (2.23 ± 1.77% year −1 ), flooded grasslands and savannas of North America (1.85 ± 0.13% year −1 ), Africa (1.17 ± 0.12% year −1 ) and South America (0.96 ± 0.19% year −1 ), tropical and subtropical grasslands, savannas, and shrublands of Oceania (1.60 ± 1.59% year −1 ) and South America (0.55 ± 0.04% year −1 ), and Mediterranean North America (0.52 ± 0.11% year −1 ) (Figure 1; Table A1).In the rest of the regions, the high severity BA constituted less than 1% of their extent.This is a direct consequence of the BA in each region (ρ = 0.71; p < 0.01), being related to a lesser extent to the mean BS (Figure A1).Relativizing the BA at high severity to the total BA in each region, similar patterns to BS were detected (ρ = 0.82; p < 0.01) (Figure A1), with around 40% of BA exhibiting a high severity in all boreal forests and non-European tundra (Figure 1; Table A1).

Temporal Trends in Burned Area and Burn Severity
World BA decreased −1.50% year −1 ; p < 0.01 which in absolute terms is equivalent to a decrease of 1.22 MKm 2 between 2003 and 2019 (25.5% less BA with respect to 2003) (Table 1).These results elucidate an intensification of the decline in BA compared to the results already reported for the period 2003-2015 [4].Analyzing BA trends by regions (Figure 2; Table A2), we found significant BA decreases in the European taiga (−6.21% year −1 ; p = 0.04), temperate grasslands, savannas and shrublands of Asia (−2.77% year −1 ; p = 0.01), Mediterranean Europe (−2.31%year −1 ; p = 0.03), tropical and subtropical dry broadleaf forests of North America (−2.11% year −1 ; p < 0.01), montane grasslands and shrublands of Africa (−2.08% year −1 ; p < 0.01) and Asia (−1.93% year −1 ; p = 0.02), and in most of tropical Africa including grasslands, savannas and shrublands (−1.64% year −1 ; p < 0.01), moist broadleaf forests (−1.57% year −1 ; p = 0.02) and mangroves (−1.37% year −1 ; p = 0.02).Generalized decreases in BA have been formerly attributed to a decrease in the number of ignitions and to a lesser extent, to decreases in fire size, driven by human activity [4].In this sense, the increase in human population, cropland area, and livestock density cause decreases in fire activity in the fire-prone open biomes [4], whereas increased efficiency in fire prevention, detection, and extinction, and abandonment of fire use in agriculture contribute to the decreasing trend in other regions [7,50].BS spatial patterns were closely linked to biomes, being more consistent among continents than BA (Figure 1, Table A1).Generally, the BS patterns detected in the present study are in line with those of fire radiative power emitted by fires detected in previous work [20], revealing at the global scale the assumed relationship between BS and fire intensity [15].In general, we found the highest mean severities in the taiga, followed by tundra, temperate coniferous forests, and Mediterranean biomes.BS was particularly high  Globally, BS exhibited a non-significant increase (0.13% year −1 ; p = 0.34) (Table 1).However, we found a large continental disparity, and regions with significant increases outnumbered those with significant decreases in BS (Figure 3, Table A2).The largest increases were found in South America, including temperate grasslands, savannas, and shrublands (2.47% year −1 ; p = 0.03) and all forest biomes (1.49 to 3.56% year −1 ; p < 0.05) except mangroves.Likewise, we detected significant increases in BS in tropical and subtropical moist forests of Australia (2.23% year −1 ; p = 0.04), Mediterranean Africa (2.16% year −1 ; p < 0.01), and Asian taiga (2.13% year −1 ; p = 0.03).Land cover changes might explain BS trends in many regions, for instance, most of Mediterranean Africa, Central and Southern Chile, and Siberian taiga have experienced long-term woody encroachment and forest expansion [51], a fuel accumulation that in these regions lead to more severe fires [3,49].Likewise, BS increases in the tropical rainforest of South America may respond to several causes.First, to increases in fuel continuity in the Eastern Brazilian coast, which is frequently burned [50]; second, to the loss of primary forests [51] that are characterized by no fire or low severity surface fires [1]; and third, to decreases in lower atmosphere moisture boosted by these forest losses that resulted in increased droughts the last years of the study period [52], making large fuel loads available to burn.Significant decreases in BS were only found in montane grasslands and shrublands of Asia (2.94% year −1 ; p = 0.04) and in the European taiga (2.56% year −1 ; p = 0.04), which inversely to the Asian, has experienced decreasing trends in fuel continuity and increasing moisture which have detrimental effects on BA [50] and likely on BS [24,26].
The analysis of BA by BS levels showed that decreases in global BA were mainly at the expense of decreasing low (−1.62% year −1 ; p < 0.01) and moderate severity BA (−1.27% year −1 ; p < 0.01) (Table 1).Globally, we found trends in BA at high severity (−0.63% year −1 ; p = 0.09) to be positively related to trends in BA (ρ = 0.67; p < 0.01) as well as to trends in BS (ρ = 0.40; p < 0.01) (Figure A2).Regionally, we detected the highest increases in BA at high severity in temperate broadleaf forests of South America (40.20% year −1 ; p = 0.04) (Figure 4A; Table A2), which can be a consequence of the expansion of exotic Pinus and Eucalyptus plantations with associated fuel load increases up to 40 Mg ha −1 between 1999 and 2006, accompanied by a large drought-driven intensification of fire activity between 2010 and 2015 [53].We observed the BA at high severity to also aggravate in Australian tropical and subtropical moist forests (28.02% year −1 ; p < 0.01), which are vulnerable to fire, as it favors alternative open biome states [44].BA at high severity also increased in Australian grasslands, savannas, and shrublands (6.85% year −1 ; p = 0.01).
In Australia, increases in fuel continuity [50] and in fire weather [31] in the last decades have been detected, which, along with the unprecedented fire events that occurred in 2019, might contribute to the detected patterns.We found significant increases in BA at high severity in the deserts and xeric shrublands and tropical and subtropical dry broadleaf forests of Asia that locally can be attributed to encroachment [49] and to increases in fuel continuity [50].The largest decreases in high severity BA were observed in European taiga (−6.03% year −1 ; p = 0.01) and in tropical and subtropical grasslands, savannas, and shrublands of North America (−5.96% year −1 ; p = 0.04), a consequence of decreases in both BA and BS.Moreover, a significant global increase in the fraction of BA that is burned at high severity was detected (0.95% year −1 ; p = 0.04), and regionally, trends in the fraction of BA burned at high severity were closer to BS (ρ = 0.51; p < 0.01) than to BA (ρ = 0.11; p = 0.38) (Figure A2), confirming the aggravation of impacts within burned areas in a vast proportion of South America, and in large parts of Northern Australia and Asia (Figure 4B, Table A2).In the tropical forests of Australia and the Amazon, it has been revealed that contemporary logging regimes and silvicultural practices exacerbate burn severity [54,55], suggesting decision makers and forest managers have a determinant role in past, present, and future fire regimes.

Relationships with Climate Variables
The climate spatial patterns (Figure A3) determined the spatial patterns of fire activity (Figures 5 and 6).Local regressions and correlation analysis (Figure 5) showed monotonic increases in BA (ρ = 0.58; p < 0.01), and monotonic decreases in BS (ρ = −0.35;p < 0.01) and in the BA at high severity (ρ = −0.54;p < 0.01) towards the warmest regions (tropics), in agreement with the results shown by the regression trees (Figure 6A-C).Moreover, the regression trees revealed mean annual temperature as the primary climate driver of regional differences in BA, BS, and BA at high severity (Figure 6A-C).Local regressions (Figure 5) and correlation analysis also showed an opposite pattern for annual temperature range, but regression trees (Figure 6A-C) suggest lower importance of annual temperature range and the other climate variables (p > 0.05), probably because of their strong interdependencies (Figure A4).
We found climate trends to be weakly related to all BA and BS trends (Figures 6D and 7), the detected trends being opposite to the assumption of increases in BS with climate warming.Thereby, we detected a significant inverse relationship between climate warming and the proportion of BA at high severity (ρ = −0.26;p = 0.04) (Figure 7), and both the local regression and the regression tree (Figure 6D) showed that the most stable regions in terms of proportion of BA burned at high severity were those experiencing the highest warming (>0.03K year −1 ).In this context, the long-term outcome of increased temperatures for fire regimes has been disclosed complex as it is caused by shifts in fire weather [31], but also by the warming effects on fuels through changes in productivity, vegetation composition [2,26], and fuel accumulation in the soil uppermost layers that are controlled by productivitydecomposition equilibriums varying at fine scales [56].Moreover, the climate shifts would have different impacts on fire regimes depending on the former climate, productivity, and induced vegetation trajectories [24,27].This complexity is expected to cause differential changes in fire regimes across regions, which some authors claim to be scarcely plausible until the coming decades [13,14].Likewise, these premises highlight the difficulty of linking climate change to generalized increases in BA and BS worldwide, as different variations are expected depending on regional intrinsic characteristics.The map in the upper part (A) shows the trends of BA at high severity calculated as absolute burned area at high severity (when calculated as percentage with respect to the total extent of the region the result is the same).The map in the bottom (B) shows the trends of BA at high severity calculated as a percentage of the burned area that burned at high severity.The color ramps are square root scaled.

Implications and Final Considerations
Our study updates BA data and constitutes the first global assessment of BS spatiotemporal patterns.The information, provided by biomes and continents, is essential in revealing increasing trends in BS in several regions and increases in the fraction of BA that is burned at high severity at the planetary scale, providing scientific support for widespread assumptions.However, we did not find evidence in line with the hypothesis that climate change is increasing the mean BS worldwide.Nevertheless, several issues should be considered when interpreting our findings.First, the used BA and BS data at 500 m spatial resolution are appropriate for global fire analysis and consistent for studying trends [4,19,31], but have limitations in terms of quantifying BA absolute values.Accordingly, it has been demonstrated that the use of finer spatial resolution imagery at regional scales can result in larger BA, particularly in those regions with a predominance of small fires [57,58].For this reason, as sensors and computational capabilities improve, we encourage future work to analyze BA and BS trends using higher spatial resolution imagery or statistically refined BA data [58], and even different methods to quantify burn severity [59].Secondly, although there are biome-fire regime relationships, it is important to note that most biome regions have been largely modified by human activities such as agriculture, livestock, or urbanization, which also have an important role on fire regimes [2,4,8,20,60].In this sense, we recommend exploring the implications of these drivers in determining the results presented in this study.Third, our study period is limited to 17 years because of no prior availability of both MODIS Aqua and Terra sensors.This could favor interferences of decadal ocean and atmospheric oscillations [61] in our results, but we assume little influence of El Niño Southern Oscillation (ENSO) due to 2-7 years oscillation of this phenomenon and the balanced events along our study period.Fourth, despite not finding strong evidence of climate change aggravation of BA and BS, we confirmed that climate is an essential variable influencing fire regimes worldwide [1,13,27], and further studies at finer scale map units of analysis are advisable to further disentangle the complex climate-fire-human interactions in the current context of climate change.

Conclusions
This study used 17 years of MODIS-derived data to analyze the spatiotemporal patterns of BA and constitutes the first BS analysis at the global scale.In addition, we showed the spatiotemporal patterns of BA and BS by regions (biomes and continents), relating them to several climate variables.
Our results updated the already known spatial patterns of BA across the globe and corroborated significant decreases in global BA.Our triple-way analyses of BS trends detected that globally (I) the mean BS exhibits non-significant increases; (II) the fraction of land affected by high BS shows non-significant decreases, as it is linked to the decreases in BA; and (III) the fraction of burned area that is affected by high BS is significantly increasing.In addition, the number of regions showing significant increases in mean BS and burned area at high BS outnumbered those with significant decreases.
We found close relationships between the spatial patterns of fire variables (BA and BS metrics) and climate variables (temperature, precipitation, and their respective interannual ranges), but our analysis by regions has not found evidence that climate warming is increasing BA nor BS, suggesting other factors as the primary drivers of change.
Given the great technical and computational advances that are currently taking place, future work may continue our analysis by using higher-resolution images, smaller analytical units, longer periods, or different methods of quantifying BS.We also encourage future research to analyze region-by-region the implications that our results may have in the field of ecology, climate regulation, and in the effectiveness and design of environmental management strategies.

Figure A1.
Synchronous relationships between the spatial patterns of mean annual burned area (BA) expressed in mean absolute values (Mkm 2 yr −1 ) and in percentage with respect to the total extent of each region (% land yr −1 ), burn severity (BS) expressed in dNBR units ranging from −2000 to 2000, and BA at high severity expressed in absolute values (Mha yr −1 ), in percentage with respect to the total extent of the land extent in each region (% land yr −1 ), and in percentage with respect the total BA in each region (% BA yr −1 ) for the period 2003-2019.Red lines in the scatterplots show the local polynomial regression fitting (loess) and shaded areas ±95% confidence intervals.Values in the panels are Spearman's correlation coefficients (ρ) between pairs of variables, and the number of asterisks denotes the level of statistical significance (p < 0.05, p < 0.01, and p < 0.001).Note that all axis for BA variables were square root scaled.Table A2.Mean trends in annual burned area (BA), in burn severity (BS), and in high severity BA (dNBR > 440) by regions (biomes × continents) for the period 2003-2019.Note that trends are normalized (value in 2003 = 100%) and expressed in percentage change per year, and therefore both, trends in mean absolute values of BA (originally measured in Mha year −1 ), and in BA with respect to the total extent of global land extent (originally measured in % land year −1 ) are the same.Slopes were calculated using the Sen slope estimator and significances with the Mann-Kendall test.Significant decreases (p < 0.05) are denoted by green boldface, significant increases by red boldface, and the number of asterisks denotes the level of statistical significance (p < 0.05, p < 0.01, and p < 0.001

Figure 1 .
Figure 1.Mean values of burned area (BA), burn severity (BS), and BA by BS levels by regions (biomes × continents) for the period 2003-2019.BA values are represented by square root scaled barsand expressed in percentage with respect to the total extent of the regional land surface.The percentage of the BA burned at each BS level is proportional to their fraction within the BA bars.BS is represented by red points and expressed in dNBR units ranging from −2000 to 2000.S.Am: South America, Oce: Oceania, N.Am: North America, Eu: Europe, Aus: Australia, As: Asia, Afr: Africa.

Figure 1 .
Figure 1.Mean values of burned area (BA), burn severity (BS), and BA by BS levels by regions (biomes × continents) for the period 2003-2019.BA values are represented by square root scaled bars and expressed in percentage with respect to the total extent of the regional land surface.The percentage of the BA burned at each BS level is proportional to their fraction within the BA bars.BS is represented by red points and expressed in dNBR units ranging from −2000 to 2000.S.Am: South America, Oce: Oceania, N.Am: North America, Eu: Europe, Aus: Australia, As: Asia, Afr: Africa.

Figure 4 .
Figure 4. Normalized trends (value in 2003 = 100%) in burned area (BA) at high severity by regions (biomes × continents) for the period 2003-2019.The map in the upper part (A) shows the trends of BA at high severity calculated as absolute burned area at high severity (when calculated as percentage with respect to the total extent of the region the result is the same).The map in the bottom (B) shows the trends of BA at high severity calculated as a percentage of the burned area that burned at high severity.The color ramps are square root scaled.

Figure 5 .
Figure 5. Scatterplots showing the synchronous relationships between the climate and fire variables for the period 2003-2019.Red lines show the local polynomial regression fitting (loess), and shaded areas ± 95% confidence intervals.Each panel also shows the Spearman's correlation coefficient (ρ), and the number of asterisks denotes the level of statistical significance (p < 0.05, p < 0.01, and p < 0.001).Note that the BA axes are square root transformed.BA: burned area, BS: burn severity, T: mean annual temperature: T range: annual temperature range, P: annual precipitation, P range: and annual precipitation range.

Figure 6 .
Figure 6.Conditional inference regression trees showing the most important climatic variables in determining the spatial patterns of burned area (BA) (A), burn severity (BS) (B), BA at high severity (C), and normalized trend (2003 = 100%) in the percentage of the BA burned at high (D).No further trees were able to grow for the other studied fire variables and trends at the selected confidence level for recursive binary partitioning (p < 0.10).R 2 values on the bottom of panels represent the variance explained by the regression tree predictions on our dataset.Red dots in the boxplots indicate the mean values.T: mean annual temperature, T range: annual temperature range, P: annual precipitation.

Figure 7 .
Figure 7. Scatterplots showing the relationship between the studied climate trends and normalized trends (value in 2003 = 100%) in fire variables for the period 2003-2019.Red lines show the local polynomial regression fitting (loess), and shaded areas ±95% confidence intervals.Each panel also shows the Spearman's correlation coefficient (ρ), and the number of asterisks denotes the level of statistical significance (p < 0.05, p < 0.01, and p < 0.001).BA: burned area, BS: burn severity, T: mean annual temperature: T range: annual temperature range, P: annual precipitation, P range: and annual precipitation range.

Figure A3 .
Figure A3.Spatial patterns (A-D) and trends (E-H) in mean annual temperature, annual temperature range, annual precipitation, and annual precipitation range for the period 2003-2019.Trends were calculated using the Sen estimator and significances with the Mann-Kendall test.

Figure A4 .
Figure A4.Synchronous relationships between the spatial patterns of mean annual temperature (T), annual temperature range (T range), annual precipitation (P) and annual precipitation range (P range) for the period 2003-2019.Red lines in the scatterplots show the local polynomial regression fitting (loess) and shaded areas ±95% confidence intervals.Values in the panels are the Spearman's correlation coefficients (ρ) between pairs of variables.the number of asterisks denotes the level of statistical significance (p < 0.05, p < 0.01, and p < 0.001).

Figure A4 .
Figure A4.Synchronous relationships between the spatial patterns of mean annual temperature (T), annual temperature range (T range), annual precipitation (P) and annual precipitation range (P range) for the period 2003-2019.Red lines in the scatterplots show the local polynomial regression fitting (loess) and shaded areas ±95% confidence intervals.Values in the panels are the Spearman's correlation coefficients (ρ) between pairs of variables.the number of asterisks denotes the level of statistical significance (p < 0.05, p < 0.01, and p < 0.001).