\chapter{The fire regimes of the Brazilian Cerrado and their changes through time (Research Paper 2)}\label{ch:paper02}


\begin{refsection}

\raggedbottom

\textit{By: Carlota Segura-Garcia, Ane A. C. Alencar, Vera L. S. Arruda, David Bauman, Wallace Silva, Dhemerson E. Conciani, Imma Oliveras Menor}

\textbf{Linking statement:} This chapter expanded on \textbf{Chapter \ref{ch:paper01} (Research Paper 1)} by quantifying the fire regimes in a more complete manner, that is, including measures of frequency, fire size, and seasonality apart from burned area. This chapter was motivated by the need to understand how the fire regimes of the Brazilian Cerrado have changed over the past few decades, both in their characteristics and in their spatial distribution across this region. This chapter used a clustering algorithm to identify the regions of the Cerrado presenting distinct fire regimes in four consecutive periods (1985-1993, 1994-2002, 2003-2011, and 20129-2020). Then, focusing on the first and last periods, it explored the anthropic and climatic factors that characterise the different fire regime regions, and studied both how these regions and their fire regimes have changed through time.

This chapter was submitted on 29$^\text{th}$ February, 2024 to a special issue on the topic \textit{Novel fire regimes under climate changes and human influences: impacts, ecosystem responses and feedbacks} in \textit{Philosophical Transactions of the Royal Society of London, Series B: Biological Sciences}, and it is currently on a second round of reviews.

\textbf{Contribution Statement:} I, Carlota Segura-Garcia, conducted the conceptualization, methodology development, data processing and data analysis, investigation, visualization of the results and writing of the manuscript. In addition, I also responded to all review comments during the first round of the peer review process. Dr David Bauman helped in the methodology, analysis and reviewing of the manuscript. Dr Ane Alencar helped in the conceptualization, investigation and reviewing of the manuscript. Vera L. S. Arruda, Dhemerson E. Conciani and Wallace Silva helped in the conceptualization and investigation of this study. Finally, Dr Imma Oliveras Menor supervised and was involved in all steps of this work.

\newpage

\section{Abstract}
The Brazilian Cerrado is a heterogeneous region of open ecosystems adapted to fire intermingled with patches of woody growth-forms, with high levels of biodiversity and endemism. In recent decades, land conversion and human activities have proliferated across the Cerrado, losing about half of its original area. These changes, coupled with climate change, are altering its fire regimes with uncertain, but possibly adverse, consequences for Cerrado ecosystems. Here, we used burned area data to characterise the fire regimes of each cell on a 30 km grid over the Cerrado, and used a spatially constrained hierarchical clustering approach to delineate the regions with different fire regimes in four consecutive 9-year periods between 1985 and 2020. Comparing the periods 1985-1993 and 2012-2020, we found substantial changes in the number and shape of the fire regime regions, and in their fire characteristics. The main factor differentiating these regions was their level of fire activity: some showed large, numerous and frequent fires, while others showed small, few and infrequent fires. We also identified a region in the north with a later peak of the fire season, characterised by small but numerous fires. Finally, we found that while the fire activity of the southern areas of the Cerrado substantially decreased, fire activity levels in the centre and north increased or remained high over time.



\textbf{\textit{Keywords -}} fire regimes, open ecosystems, disturbance ecology, land use effects, climate change

\section{Introduction}

Across the world, areas with different types of vegetation, climatic ranges and socio-economic characteristics present different patterns of fire occurrence, referred to as fire regimes \citep{chuvieco_global_2008, archibald_defining_2013, cattau_modern_2022}. The relations between characteristics of fire like size, frequency, intensity or seasonality, and the factors that affect them are complex, often forming feedbacks by which one factor influences the others and vice versa \citep{archibald_biological_2018}. Further, the effects of different drivers of fire may become apparent only at certain scales - landscape, regional or global \citep{krawchuk_constraints_2011, cattau_modern_2022}. Fire regimes are a useful framework to study how the characteristics of fire vary geographically in relation to these factors \citep{cattau_modern_2022}, understand how fire may contribute to shaping different landscapes \citep{bond_global_2005, bond_fire_2005}, weigh its contribution to carbon emissions and other biogeochemical cycles \citep{van_der_werf_global_2010}, and foresee how fire may respond to changes in the climate. At regional to local scales, characterising fire regimes can also contribute towards the conceptualisation of fire management strategies and conservation plans that are adequate to specific contexts.


Certain open ecosystems, like savannas and grasslands, are adapted to fire as their vegetation has co-evolved with this phenomenon over millions of years \citep{pausas_burning_2009, bond_open_2019}. Therefore, given a certain vegetation composition, it is argued that there exists an underlying fire regime the vegetation is specifically adapted to \citep{pausas_burning_2009}. However, biomes of open ecosystems are usually very heterogeneous in their vegetation composition, topographic features, and climatic characteristics and, as a consequence, can present different fire regimes across their area \citep{archibald_defining_2013}. Furthermore, the expansion of modern human land use activities and anthropic climate change are rapidly altering the fire regimes of these regions \citep{archibald_managing_2016}. Indeed, in tropical savannas, humans shape fire regimes through landscape fragmentation \citep{rosan_fragmentation-driven_2022}, shifting the fire season \citep{ramos-neto_lightning_2000, archibald_managing_2016} by igniting fires for agricultural and land conversion purposes \citep{klink_conservation_2005, pivello_use_2011}, or actively preventing or suppressing fires - which may lead to more infrequent but larger and more intense fires as a result of fuel accumulation \citep{stevens_savanna_2017, stevens_trillion_2023}. Regarding climate, even though there is established evidence that climate change is altering savannas' fire regimes \citep{jolly_climate-induced_2015}, and that it will continue to do so \citep{silva_future_2016, silva_impacts_2019}, the specific long-term effects are still uncertain. While climate change is associated with surges in the number, size and intensity of fires in relation to increasing temperatures and the frequency of climate anomalies \citep{van_der_werf_climate_2008}, or the lengthening of the dry season \citep{jolly_climate-induced_2015}; in the long run, the effects on fire regimes will also depend on its effects on vegetation productivity during the growing season, which ultimately influences the amount of available fuel \citep{pausas_global_2013, krawchuk_constraints_2011}. In turn, the alteration of savannas' fire regimes as a result of the expansion of human activities and climate change may further couple with these two processes to affect the vegetation of these open ecosystems \citep{pivello2021understanding}, as well as having broader effects on the planet's carbon cycle \citep{beer_terrestrial_2010}. 

The Brazilian Cerrado is a clear example of large and heterogeneous mosaic of open ecosystems where both its landscape and fire regimes are undergoing rapid and profound changes in their configuration \citep{silva_putting_2021, mataveli_satellite_2018}. Indeed, despite its high levels of endemism and rich biodiversity \citep{myers_biodiversity_2000, klink_conservation_2005}, around 50\% of the Cerrado's 2 million km$^2$ have been converted to agricultural land uses in recent decades \citep{alencar_mapping_2020} roughly expanding from the south towards the Matopiba region in the north, the Cerrado's current agricultural frontier \citep{spera_land-use_2016}. Simultaneously, this region has also become warmer and drier \citep{hofmann_brazilian_2021} under climate change. Apart from its heterogeneous vegetation, the Cerrado presents a diverse topography, climate and soil composition, as well as varying levels of landscape fragmentation and different anthropic land uses, giving rise to different ecoregions \citep{sano_cerrado_2019}. All of these elements influence its fire regimes \citep{libonati_algorithm_2015, bowman_vegetation_2020, pausas_global_2013}. Indeed, studies across the Cerrado and its ecoregions have observed a great heterogeneity in the patterns of fire, with overall low burned area levels in the south of the Cerrado and higher levels in the centre of this region  \citep{silva_putting_2021, mataveli_satellite_2018}, as well as heterogeneous changes in burned area through time. Decreases in burned area are generally associated with increasing landscape fragmentation as a consequence of agricultural expansion, while increases in burned area are related to both climate change and the growing presence of human activities in the least converted parts of the Cerrado \citep{silva_putting_2021, da2024assessing, segura-garcia_human_2024}, as humans cause ignitions in relation to agricultural practices, accidents, fire management or in the conversion of certain vegetation forms \citep{pivello_use_2011, ribeiro_time_2024}. However, a clear delineation of the different regions of the Cerrado presenting different fire regimes and how these may have changed over time is still lacking.

In this study, we analysed the fire regimes of the Brazilian Cerrado in the period 1985-2020. First, we aimed to identify the different regions of the Cerrado presenting distinct characteristics of fire in four 9-year periods between 1985 and 2020, with a particular focus on the first (1985-1993) and last (2012-2020) periods. For this purpose, we combined monthly burned area data with 30 m spatial resolution with advanced spatially constrained clustering methods to delineate the regions of the Cerrado with different fire regimes based on a 30 km grid. Then, we studied the characteristics of these regions in terms of various factors related to the presence of human activities in the landscape and in terms of various climatic variables that may explain the differences in the fire regimes of the various regions identified. Finally, we analysed how these regions have changed over time, both geographically and in terms of their fire characteristics, with the ultimate aim of understanding how the rapid changes the Cerrado has experienced over the past four decades have shaped its fire regimes. 

\section{Methods}

\subsection{Study area}

The Brazilian Cerrado is formed by a range of vegetation growth-forms \citep{ribeiro_cerrado_2008}, predominantly savannas and grasslands. The Cerrado occupies an area over 2 million km$^2$ spanning eleven Brazilian states, and presents varying levels of anthropic land uses \citep{alencar_mapping_2020}. In the 1970s, large-scale agricultural expansion rapidly converted the Cerrado to anthropic land uses \citep{de_albuquerque_prado_expansao_2012, rocha_deteccao_2012}, affecting around 50\% of its native vegetation area \citep{alencar_mapping_2020}.

The Cerrado presents a latitudinal gradient of temperatures increasing towards the Equator, ranging from around 18$^\text{o}$C to 28$^\text{o}$C \citep{nascimento_clima_2020}. The climate of this biome is strongly seasonal, with a wet and a dry season (roughly May to October). In addition, there is a gradient of precipitation increasing from east to west varying between 800 mm and over 1,800 mm \citep{nascimento_clima_2020}.

To characterise the fire regimes of these savannas, we used a 30x30 km grid laid over the Cerrado and divided the study period in four 9-year periods: 1985-1993, 1994-2002, 2003-2011, 2012-2020. The size of the cells was chosen to have enough spatial resolution, while being large enough to simultaneously limit the number of fire polygons larger than one cell (see next Section), as well as minimising the number of cells with little to no fire events in a period.

\subsection{Data sources and processing}

\subsubsection{Burned area data and fire polygons}

We used burned area data from MapBiomas Fogo (Collection 2.0, \citep{alencar_long-term_2022}), a publicly available dataset covering the period 1985-2022 based on Landsat images (30 m spatial resolution with a 16-days interval) and created using machine learning algorithms in the Google Earth Engine platform \citep{gorelick_google_2017}. This dataset provides annual maps of burned area where pixels identified as burned are labelled with their mapping month. We downloaded the data using the MapBiomas toolkit for the aforementioned Google Earth Engine platform.

We processed this raster data using a custom algorithm to create fire polygons in vector format so that we could have an estimation of the size of the different fire events. Considering the limitations of Landsat data regarding the temporal resolution, as well as those possibly introduced by the mapping methodology, we mapped together as part of the same fire event (i.e. fire polygon) those spatially contiguous pixels belonging to consecutive months, as well as those polygons separated by up to 150 m from border to border. The reason was three-fold: first, two physically separated fires occurring close to one another in space and time may be part of the same fire event \citep{martin_fires_1992}; second, Landsat's spectral resolution and the mapping algorithm may fail to detect certain burned pixels within a fire event; and third, when simply mapping the polygons without this last condition, we found an excess number of small fires compared to the expected heavy-tail fire size distribution \citep{hantson_status_2016}. Finally, we kept only those fire polygons larger than 3 ha.

We quantified the fire regimes of a certain 30 km cell in a specific 9-year period by focusing our analyses on six fire variables: (1) \textbf{number of fires}, the number of polygons normalised to flammable area (see next Section); (2) the \textbf{50th quantile} and (3) \textbf{99th quantile} of the polygons' size distribution; (4) \textbf{burned area}, the sum of burned pixels in the period normalised to flammable area; (5) \textbf{fire frequency} as the 99th quantile of the smoothed Empirical Cumulative Distribution Function (sECDF) of the number of times each cell's pixel burned in the period (between 0 and 9); and (6) \textbf{fire seasonality}, as the peak of the fire season using the median value of the sECDF of monthly burned area in a fire-year. We established the fire-year for each cell as the 12-month period centered at the month with highest burned area. This consideration was important given the latitudinal gradient observed in the peak of the fire season . Henceforth, we will refer to the fire regime of a grid cell or region in a period as the combination of these six quantities.

Finally, a few grid cells in each 9-year period presented no burned area, indicating absence of fire. However, the clustering algorithm requires each cell to have a value for each fire variable. Therefore, we assigned these cells the minimum fire size (3 ha) both for the median and 99th percentile fire size and for burned area, and a frequency and number of fires of 1. We used the average peak of the fire season from the immediate 8 neighbouring cells. There were few cells with no fire occurrences, and these were generally spread out across the south of the Cerrado, with the exception of a small group of contiguous cells presenting no fires in the western-most tip of the Cerrado in the period (1985-1993), which were classified as part of cluster wLFA\_P1 (see Results and Figure \ref{fig:2maps_PCAs} for a definition of the clusters and their naming).

Both the creation of fire polygons and the quantification of the six different fire variables in each 30 km grid cell and period were carried out in Python (Python Software Foundation, \url{https://www.python.org/}), using the packages GeoPandas \citep{kelsey2020geopandas}, Pandas \citep{reback2020pandas} and Rasterio \citep{gillies2019rasterio}, as well as the QGIS \citep{QGIS_software} Python API.

\subsubsection{Land use area data and fragmentation indices}

We used annual maps of land use and land-cover from the MapBiomas project (Collection 8.0, \citep{souza_reconstructing_2020}). For each grid cell in the first (1985-1993) and last (20120-2020) periods, we calculated the average percentage of anthropic area (comprising "Urban Area", "Mining" and all land use categories under "Farming"), "Forest Formation", "Savanna Formation", "Grassland", and "Wetland". Additionally, we calculated the change in anthropic area within each period as the difference in anthropic area between the last and the first year in the period. Finally, we calculated the average flammable area as the sum of all those categories describing landscapes with enough vegetation to sustain a fire (that is, all categories except "Hypersaline Tidal Flats", "Rocky Outcrops", "Non Vegetated Area" and "Water"). We downloaded the land use and land-cover data using the MapBiomas toolkit for Google Earth Engine, and processed it using various Python packages (Python Software Foundation, \url{https://www.python.org/}).

We also calculated different types of landscape fragmentation indices using the Python package "PyLandStats" \citep{bosch_pylandstats_2019}. However, these indices were not used in the final analysis after inspection of collinearity and variable importance using Random Forests algorithm (see below).

\subsubsection{Climate reanalysis data}

For each 30 km grid cell in the first (1985-1993) and last (2012-2020) periods, we calculated the average annual and dry season temperature and precipitation, and the peak of the dry season. We used the 2 m temperature variable from ERA5-Land Monthly Averaged climate reanalysis data \citep{munoz-sabater_era5-land_2021} (9 km spatial resolution) from the European Centre for Medium-range Weather Forecasts (ECMWF). For precipitation, we used the 5-day precipitation accumulation in mm from the Climate Hazards Group InfraRed Precipitation with Station data \citep{funk_climate_2015} (CHIRPS, version 2.0), a quasi-global rainfall dataset with 0.05$^\text{o}$ spatial resolution. We chose the ERA5-Land dataset for temperature because the ECMWF data products have been identified as the most reliable air temperature datasets \citep{burton2018inter}. We used CHIRPS v2.0 for precipitation, as it is considered the best data source for precipitation across the tropics \citep{burton2018inter} - along with the Tropical Rainfall Measuring Mission (TRMM) \citep{liu2012tropical}. Finally, we used the potential evapotranspiration (PET) data provided by TerraClimate \citep{abatzoglou_terraclimate_2018}, a global dataset of climatic data with 0.5$^\text{o}$ spatial resolution. 

First of all, we calculated the average value of each climatic variable - temperature, precipitation and PET - for each 30 km grid cell, year and month as the average of all pixels of the corresponding data products intersecting the cell in the corresponding year and month (weighting each pixel's contribution to the average depending on the pixel's area within the cell, using the R package Exactextractr \citep{baston2020exactextractr}).

For each period and grid cell, we considered the dry season to be those months in which, on average, PET exceeded precipitation. We identified the peak of the dry season as the month with the most negative Climatological Water Deficit (CWD) on average. We calculated the CWD based on Malhi et al. 2009 \citep{malhi_exploring_2009}, as the cumulative difference between precipitation and PET ($\text{CWD}_m = \text{CWD}_{m-1} + \text{Prec.}_m - \text{PET}_m$, $\text{Max}(\text{CWD}_m) = 0$, where $m$ represents a certain month) in every 12-month cycle starting on the wettest month on average of each grid cell and period. Finally, we calculated the average monthly CWD for each period and grid cell, and identified the month with the most negative value.

Finally, for each 30 km grid cell and period, we calculated the average annual temperature as the average temperature across all months in the period, the dry season temperature as the average temperature across dry season months, the annual precipitation as the average total precipitation in a year, and the dry season precipitation as the average total precipitation across dry season months.

We downloaded the climatic data from the Google Earth Engine platform using the Geemap Python package \citep{wu2020geemap} and Jupyter Notebook \citep{kluyver_jupyter_2016}. In this case, we processed the data in R \citep{r2020r}, using the package Exactextractr \citep{baston2020exactextractr}.


\subsubsection{Population and livestock density data}

For each grid cell in the first (1985-1993) and last (2012-2020) periods, we calculated the median annual population density and livestock density. The density data were obtained from the Instituto Brasileiro de Geografia e Estatística (IBGE) \citep{IBGE_populacao_2023, IBGE_livestock_2023} at municipality level. We resampled these two sets of data from the municipality level to the grid level using various Python packages.


\subsection{Analysis}

Here, we used a combination of ordination and clustering techniques, as well as linear regressions, to address our research questions. As described with more detail below, we first summarised the main axes of variation of the six fire variables (burned area, frequency, largest and median fire sizes, peak of the fire season, and number of fires) using Principal Component Analyses (PCA) to characterise the fire regimes at the grid cell level in each 9-year period. We then used a spatially constrained clustering algorithm to identify the fire regime regions, that is, the geographic areas formed by each cluster of grid cells with similar combinations of fire variables (similar fire regimes). Next, we focused on the first (1985-1993) and last (2012-2020) periods as we were mostly interested in understanding how the fire regime regions and their climatic and anthropic characteristics may have changed between 1985 and 2020. For each of these two periods, we used PCAs and redundancy analysis to investigate the characteristics of the different fire regime regions in terms of the climate and the presence of human activities, and their association with the characteristics of fire. Finally, we explored how these fire regime regions changed through time geographically and also how their characteristics of fire evolved using linear regression on time series data.

\subsubsection{Clustering algorithm}

In preparation for the clustering algorithm, we calculated the fire regimes for each cell and 9-year period, transformed seeking normality by taking the natural logarithm of the median and 99th percentile fire size, burned area, number of fires, and fire frequency (after adding a unit). Then, for each 9-year period, we performed a PCA on the correlation matrix of the set of standardised fire variables \citep{borcard_numerical_2011} (henceforth PCA\_fire). We ran the clustering algorithm on the first three axes of each period's PCA as they described more than 95\% of the variation in the data.

To determine the different fire regime regions of the Cerrado, we used a hierarchical clustering algorithm with a spatial contiguity constraint \citep{guenard_hierarchical_2022} because we were interested in delineating spatially contiguous regions where fire regimes have similar features in a certain period, and investigate the environmental characteristics that may explain such geographical patterns. As a result of such constraint, areas with similar fire regimes but geographically disconnected would be identified as different clusters - i.e. different fire regime regions. We explored the fire regime characteristics of each cluster using the first two components of the PCA\_fire (PC1\_fire and PC2\_fire, Supplementary Information, Table S1).

Because the algorithm we used is an unsupervised method \citep{guenard_hierarchical_2022} - that is, the classification does not use any \textit{a priori} knowledge on the clustering of the data -, we explored a few metrics to identify an optimal clustering solution (i.e. the number of clusters) delineating compact and dissimilar clusters. We finally based the decision on the Bayesian Information Criteria (BIC) (Supporting Information, Figure S1), although also consulting the Dunn Index and a Silhouette plot.

Finally, we repeated the clustering procedure on the first (1985-1993) and last (2012-2020) periods excluding the peak of the fire season. The objective was to understand the distinctiveness of certain areas in absence of this fire attribute (see Results and Supplementary Information, Figure S9).

We carried out this part of the analysis using R \citep{r2020r}. In particular, we used the package Vegan \citep{r_vegan} to run the PCA and the function constr.hclust() \citep{guenard_hierarchical_2022} distributed with the package adespatial \citep{adespatial} to run the clustering algorithm.


\subsubsection{Environmental characteristics of the clusters}

Focusing on the first (1985-1993) and last (2012-2020) periods, we analysed the fire regime regions in terms of two sets of explanatory variables related to the level of anthropic activities in the landscape (percentage of anthropic area, the change in anthropic area percentage within the period, and the population and livestock densities), and to the climatic characteristics (average annual and dry season temperature, average total annual and dry season precipitation, and the peak of the dry season) of each fire regime region. For both groups of variables, we first assessed collinearity between variables and performed transformations as appropriate seeking to improve normality (Supplementary Information, Figures S3 and S4 and Table S8).

For each group of explanatory variables separately and together, we first ran a Random Forests algorithm \citep{breiman_random_2001} to assess the classification power of each set and evaluate the variables' predictive importance. Then, we performed a PCA on each set of variables in each of the two periods to analyse the climatic and anthropic characteristics of each region (henceforth PCA\_clim and PCA\_anth). Finally, we analysed the relation between each of the first components of the anthropic (PC1\_ant) and climatic (PC1\_clim) variables and the PC1\_fire. These post-hoc assessments of PC1 components and environmental variables were further tested using constrained ordinations - redundancy analyses (RDA), a method that allows to use multiple explanatory variables to model simultaneously a set of response variables \citep{borcard_numerical_2011} - to test for the statistical significance and explanatory power of the environmental variables sets to explain the variation in the multivariate space of fire variables.

To run the Random Forests algorithm, we used the R package randomForest \citep{randomForest}. For the various PCAs and RDAs, we again used the package Vegan \citep{r_vegan}.


\subsubsection{Temporal changes of the fire regimes}

We studied temporal changes in each of the six fire variables and the PC1\_fire. To make the periods comparable, we projected the PC1\_fire of each 9-year period on the PC1\_fire of the first period (Figure \ref{fig:map_changePCA_periods1_proj4_fixed}). Then, we calculated the temporal changes of PC1\_fire and of each fire variable by fitting a linear regression model to each cell using either the time period (for PC1\_fire and fire frequency) or the year (for the other fire variables) as the predictor (Figure \ref{fig:maps_change_fireVars} and Supplementary Information, Figure S7). Hence, we fitted
\begin{align}\label{1.1}
\begin{split}
\text{x}_{f} = \alpha + \beta \text{t},
\end{split}
\end{align}
where x$_{f}$ represents PC1\_fire or each of the fire variables, t is the time variable, $\alpha$ represents the average value of the grid cell at the beginning of each time series, and $\beta$ is the rate of change of x$_{f}$ over time. We fitted the model using Ordinary Least Squares (OLS) in Python with the statsmodels library \citep{seabold2010statsmodels}. We used an 80\% confidence interval to assess the statistical significance of the trends.

Finally, because the PC1\_fire - capturing the largest variation in the data - had a meaning in terms of the characterisation of the fire regimes of the various regions (see Results), we further compared the PC1\_fire value of each cell in the first (1985-1993) and last (2012-2020) periods to investigate whether this PCA\_fire score had increased or decreased over time.

\section{Results}

\subsection{The fire regime regions of the Cerrado}

\begin{figure}[t!]%
\centering
\includegraphics[width=\textwidth]{_chapter_02_revisions/Figure01.png}
\caption{Fire regime regions obtained for (a) the first period (1985-1993) and (c) the last period (2012-2020) corresponding to the different clusters identified by the spatially constrained clustering algorithm performed on a 30 km grid over the Cerrado. Plots (b) and (d) present, for each period respectively, the first two axes (PC1\_fire and PC2\_fire) resulting from the Principal Component Analysis performed on the six fire variables (PCA\_fire): the median and 99th percentile fire sizes, the burned area and the number of fires, the 99th percentile of the fire frequency distribution, and the median (peak) of the fire season. Each point in the PCAs represents a 30 km grid cell in the period coloured depending on the fire regime region it was classified into. Overall, in both periods, the PC1\_fire captures a gradient of fire activity (see Results) and the fire regime regions appear ordered along this axis. In the first period, five different fire regime regions were identified: one in the north characterised by small numerous fires and a later peak of the fire season (LS\_P1), a region in the centre-west with comparatively higher fire activity (HFA\_P1), another region in the centre-east with intermediate fire activity (IFA\_P1), and two regions with lower fire activity, one in the south (sLFA\_P1) and another one in the western-most tip of the Cerrado (wLFA\_P1). In the fourth and last period, three regions were identified: one in the north presenting a later peak of the fire season (LS\_P4), one in the centre with higher fire activity (HFA\_P4), and another one in the south with comparatively lower fire activity (LFA\_P4). Note that even though the same colour scheme was used in both periods to denote analogous fire regime regions, neither the geographical limits nor the overall fire characteristics of the equally-coloured fire regime regions between periods are necessarily the same.}
\label{fig:2maps_PCAs}
\end{figure}

Using the hierarchical clustering algorithm with a spatial contiguity constraint and based on the BIC criteria, we identified a different optimal number of fire regime regions (the different clusters of grid cells located in a certain area of the Cerrado) for each 9-year period: five regions for the first period (1985-1993), and three regions for the following three periods (1994-2002, 2003-2011, 2012-2020) (Figure \ref{fig:2maps_PCAs} and Supplementary Information, Figure S8).

Focusing on the first (1985-1993) and last (2012-2020) periods, the PCA\_fire in these two periods (Figures \ref{fig:2maps_PCAs} b and d, and Supplementary Information, Table S1) showed a similar composition of the first and second axes (PC1\_fire and PC2\_fire respectively) in terms of fire variables, and a similar proportion of the variance described (88.9\% and 87.9 in the first and last periods respectively). In both periods, burned area, fire frequency, largest fire size, and the number of fires contributed the most, and in the same direction, to PC1\_fire. Regarding PC2\_fire, it was the peak of the fire season and the number of fires in one direction, and the median fire size in the opposite direction, that contributed the most. Further, even though the contribution of the 99th percentile fire size to PC2\_fire was weaker than the three aforementioned variables, it did so in the same direction of the median fire sizes, and opposite direction to the number of fires. This opposite direction between the size and the number of fires indicated a trade-off between these two variables of fires, even though we still observed some cells with simultaneously large and rather numerous fires.

In both periods, the clusters appeared ordered along PC1\_fire indicating that, overall, cells with high levels of burned area, also had high levels of fire frequency, fire size or number of fires, and vice versa for cells with low and intermediate values (Figure \ref{fig:2maps_PCAs}). That is, cells with similar overall low, intermediate or high levels of the various variables were generally classified as part of the same fire regime region. Therefore, the PC1\_fire captured a gradient of increasing fire activity, and the level of fire activity (PC1\_fire) was the main factor differentiating the different clusters. 

In the first period (1985-1993), the clustering algorithm delineated a fire regime region whose grid cells presented the highest fire activity levels (henceforth HFA\_P1) overall when compared to the other four regions. This HFA\_P1 cluster was the largest region delineated in the first period (occupying 38\% of the Cerrado), located in the west of the Cerrado and spanning a large latitudinal range. East of this region, the algorithm identified a fire regime region that also covered a large latitudinal variation and presented intermediate values for all fire characteristics (intermediate fire activity, IFA\_P1). The algorithm identified two clusters with low fire activity when compared to the other fire regime regions, one located in the southern-most part of the Cerrado (sLFA\_P1) and another located in the western-most tip of the Cerrado (wLFA\_P1) - the smallest region identified in the first period. Finally, the last fire regime region delineated was located in the northern-most part of the Cerrado, standing out for presenting small but numerous fires and, particularly, a later peak of the fire season - the month when most fires occur - compared to the other regions (late season region, LS\_P1). Particularly, most fires in the LS\_P1 fire regime region took place between October and November, while the peak of the fire season for the other clusters took place between August and September (Supplementary Figure S2 f).

In the fourth and last period (2012-2020), the algorithm identified again a fire regime region with overall higher fire activity (HFA\_P4) when compared to the other fire regime regions, located in the central part of the Cerrado. In addition, the algorithm identified a region in the south with overall lower levels of fire activity (LFA\_P4), except for a set of cells showing some of the largest fires (highlighted in blue in Figures \ref{fig:2maps_PCAs} c and d). In this last period, however, this low fire activity region was much larger than in the first period, spanning 59.5\% of the Cerrado. Finally, there was only one fire regime region that persisted from the first period presenting, again, overall high number of fires and a later peak of the fire season (late season region, LS\_P4). In this last period, however, the size of the largest fires observed in LS\_P4 were larger than those in LS\_P1, and comparable in size to HFA\_P4 (Figure \ref{fig:2maps_PCAs} d). Again, while the peak of the fire season of grid cells in LS\_P4 took place in October and November; the peak of the fire season for the other two regions was between August and September (Supplementary Figure S2 f).

Finally, we also explored the relative importance of the peak of the fire season when delineating the fire regime regions by excluding this variable from the clustering. Even when not including the peak of the fire season, the fire regime region in the north was still delineated in the first period (Supplementary Information, Figure S9 a and b), while in the last period it was classified as part of the high fire activity region but still presented a distinctive combination of fire characteristics in comparison to other grid cells in the same fire regime region (Supplementary Information, Figure S9 c and d).

\subsection{Environmental characteristics of the fire regime regions}

\begin{figure*}[!p]%
\centering
\includegraphics[width=\textwidth]{_chapter_02_revisions/Figure02.png}
\caption{The two major axes of variation obtained from the PCA performed on each set of explanatory variables, anthropic and climatic variables. The anthropic variables (PCA\_anth) used in the first (a) and last periods (b) were anthropic area percentage, livestock density, population density, and change in anthropic area percentage within the period. The climatic variables (PCA\_clim) used in the first (c) and last (d) periods were the total annual precipitation, total precipitation in the dry season, average temperature in the dry season, and peak of the dry season. Points represent a 30 km grid cell in the corresponding period coloured according to the fire regime region they were classified into (see Figure \ref{fig:2maps_PCAs}). Note that even though the same colour scheme was used in both periods to denote analogous fire regime regions, neither the geographical limits nor the overall nor the fire, anthropic or climatic characteristics of the equally-coloured fire regime regions are necessarily the same.}
\label{fig:map_PCA_expVars}
\end{figure*}

\begin{figure*}[!p]%
\centering
\includegraphics[width=0.9\textwidth]{_chapter_02_revisions/Figure03.png}
\caption{Relation between the first component of the PCA performed on the anthropic variables (PC1\_anth), (a) and (b) in the first and last periods respectively, or the climatic variables (PC1\_clim), (c) and (d) in each period respectively, and the first component of the PCA performed on the fire variables (PC1\_fire). PC1\_fire in each period represents a gradient of increase in burned area, frequency, and the size and number of fires (Figure \ref{fig:2maps_PCAs} a, Supplementary Information, Table S1), capturing a gradient of fire activity. PC1\_anth in both periods represents a gradient of increase in anthropic area, population and livestock density (Figures \ref{fig:map_PCA_expVars} a and b, Supplementary Information, Table S4), that is, an increase in the presence of human activities and land use area in the landscape. Finally, PC1\_clim represents a gradient of flammability conditions, with an increase in dry season temperature and in annual precipitation (contributing to vegetation growth \citep{krawchuk_constraints_2011, pausas_global_2013}), and a decrease in dry season precipitation, when most fires occur \citep{silva_impacts_2019} (Figures \ref{fig:map_PCA_expVars} c and d, Supplementary Information, Table S6). Points represent a 30 km grid cell in the corresponding period coloured according to the fire regime region they were classified into (see Figure \ref{fig:2maps_PCAs}). Note that even though the same colour scheme was used in both periods to denote analogous fire regime regions, neither the geographical limits, nor the fire, anthropic or climatic characteristics of the equally-coloured fire regime regions are necessarily the same.}
\label{fig:expVars_relation}
\end{figure*}

Exploring the climatic characteristics and the level of human activities of the different regions, both PCA\_clim and PCA\_anth revealed that the fire regime regions obtained independently through the clustering algorithm appeared well separated by each set of explanatory variables in each period (Figure \ref{fig:map_PCA_expVars}), indicating that the variation in these variables shared geographical structure with that of the fire variables. Further, we evaluated the ability of the explanatory variables to classify the fire clusters running a Random Forests algorithm. The classification accuracy was very high in both periods for both the sets of anthropic (81.7 [78.1 84.9]\% and 82.3 [79.3, 85.9]\%, both with 95\% confidence interval, CI) and climatic (82.3 [78.7, 85.4]\% and 90.0 [87.1, 92.4]\%, 95\% CI) variables (Supplementary Information, Table S2). This accuracy was even higher when using both sets of variables together (88.7 [85.9, 91.2] \% and 94.0 [91.8 and 95.8] \% in each period respectively, 95\% CI). 

Regarding the PCA\_anth, in both periods PC1\_anth was contributed primarily by the percentage of anthropic area, and the livestock and human population densities, all pointing in the same direction (Figures \ref{fig:map_PCA_expVars} a and b, and Supplementary Information, Table S4). Therefore, this axis seemed to capture the variation in the level of anthropic presence and activities in the landscape. Moreover, the clusters appeared aligned along this gradient. In the first period, region sLFA\_P1 in the south presented the highest levels of anthropic activity, followed by regions IFA\_P1 in the centre-east, LS\_P1 in the north, and part of cluster HFA\_P1 in the centre-west; while region wLFA\_P1 in the western-most tip of the Cerrado and most of region HFA\_P1 presented the lowest levels of anthropic presence (Figure \ref{fig:map_PCA_expVars} a). In the last period, the region LFA\_P4 in the south had the highest levels of anthropic presence, followed by regions LS\_P4 in the north and HFA\_P4 in the centre, in opposite order to the gradient of fire activity of these regions (Figure \ref{fig:map_PCA_expVars} b). In both periods, the change in anthropic area within the period was the variable that contributed the most to PC2\_anth. 

Considering the PCA\_clim, both in the first and last periods the dry season temperature and mean annual precipitation in one direction, and the dry season precipitation in the opposite, contributed the most to PC1\_clim, while the peak of the dry season contributed the most to PC2\_clim (Figures \ref{fig:map_PCA_expVars} c and d, and Supplementary Information, Table S6). In this case, PC1\_clim seemed to capture a gradient of availability of fuel and suitability of conditions to burn, that is, the flammability: hotter and drier climate in the dry season, when most of the fires occur \citep{silva_impacts_2019}; and amount of annual precipitation, which contributes to vegetation growth \citep{krawchuk_constraints_2011, pausas_global_2013}. In this case, however, the clusters aligned along the direction of dry season temperature - hence following the latitudinal gradient -, and less perfectly along the flammability gradient (PC1\_clim). Finally, in both periods, the northern regions LS\_P1 and LS\_P4 presented a later peak of the dry season than the other regions.

Lastly, since the PC1 axes represented environmental gradients, we looked at their relation with the fire activity gradient (Figure \ref{fig:expVars_relation}). On the one hand, PC1\_anth and PC1\_fire were inversely related: as the presence of humans and their activities in the landscape decreased (lower PC1\_anth), fire activity increased (higher PC1\_fire) for all clusters in both periods (Figures \ref{fig:expVars_relation} a and b). The only exception to this relation were region wLFA\_P1 and some cells in the northern region LS\_P1 that presented a lower fire activity than other cells with the same level of human presence. On the other hand, PC1\_clim and PC1\_fire were directly related: as the flammability increased (higher PC1\_clim), fire activity also increased in both periods (Figures \ref{fig:expVars_relation} c and d). Again, most cells in regions wLFA\_P1 and LS\_P1 presented lower fire activity than those cells in region HFA\_P1 with similar flammability scores.

These negative and positive correlations with fire activity of human presence and flammability respectively were further confirmed when modelling the correlations of each set of explanatory variables with the fire variables using a redundancy analysis approach (r$^{2}_{\text{adj}}$=0.316, 0.305 for anthropic and r$^{2}_{\text{adj}}$=0.242, 0.305 for climatic variables in each period respectively, and p-value<0.001 for all variables and cases; Supplementary Information, Figure S6 and Tables S5 and S7). These models further revealed that precipitation was only weakly correlated with the fire variables, and that there was a strong correlation between the peak of the dry season and that of the fire season (Supplementary Information, Figures S6 c and d, and Table S7).


\subsection{Changes of fire regime regions over time}

\begin{figure*}[t!]%
\centering
\includegraphics[width=\textwidth]{_chapter_02_revisions/Figure04.png}
\caption{Bivariate maps showing the PCA\_fire coordinates of each 30 km grid cell in the first and last periods projected on the PCA\_fire-space obtained for the first period PCA\_fire. Map (a) shows the PCA\_fire scores directly obtained from the fire variables in the first period (1985-1993). Map (b) shows the projection of the fire variables of the last period (2012-2020) on the PCA\_fire-space of the first period. In this manner, the fire regimes (6-dimensional) are summarised into two dimensions, and the values of the grid cells in each period can be compared. PC1\_fire (x-axis of the colour map) represents an increase in fire activity (larger burned area, frequency, and larger fire sizes) from white to blue, while PC2\_fire (y-axis) represents the direction of a later fire season, an increase in the number of fires and a decrease in median fire size from white to yellow.}
\label{fig:map_changePCA_periods1_proj4_fixed}
\end{figure*}

\begin{figure*}[t!]%
\centering
\centerline{\includegraphics[width=1\textwidth]{_chapter_02_revisions/Figure05.png}}
\caption{Relation between PC1\_fire obtained in the first period, 1985-1993, (x-axis), and the projection of PC1\_fire obtained for the last period, 2012-2020, projected on the PCA\_fire of the first period (y-axis). (a) subset of 30 km grid cells classified into the northern fire regime region in the last period with small numerous fires and a later peak of the fire season, LS\_P4; (b) subset of grid cells classified into the high fire activity region of the last period, in the centre of the Cerrado, HFA\_P4; and (c) cells classified into the lower fire activity region of the last period, in the south, LFA\_P4. Points represent 30 km grid cells and are coloured depending on the fire regime region they were classified into in the first period if they present clear trends in time in PC1\_fire (with 80\% confidence interval), or coloured in grey if they do not present clear trends in time. The dashed black lines represent the one-on-one correspondence between PC1\_fire in the first period and the projection of PC1\_fire in the last. Thus, points above the dashed black line represent an increase in fire activity over time, while data points below the line represent a decrease. Dashed gray lines represent the 0 values of PC1\_fire and are included for reference.}
\label{fig:change_PC1}
\end{figure*}

\begin{figure*}[htbp!]%
\centering
\includegraphics[width=\textwidth]{_chapter_02_revisions/Figure06.png}
\caption{Trends in time (slope coefficients) obtained by fitting a linear regression (using Ordinary Least Squares) to the time series of each of the six fire variables, for each 30 km grid cell over the Cerrado. Plots show the median fire size (a), the 99th percentile fire size (b), burned area (c), the number of fires (d), and the median fire season (f). The time series for these five variables comprised the annual values for each year in the period 1985-2020. The time series for the 99th percentile of the fire frequency distribution (e) comprised the frequency value in each of the four periods, since this value expresses the largest number of times the pixels in a cell burned in 9 years. Pink indicates decreases in the fire variable over time, green indicates increases, and grey indicates no clear trends in time with 80\% confidence interval.}
\label{fig:maps_change_fireVars}
\end{figure*}


Between the first (1985-1993) and last (2012-2020) periods, there was a reduction in the number and a change in the shape of the fire regime regions (Figures \ref{fig:2maps_PCAs} and \ref{fig:map_changePCA_periods1_proj4_fixed}). Moreover, the last period clusters were formed by grid cells coming from a mix of different first-period regions (Figure \ref{fig:map_changePCA_periods1_proj4_fixed}).

Apart from these geographical changes, the fire characteristics of the individual grid cells also changed over time. In particular, there were increases and decreases in fire activity heterogeneously distributed across the Cerrado (the PC1\_fire score, Figure \ref{fig:change_PC1} and Supplementary Information, Figure S7). Thus, even though the fire regime regions were organised along a gradient of fire activity in both periods, the levels of high fire activity of region HFA\_P1 and low levels of regions sLFA\_P1 and wLFA\_P1 in the first period, did not perfectly correspond to the levels of high fire activity of region HFA\_P4 and low fire activity of region LFA\_P4 in the last period (Figures \ref{fig:change_PC1} b and c). Furthermore, the last period clusters were not the result of unidirectional changes, but rather the result of changes in different directions of each of the different fire variables (Figure \ref{fig:maps_change_fireVars}).

Clusters LS\_P1 and LS\_P4 in the north comprise the only region whose geographical shape and location did not substantially change between periods, and maintained its unique attribute of a later fire season peak. Regarding fire activity, however, 29\% of its area experienced an increase in fire activity over the study period (Figure \ref{fig:change_PC1} a). On the other hand, region HFA\_P4 in the centre of the Cerrado - with the highest levels of fire activity in the last period - was composed of cells showing increases (14.3\% of its area) and decreases (15.9\%) in fire activity over time (Figure \ref{fig:change_PC1} b). However, these decreases were limited and remained mostly within the range of activity of cluster HFA\_P1 (around a minimum value of 0 in the PC1\_fire, Figure \ref{fig:change_PC1} b). Regarding region LFA\_P4 in the south - the cluster with lowest fire activity in the last period -, 27.1\% of its cells experienced decreases in fire activity (Figure \ref{fig:change_PC1} c) over the study period. Most of its cells transitioned from clusters HFA\_P1 and IFA\_P1 - initially with high and intermediate fire activity - experiencing mild to strong decreases in fire activity (47.9\% and 22.8\% of each cluster respectively), while 16.9\% of the cells transitioning from region sLFA\_P1 - already with low levels of fire activity - experienced further reductions. 

\section{Discussion}

In this study, we identified the regions of the Cerrado presenting different fire regimes, the fire regime regions, using a hierarchical clustering algorithm with a spatial contiguity constraint in each of four 9-year periods between 1985 and 2020. We identified five fire regime regions in the first period, 1985 to 1993, and three fire regime regions in the last period, 2012 to 2020, and in the two periods in-between (Figure \ref{fig:2maps_PCAs} and Supplementary Information, Figure S8). These regions reflect the spatial patterns of fire across the Cerrado and allowed us to explore how its fire regimes have changed over time. We found considerable differences in the geographical shape of the fire regime regions between the first and last periods, as well as in the fire characteristics of the individual 30 km grid cells that form the different regions (Figures \ref{fig:map_changePCA_periods1_proj4_fixed} and \ref{fig:change_PC1}).

In each period, the fires regime regions of the Cerrado were primarily differentiated by their level of fire activity, from low fire activity - that is, with overall infrequent, small, and non-numerous fires - to high fire activity - with overall frequent, large, and numerous fires - (Figures \ref{fig:2maps_PCAs} b and c). This gradient of fire activity was captured by the first component of the PCA performed on the six fire variables (PCA\_fire), along which the fire regime regions appear ordered (Figure \ref{fig:2maps_PCAs}). Therefore, the primary feature differentiating the fire regimes of the Cerrado is their overall level of fire activity, prevailing over other differences in fire characteristics, in agreement with previous studies characterising the fire regimes of the Cerrado \citep{silva_putting_2021}. Still, we observed a certain trade-off between the size and the number of fires (Figures \ref{fig:2maps_PCAs} b and d), indicating that the grid cells with the largest fires did not experience the highest numbers of fires, and vice versa. The trade-off between these two variables of fire is related to the fact that when a large fire occurs, it consumes a large portion of the landscape's fuel bed, thus inevitably hampering the occurrence of other fires in the same period and limiting the overall number of fires, as reported in previous studies \citep{conciani_human-climate_2021, batista_evaluation_2018, alvarado_management_2018, fidelis_year_2018}.

The Cerrado is very heterogeneous in terms of vegetation \citep{ribeiro_cerrado_2008, klink_conservation_2005}, topography and soil composition \citep{sano_cerrado_2019}, climate \citep{nascimento_clima_2020}, hydrology \citep{oliveira_trends_2014}, and land use area \citep{souza_reconstructing_2020}. Based on these factors, different ecoregions have been described in the Cerrado \citep{sano_cerrado_2019}, and their patterns of fire subsequently quantified \citep{silva_putting_2021}. Because these factors affect the patterns of fire \citep{archibald_what_2009}, we could expect some similarities between the fire regimes of these ecoregions \citep{silva_putting_2021} and the fire regime regions identified here. Even though neither the number nor the overall geographical shape of the ecoregions (Supplementary Information, Figure S10) corresponds to those of the fire regime regions identified (Figure \ref{fig:2maps_PCAs}), we did find some similarities. For example, the ecoregions located in the centre of the Cerrado were found to present large burned area and fire sizes overall \citep{silva_putting_2021}, in agreement with the high and intermediate fire activity regions identified here (Figure \ref{fig:2maps_PCAs}). In particular, the Bananal ecoregion in the centre-west of the Cerrado (Supplementary Information, Figure S10, \citep{sano_cerrado_2019}) appears clearly delineated as part of the high fire activity region of the last period (HFA\_P4), presenting some of the highest levels of fire activity (Figure \ref{fig:map_changePCA_periods1_proj4_fixed}), in accordance with this ecoregion's fire regimes \citep{silva_putting_2021, da2024assessing}. Similarly, the ecoregion Chapadão do São Francisco (Supplementary Information, Figure S10, \citep{sano_cerrado_2019}), characterised with large burned areas and fire sizes \citep{silva_putting_2021}, can also be recognised as part of the high fire activity region of the first period (HFA\_P1).

We did find, however, a remarkable difference between our fire regime regions and the ecoregions located in the south-west of the Cerrado, comprising the states of Mato Grosso and Mato Grosso do Sul (Supplementary Information, Figure S10). The ecoregion covering this area, Paraná Guimarães  (Supplementary Information, Figure S10, \citep{sano_cerrado_2019}), was classified as having intermediate levels of burned area and big fires in the period 2001-2019 \citep{silva_putting_2021}. However, we observed the lowest levels of fire activity in the period 2012-2020 precisely within this area (LFA\_P4, Figure \ref{fig:map_changePCA_periods1_proj4_fixed}). This difference in the level of fire activity could be related to the fact that the Paraná Guimarães ecoregion includes part of the high fire activity region (HFA\_P4) identified here, which could be contributing to the higher burned area and larger fire sizes reported for the whole ecoregion \citep{silva_putting_2021}. This highlights the relevance of exploring patterns of fire at different scales, as it can reveal different processes and dynamics \citep{krawchuk_constraints_2011}. 

The fire regime regions we identified are overall large and somewhat heterogeneous in the fire characteristics of their grid cells (Figure \ref{fig:map_changePCA_periods1_proj4_fixed}). For example, we observed the subregion of large fires located within the low fire activity region of the last period (LFA\_P4, highlighted in blue in Figures \ref{fig:2maps_PCAs} c and d). Hierarchical clustering algorithms with a spatial contiguity constraint generate solutions (i.e. classifications for different numbers of clusters) by grouping together grid cells and clusters consecutively only if they are spatially contiguous \citep{borcard_numerical_2011, guenard_hierarchical_2022}. Here, we focused on the first optimal clustering classification of the data, which we found to discriminate generally large fire regime regions primarily in relation to their level of fire activity. However, there could be further statistically relevant clustering solutions delineating clusters at smaller scales \citep{lambiotte_random_2014}, which could uncover further patterns of fire regimes perhaps recovering some of the Cerrado's ecoregions, opening an avenue for future studies of the Cerrado's fire regimes.

In terms of the two major drivers of fire, humans and the climate \citep{archibald_managing_2016, archibald_southern_2010, hantson_global_2015}, both the anthropic and the climatic variables have high classification accuracy of the fire characteristics (Supplementary Information, Table S2), as expected \citep{bowman_fire_2009, archibald_what_2009}. We found that the fire regime regions appear distributed along the first axis of the PCA performed on the anthropic variables (PC1\_anth, Figures \ref{fig:map_PCA_expVars} a and b), which indicates a gradient of increasing presence of human activities. This anthropic presence is particularly associated with agriculture - land use area and livestock density - and, to a lesser extent, with population density (Figure \ref{fig:map_PCA_expVars} a). Importantly, we observe that as anthropic presence increases, fire activity decreases (Figure \ref{fig:expVars_relation}), in agreement with previous studies relating the expansion of land use area and fragmentation in savannas with decreases in burned area \citep{silva_putting_2021, segura-garcia_human_2024, da2024assessing, archibald_southern_2010, andela_human-driven_2017, andela_recent_2014}. 

Indeed, both land use area and livestock grazing are associated with smaller fires and reduced burned area as they create fuel discontinuities through landscape fragmentation \citep{hantson_global_2015, archibald_what_2009} and diminishing fuel loads \citep{archibald_evolution_2012, eloy_seasonal_2019}. Population density, however, has a more intricate relation with fire activity as humans can ignite fires in higher or smaller numbers \citep{hantson_global_2015, archibald_what_2009}, but they can also limit fire activity through fire suppression efforts \citep{archibald_managing_2016, archibald_evolution_2012} - a widespread practice in the Cerrado \citep{durigan_zero-fire_2020, durigan_need_2016, fidelis_year_2018}. On the other hand, we observed the highest level of fire activity in the areas of the Cerrado with larger remnants of native vegetation (Supplementary Information, Figures S2 and S3), in agreement with previous studies \citep{mataveli_satellite_2018, silva_putting_2021, segura-garcia_human_2024, da2024assessing}. These dynamics in the more intact areas of the Cerrado may result from a confluence of factors comprising human-started fires escaping from nearby agricultural lands into native vegetation \citep{pivello_use_2011} and the timing of ignitions in the dry season \citep{pivello_use_2011}, when the fuel bed of native vegetation is more flammable and effectively more continuous \citep{zanzarini2022flammability}. 

The western fire regime region in the first period ( wLFA\_P1) is, however, an exception to the inverse relation between the gradient of human presence and that of fire activity (Figure \ref{fig:2maps_PCAs} a and b). This fire regime region simultaneously presents both the lowest fire activity levels and the lowest presence of human activities in this period (Figures \ref{fig:2maps_PCAs} b and \ref{fig:map_PCA_expVars} a). The low fire activity of this region could be the result of a vegetation composition closer to that of the Amazon tropical rainforest, where natural fires are extremely rare \citep{barlow_clarifying_2020, rosan_fragmentation-driven_2022} due to moister conditions and a closed canopy with scant surface fuels, particularly when human presence is low. The contrast in terms of human presence between this region and the other low fire activity region in the same period (sLFA\_P1) also shows how areas with different landscape characteristics can present similar fire regimes \citep{archibald_defining_2013}.

The distribution of the fire regime regions on the PCA\_clim seems to reflect their geographical location, as they appear ordered along the latitudinal gradient of temperature \citep{nascimento_clima_2020}, with the northern regions (LS\_P1 and LS\_P4) presenting the hottest values, and the southern regions (sLFA\_P1 and LFA\_P4) presenting the lowest values (Figures \ref{fig:map_PCA_expVars} c and d). Perpendicularly, the fire regime regions appear organised along the direction of increase in annual precipitation from east to west \citep{nascimento_clima_2020}. For example, in the first period, the intermediate fire activity (IFA\_P1) in the east presents lower annual precipitation than the high fire activity (HFA\_P1) and the western low fire activity (wLFA\_P1) regions. Still, the first axis of the PCA\_clim captures a gradient of climatic flammability conditions (Figures \ref{fig:map_PCA_expVars} c and d), which is positively related with fire activity (Figures \ref{fig:expVars_relation} c and d). That is, fire activity increases with drier and hotter conditions in the dry season as expected - when most fires occur \citep{mataveli_satellite_2018, silva_impacts_2019} -\citep{silva_drivers_2020, silva_impacts_2019, westerling_warming_2006, flannigan_global_2013}, as well as with greater annual precipitation, which contributes to productivity during the growing season \citep{pausas_global_2013} and fuel load \citep{krawchuk_constraints_2011}. However, we observed that, comparing cells with similar high climatic flammability values, some cells with high presence of anthropic activities in the last period (region LFA\_P4) have lower fire activity than cells with low anthropic activity levels (region HFA\_P4), suggesting an interaction between humans and climate. Indeed, it has been demonstrated that the effect of the climate on fire is subject to the level of anthropic presence in savannas \citep{conciani_human-climate_2021, segura-garcia_human_2024, andela_recent_2014}.

Only one fire regime region was consistently identified in each period (LS\_P1 and LS\_P4, Figure \ref{fig:2maps_PCAs} and Supplementary Information, Figure S8) presenting a later peak of the fire season. In addition, this area presents small but numerous fires, which makes this region stand out even when excluding the peak of the fire season from the clustering analysis (Supplementary Information, Figure S9). That is, this region's fire regimes are distinct from those of the other fire regime regions and not merely shifted in time. A number of factors may explain the particularities of this region's fire regimes. First, this area comprises the Floresta de Cocais and Costeiro ecoregions (Supplementary Information, Figure S10, \citep{sano_cerrado_2019}), which is noteworthy for being an ecological transition region between tropical dry forest, rainforest and savanna \citep{sano_land_2010} composed primarily of broadleaf tree species \citep{oliveira-filho_neotroptree_2017} and palm trees \citep{sano_cerrado_2019}. In fact, in early versions of the Cerrado's geographical limits, this area was not fully included and was classified instead as either a mix of Amazon and Cerrado \citep{ibge_mapa_1993} or as dry broadleaf forest \citep{olson_terrestrial_2001}. Second, this area also presents unique climatic characteristics as it is affected by the Intertropical Convergence Zone, which carries air moisture from the Atlantic towards the Amazon rainforest \citep{nascimento_clima_2020}, and presents an overall later peak of the dry season (Figure \ref{fig:map_PCA_expVars}). Third, in terms of human presence, this area presents a number of municipalities with high densities of \textit{quilombola} population \citep{censo2022quilombola}, communities of people with African descent \citep{censo2022quilombola} who traditionally use fire to renew native pastures or as part of shifting cultivation practices, among other uses \citep{eloy_seasonal_2019}. Finally, this area is found within the Matopiba region, considered the Cerrado's current agricultural frontier and where land conversion has been more widespread in recent years  \citep{spera_land-use_2016, silva_drivers_2020, teixeira_pitta_impacts_2017}.

We observed opposing trends over time in fire activity across the Cerrado (Figure \ref{fig:change_PC1}), in agreement with previous studies \citep{silva_putting_2021, segura-garcia_human_2024}. In the northern half of the Cerrado, comprising the Matopiba - encompassing regions LS\_P4 and HFA\_P4 -, we observed the largest concentration of increases in fire activity. The surges in fire activity of these regions could be the result of the use of fire in land conversion and as part of modern agricultural practices \citep{klink_conservation_2005, pivello_use_2011} alongside large native vegetation remnants with little fuel discontinuities, and in the context of climate change. On the other hand, the southern half of the Cerrado mostly experienced decreases in fire activity. Such reduction was particularly striking in the state of Mato Grosso (Supplementary Information, Figures S7 and S10), in agreement with previous studies \citep{silva_putting_2021, da2024assessing}, which has undergone extensive Cerrado vegetation fragmentation due to conversion to pastures and croplands in the 1990s and early 2000s \citep{rocha_deteccao_2012, alencar_mapping_2020}. Further, the shape of the low fire activity region in the southern half of the Cerrado in the last period (LFA\_P4) clearly reflects the geographical patterns of anthropic land uses \citep{alencar_mapping_2020}, highlighting the depth of the footprint of land use area on the Cerrado's fire regimes.

Finally, between 1985 and 2020, many aspects of the fire management policies across Brazil have changed. Indeed, at the beginning of the study period, the use of fire had been banned by successive laws since colonial times \citep{pivello2021understanding}, further supported by misunderstandings about the ecology of fire-adapted ecosystems and the occurrence of destructive fires in iconic fire-sensitive ecosystems. In the late 1990s, however, there was a shift in the discussion as fire started to be recognised as a disturbance agent that could promote and sustain Cerrado biodiversity \citep{durigan_zero-fire_2020}. There was hence a push to (re)introduce fires into protected areas, which eventually translated into fire policy changes \citep{pivello2021understanding}. Indeed, the 2012 New Forest Code contemplated the use of fire in fire-adapted ecosystems for ecological and cultural purposes. Therefore, part of the increase in fire activity that we observe in certain areas of the Cerrado could be related to the use of prescribed fires as a management tool. However, prescribed burns are usually performed in the early dry season \citep{schmidt_fire_2020, batista_evaluation_2018}, while most Cerrado fires currently occur later in the dry season \citep{da2024assessing}, indicating that other uses of fire may still prevail. In this sense, the recently approved Brazilian National Policy on Integrated Fire Management (IFM) is a promising avenue for Cerrado conservation, as these practices are demonstrating to be an effective tool in Cerrado conservation and preventing more extreme fire behaviour \citep{schmidt_fire_2020, fidelis_year_2018}.


\section{Conclusion}\label{sec12}

Our findings indicate that the main factor differentiating Cerrado's fire regime regions is the level of fire activity, which is inversely related to the extent of modern human presence in the landscape. On the one hand, we report that recent agricultural expansion and fragmentation have profound effects on the Cerrado's fire regimes by primarily diminishing fire activity. Indeed, high livestock densities, modern mechanised agricultural practices and croplands such as soybean plantations where fire use is limited now occupy large areas of the Cerrado. On the other hand, climate change could be setting the stage for increased fire sizes and frequency in areas where most native vegetation remains. Therefore, both agricultural expansion and climate change could impact Cerrado's vegetation and contribute to its homogenization and loss of biodiversity, highlighting the profound multi-faceted effects of large-scale Cerrado native vegetation loss. Future research efforts should be allocated to decipher the variation in fire patterns at finer spatial scales, as well as establishing in more detail the causal relations between the different drivers and specific fire characteristics. This knowledge about the Cerrado's fire regimes is essential to develop meaningful policies aimed at conserving this unique biome in the context of global change.

\section{Data and code availability statement}

\sloppy The climate data used in this study are publicly available from ERA5-Land Monthly Averaged at \url{https://doi.org/10.5194/essd-13-4349-2021}, the Climate Hazards Group InfraRed Precipitation with Station data version 2.0 at \url{https://doi.org/10.1038/sdata.2015.66}, and TerraClimate at \url{https://doi.org/10.1038/sdata.2017.191}. The land use data are publicly available from MapBiomas collection 8.0 at \url{https://doi.org/10.3390/rs12172735}. The burned area data are publicly available from MapBiomas Fogo Collection 2.0 at \url{https://doi.org/10.3390/rs14112510}. Livestock density data for Brazil is available at \url{https://sidra.ibge.gov.br/tabela/3939} and population density of Brazil \url{https://www.ibge.gov.br/estatisticas/sociais/populacao/9103-
estimativas-de-populacao.html}. The processed data (including the custom fire polygons created from MapBiomas Fogo Collection 8.0) and custom code are available at \url{https://doi.org/10.5281/zenodo.14025420}.


\section{Acknowledgements}

We thank the Ecosystems Lab Group at the Environmental Change Institute, School of Geography and the Environment, University of Oxford, and the Science Team at the Instituto de Pesquisa Ambiental da Amazônia (IPAM), Brazil, for helpful feedback during discussions. We thank Prof. Emilio Chuvieco, Dr Blas Kolic and Francisco Navarro-Rosales for their comments and suggestions. We would like to thank the reviewers for providing useful feedback to improve the readability and content of the manuscript.


\begingroup
\setstretch{1.2}
% \setlength\bibitemsep{0pt}
\printbibliography[heading=subbibnumbered]
\endgroup


\end{refsection}

\clearpage
\newpage


\begin{refsection}
    


\supplementarysection

\section{Supplementary Information}


\begin{figure*}[htbp]%
\centering
\includegraphics[width=0.8\textwidth]{_chapter_02_revisions/BIC_grid30km_period1&4.pdf}
\caption{Bayesian Information Criteria index calculated for each clustering solution with an increasing number of clusters between 2 and 25 for the first (a) and last (b) periods.}
\label{fig:BIC_period1_4}
\end{figure*}


\begin{figure*}[htbp]%
\centering
\includegraphics[width=\textwidth]{_chapter_02_revisions/boxplot_fireVars_bothPeriods_review_newLabels.png}
\caption{Distribution of each of the six fire variables for each fire regime region in the first (1985-1993, shaded in blue) and last (2012-2020, shaded in orange) periods. Plots show the median fire size (a), the 99th percentile fire size (b), the burned area (c), the number of fires (d), the 99th percentile fire frequency (e), and median fire season (f), expressed in months.}
\label{fig:boxplot_fireVars_bothPeriods}
\end{figure*}


\begin{figure*}[htbp]%
\centering
\includegraphics[width=\textwidth]{_chapter_02_revisions/boxplot_expVars_anthropic_bothPeriods_reviews_newLabels.png}
\caption{Distribution of each anthropic variable for each fire regime region in the first (1985-1993, shaded in blue) and last (2012-2020, shaded in orange) periods. Plots show the anthropic area percentage (a), the change in anthropic area percentage within the period (b), the patch density (c), the population density (d), and the livestock density (e).}
\label{fig:boxplot_expVars_anthropic_bothPeriods}
\end{figure*}

\begin{figure*}[htbp]%
\centering
\includegraphics[width=\textwidth]{_chapter_02_revisions/boxplot_expVars_climatic_bothPeriods_reviews_newLabels.png}
\caption{Distribution of each climatic variable for each fire regime region in the first (1985-1993, shaded in blue) and last (2012-2020, shaded in orange) periods. Plots show the mean annual precipitation (a), mean annual temperature (b), dry-season length (c), precipitation in the dry season (d), the mean annual temperature in the dry season (e), and the peak of the dry season (f).}
\label{fig:boxplot_expVars_climatic_bothPeriods}
\end{figure*}


\begin{figure*}[htbp]%
\centering
\includegraphics[width=\textwidth]{_chapter_02_revisions/boxplot_expVars_vegType_bothPeriods_reviews_newLabels.png}
\caption{Distribution of the percentage area covered by different vegetation and land use categories (as per MapBiomas Collection 8.0 categories, see Methods) for each fire regime region in the first (1985-1993, shaded in blue) and last (2012-2020, shaded in orange) periods. Plots show savanna (a), forest (b), grassland (c), wetland (seasonally flooded savanna) (d), pasture (e), and cropland (f) land use categories.}
\label{fig:boxplot_expVars_vegType_bothPeriods}
\end{figure*}


\begin{figure*}[htbp]%
\centering
\includegraphics[width=0.9\textwidth]{_chapter_02_revisions/map_RDA_periods1_4_varsSubsetRF_review_newLabels.png}
\caption{Results of the redundancy analysis (RDA) performed using each set of anthropic (top row) and climatic parameters (bottom row) as explanatory variables and the six fire characteristics as response variables, for the first (1985-1993, left column) and last (2012-2020, right column) periods. The anthropic variables comprise anthropic area percentage, livestock density, population density, and change in anthropic area percentage within the period. The climatic variables comprise the total annual precipitation, total precipitation in the dry season, mean temperature in the dry season and the peak of the dry season. The projection used to represent the RDA results is such that the angles between the explanatory and response variables represent the correlation between the variables. The data points are projected in the space of the response variables and coloured according to the fire regime region they were classified in for each period (see Methods). Note that even though the same colour scheme was used in both periods to denote analogous fire regime regions, neither the geographical limits, nor the fire, anthropic or climatic characteristics of the equally-coloured fire regime regions are necessarily the same. All variables in all cases had a p-value $<$ 0.001. }
\label{fig:map_RDA_expVars}
\end{figure*}



\begin{figure*}[htbp]%
\centering
\includegraphics[width=0.8\textwidth]{_chapter_02_revisions/maps_change_PC1_sig_review.png}
\caption{Trends in time (slope coefficients) over the four 9-year periods of fire activity (see Methods) for each 30 km grid cell in the Cerrado. The time series of fire activity comprised the PC1\_fire values for each period projected on the PCA\_fire space of the first period (see main, Figure 1) so that these consecutive values are comparable. The time series was calculated as a linear regression and fitted using Ordinary Least Squares. Pink indicates decreases in fire activity over time, green indicates increases, and grey indicates no clear trends in time with 80\% confidence interval.}
\label{fig:maps_change_PCA_sig}
\end{figure*}




\begin{figure*}[htbp]%
\centering
\includegraphics[width=0.8\textwidth]{_chapter_02_revisions/map_PCA_allPeriods_reviews_newLabels.png}
\caption{Maps showing the fire regime regions identified by the spatially constrained hierarchical algorithm (see Methods) accompanied by plots of the two main axes of variation of the PCA\_fire for the period 1985-1993 (a-b), 1994-2002 (c-d), 2003-2011 (e-f), and 2012-2020 (g-h). Percentage values indicate the explanatory power of each PCA axis.}
\label{fig:map_PCA_allPeriods_fixed}
\end{figure*}

\begin{figure*}[htbp]%
\centering
\includegraphics[width=\textwidth]{_chapter_02_revisions/noSeason_map_PCA_periods1_4_review_newLabels.png}
\caption{Maps (a and c) showing the alternative fire regime regions obtained in the first and last periods respectively when the peak of the fire season is not included as a fire variable when running the spatially constrained clustering algorithm. Figures (b) and (d) show the first two axes of the PCA\_fire used to classify the grid cells. These maps and PCAs were generated to explore whether regions LS\_P1 and LS\_P4 in each period would present any distinct characteristic in addition to the later peak of the fire season, and would hence be delineated when running the clustering algorithm without this variable. In this case, the Bayesian Information Criteria (BIC) suggests that four regions in the first period and two regions in the last period are the clustering solutions optimising data partition. In the first period, the four-clusters solution does identify a region equivalent to LS\_P1, confirming the distinct fire regime of this region compared to HFA\_P1 (and its HFA\_P1$_\text{alt.}$ equivalent alternative), with an overall lower fire activity and overall small and numerous fires. In the last period, the equivalent LS\_P4 region is not identified in the two-cluster solution suggested by the BIC, suggesting that the increases in fire activity experienced by this area over time resulted in fire regimes closer to those of region HFA\_P4 (and its HFA\_P4$_\text{alt.}$ equivalent). Nevertheless, this northern area is identified in the four-clusters solution for the last period when leaving the peak of the fire season out (outlined in green in c and d, LS\_P4$_\text{alt.}$), suggesting this area presents characteristics relevant to data classification. Moreover, despite the fire characteristics of the northern region being closer to those of HFA\_P4$_\text{alt.}$, we can observe these cells are still clearly clustered in the PCA\_fire with numerous and medium-sized fires.}
\label{fig:map_PCA_allPeriods_fixed}
\end{figure*}

\begin{figure*}[htbp]%
\centering
\includegraphics[width=0.7\textwidth]{_chapter_02_revisions/ecoregions_states_cerrado_labels.png}
\caption{Map showing the limits of the Cerrado (orange), the limits of the different Cerrado ecoregions (green), and the different Brazilian states (grey) that intersect the Cerrado. Ecoregions: 1, Chapada dos Parecis; 2, Depressão Cuiabana; 3, Paraná Guimarães; 4, Complexo Bodoquena; 5, Bananal; 6, Araguaia Tocantins; 7, Planalto Central; 8, Basaltos do Paraná; 9, Alto São Francisco; 10, Vão do Paraná; 11, Paracatu; 12, Chapadão do São Francisco; 13, Bico do Papagaio; 14, Alto Parnaíba; 15, Floresta de Cocais; 16, Centro-norte Piauiense; 17, Costeiro; 18, Parnaguá; 19, Jequitinhona; 20, Depressão Cárstica do São Francisco. PA, Pará; MT, Mato Grosso; MS, Mato Grosso do Sul; GO, Goiás; DF, Distrito Federal; PR, Paraná; SP, São Paulo; MG, Minas Gerais; MA, Maranhão; TO, Tocantins; PI, Piauí; BA, Bahía.}
\label{fig:ecoregions}
\end{figure*}


\newpage

\clearpage

\begin{table}[htbp]
\caption[]{Components of each fire variable on the six axes of the PCA\_fire in the first (top) and in the last (bottom) periods.}
\centering
\hspace*{-0.6cm}\begin{tabular}{lcccccc}
                            & \multicolumn{6}{c}{\textbf{Period 1985-1993}}                            \\
                            \cmidrule(rl){2-7}
\textbf{Fire variables}     & PC1\_fire      & PC2\_fire      & PC3\_fire      & PC4\_fire      & PC5\_fire      & PC6\_fire      \\
\midrule
Peak fire   season          & -2.07 & 8.11  & 7.46  & -0.33  & 0.51  & -0.03 \\
Median   fire size          & 4.15  & -4.70 & 5.89  & -7.12  & 0.76  & -0.89 \\
99th   percentile fire size & 5.03  & -1.77 & 3.92  & 7.34   & -4.12 & -3.40 \\
Log.   burned area          & 5.51  & 2.12  & -0.57 & -0.08  & -2.45 & 9.21  \\
Log.   number fires         & 4.48  & 5.18  & -4.33 & -4.16  & -3.96 & -5.23 \\
Log. 99th   perc. Frequency & 5.37  & 1.92  & -1.15 & 2.06   & 9.30  & -1.23 \\
\midrule
                            & \multicolumn{6}{c}{\textbf{Period 2012-2020}}                            \\
                            \cmidrule(rl){2-7}
\textbf{Fire variables}     & PC1\_fire      & PC2\_fire      & PC3\_fire      & PC4\_fire      & PC5\_fire      & PC6\_fire      \\
\midrule
Peak fire   season          & -1.92 & 6.71  & -8.71 & -1.17  & 0.38  & 0.01  \\
Median   fire size          & -2.58 & -8.03 & -4.85 & -5.52  & 0.58  & 0.82  \\
99th   percentile fire size & -5.07 & -2.60 & -2.13 & 8.09   & -3.69 & 3.18  \\
Log.   burned area          & -5.81 & 0.29  & 1.48  & -0.63  & -2.36 & -9.17 \\
Log.   number fires         & -5.03 & 3.06  & 3.98  & -5.24  & -4.35 & 5.41  \\
Log. 99th   perc. Frequency & -5.57 & 0.69  & 2.04  & 0.99   & 9.35  & 1.41  \\
\bottomrule
\end{tabular}
\label{tab:pca_fire}
\end{table}


\begin{table}[htbp]
\centering
\caption{Outputs from the Random Forest algorithm, decrease in the GINI index (top) and classification accuracy of the resulting model (bottom). The algorithm was run for each set of variables - anthropic and climatic - separately and together, and per period. Parenthesis next to classification accuracy express the 95\% confidence interval. Highlighted in bold are the four most important variables for classification according to the GINI index.}
\centerline{\begin{tabular}{lcccc}
                                          & \multicolumn{2}{c}{\textbf{Each set separately}} & \multicolumn{2}{c}{\textbf{All variables together}} \\
                                          \cmidrule(rl){2-3} \cmidrule(rl){4-5}
\textbf{Mean GINI decrease}               & 1985-1993              & 2012-2020             & 1985-1993              & 2012-2020             \\
\midrule
Log. population density                   &  441.6                 &  247.7               &  \textbf{237.5}        &   \textbf{117.1}       \\
Log. livestock density                    &  338.2                 &  217.1               &  \textbf{225.6}        &   109.7       \\
Sqrt.   anthropic area perc.              &  403.3                 &  388.2               &  \textbf{215.5}        &   \textbf{170.5}       \\
Log.   change in anthropic area perc.     &  205.5                 &  265.2               &   89.6        &   104.6       \\
\cmidrule(rl){1-3}
Average annual precipitation              &  353.4                 &  214.9               &  147.5         &  94.3       \\
Average temperature in the dry season     &  427.2                 &  503.9               &  \textbf{189.7}         &  \textbf{273.6}       \\
Average precipitation in the dry season   &  410.7                 &  184.6               &  164.2         &  94.9       \\
Peak of the dry season                    &  197.6                 &  214.1               &   81.5         &  \textbf{123.6}       \\
\midrule
                                          & \multicolumn{2}{c}{\textbf{Each set separately}} & \multicolumn{2}{c}{\textbf{All variables together}} \\
                                          \cmidrule(rl){2-3} \cmidrule(rl){4-5}
\textbf{Classification accuracy (\%)}               & 1985-1993              & 2012-2020             & 1985-1993              & 2012-2020             \\
\midrule
Anthropic variables                   & 81.7 (78.1, 84.9)          & 82.8 (79.3, 85.9)         &   \multirow{2}{*}{88.8 (85.9, 91.2)} &  \multirow{2}{*}{94.0 (91.8, 95.8)}       \\
Climatic variables                    & 82.3 (78.7, 85.4)          & 90.0 (87.1, 92.4)         &                                     \\

\bottomrule
\end{tabular}}
\label{tab:RDA_climatic}
\end{table}



\begin{table}[htbp]
\centering
\caption[]{Percentage area of each fire regime region in the last period (2012-2020) that experienced increasing ($>0$), decreasing ($<0$), or no clear trends ($==$) in time for each one of the fire variables and each of the PC1\_fire with 80\% confidence level. The trends in burned area, number of fires, median fire size, 99th percentile fire size, and median fire size were calculated by fitting a linear regression on the annual data for each variable in the period 1985-2020. The trends in the 99th percentile of fire frequency, PC1\_fire and PC2\_fire were calculated in the same manner but on a 9-year interval time series, since the frequency was estimated as the number of times a pixel burned in each of the four time intervals (see Methods).}
\begin{tabular}{rccccccccc}
%\toprule
\multicolumn{1}{c}{}  & \multicolumn{3}{c}{\textbf{LS\_P4}}             & \multicolumn{3}{c}{\textbf{HFA\_P4}}             & \multicolumn{3}{c}{\textbf{LFA\_P4}}             \\
\cmidrule(rl){2-4} \cmidrule(rl){5-7} \cmidrule(rl){8-10}
\textbf{Fire variables}  & \textless 0 & ==   & \textgreater 0 & \textless 0 & ==   & \textgreater 0 & \textless 0 & ==   & \textgreater 0 \\
\midrule
Log. 99th frequency   & 2.2         & 69.2 & 28.6           & 14          & 74.4 & 11.6           & 27.5        & 69   & 3.5            \\
Log. burned area      & 8.2         & 44   & 47.8           & 30.8        & 53.3 & 15.9           & 53.4        & 36.4 & 10.2           \\
Log. number fires     & 15.9        & 53.3 & 30.8           & 62.7        & 32   & 5.3            & 63.8        & 27.7 & 8.6            \\
Log. median fire size & 3.3         & 57.7 & 39             & 7.5         & 66.9 & 25.6           & 11.4        & 70.7 & 17.9           \\
Log. 99th fire size   & 3.3         & 37.9 & 58.8           & 9.3         & 55.5 & 35.2           & 36.9        & 50.3 & 12.8           \\
Median fire season    & 7.7         & 88.5 & 3.8            & 4.8         & 89.4 & 5.8            & 4.5         & 91.8 & 3.7            \\
PC1\_fire             & 2.2         & 68.7 & 29.1           & 15.9        & 69.8 & 14.3           & 27.1        & 67.7 & 5.2             \\
PC2\_fire             & 7.7         & 85.7 & 6.6            & 5.7         & 86.6 & 7.7            & 16.8        & 79.2 & 4.1             \\
\bottomrule
\end{tabular}
\label{tab:perChange_perCluster_p4}
\end{table}

\begin{table}[htbp]
\centering
\caption[]{Components of each anthropic variable on the four axes of the PCA\_anth in the first (top) and in the last (bottom) periods.}
\begin{tabular}{lccccc}
                                      & \multicolumn{4}{c}{\textbf{Period 1985-1993}}                 \\
                                      \cmidrule(rl){2-5}
\textbf{Anthropic variables}          & PC1\_anth      & PC2\_anth      & PC3\_anth      & PC4\_anth     \\
\midrule
Log.   population density             & 5.34 & -4.21 & 6.01  & -4.52  \\
Log.   livestock density              & 6.31 & -0.13 & 0.27  & 7.94   \\
Sqrt.   anthropic area perc.          & 5.78 & 2.39  & -6.72 & -6.72  \\
Log.   change in anthropic area perc. & 1.06 & 8.91  & 4.65  & 4.65   \\
\midrule
                                      & \multicolumn{4}{c}{\textbf{Period 2012-2020}}                 \\
                                      \cmidrule(rl){2-5}
\textbf{Anthropic variables}          & PC1\_anth      & PC2\_anth      & PC3\_anth      & PC4\_anth     \\
\midrule
Log.   population density             & 4.32  & 5.25 & 7.53 & -0.16    \\
Log.   livestock density              & 5.91  & 0.66 & -3.70 & 7.34    \\
Sqrt.   anthropic area perc.          & 5.78  & 1.31 & -4.38 & -6.98   \\
Log.   change in anthropic area perc. & -3.99   & 8.56  & -3.66 & 0.60   \\
\bottomrule
\end{tabular}
\label{tab:pca_anthropic}
\end{table}

\begin{table}[htbp]
\centering
\caption[]{Outputs from the Redundancy Analysis performed on the fire variables and the anthropic variables. Showing the first two components of each fire variable (top) and the first two coefficients of each anthropic variable (bottom) for the first (left) and last (right) periods.}
\hspace*{-0.5cm}\begin{tabular}{lcccc}
                                      & \multicolumn{2}{c}{\textbf{Period 1985-1993}} & \multicolumn{2}{c}{\textbf{Period 2012-2020}} \\
                                      \cmidrule(rl){2-3} \cmidrule(rl){4-5}
\textbf{Fire variable}                & RDA1\_fire              & RDA2\_fire             & RDA1\_fire              & RDA2\_fire             \\
\midrule
Peak fire   season                    & -0.51         & 1.64         & 1.83          & 0.86         \\
99th   percentile fire size           & 2.93          & -0.25        & 2.40          & -1.17         \\
Log.   burned area                    & 2.62          & 0.55         & 2.81          & 0.10          \\
Log.   number fires                   & 1.60          & 1.24         & 2.38          & 1.13         \\
Median   fire size                    & 2.66          & -1.21        & 0.99          & -1.52         \\
Log. 99th   perc. Frequency           & 2.77          & 0.61         & 2.94          & -0.07          \\
\midrule
                                      & \multicolumn{2}{c}{\textbf{Period 1985-1993}} & \multicolumn{2}{c}{\textbf{Period 2012-2020}} \\
                                      \cmidrule(rl){2-3} \cmidrule(rl){4-5}
\textbf{Anthropic variable}           & RDA1\_anth              & RDA2\_anth             & RDA1\_anth              & RDA2\_anth             \\
\midrule
Log.   population density             & -0.86          &  0.33         & -0.41          &  0.90         \\
Log.   livestock density              & -0.84          & -0.03         & -0.77          &  0.14         \\
Sqrt.   anthropic area perc.          & -0.74          & -0.65         & -0.910         & -0.02         \\
Log.   change in anthropic area perc. &  0.26          & -0.48         &  0.66          &  0.13         \\
\bottomrule
\end{tabular}
\label{tab:RDA_anthropic}
\end{table}


\begin{table}[htbp]
\centering
\caption[]{Components of each climatic variable on the four axes of the PCA\_clim in the first (top) and in the last (bottom) periods.}
\hspace*{-0.275cm}\begin{tabular}{lcccc}
                                          & \multicolumn{4}{c}{\textbf{Period 1985-1993}}      \\
                                          \cmidrule(rl){2-5}
\textbf{Climatic variables}               & PC1\_clim      & PC2\_clim     & PC3\_clim      & PC4\_clim      \\
\midrule
Average annual precipitation            & 4.52 & 5.75  & -6.23 & -3.22 \\
Average temperature in the dry season   & 5.80 & -3.74 & 3.99  & -6.27 \\
Average precipitation in the dry season & -5.88  & -3.50 & -6.11 & -6.11 \\
Peak of the dry season                  & 3.77 & -6.60 & 3.99  & 3.99 \\
\midrule
                                          & \multicolumn{4}{c}{\textbf{Period 2012-2020}}      \\
                                          \cmidrule(rl){2-5}
\textbf{Climatic variables}               & PC1\_clim      & PC2\_clim      & PC3\_clim      & PC4\_clim      \\
\midrule
Average annual precipitation             & 4.56 & 6.28 & -0.63  & -6.51 \\
Average temperature in the dry season    & 5.43 & -3.39 & 7.86  & -0.22  \\
Average precipitation in the dry season  & -6.41  & -2.06 & 3.35  & -6.81 \\
Peak of the dry season                   & 3.39 & -6.91 & -5.43 & -3.77 \\
\bottomrule
\end{tabular}
\label{tab:pca_climatic}
\end{table}

\begin{table}[htbp]
\centering
\caption[]{Outputs from the Redundancy Analysis performed on the fire variables and the climatic variables. Showing the first two components of each fire variable (top) and the first two coefficients of each climatic variable (bottom) for the first (left) and last (right) periods.}
\hspace*{-0.5cm}\begin{tabular}{lcccc}
                                          & \multicolumn{2}{c}{\textbf{Period 1985-1993}} & \multicolumn{2}{c}{\textbf{Period 2012-2020}} \\
                                          \cmidrule(rl){2-3} \cmidrule(rl){4-5}
\textbf{Fire variables}                    & RDA1\_fire              & RDA2\_fire             & RDA1\_fire              & RDA2\_fire             \\
\midrule
Peak fire   season                        & 0.51          & -1.84         & 2.24          & -1.18         \\
99th   percentile fire size               & 1.78          & 0.58          & 2.38          & 0.65         \\
Log.   burned area                        & 2.55          & -0.02         & 2.89          & 0.13         \\
Log.   number fires                       & 2.56          & -0.65         & 2.69          & -0.27         \\
Median   fire size                        & 0.92          & 0.99          & 0.75          & 0.67         \\
Log. 99th   perc. Frequency               & 2.76          & 0.25          & 2.95          & 0.33         \\
\midrule
                                          & \multicolumn{2}{c}{\textbf{Period 1985-1993}} & \multicolumn{2}{c}{\textbf{Period 2012-2020}} \\
                                          \cmidrule(rl){2-3} \cmidrule(rl){4-5}
\textbf{Climatic variables}               & RDA1\_clim              & RDA2\_clim             & RDA1\_clim              & RDA2\_clim             \\
\midrule
Average annual precipitation            & 0.37          & 0.62          & 0.05          & 0.76         \\
Average temperature in the dry season   & 0.85          & -0.39         & -0.82           & 0.14          \\
Average precipitation in the dry season & -0.79          & -0.37        & 0.61          & -0.65         \\
Peak of the dry season                  & 0.28          & -0.78         & -0.65           & -0.53         \\
\bottomrule
\end{tabular}
\label{tab:RDA_climatic}
\end{table}



\begin{table}[htbp]
\centering
\caption[]{Pairwise correlation coefficients between all the explanatory variables used in the study. Coefficients were obtained using together the first and last periods' data.}
\hspace*{-1.25cm}\begin{tabular}{lcccccccc}
                                                 & \textbf{(1)} & \textbf{(2)} & \textbf{(3)} & \textbf{(4)} & \textbf{(5)} & \textbf{(6)} & \textbf{(7)} & \textbf{(8)} \\
\midrule
\textbf{Av. annual precipitation (1)}            & 1.00         &              &              &              &              &              &              &              \\
\textbf{Av. temperature in the dry season (2)}   & 0.01         & 1.00         &              &              &              &              &              &              \\
\textbf{Av. precipitation in the dry season (3)} & -0.54        & -0.20        & 1.00         &              &              &              &              &              \\
\textbf{Peak of the dry season (4)}              & -0.12        & 0.32         & -0.06        & 1.00         &              &              &              &              \\
\textbf{Log. population density (5)}             & -0.31        & -0.27        & -0.15        & 0.19         & 1.00         &              &              &              \\
\textbf{Log. livestock density (6)}              & 0.26         & 0.24         & 0.26         & 0.02         & 0.44         & 1.00         &              &              \\
\textbf{Sqrt. anthropic area perc. (7)}          & 0.07         & -0.39        & 0.01         & -0.15        & -0.01        & 0.60         & 1.00         &              \\
\textbf{Log. change in anthropic area perc. (8)} & 0.16         & -0.29        & 0.08         & 0.24         & -0.17        & -0.07        & -0.10        & 1.00         \\
\bottomrule
\end{tabular}
\end{table}


\end{refsection}
