\chapter{Human land occupation regulates the effect of the climate on the burned area of the Brazilian Cerrado (Research Paper 1)}\label{ch:paper01}

\begin{refsection}


\raggedbottom

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

\textbf{Linking statement:} This first chapter was motivated by the need to understand whether and how the opposing trends in burned area observed across the Brazilian Cerrado were explained by the generalised increases in temperature and in land use area in this region over the past four decades. Some studies in the Cerrado had suggested that decreases in burned area could be related to increases in land use area, while climate change was related to increases in burned area. However, no studies thus far had provided a clear explanation for how these two drivers of fire, land use area and the climate, combined to affect the burned area of the Cerrado. Therefore, as a first step to understanding how the dynamics of fire of the Cerrado have been changing over the past few decades, this chapter focused on characterising the interaction between land use area and the climate as drivers of burned area.

This chapter was published on 6$^\text{th}$ July, 2024 in \textit{Communications Earth \& Environment} (\url{https://doi.org/10.1038/s43247-024-01521-5 }). It has been included in this thesis in the same structure in which it was published, with the Methods section after the Discussion.

Results were presented as an oral presentation at the 59$^\text{th}$ Annual Meeting of the Association for Tropical Biology and Conservation in Coimbatore, India, and at the General Assembly of the European Geosciences Union in Vienna, Austria.

\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 responded to all review comments during the peer review process. Dr David Bauman helped in the methodology, analysis and reviewing of the manuscript and of the responses to the reviewer's comments. Vera L. S. Arruda and Dr Ane Alencar helped in the conceptualization, investigation and reviewing of the manuscript. Finally, Dr Imma Oliveras Menor supervised and was involved in all steps of this work.

\newpage

\section{Abstract}

Human activities and climate change are transforming fire regimes globally. The interaction between these two drivers is poorly understood, yet critical if we aim at predicting how biomes will respond to novel fire regimes. In the Brazilian Cerrado, altered fire regimes are threatening its unique biodiversity and ecosystem functioning. Here, using geospatial data for the period 1985-2020 and a causal inference framework to design Bayesian statistical models, we demonstrate that a larger human presence in the landscape ($\geq$40\% land use area) reduces the Cerrado's burned area and hinders its responsiveness to climate; while climatic effects only become apparent in landscapes with little human presence, where hotter and drier conditions increase burned area. Finally, we find spatially heterogeneous burned area trends over time, with increases associated with climate change in landscapes that have remained mostly intact, and decreases caused by anthropic expansion. Both diverging trends have important implications for the conservation of the Cerrado as land use expansion and climate change continue to unfold.

\textbf{\textit{Keywords -}} fire, climate change, climate-human interactions, neotropical savannas

% \newpage

\section{Introduction}

Fire is a complex process that involves the interplay between anthropic activities, the climate, and vegetation structure that both shape and are in turn affected by this phenomenon \citep{archibald_biological_2018}. Understanding the drivers of fire and how they interact is essential to tackling the many challenges which fire poses globally, from human health and the economy \citep{de_mendonca_economic_2004, xu_wildfires_2020} to its influences on the Earth's biogeochemical cycles \citep{archibald_biological_2018, wu_reduced_2022}, and the conservation of both fire-adapted and fire-sensitive ecosystems \citep{pivello_use_2011}. Indeed, fire is a natural disturbance process playing a crucial role in shaping ecosystems \citep{bond_global_2005}, but human activities are altering the Earth's fire regimes \citep{kelly_fire_2020, archibald_managing_2016}. Over the last two decades, global burned area has declined in association with agricultural expansion and intensification, with greater declines concentrated in savannas \citep{andela_human-driven_2017}. In other regions, modern human activities have escalated the occurrence of fire \citep{barlow_clarifying_2020} indicating that anthropic effects vary depending on the context \citep{andela_human-driven_2017}. Simultaneously, climate change is intensifying the risk of fire weather conditions \citep{jolly_climate-induced_2015, abatzoglou_global_2019} and the frequency of extreme weather events \citep{stott_how_2016}, exacerbating fire activity and burned area. Therefore, globally, the two main drivers of fire, humans and the climate, could be pushing burned area in opposite directions, but there is still considerable uncertainty in how they interact \citep{archibald_what_2009, aldersley_global_2011}. Quantifying such interactions is fundamental to improve our understanding of feedbacks between the climate, fires and vegetation \citep{wu_historical_2021, wu_reduced_2022}. %predictive climate-vegetation models and 

Savanna ecosystems stand out as a paradigmatic case of these apparently contrasting processes \citep{archibald_southern_2010}. Savannas are a predominant biome covering around 20\% of the terrestrial surface with extensive presence of humans. These ecosystems play a relevant role on the Earth's carbon cycle contributing to around 25\% of the total gross primary productivity \citep{beer_terrestrial_2010} and have an - often neglected - high biodiversity and ecological value \citep{parr_tropical_2014, veldman_where_2015, veldman_toward_2015}. Savannas are composed of fire-adapted vegetation that depends on this natural disturbance for their continuity and the preservation of their biodiversity \citep{he_fire_2019}. Therefore, departures from the historical fire regimes in either direction can have important ecological consequences potentially resulting in ecosystem shifts to other forms \citep{staver_global_2011, oliveras_many_2016}, as well as effects on the Earth's carbon cycle that can feed back to the climate \citep{wu_reduced_2022}. In these regions, the observed declines in burned area \citep{archibald_southern_2010, silva_drivers_2020} are often linked to landscape fragmentation \citep{archibald_southern_2010, rosan_fragmentation-driven_2022}, population density, certain agricultural practices or increasing grazing pressure that disrupt fuel continuity impairing fire spread. Active fire suppression is also a relevant cause of burned area reduction, known to lead to woody encroachment causing biodiversity losses \citep{stevens_savanna_2017, rosan_extensive_2019} and increasing the risk of more severe wildfire events \citep{fidelis_year_2018}.  Regarding the climate, some studies forecast increases in savannas' burned area associated with warming and the lengthening of the dry season \citep{wu_historical_2021}. However, it has also been suggested that human controls on savannas' burned area may limit its responsiveness to climatic conditions \citep{archibald_climate_2010, conciani_human-climate_2021}, but it is still unclear how much human activities regulate the climate effects on fire.

Here, we focus on the Cerrado, a Brazilian ecoregion comprising several savanna ecosystems and considered to be a biodiversity hotspot \citep{myers_biodiversity_2000}. The Cerrado occupies over 2 million km$^2$, with around 50\% of its area already converted to anthropic land uses \citep{souza_reconstructing_2020}, and undergoing a warming and drying process \citep{hofmann_brazilian_2021} associated with climate change \citep{li_anthropogenic_2021}. These savannas present spatially heterogeneous patterns of burned area \citep{silva_putting_2021} in accordance with its mosaic vegetation and widespread agricultural areas, with both increasing and decreasing burned area trends \citep{silva_drivers_2020, rosan_fragmentation-driven_2022, mataveli_satellite_2018}. Much research has been devoted to studying climate effects on the intra-annual \citep{libonati_spatio-temporal_2015, moreira_de_araujo_distribution_2012, nogueira_spatial_2017} and inter-annual \citep{libonati_spatio-temporal_2015, li_anthropogenic_2021} variability of burned area, identifying increasing burned area trends in relation to climate change \citep{silva_drivers_2020, silva_impacts_2019}. Regarding anthropic area and landscape fragmentation, these seem to have a more complex association with fire activity \citep{silva_drivers_2020, rosan_fragmentation-driven_2022, conciani_human-climate_2021}, while burned area is concentrated in regions with larger natural remnants \citep{moreira_de_araujo_distribution_2012, silva_putting_2021}. Disentangling the effects of climate change and human presence is essential to devise effective fire management policies and conservation plans that are targeted to the characteristics of the landscape and to its vegetation composition.

Using a causal inference framework and a Bayesian approach, we aim to study whether and how the level of human presence in the landscape regulates the effects of the climate on the burned area of the Cerrado. We use the percentage of anthropic area - comprising agricultural, urban and mining land uses - as a proxy for the level of human presence in the landscape to, first, determine its effects on burned natural area - measured as the percentage area of Cerrado native vegetation that burned. Then, we quantify the interaction between the climate and human presence by modelling the effects of the climate on burned area in landscapes with different levels of human presence using Bayesian multilevel regression models. Finally, we explore the spatial distribution of the temporal trends in burned area, climate and anthropic expansion and build a model to determine the causal effects of the last two as predictors of the former.


\section{Results}

\subsection{Limited burned natural area in human-dominated landscapes}


First, we modelled the causal effect of anthropic area percentage on the burned natural area percentage of the Cerrado. With this aim, we combined annual land use and burned area data on a 0.2$^\text{o}$ grid over the Cerrado for the period 1985 to 2020 (see Methods). As we were interested in studying the causal effect of the predictor - anthropic area -, we used a structural causal modelling framework \citep{arif_applying_2023} which, based on assumptions and \textit{a priori} knowledge about our study system, provides a set of rules - known as the backdoor criterion \citep{pearl_causal_2016, arif_applying_2023} - to design models with the sufficient statistical adjustments to limit bias and confounding, thus enabling causal inference.

Therefore, we first designed a Directed Acyclic Graph (DAG), DAG1 (Supplementary Fig. \ref{fig:dag01}), a form of heuristic causal model used to visually represent the causal relations between the different variables of our system and inform covariate selection (see Methods and Supplementary Methods). This DAG was designed based on author's expert knowledge about savanna fire and its drivers, the literature and exploratory data analysis, and reflects our assumptions about how the different variables affect each other (see Methods and Supplementary Methods). Hence, based on DAG1, we built model M1 (Equation \ref{model1}). This model included anthropic area percentage as the explanatory variable of interest for which we wanted to determine its total causal effects on burned natural area, along with protected area and latitude as control variables. These last two covariates were needed according to the backdoor criterion to enable the causal interpretation of anthropic area effects, but it is important to stress that their own modelled effects cannot be interpreted causally (see Methods).

With model M1, we found that the higher the presence of humans in the landscape in terms of area occupied, the smaller the fraction of burned natural area (Supplementary Fig. \ref{fig:anthropic_brms_summary}, slope coefficient: -0.130, [-0.135, -0.126], 95\% Highest Posterior Density Interval, HPDI, \% burned natural area \% anthropic area$^{-1}$), implying that human dominance of a landscape acts to limit the fraction of burned area in Cerrado's native vegetation as hypothesised.

Finally, we calculated the R-squared coefficient adjusted for Bayesian inference \citep{gelman_r-squared_2019} of model M1 (0.544 [0.542, 0.547] 95\% confidence interval, Supplementary Table \ref{tab:ch1_01}). However, it is important to bear in mind that the aim of this study was to determine total causal effects, not to maximise model predictability. Therefore, in this and subsequent models, the R-squared coefficient was calculated just as an indication of the proportion of variance explained by the model, but it should not be interpreted as a measure of the causal relation between the predictors and the response variable.

\subsection{Climate effects on burned natural area are conditional on human landscape dominance}

\begin{figure*}[t!]
\centering
\includegraphics[width =\textwidth]{_chapter_01/figure_01.pdf}
\caption{\textbf{Posterior distributions of the covariate slopes from multilevel regression models M2 and M3.} These models use a focal climatic variable of interest as predictor to explain percentage of burned natural area, while allowing the climate effect (x-axis) to vary by decile of anthropic area percentage (y-axis). Predictor variables: \textbf{a} temperature, \textbf{b} vapour pressure deficit (VPD), \textbf{c} same year dry season precipitation, and \textbf{d} previous-year total precipitation. As calculated, natural area percentage and anthropic area percentage add to 100\%. Green shaded areas indicate the regions of the posterior distributions falling outside the 90\% HPDI, which we use as credibility interval. Red dashed lines are drawn as a reference to indicate whether 0 (i.e. no clear effect) is a likely outcome of the posterior distributions. }
\label{fig:rq01_option04_aa}
\end{figure*}

Then, we assessed the causal effects of climatic controls on the burned natural area of the Cerrado, and the extent to which these controls depend on anthropic factors such as anthropic area. Thus, we quantified the causal effects of temperature, vapour pressure deficit (VPD), same and previous year total precipitation and dry season precipitation (Supplementary Fig. \ref{fig:climate_maps}) on the observed percentage of burned natural area (Supplementary Fig. \ref{fig:ba_maps}) using climate data from different sources (see Methods). 

Using the same DAG1 as before, we built two Bayesian multilevel models (M2 and M3, Equations \ref{model2} and \ref{model3}, see Methods) where we allowed the climate effect coefficients to vary by decile of anthropic area percentage (Supplementary Fig. \ref{fig:anthropic_maps}) to determine the interaction of the climate covariates with the level of human occupation (Fig. \ref{fig:rq01_option04_aa}, see R-squared coefficients in Supplementary Table \ref{tab:ch1_02}). Using the backdoor criterion, we included latitude as a control variable to avoid confounding in model M2 for temperature and VPD (see Methods), while no covariates were needed in model M3 for the different precipitation variables.

We found that the slope coefficients for all climate predictors showed a non-linear change in behaviour with the level of human presence in the landscape (Fig. \ref{fig:rq01_option04_aa}). In cells with a large area of natural vegetation - around 40\% or less of anthropic area -, we observed the expected effect of the climate on burned natural area: increases in temperature and VPD cause increases in burned natural area; the lower the dry season precipitation, the larger the burned area; and higher precipitation in the previous year brought about larger burned area in the following year. We found no strong causal effect of same year total precipitation (Supplementary Fig. \ref{fig:option_08_samePrec}a). In contrast, where anthropic land uses dominate the landscape - around 40 \% or more of the cell's area -, we observed no effect of these climatic variables on the percentage of burned natural area, suggesting strong controls of land use activities.

Those observations in the [90-100]\% decile of anthropic area showed negligible effects of any climatic variables on burned natural area (Fig. \ref{fig:rq01_option04_aa}). We used these results as a reference to test the presence of thresholds from which the effect of the climate in following deciles was different compared to this extreme case. With this aim, for each climate covariate, we generated pairwise contrasts by subtracting the posterior distribution of the slope coefficient obtained for the [90, 100] \% decile from the corresponding posterior of every other decile (Supplementary Fig. \ref{fig:option_07}). We considered that if the 90\% HPDI of the pairwise contrast (the difference between the posteriors) did not encompass 0, then there was an effect of the climate on the burned natural area in that decile compared to the no-effect [90, 100] \% reference. Indeed, we found that below some anthropic area percentage, temperature (50\% anthropic area), VPD (40\%), previous year precipitation (40\%), and dry season precipitation (20\%) had an effect on burned natural area. Hence, for each climatic variable there was a certain threshold value of anthropic area percentage at which burned natural area transitioned from being not responsive to becoming increasingly affected by climate. This transition was sudden, that is, non-linear.

Among these climatic variables, VPD and temperature seemed to have a more marked effect on burned area than dry season and previous year precipitation (Fig. \ref{fig:rq01_option04_aa}), suggesting the former two may emerge as stronger climatic controls of burned natural area. We explored this further by comparing the predicted effect on burned natural area of each climate variable at their average value to their predicted effect at their average value plus (minus in the case of dry season precipitation) 2 standard deviations in each anthropic area decile (Supplementary Fig. \ref{fig:counterfactual_delta}). We observed that, in the less anthropic landscapes (particularly the [0, 10)\% decile, see Supplementary Fig. \ref{fig:counter_delta_0-10AA}), when increasing (decreasing) each climatic variable's value by 2 standard deviations, VPD was the variable with the largest increase in burned natural area, followed by temperature; while the change in effect of dry season and previous year precipitation was more limited.


\subsection{Trends in burned natural area, climate change and anthropic expansion}

\begin{figure*}[t!]
\centering
\includegraphics[width=\textwidth]{_chapter_01/figure_02.pdf}
\caption{\textbf{Spatial distribution of the trends in time (slope coefficients) obtained from model M4 for each cell in a 0.2$^\text{o}$ grid over the Cerrado using the year as the predictor over the period 1985-2020.} Trends for \textbf{a} burned natural area percentage, \textbf{b} anthropic area percentage, and \textbf{c} temperature. Grey indicates cells with no evidence of clear trends - the 95\% HPDI encompasses 0.}
\label{fig:rq02_trendMaps}
\end{figure*}

\begin{figure*}[t!]
\centering
\includegraphics[width=\textwidth]{_chapter_01/figure_03.pdf}
\caption{\textbf{Predicted temporal trends in burned natural area (tBA) from model M5 for those cells that were mostly intact in 1985, with 0 - 20\% of anthropic area (80 - 100\% of natural area).} These predictions are generated by fixing one predictor and exploring the effects of varying the other. \textbf{a} The effects of trends in temperature (tT) on tBA for two fixed values of trends in anthropic area (tA) corresponding to the 10th percentile (0.00 \% year$^{-1}$) and 90th percentile (1.36 \% year$^{-1}$, or 48.96 \% in the period) of the range of tA values observed in this group of cells. \textbf{b} The effects of tA  after fixing tT at the 10th (0.027 $^\text{o}$C year $^{-1}$, or 0.97 $^\text{o}$C in the period) and 90th (0.039 $^\text{o}$C year $^{-1}$, or 1.40 $^\text{o}$C in the period) percentiles of its observed data values for this group of cells. Grey shaded areas correspond to the 80\% HPDI.} 
\label{fig:rq02_g05counterfactual}
\end{figure*}

We explored the spatial distribution of the temporal trends in burned natural area, the climate and anthropic area (Fig. \ref{fig:rq02_trendMaps}) using a Bayesian linear regression model with year as the independent variable (model M4, Equation \ref{model4}). We observed a heterogeneous spatial structure of increases and decreases in burned natural area percentage over the period 1985-2020 (Fig. \ref{fig:rq02_trendMaps}a). We also detected a widespread expansion of human land uses (Fig. \ref{fig:rq02_trendMaps}b) particularly acute in some regions that were mostly natural at the beginning of the study period (Supplementary Fig. \ref{fig:anthropic_maps}a). Regarding the climate, there was a generalised increase in temperature in the Cerrado (Fig. \ref{fig:rq02_trendMaps}c), but we found no evidence of clear trends in precipitation, in agreement with previous studies \citep{hofmann_brazilian_2021}.

Then, we studied whether these trends in temperature (tT) and in anthropic area (tA) had a causal effect on the temporal changes in burned natural area (tBA), and whether these variables mediated each other's effects. With this aim, we built a second DAG (DAG2, Supplementary Fig. \ref{fig:dag02}) representing the relations between these different trends in time. Since the extent of anthropic area expansion was contingent on the level of human occupation at the beginning of the study period, we grouped the cells in five intervals of 20\% anthropic area in 1985. Hence, we modelled the causal effects of tT and tA allowing the slope coefficients to vary by interval using a Bayesian multilevel model (model M5, Equation \ref{model5}) based on DAG2. In addition, we included an interaction term between tT and tA to allow the effect of climate change severity to depend on the level of change in anthropic area. Finally, we explored the implications of model M5 looking at the implied posterior predictions of tT and tA on tBA in the study period (Fig. \ref{fig:rq02_g05counterfactual} and Supplementary Figs. \ref{fig:all_groups_tT_corrected} and \ref{fig:all_groups_tAA_corrected}). That is, we explored the values of tBA predicted by model M5 for different combinations of tA and tT.

We found that, indeed, the effects on tBA exerted by tT and tA varied depending on the pre-existing levels of anthropic area in 1985 (Supplementary Figs. \ref{fig:all_groups_tT_corrected} and \ref{fig:all_groups_tAA_corrected}, see R-squared coefficients in Supplementary Table \ref{tab:ch1_03}). For those cells that were mostly anthropic in 1985 (more than 40\% of anthropic area), we found that both predictors had limited to no effects on tBA. However, for those cells that were mostly natural at the beginning of the period (less than 40\% of anthropic area), model M5 predicted a different effect for each predictor. These differing effects were even more evident in cells with $<$20\% of anthropic area in 1985 (Fig. \ref{fig:rq02_g05counterfactual}). 

In these cells, the effects of the rate of temperature increase seemed to be mediated by the extent of anthropic expansion (Fig. \ref{fig:rq02_g05counterfactual}a). Where human presence in the form of land uses markedly expanded (tA value of 1.36\% per year, or 48.96\% in the period, red line in Fig. \ref{fig:rq02_g05counterfactual}a), burned natural area decreased over time (negative tBA) regardless of tT. In contrast, where anthropic area remained low and roughly at the same level (tA of 0\%, blue line in Fig. \ref{fig:rq02_g05counterfactual}a), we found that warming caused increases in burned natural area over time (positive tBA), and the more acute the temperature increase, the greater the increase in burned natural area over time.

This behaviour became even clearer when we fixed tT at a large (0.039 $^{\text{o}}$C year$^{-1}$, or 1.40 $^{\text{o}}$C in 36 years, red line in Fig. \ref{fig:rq02_g05counterfactual}b) and small (0.027 $^{\text{o}}$C year$^{-1}$, or 0.97 $^{\text{o}}$C in 36 years, blue line in Fig. \ref{fig:rq02_g05counterfactual}b) value and looked at the effect on tBA of varying tA (Fig. \ref{fig:rq02_g05counterfactual}b). In both cases, when tA is small (limited anthropic expansion), burned natural area increases over time (positive tBA), and the larger the increase in temperature (tT), the larger the tBA (larger tBA value in the red line than in the blue line in Fig. \ref{fig:rq02_g05counterfactual}b). However, when tA is large (large anthropic expansion) tBA becomes negative (burned area decreases over time) independently of tT value (similar tBA values in the red and blue predicted values).


\section{Discussion}\label{discussion}

This study demonstrates that anthropic land use activities limit the responsiveness of burned natural area to climatic conditions in areas of the Brazilian savannas where humans have a substantial presence (Fig. \ref{fig:rq01_option04_aa}). Moreover, we find that the interaction between these two drivers is non-linear: humans completely limit burned natural area responses to climate up to a certain level of anthropic area percentage beneath which the expected relation between the climate and burned natural area is recovered. These mitigating effects extend to areas where anthropic land uses are not the major land cover (up to 30-40\% of the area, Supplementary Fig. \ref{fig:option_07}), indicating that human controls on fire stretch to landscapes where considerable natural vegetation remains. In addition, we find that a larger anthropic area is associated with an overall reduced percentage of burned area in natural land cover areas (Supplementary Fig. \ref{fig:anthropic_brms_summary}).

We use anthropic area percentage as an indicator of the presence of human activities and their influence in the landscape, as it encompasses a series of direct anthropic effects on burned natural area. Most anthropic area in the Cerrado consists of agricultural land \citep{souza_reconstructing_2020} conforming fragmented landscapes \citep{rosan_fragmentation-driven_2022}, which limits fire spread through fuel discontinuities \citep{archibald_what_2009}. Another cause of burned area limitation is active fire suppression, a common practice in the Cerrado enforced by policies \citep{durigan_need_2016}, although there have been substantial efforts to change these approaches \citep{durigan_need_2016, durigan_zero-fire_2020, schmidt_fire_2020}. Among other consequences, this reduction in burned natural area contributes to woody encroachment \citep{rosan_extensive_2019, stevens_savanna_2017}, in turn leading to the loss of biodiversity \citep{abreu_biodiversity_2017} and potential ecosystem shifts \citep{oliveras_many_2016}. 

On the other hand, where anthropic area is much smaller, landscape fragmentation is reduced and, consequently, fuel bed continuity may no longer be a limiting factor to fire spread \citep{archibald_what_2009}. Furthermore, access to these areas may be more challenging, reducing fire suppression opportunities. Accordingly, in mostly natural land cover areas fuel load and its flammability may become more relevant determinants of burned area \citep{rego_fire_2021}, both of which are influenced by climate. Indeed, our results demonstrate that in these landscapes climate is a driver of burned natural area (Fig. \ref{fig:rq01_option04_aa}). We find that a higher temperature and VPD causes an increase in burned natural area, as expected, since temperature increases the moisture holding capacity of the air, drying the fuel and making it more flammable \citep{westerling_warming_2006}. In agreement with earlier studies, we observe that lower precipitation in the dry season contributes to larger burned areas as it also influences fuel flammability \citep{flannigan_global_2013, conciani_human-climate_2021}, while higher precipitation in the previous year causes larger burned areas as it affects vegetation productivity and, hence, fuel load \citep{pausas_global_2013}. However, although these fires are occurring in mostly non-anthropic land cover areas, they are not necessarily naturally ignited (e.g. by lightning), as humans are still the main ignition source \citep{pivello_use_2011}. This fact does not conflict with the low presence of humans in these landscapes because most burned area is produced by few very large fires \citep{hantson_status_2016}, so a small number of ignitions can result in large burned areas \citep{archibald_what_2009}, especially under hot and dry conditions. 

Regarding the trends in time, we find a heterogeneous pattern of increases and decreases in burned area (Fig. \ref{fig:rq02_trendMaps}), in line with previous studies \citep{mataveli_satellite_2018, silva_drivers_2020, silva_putting_2021}. When exploring the relation between trends in burned natural area, anthropic expansion and climate change, model predictions reinforce our results presented above. First, changes in burned natural area over time are conditional on the extent of anthropic area at the beginning of the period. Landscapes that were mostly anthropic in 1985 do not present strong trends because burned natural area was already limited and not responsive to climate (Supplementary Figs. \ref{fig:all_groups_tT_corrected} and \ref{fig:all_groups_tAA_corrected}). As such, significant trends in burned natural area are concentrated in cells that were mostly natural in 1985. We find that, in these areas, trends are firstly subject to the extent of anthropic area increases (Fig. \ref{fig:rq02_g05counterfactual}), which brings about a reduction in burned  natural area - adding to 31.1\% of the cells -, probably in relation to increasing fragmentation and human controls on fire. In contrast, those cells that have remained more intact are the ones where burned natural area increased over time, totalling 9.7\% of the cells. Our model indicates that these increases are associated with the warming and drying of the Cerrado due to climate change. 

Indeed, both VPD and temperature, which we identify as the strongest climatic drivers of burned natural area change out of those considered (Supplementary Figs. \ref{fig:counterfactual_delta} and \ref{fig:counter_delta_0-10AA}), have increased over the Cerrado in the past few decades and are expected to continue rising \citep{hofmann_brazilian_2021}. Such increases have been particularly acute towards the end of the dry season \citep{hofmann_brazilian_2021} when climatic conditions are most favourable for fire spread and when most human-started fires occur \citep{silva_putting_2021}. Our findings indicate that as the warming and drying conditions continue lengthening and exacerbating under climate change, we can expect these increasing trends in the burned areas of the more intact regions of the Cerrado to continue. Therefore, we identify those areas with larger natural vegetation remnants as being more vulnerable to climate change. Among these, the region of the Parque Estadual do Araguaia in the centre-west of the Cerrado stands out for its considerable increase in burned natural area over the period. This area is formed by seasonally flooded savannas particularly vulnerable to droughts \citep{de_souza_neves_dinamica_2018}, suggesting that the type of vegetation may be important in the evolution of fire activity. Although savannas are fire-adapted ecosystems, changes in fire regimes towards larger, more frequent, and intense events under climate change may have degradation effects on these ecosystems, which could hamper their recovery following disturbances in the future \citep{pivello_understanding_2021}. In this sense, practices such as integrated fire management are proving to be effective strategies to prevent large and severe fires while contributing to preserve the Cerrado's vegetation mosaic \citep{fidelis_year_2018, santos_prescribed_2021, franke_prescribed_2024}.

In conclusion, climate change and human land use activities in the Cerrado are altering the fire regimes of the region, both with important ecological consequences. Importantly, we show that human controls over the Cerrado's fire dynamics are substantial and act to limit burned area and to regulate the effect of climate and climate change, contributing to the growing recognition of the role of humans on shaping fire patterns \citep{andela_human-driven_2017, bowman_fire_2009, archibald_what_2009, conciani_human-climate_2021} and emphasising the need to improve their representation in projections of fire activity and in climate-vegetation models \citep{wu_historical_2021, wu_reduced_2022}. We find that the greatest changes in burned natural areas over time occurred in the most intact areas, which experienced increases in burned area due to climate change, and in areas with more severe agricultural expansion, which brought about decreases in the proportion of burned natural area. Both directions of burned area change can have important deleterious effects on Cerrado vegetation. Future research could build upon these results by exploring whether these responses to human activities and the climate differ between vegetation formations or land use type - mainly croplands or pastures \citep{rosan_fragmentation-driven_2022}. The widespread changes of fire activity in either direction, in combination with other consequences of large-scale agricultural expansion and climate change, highlight the multiple pressures the Brazilian savannas are subject to. Considering the importance of the Cerrado in providing local livelihoods, its major role providing and distributing regional water resources \citep{oliveira_trends_2014}, and its rich biodiversity and high levels of endemism \citep{myers_biodiversity_2000}, our findings stress the urgent need to extend conservation practices beyond protected areas and to deploy appropriate fire management policies \citep{durigan_need_2016, durigan_zero-fire_2020}. 


\section{Methods}

\subsection{Study area}
In this study we focus on the Cerrado, a region of over 2 million km$^2$ comprising the Brazilian savannas. These savannas are very biodiverse \citep{myers_biodiversity_2000} and spatially heterogeneous, with vegetation ranging from grasslands to tracts of gallery forests and other woody forms \citep{ribeiro_cerrado_2008}. The climate of the Cerrado is characterised by a strong seasonality, with a wet and a dry season, when most burned area occurs \citep{silva_putting_2021}. The average temperatures range between 18$^\text{o}$C and 28$^\text{o}$C increasing about a latitudinal gradient towards the Equator (Supplementary Fig. \ref{fig:climate_maps}a). Precipitation varies from about 800 mm to 2000 mm of total annual precipitation, roughly increasing from east to west towards the Amazon (Supplementary Fig. \ref{fig:climate_maps}b).

Over the last four decades, the Cerrado has experienced a rapid agricultural expansion affecting around 45\% of its area \citep{alencar_mapping_2020}. In this study, we focus on the 36 year period from 1985 to 2020. It is important to note that already in 1985 large areas in the south of the Cerrado had been converted to pastures and croplands (Supplementary Fig. \ref{fig:anthropic_maps}a).

As we use data from different sources with various spatial resolutions, we divide the Cerrado into a 0.2$^{\text{o}}$ grid (about 22.2 km$^2$) and summarise the different variables at the cell level (see details below). In the analysis, we only include cells with more than 90\% of their area within the Cerrado.


\subsection{Land use data}
We use annual maps of land use and land cover from MapBiomas (Collection 6.0, \citep{souza_reconstructing_2020}) for the period 1985 to 2020. This open data and open source project used a random forest algorithm to classify pixels from annual mosaics of Landsat images at 30 m spatial resolution into several categories. This dataset was developed using Google Earth Engine \citep{gorelick_google_2017} and is available through the MapBiomas toolkit on the same platform.

In this study, we first simplify the MapBiomas categorisation to either natural land cover pixels (herein natural pixels), anthropic land use pixels (anthropic pixels), or other pixels based on the legend codes of Collection 6.0 \citep{souza_reconstructing_2020}. Under natural land cover, we group all subcategories under "Forest", those of "Non Forest Natural Formation", and "Beach, Dune and Sand Spot". As anthropic land uses, we group all subcategories under "Farming", along with "Urban Area" and "Mining". The remaining categories - "Other non Vegetated Area", "Water" and "Non Observed" - are grouped as "other". 

For each cell and year, we calculate the percentage of anthropic area as
\begin{equation}
    \text{\% anthropic area} = 100 \cdot \frac{\text{N anthropic pixels}}{\text{N anthropic pixels} + \text{N natural pixels}},
\end{equation}
where N stands for number. Therefore, $\text{\% anthropic area} = 100 - \text{\% natural area}$, as we are leaving out of the calculation the "other" categories because they are not susceptible to burning.

\subsection{Burned area data}
We obtain the burned area data from MapBiomas Fogo (Collection 1.0, \citep{alencar_long-term_2022}) for the study period, 1985-2020. This dataset consists of maps where pixels are identified as either burned or not burned, along with the month in which the pixel was mapped as such. These maps were created from mosaics of Landsat images at 30 m spatial resolution with a 16-day interval.

For each cell and year, we use the number of natural pixels classified as burned to calculate our main variable of interest, the percentage of natural area that burned (herein burned natural area),
\begin{equation}
    \text{\% burned natural area} = 100 \cdot \frac{\text{N burned natural pixels}}{\text{N natural pixels}}.
\end{equation}

\subsection{Climate reanalysis data}
We use various climatic variables from different data products because of their relation to burned area. We use the 2 metre temperature parameter from ERA5-Land Monthly Averaged climate reanalysis data at 9 km spatial resolution \citep{munoz-sabater_era5-land_2021}. We use this variable along with the 2 metre dewpoint temperature from the same dataset to calculate vapour pressure deficit (VPD), a direct measure of the atmospheric demand for water. We calculate VPD as the difference between saturation pressure (e$_\text{s}$) and actual vapour pressure (e$_\text{a}$) following the procedure detailed in \citep{allen_crop_1998}. Finally, we calculate the annual average temperature and VPD for each cell and year.

We use precipitation data from the Climate Hazards Group InfraRed Precipitation with Station data (CHIRPS, version 2.0), a quasi-global rainfall dataset with 0.05$^\text{o}$ spatial resolution covering more than 35 years \citep{funk_climate_2015}. This data was acquired in the Pentad version (5-day precipitation accumulation in mm). To obtain the total precipitation per year and cell, we first calculate the monthly total precipitation per pixel, then we calculate the monthly precipitation per cell taking the spatial average, finally adding up the monthly averages to the annual precipitation.

For each calendar year, we calculate the total dry season precipitation as the sum of the monthly precipitation of the dry season months (usually around May to September). We identify the dry season months for each year as those in which precipitation is smaller than the potential evapotranspiration (PET). We obtain the PET from TerraClimate \citep{abatzoglou_terraclimate_2018}, a dataset of monthly climate data at the global scale and 0.5$^\text{o}$ resolution, and calculate the cell's spatial average PET.

All climate data from the different products were downloaded through Google Earth Engine \citep{gorelick_google_2017}.

\subsection{Protected area data}

In the first model, we include the percentage of protected area in a cell and year as a covariate. We obtained the data on the conservation units and indigenous lands of Brazil from the Ministério de Meio Ambiente do Brasil (MMA) \citep{MMA_protected-areas_2019}. We calculated the percentage area of each cell and year occupied by those protected area established on or before the year.

\subsection{Data analysis}

To tackle the different research questions, we use various Bayesian linear regression models and Bayesian multilevel models (M1 through to M5), whose structure and variable inclusions were guided to allow causal inference. With this aim, we first draw two Directed Acyclic Graphs (DAGs) describing our system (DAG1 and DAG2, Supplementary Figs. \ref{fig:dag01} and \ref{fig:dag02}). Based on these DAGs, and depending on the causal effects of interest, we applied the backdoor criterion - a test that allows to determine whether and how we can compute a causal effect - to inform our modelling strategy and the selection of covariates to explore causal implications \citep{mcelreath_statistical_2018, pearl_causal_2016, arif_applying_2023, addicott_toward_2022, laubach_biologists_2021}. 

\subsubsection{Creation of a causal inference framework}
Directed Acyclic Graphs (DAGs) are graphical representations of the causal relations between the different variables of interest, along with other variables that may bias or alter the causal association between a given explanatory variable of interest and the response variables \citep{pearl_causal_2016, arif_applying_2023}. We build DAG1 (Supplementary Fig. \ref{fig:dag01}) to represent the relations between the variables that may affect the percentage of burned natural area, and DAG2 (Supplementary Fig. \ref{fig:dag02}) to depict the relations between the different variables that may affect trends in burned natural area with time. The design of both DAGs is based on the authors' expertise, on supporting literature, and exploratory data analysis of the variables that may play a role in the system and their relations. 

Fire, and burned area in specific, is a complex phenomenon that involves many different drivers and actors interacting in an intricate manner, and while fully describing complex multifactorial systems such as the one studied here can be challenging, we explicitly acknowledged and accounted for what we consider to be all the major nodes and causal links relevant to the questions of the present work (details in Supplementary Methods). These DAGs therefore provided a rational basis upon which we conditioned causal inference using statistical models, while making the ecological assumptions underlying our models and causal inferences transparent. DAGs entail testable implications of relations between variables, like conditional independencies \citep{mcelreath_statistical_2018, addicott_toward_2022}, that we have checked based on the associations and conditional associations in the data to validate the DAG structure, wherever possible. See Supplementary Methods for an in-depth explanation of the DAGs designed.


The main interest of this study is to explore causal effects of different explanatory variables of interest to help understand the specific mechanisms by which human land occupation and the climate affect burned area, and how they interact. Thus, each model we build (see following sections) is aimed at quantifying the total causal effects of a variable of interest, notwithstanding the predictive power of the model. For each explanatory variable of interest (anthropic area percentage, climatic variables, or trends in these variables), the corresponding DAG shows the direct and indirect causal paths between this causal variable and the response variable (burned natural area percentage, or its trends), and helps identify the covariates that we need to control for to avoid confounding: that is, those variables that are necessary to close the non-causal paths between the explanatory and response variables that could introduce bias. These control variables (protected area and latitude in model M1, and latitude in model M2, see following sections) are identified using the backdoor criterion. Furthermore, this criterion also indicates those variables that must not be included to avoid opening non-causal paths. As the non-causal paths may be different for each variable of interest, the backdoor criterion may indicate different choices of covariates for each causal question and, thus, imply building different models. It is important to emphasise that, as a consequence of building our models in this manner, only the effects of the variables of interest can be interpreted causally, while the model outcomes for the control covariates have no causal interpretation as they could be confounded - and thus their specific quantification is not relevant to our research questions.

Finally, we would like to remark that, as a consequence of the explained above, our models include only those variables needed to address our research questions about causal relations. Therefore, there may be of course other variables that also affect burned area causally, or that are good predictors of this variable, but that we have not included in our models because they are neither the variables of interest to the research questions nor needed for control. In addition, because the models are not designed to maximise their predictability, their overall explanatory power may be limited. Instead, the results from this work should be assessed in the light of the uncertainty of the model parameters, which, as we are using a Bayesian approach, stems from the posterior probability densities - which we summarise using the High Posterior Density Interval, the shortest interval in a posterior probability representing a certain confidence interval.

\subsubsection{Effect of anthropic area on burned natural area}
Model M1 is a Bayesian linear regression model with anthropic area percentage as the predictive variable of interest and the percentage of burned natural area as the response variable. We run this model for all data points, consisting of each grid cell per year. We construct the model based on DAG1 (Supplementary Fig. \ref{fig:dag01}). In our system, the causal effect of anthropic area percentage on natural area burned may be confounded by the extent of protected area and by the latitude (see Supplementary Methods for further explanation), since both are common causes of the explanatory and the outcome variables (backdoor paths). Therefore, we include these two variables as controls to block the non-causal paths, so that the slope coefficient of the anthropic area covariate will be the total causal effect of this variable on burned natural area (conditional on our DAG) \citep{mcelreath_statistical_2018, laubach_biologists_2021}. We therefore model burned natural area as
\begin{align}
    & \text{BA}_i \sim \text{N}\left( \mu_i, \sigma\right) \nonumber \\
    & \mu_i = \alpha + \beta \text{A}_i + \eta \text{L}_i + \psi \text{PA}_i + \gamma_{\text{cell}} + \theta_t, \label{model1}
\end{align}
where BA represents percentage of burned natural area; $\alpha$ is the expected average BA when all covariates are zero (intercept); $\beta$, $\eta$, and $\psi$ represent the burned natural area response to anthropic area percentage (A), latitude (L) and percentage of protected area (PA). It is worth stressing here that in this model only $\beta$ is to be interpreted causally, while $\eta$, and $\psi$ are only there to block non-causal paths (see above). We use parameters $\gamma_\text{cell}$ and $\theta_t$ as varying intercepts to capture the residual variation in expected burned natural area among grid cells and years respectively. Full model specification is in Supplementary Methods.

\subsubsection{Climate effects on burned natural area and mediation by anthropic area percentage}

In models M2 and M3, we use data for the cell per year observations to test the effect of the climate on the percentage of burned natural area for different levels of anthropic area percentage. However, burned area is a stochastic process that usually follows a distribution within the family of heavy-tail distributions \citep{hantson_status_2016}: small burned areas are considerably more frequent than extreme events (larger burned areas), which occur only rarely. In our data, we observe this structure for each interval of each climatic variable (Supplementary Fig. \ref{fig:rq01_boxplots}). We notice that it is the likelihood of the larger burned areas - the tail of the power law - that changes with the climate (Supplementary Fig. \ref{fig:rq01_boxplots}). Therefore, we build these two models using a subset of the data comprising the larger events per climatic interval. That is, we divide the values of a certain climatic variable into regular intervals, and for each interval we subset those data points corresponding to the upper 95th percentile of burned natural area.

Additionally, we observe that the cell-year observations experiencing the highest temperatures in the study period (1985-2020) do not follow the same behaviour than the rest of the data points (Supplementary Fig. \ref{fig:rq01_hotcells}a): all the burned natural area events for these temperatures are far smaller than the events with temperature immediately lower. When mapping the number of years each cell has experienced an average temperature equal or larger than 28.5 $^\text{o}$C, we find that these events occurred in a small number of cells, which are clustered together and located at the border with the Caatinga biome (Supplementary Fig. \ref{fig:rq01_hotcells}b). This difference in behaviour could hence be related to this area's vegetation being more similar to that of the Caatinga (xeric shrubland and thorn forest), with distinct fuel types and loads. We briefly explored the climate of these cells and found that they also present some of the more extreme maximum climatological water deficit (MCWD) values in the Cerrado (Supplementary Fig. \ref{fig:rq01_hotcells}c), but we do not investigate this matter further as it is beyond the scope of the study. Since these cells are clustered together and jointly present a markedly different behaviour, we excluded them from this part of the analysis.

Based on DAG1 (Supplementary Fig. \ref{fig:dag01}), we use two different models to quantify the causal relations of the climate on the burned natural area. When studying the effect of temperature and VPD, our model must include latitude as a predictor (model M2). This is because latitude is a parent cause to temperature (and therefore to VPD) and to burned natural area, hence expected to confound the causal effect of temperature on burned natural area if not controlled for. To explore the causal effects of the different precipitation-related variables (same year total precipitation, same year dry season precipitation, and previous year total precipitation), there is no need to control for any other variable (model M3), as there are no alternative paths going through a common cause between the predictor and the response variable (Supplementary Fig. \ref{fig:dag01}).

Models M2 and M3 have a two-level hierarchical structure where we use a linear regression to predict percentage of burned natural area (BA) as a function of a climatic variable (C) while allowing both the slope and intercept coefficients to vary by decile of the anthropic area percentage factor, with 10 levels. That is,
\begin{align}
    & \text{BA}_i \sim \text{N} \left(\mu_i, \sigma \right) \nonumber \\
    & \mu_i = \alpha_{\text{[A]}} + \beta_{\text{[A]}} \text{C}_i + \gamma_{\text{cell}} + \theta_t + \eta \text{L}_i, \label{model2}
\end{align}
and
\begin{align}
    & \text{BA}_i \sim \text{N} \left(\mu_i, \sigma \right) \nonumber \\
    & \mu_i = \alpha_{\text{[A]}} + \beta_{\text{[A]}} \text{C}_i + \gamma_{\text{cell}} + \theta_t. \label{model3}
\end{align}
\indent In the models, [A] indicates that the corresponding coefficient is allowed to vary by decile of anthropic area, $\alpha_{\text{[A]}}$ are the expected average burned natural area for each decile (intercepts), and $\beta_{\text{[A]}}$ correspond to the effect of each climatic variable on the percentage of burned natural area in each decile. Parameters $\gamma_\text{cell}$ and $\theta_t$ are varying intercepts used to capture the residual variation in expected burned natural area among cells and years respectively. Lastly, $\eta$ represents the relation between burned natural area and latitude (L). Again, only the $\beta_\text{[A]}$ coefficients are to be interpreted causally. The full model specifications including group-level covariance matrices, prior and hyperprior specifications can be found in the Supplementary Methods. Model M2 was run separately for T and VPD, and model M3 was run separately for same year total precipitation, same year dry season precipitation, and previous year total precipitation.

Finally, we compare the effects of the different climatic variables to explore whether any of them may have a stronger effect on burned natural area than the others. With this aim, we explore the change in burned natural area when increasing (decreasing in the case of dry season precipitation) the value of each climatic variable by 2 standard deviations from their average observed value for each decile of anthropic area. That is, we use the joint posterior of the fitted models to generate burned natural area distributions under a "mean value" and under a "mean value + (or -) 2 standard deviations" of the climatic variables. This results in two predictive posteriors for each anthropic area decile group. We then calculate the "contrasts", that is, we subtract the two group-level posteriors for each variable and decile. In this manner, we can compare the change in burned natural area when varying all climatic variables by the same amount (within each variable's distribution).

\subsubsection{Trends in burned natural area, climate change and anthropic expansion}
For each cell in the 0.2$^\text{o}$ grid over the Cerrado, we estimated the annual rate of change (temporal trends) over the period 1985-2020 of burned natural area percentage, anthropic area percentage and various climatic variables. We calculated the trends using Bayesian linear regression,
\begin{align}
    & \text{V}_i \sim \text{N} \left(\mu_i, \sigma \right) \nonumber \\
    & \mu_i = \alpha + \beta \text{t}, \label{model4}
\end{align}
where V generically represents the  response variables of interest for a cell $i$ in a certain year, t represents time measured in years, $\alpha$ characterises the average value V of a cell in 1985, and $\beta$ is the annual rate of change of V - the parameter of interest. We ran this model separately for burned natural area, anthropic area, temperature, VPD, and precipitation using uninformative priors (see Supplementary Methods). We consider that a certain cell does not show clear trends in time if the its 95\% HPDI of the posterior distribution for $\beta$ encompasses 0. Finally, we map the temporal trends obtained for each variable (Fig. \ref{fig:rq02_trendMaps}) over the Cerrado.

We then proceed with a model (M5) aiming to test whether and how the temporal trends in anthropic area (tA) expansion and in temperature (tT) are responsible for the spatial variation in the temporal change in burned natural area (tBA). These temporal trends are the $\beta$ parameters obtained in model M4. In this case, we only use those cells that show clear trends in slopes for tA, tT and tBA (95\% HPDI). We built model M5 based on DAG2 (Supplementary Fig. \ref{fig:dag02}), which we designed using existing knowledge of similar systems and testing the associations and conditional associations implied by the DAG among the variables to validate its statistical plausibility. In the model, we include a linear interaction term between tT and tA to allow each predictor to regulate the effect of the other. However, the extent of anthropic expansion in the study period is contingent on the percentage of anthropic area in each cell in the year 1985: cells that were mostly intact at the beginning of the period could experience any degree of expansion in anthropic area, large or small, and they consequently show a wide range of tA values. By contrast, cells that were already largely converted in 1985 could only experience small tA values. Hence, we can expect different effects of the predictor variables depending on the anthropic area extent of the cells at the start of the period. Therefore, we classify cells in five intervals of 20\% of anthropic area percentage in 1985, and allow the intercept and slope coefficients of tT and tA to vary by group by building M5 as a Bayesian multilevel model,
\begin{align}
    & \text{tBA}_i \sim \text{N}\left( \mu_i, \sigma \right) \nonumber \\
    & \mu_i = \alpha_{[\text{A}_{85}]} + \beta_{1,[\text{A}_{85}]} \text{tA}_i + \beta_{2,[\text{A}_{85}]} \text{tT}_i + \beta_{3,[\text{A}_{85}]} \text{tA}_i \text{tT}_i, \label{model5}
\end{align}\\
where $\alpha_{[\text{A}_{85}]}$ characterises the average tBA of cells with a certain level of anthropic area in 1985 [A$_{85}$], and $\beta_{1,[\text{A}_{85}]}$, $\beta_{2,[\text{A}_{85}]}$, and $\beta_{3,[\text{A}_{85}]}$ characterise the response of tBA to tA, tT and their interaction, respectively, for each [A$_{85}$] group. The full model specification is in the Supplementary Methods.

Finally, we explore the implications of model M5 generating predictive posterior draws for each A$_{85}$ group by fixing one of the predictors (at its 90th and 10th percentile for each A$_{85}$ group) and exploring the effects of the other predictor on tBA throughout its observed value range (Fig. \ref{fig:rq02_g05counterfactual} and Supplementary Figs. \ref{fig:all_groups_tT_corrected} and \ref{fig:all_groups_tAA_corrected}).


\section{Data availability}

\sloppy The climate data used in this study are publicly available from ERA5-Land Monthly Averaged at \url{https://www.ecmwf.int/en/forecasts/dataset/ecmwf-reanalysis-v5-land}, Climate Hazards Group InfraRed Precipitation with Station data version 2.0 at \url{https://data.chc.ucsb.edu/products/CHIRPS-2.0/}, and TerraClimate at \url{https://www.climatologylab.org/terraclimate.html}. The land use data are publicly available from MapBiomas collection 6.0 at \url{https://brasil.mapbiomas.org/en/colecoes-mapbiomas/}. The burned area data are publicly available from MapBiomas Fogo Collection 1.0 at \url{https://brasil.mapbiomas.org/en/colecoes-mapbiomas/}. Finally, the shapefiles used for protected area data (conservation units and indigenous lands) are publicly available from the Ministério de Meio Ambiente do Brasil \url{http://www.mma.gov.br/areas-protegidas/cadastro-nacional-de-ucs}.

\section{Code availability}

All code used in the study for data processing, analysis and modelling is publicly available at: \url{https://doi.org/10.5281/zenodo.10908497}.

\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 Michael Collins for feedback on the manuscript, and Dr Fabian Aguirre-Lopez and Dr Blas Kolic for their suggestions on early stages of the work. We would like to thank the reviewers for providing useful feedback to improve the readability and content of the manuscript. CSG was supported by the UK Natural Environment Research Council under grant number NE/S007474/1. DB was supported by the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement no. 895799.

\section{Author contributions}

C. S. G. contributed to the conceptualization, data curation, methodology, formal analysis, investigation, visualization, writing original draft, and reviewing and editing. D. B. contributed to the methodology, formal analysis, investigation, visualization, and reviewing and editing; V. A. and A. A. contributed to the methodology, investigation, and reviewing and editing; I. O. M. contributed to the conceptualization, methodology, formal analysis, investigation, reviewing and editing, and the supervision.



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

\end{refsection}

\clearpage
\newpage

\begin{refsection}

\supplementarysection

\section{Supplementary Information}

\raggedbottom

\subsection{Supplementary Methods}

\subsubsection{Directed Acyclic Graphs}

We used Directed Acyclic Graphs to inform the covariate selection in our models in order to perform causal inference. We used two different DAGs: DAG1 represented those variables that may affect the percentage of burned natural area, and we used it to design models M1, M2 and M3; DAG2 represented those variables that may affect the temporal trends of change in the percentage of natural burned area, and we used it to inform model M5.

In this section, we explain the reasoning behind each of the DAGs, discuss the implications of the main causal relations, as well as the relative importance of the presence and direction of some causal relations in designing the models and making causal inference.

\textbf{DAG1}

DAG1 (Supplementary Fig. \ref{fig:dag01}) represents the various variables that are - or may be - involved in causal relations with burned natural area, either directly or indirectly. The main variables of interest for this study were the percentage of burned natural area as the outcome, and anthropic area percentage, temperature, VPD, dry season precipitation, and same and previous year precipitation as separate predictors. 

The different climatic variables affect the effective availability of fuel and its continuity - its load and flammability - which directly affect burned area. Temperature, VPD and the amount of precipitation in the dry season - when most fires occur \citep{silva_putting_2021} - influence fuel flammability, as they directly affect how dry the fuel is \citep{westerling_warming_2006, flannigan_global_2013}. In turn, the flammability state of the fuel influences the amount of combustible that will readily burn - the effective fuel load \citep{rego_fire_2021}. As well, precipitation in the growing season affects vegetation productivity and, hence, fuel load. Both fuel load and flammability indirectly affect burned area through ignitions, because they affect the success of the fire onset \citep{rego_fire_2021}.

On the other hand, we are using anthropic area percentage as an indicator of the presence and extent of human land use activities in the landscape. These activities bring about landscape fragmentation \citep{rosan_fragmentation-driven_2022}, ignition \citep{pivello_use_2011}, but also active fire suppression \citep{durigan_need_2016, schmidt_fire_2020}, and influence road density. We calculated anthropic area as the corresponding percentage of the sum of natural and anthropic area (see Methods). However, the presence and extent of anthropic area in any cell could be argued to being conditional on the presence of areas demarcated as protected, so there would be a causal link between the extent of protected areas and those of human land uses, as supported by a statistical association in the data (not shown). But protected areas also shape the level of fragmentation and the road density, they would also have their own fire suppression or prevention policies \citep{durigan_need_2016}, and would also probably be associated with a particular ignition regime. The potential existence of these causal relations that involve protected area percentage as a cause of burned natural area both through anthropic area - the variable of interest - and through an alternative \textit{backdoor causal path} \citep{pearl_causal_2016, laubach_biologists_2021}, imply that protected area must be included when studying the effects of anthropic area in model M1, based on the backdoor criterion.

Another aspect to highlight of DAG1 is the inclusion of \textit{latitude} and \textit{history}, which has important implications for the models' design. In addition to the well known causal relation between latitude and temperature due to the angle of solar radiation, there is a path between latitude and anthropic area as a result of the fact that human expansion in the Cerrado primarily developed from the southern regions towards the northern areas \citep{jepson_disappearing_2005}. This link is symbolised through the inclusion of \textit{history} as an unobserved (or latent) variable, and implies that latitude is simultaneously a causal parent (i.e. a cause) of anthropic area and of temperature and VPD. For this reason, models M1 and M2, studying the effects of anthropic area and temperature or VPD respectively, must include latitude as a covariate to allow causal inference of the variables of interest, by blocking the non-causal path going to percentage of burned natural area through latitude (see Supplementary Fig. \ref{fig:dag01}).

Finally, it is worth mentioning that the conditional independencies implied by DAG1 and including unobserved variables (grey circles, Supplementary Fig. \ref{fig:dag01}) could not be assessed from the raw data, given the unobserved nature of such variables. However, regarding the variables of interest and the overall causal structure of the DAG that concerns them, nor the presence/absence of most of the arrows related to unobserved variables nor their direction would change the models' design nor our ability to draw causal inference - conditional on the rest of the DAG being correct, of course. This is because the potential backdoor paths that these arrows create between these unobserved variables include causal structures called inverted forks or colliders, that is, non-causal paths that are already closed by default, and would only be opened by conditioning on the collider (see details in \citep{mcelreath_statistical_2018, elwert_endogenous_2014}).


\textbf{DAG2}

DAG2 (Supplementary Fig. \ref{fig:dag02}) represented the relations between the temporal trends in burned natural area percentage (tBA) as the outcome, and temporal trends in anthropic area percentage (tA) and in temperature (tT) as predictors. We included the percentage of anthropic area in 1985 (A$_{85}$) because it constrains the extent of further land use expansion. In this case, we validated the DAG by using the data to test the conditional independencies implied by the DAG, i.e. the correlations and conditional correlations among the represented variables. This step also allowed us to discard alternative DAGs and additional arrows between the observed variables \citep{mcelreath_statistical_2018, laubach_biologists_2021}.


\subsubsection{Models Specifications}

\vspace{3pt}

\textbf{Model M1}
\begin{align*}
    &\text{BA}_i \sim \text{N}\left( \mu_i, \sigma\right)& \text{[Likelihood]}&\\
    &\mu_i = \alpha + \beta \text{A}_i + \eta \text{L}_i + \psi \text{PA}_i + \gamma_{\text{cell}} + \theta_t, & \text{[Linear model]}&\\
    \\
    &\alpha, \beta, \eta, \psi \sim \text{N} (0, \sigma_{\alpha, \beta, \eta, \psi}) & \text{[Priors for intercept and slopes]}&\\
    &\gamma_{\text{cell}} \sim \text{N} \left(0, \sigma_{\gamma} \right) & \text{[Adaptive priors for cells]}&\\
    &\theta_t \sim \text{N}\left(0, \sigma_t\right) & \text{[Adaptive priors for the \textit{t} years]}&\\
    \\
    &\sigma_{\alpha}, \sigma_{\beta}, \sigma_{\gamma}, \sigma_{\psi}, \sigma_t, \sigma_{\text{lat}}, \sigma \sim \text{Exp}(1) & \text{[Priors for standard deviation]}&\\
\end{align*}

\noindent \textbf{Markov chain Monte Carlo settings:} Model was run using R package brms \citep{brms} with 2000 total iterations per chain, 1000 warmup iterations, 8 Markov chains, a 0.98 threshold to accept a sampling step of the NUTS (no-U-turn sampler), a seed value of 9.

\textbf{Model M2}
\begin{align*}
    &\text{BA}_{i} \sim \text{N}\left(\mu_{i}, \sigma\right) & \text{[Likelihood]}&\\
    &\mu_i = \alpha_{\text{[A]}} + \beta_{\text{[A]}} \text{C}_i + \gamma_{\text{cell}} + \theta_t + \eta L_i & \text{[Linear Model]}&\\
    \\
    &\left( \begin{array}{c} \alpha_{\text{[A]}} \\ \beta_{\text{[A]}} \end{array} \right) \sim \text{MVNormal}\left(
    \left( \begin{array}{c} \alpha \\ \beta \end{array} \right), \textbf{S}\right) & \text{[Adaptive priors for each A]}&\\
    &\textbf{S} = \left( \begin{array}{cc} \sigma_{\alpha} & 0 \\ 0 & \sigma_{\beta} \end{array}\right) \textbf{R} \left( \begin{array}{cc} \sigma_{\alpha} & 0 \\ 0 & \sigma_{\beta} \end{array}\right) & \text{[Covariance matrix]}&
    \\
    & \textbf{R} = \left( \begin{array}{cc} 1 & \rho_{\alpha, \beta} \\ \rho_{\beta, \alpha} & 1 \end{array} \right) & \text{[Correlation matrix for A]} \\
    &\textbf{R} \sim \text{LKJcorr}(2) & \text{[Prior for correlation matrix]}&\\
    \\
    &\alpha, \beta \sim \text{N} (0, 10) & \text{[Priors for intercept and slope]}&\\
    &\gamma_{\text{cell}} \sim \text{N} \left(0, \sigma_{\gamma} \right) & \text{[Adaptive priors for cells]}&\\
    &\theta_t \sim \text{N}\left(0, \sigma_t\right) & \text{[Adaptive priors for the \textit{t} years]}&\\
    &\eta \sim \text{N}\left(0, \sigma_{\text{lat}}\right) & \text{[Adaptive priors for \textit{lat} latitude]}&\\
    \\
    &\sigma_{\alpha}, \sigma_{\beta}, \sigma_{\gamma}, \sigma_t, \sigma_{\text{lat}}, \sigma \sim \text{Exp}(1) & \text{[Priors for standard deviation]}&\\
\end{align*}

\noindent \textbf{Markov chain Monte Carlo settings:} Model was run using R package brms \citep{brms} with 2000 total iterations per chain, 1000 warmup iterations, 8 Markov chains, a 0.98 threshold to accept a sampling step of the NUTS (no-U-turn sampler), a seed value of 9.

\textbf{Model M3}
\begin{align*}
    &\text{BA}_{i} \sim \text{N}\left(\mu_{i}, \sigma\right) & \text{[Likelihood]}&\\
    &\mu_i = \alpha_{\text{[A]}} + \beta_{\text{[A]}} \text{C}_i + \gamma_{\text{cell}} + \theta_t & \text{[Linear Model]}&\\
    \\
    &\left( \begin{array}{c} \alpha_{\text{[A]}} \\ \beta_{\text{[A]}} \end{array} \right) \sim \text{MVNormal}\left(
    \left( \begin{array}{c} \alpha \\ \beta \end{array} \right), \textbf{S}\right) & \text{[Adaptive priors for each A]}&\\
    &\textbf{S} = \left( \begin{array}{cc} \sigma_{\alpha} & 0 \\ 0 & \sigma_{\beta} \end{array}\right) \textbf{R} \left( \begin{array}{cc} \sigma_{\alpha} & 0 \\ 0 & \sigma_{\beta} \end{array}\right) & \text{[Covariance matrix]}&
    \\
    & \textbf{R} = \left( \begin{array}{cc} 1 & \rho_{\alpha, \beta} \\ \rho_{\beta, \alpha} & 1 \end{array} \right) & \text{[Correlation matrix for A]} \\
    &\textbf{R} \sim \text{LKJcorr}(2) & \text{[Prior for correlation matrix]}&\\
    \\
    &\alpha, \beta \sim \text{N} (0, 10) & \text{[Priors for intercept and slope]}&\\
    &\gamma_{\text{cell}} \sim \text{N} \left(0, \sigma_{\gamma} \right) & \text{[Adaptive priors for cells]}&\\
    &\theta_t \sim \text{N}\left(0, \sigma_t\right) & \text{[Adaptive priors for the \textit{t} years]}&\\
    \\
    &\sigma_{\alpha}, \sigma_{\beta}, \sigma_{\gamma}, \sigma_t, \sigma \sim \text{Exp}(1) & \text{[Priors for standard deviation]}&\\
\end{align*}

\noindent \textbf{Markov chain Monte Carlo settings:} Model was run using R package brms \citep{brms} with 2000 total iterations per chain, 1000 warmup iterations, 8 Markov chains, a 0.98 threshold to accept a sampling step of the NUTS (no-U-turn sampler), a seed value of 9.

\textbf{Model M4}
\begin{align*}
    &\text{V}_i \sim \text{N} \left(\mu_i, \sigma \right) & \text{[Likelihood]}\\
    &\mu_i = \alpha + \beta \text{t} & \text{[Linear Model]}\\
    &\alpha, \beta \sim \text{N} (0, 10) & \text{[Priors for coefficients]}\\
    &\sigma \sim \text{Exp}(1) & \text{[Prior for st. deviation]}\\
\end{align*}

\noindent \textbf{Markov chain Monte Carlo settings:} Model was run using R package brms \citep{brms} with 2000 total iterations per chain, 1000 warmup iterations, 2 Markov chains, a 0.98 threshold to accept a sampling step of the NUTS (no-U-turn sampler), a seed value of 9.

\textbf{Model M5}
\begin{align*}
    &\text{tBA}_{i} \sim \text{N}\left(\mu_{i}, \sigma\right) & \text{[Likelihood]}\\
    &\mu_i = \alpha_{[\text{A}_{85}]} + \beta_{1,[\text{A}_{85}]} \text{tA}_i + \beta_{2,[\text{A}_{85}]} \text{tT}_i + \beta_{3,[\text{A}_{85}]} \text{tA}_i \text{tT}_i &\text{[Linear Model]}\\
    \\
    &\left( \begin{array}{c} \alpha_{[\text{A}_{85}]} \\ \beta_{1,[\text{A}_{85}]} \\ \beta_{2,[\text{A}_{85}]} \\ \beta_{3,[\text{A}_{85}]} \end{array} \right) \sim \text{MVNormal}\left(
    \left( \begin{array}{c} \alpha \\ \beta_1 \\ \beta_2 \\ \beta_3 \end{array} \right), \textbf{S}\right) & [\text{Adaptive priors for A}_{85}]\\
    \\
    &\textbf{S} = \left( 
    \begin{array}{cccc} 
    \sigma_{\alpha} & 0 & 0 & 0 \\ 
    0 & \sigma_{\beta_1} & 0 & 0 \\
    0 & 0 & \sigma_{\beta_2} & 0 \\
    0 & 0 & 0 & \sigma_{\beta_{3}} \end{array}\right) 
    \textbf{R} \left( 
    \begin{array}{cccc} 
    \sigma_{\alpha} & 0 & 0 & 0 \\ 
    0 & \sigma_{\beta_1} & 0 & 0 \\
    0 & 0 & \sigma_{\beta_2} & 0 \\
    0 & 0 & 0 & \sigma_{\beta_{3}} \end{array}\right) & \text{[Covariance matrix]}\\
    \\
    & \textbf{R} = \left( 
    \begin{array}{cccc} 
        1 & \rho_{\alpha, \beta_1} & \rho_{\alpha, \beta_2} & \rho_{\alpha, \beta_3}\\ 
        \rho_{\beta_1, \alpha} & 1 & \rho_{\beta_1, \beta_2} & \rho_{\beta_1, \beta_3} \\
        \rho_{\beta_2, \alpha}& \rho_{\beta_2, \beta_1} & 1  & \rho_{\beta_2, \beta_3} \\
        \rho_{\beta_3, \alpha}& \rho_{\beta_3, \beta_1} & \rho_{\beta_3, \beta_2} & 1  
    \end{array} \right) & \text{[Correlation matrix]}\\
    &\textbf{R} \sim \text{LKJcorr}(2) & \text{[Prior for corr. matrix]}\\
    \\
    &\alpha, \beta_1, \beta_2, \beta_3 \sim \text{N} (0, 10) & \text{[Priors for coefficients]}\\
    &\gamma_{\text{cell}} \sim \text{N} \left(0, \sigma_{\gamma} \right) & \text{[Adaptive priors for cells]}\\
    &\theta_t \sim \text{N}\left(0, \sigma_t\right) & \text{[Adaptive prior for \textit{t}]}\\
    &\eta \sim \text{N}\left(0, \sigma_{\text{lat}}\right)&  \text{[Adaptive priors for \textit{lat}]}\\
    &\sigma_{\alpha}, \sigma_{\beta_1}, \sigma_{\beta_2}, \sigma_{\beta_3}, \sigma_{\gamma}, \sigma_t, \sigma_{\text{lat}}, \sigma \sim \text{Exp}(1) & \text{[Priors for st. deviation]}\\
\end{align*}

\noindent \textbf{Markov chain Monte Carlo settings:} Model was run using R package brms \citep{brms} with 2000 total iterations per chain, 1000 warmup iterations, 8 Markov chains, a 0.98 threshold to accept a sampling step of the NUTS (no-U-turn sampler), a seed value of 9.




\newpage
\subsection{Supplementary Figures}
\begin{figure}[h!]
\centering
\includegraphics[width=\textwidth]{_chapter_01/DAG1.pdf}
\caption[]{\textbf{Directed Acyclic Graph, DAG1, describing the set of variables that influence directly or indirectly the burned natural area percentage and the causal relations between them.} We used this DAG in combination with the backdoor criterion (see \textit{Materials and Methods}) to inform the design of models M1, M2 and M3. Green circles indicate those variables that we have included in the study either directly as covariates or in the calculation of other covariates (like in the case of dew temperature, used to calculate VPD), the red circle indicates the outcome variable, and grey circles indicate unobserved variables. Unobserved variables comprise variables that could be measured but were not included in the study, and variables that could not be directly measured. The variable "history" serves to indicate the link between latitude and anthropic area expansion in the Cerrado, as land use conversion mostly started developing in the southern regions towards the northern areas (see Supplementary Fig. \ref{fig:anthropic_maps}).}
\label{fig:dag01}
\end{figure}

\begin{figure}
\centering
\includegraphics[width=0.75\textwidth]{_chapter_01/anthropic_brms_summary_114.pdf}
\caption[]{\textbf{Posterior distribution of the slope coefficients obtained for the different covariates in model M1}. Effect of anthropic area percentage (top), the explanatory variable of interest and the only covariate whose effects can be interpreted causally. Slope coefficients of control variables latitude (middle) and protected area percentage (bottom), which were included in the model to avoid confounding, and whose effects cannot be interpreted causally (see Methods). Green coloured areas indicate the region of the posterior distributions falling outside the 90\% HPDI, which we use as credibility interval. The red dashed line is drawn as a reference to indicate whether 0 is a likely outcome of the posterior distributions.}\label{fig:anthropic_brms_summary}
\end{figure}

\begin{figure}%
\centering
\includegraphics[width=\textwidth]{_chapter_01/average_climate_maps_fourVars_NCC_latlon.pdf}
\caption[]{\textbf{Maps showing the average values of the different climatic variables in the period 1985-2020}. \textbf{a} Annual temperature, \textbf{b} average total annual precipitation, \textbf{c} average VPD, and \textbf{d} average dry season precipitation. We can see a latitudinal gradient of increasing average temperature, and an east-west gradient of increase in precipitation.}\label{fig:climate_maps}
\end{figure}

\begin{figure}%
\centering
\includegraphics[width=0.5\textwidth]{_chapter_01/average_ba_maps_NCC_latlon.pdf}
\caption[]{\textbf{Map of the logarithm average burned area in the period 1985-2020.}}\label{fig:ba_maps}
\end{figure}

\begin{figure}
\centering
\includegraphics[width=\textwidth]{_chapter_01/anthropic_maps_categorical_NCC_latlon.pdf}
\caption[]{\textbf{Maps of anthropic area percentage at the beginning and end of the study period}. Anthropic area in \textbf{a} 1985 and in \textbf{b} 2020, ranging from 0\% (green) to 100\% (orange).}
\label{fig:anthropic_maps}
\end{figure}

\begin{figure}[h]
\centering
\includegraphics[width=\textwidth]{_chapter_01/option_08_samePrec_NCC.pdf}
\caption[]{\textbf{Posterior distributions of same year precipitation effects on burned natural area.} \textbf{a} Posterior distribution of the covariate slopes from multilevel model M3 using same year precipitation as the variable of interest to explain percentage of burned natural area. We use a multilevel model to allow the climate effect (x-axis) to vary by decile of anthropic area percentage (y-axis). \textbf{b} Pairwise posterior distribution contrasts showing the difference in the effect of same year precipitation between each decile of anthropic area and the case where human land uses occupy most of the landscape (90-100\% of the cell's area) - see Supplementary Fig. \ref{fig:counterfactual_delta} for a more detailed explanation. Green and red coloured areas indicate the region of the posterior distributions falling outside the 90\% HPDI, which we use as credibility interval. Red dashed lines are drawn as a reference to indicate whether 0 is a likely outcome of the posterior distributions.}\label{fig:option_08_samePrec}
\end{figure}


\begin{figure}
\centering
\includegraphics[width=\textwidth]{_chapter_01/option_07_NCC.pdf}
\caption[]{\textbf{Pairwise posterior distribution contrasts (counterfactual plots) showing the difference in the effect of each climatic variable between each decile of anthropic area and the case where human land uses occupy most of the landscape (90-100\% of the cell's area).} We used this decile as a reference because no climatic variables show an effect on burned natural area in this extreme case (Main Fig. 1). Plots show the difference in effects of \textbf{a} temperature, \textbf{b} VPD, \textbf{c} dry season precipitation, and \textbf{d} previous year precipitation. The y-axis labels indicate explicitly the anthropic area deciles used in each pairwise contrast calculation. That is, for each climatic variable and decile (e. g. the (20, 30] \% anthropic area decile), the posterior contrasts (x-axis) were obtained by calculating the difference between the draws from the posterior distribution of the slope coefficients obtained in models M2 or M3 for this decile and the draws obtained for the slope coefficients of the [90, 100]\% decile (and hence it is indicated as "[20, 30) - [90, 100]" on the y-axis). Red coloured areas indicate the region of the posterior distributions falling outside the 90\% HPDI, which we use as credibility interval. Red dashed lines are drawn as a reference to indicate whether 0 is a likely outcome of the posterior distributions.}\label{fig:option_07}
\end{figure}

\begin{figure}[h]
\centering
\includegraphics[width=\textwidth]{_chapter_01/counterfactual_delta_allAA_sd2.pdf}
\caption[]{\textbf{Pairwise posterior distribution contrasts of burned area modelled when each climatic variable is increased (temperature, VPD and previous year precipitation) or decreased (dry season precipitation) by 2 standard deviations from its average value in each anthropic area decile.} Pairwise contrasts for \textbf{a} temperature, \textbf{b} VPD, \textbf{c} previous year precipitation, and \textbf{d} dry season precipitation. To obtain these pairwise contrasts, we used models M2 (for temperature and VPD) and M3 (for previous year and dry season precipitation) to model burned natural area for each climatic variable at its average value in each decile, and subtracted it from the modelled burned area at the climatic variable's average value displaced by 1.5 standard deviations. The aim was to ascertain what climatic variable may have a more prominent effect on burned area when their values change by the same standard deviation.}\label{fig:counterfactual_delta}
\end{figure}


\begin{figure}
\centering
\includegraphics[width=\textwidth]{_chapter_01/counterfactual_delta_0-10AA_sd2.pdf}
\caption[]{\textbf{Pairwise posterior distribution contrasts (counterfactual plots) of burned area modelled when each climatic variable increases (temperature, VPD, and previous year precipitation) or decreases (dry season precipitation) by 2 standard deviations from its average value in landscapes with least anthropic presence ([0-10) \% anthropic area).} To obtain these pairwise contrasts, we used models M2 (for temperature and VPD) and M3 (for previous year and dry season precipitation) to predict burned natural area for each climatic variable at its average value in this anthropic area decile, and subtracted it from the modelled burned area at the climatic variable's average value displaced by 2 standard deviations. The aim was to ascertain what climatic variable may have a more prominent effect on burned area when their values change by the same amount of standard deviations.}\label{fig:counter_delta_0-10AA}
\end{figure}

\begin{figure}
\centering
\includegraphics[width=0.75\textwidth]{_chapter_01/DAG2_noUnobserved.pdf}
\caption[]{\textbf{Directed Acyclic Graph, DAG2, describing the set of variables that may influence trends in burned natural area that were included in this study.} We used this DAG to inform the design of model M4. All of the conditional independencies \citep{laubach_biologists_2021} implied by this and alternative DAGs that we considered (comprising the same variables as in DAG2) were tested and validated. Green circles indicate those variables that we have included in the study either directly as covariates or in the calculation of other covariates and the red circle indicates the outcome variable.}
\label{fig:dag02}
\end{figure}

\begin{figure}
\centering
\includegraphics[width=\textwidth]{_chapter_01/all_groups_tT_brmsTrends_corrected.pdf}
\caption[]{\textbf{Model predictions of trends in burned area from model M4 using trends in temperature (tT) as the predictor after fixing the value of trends in anthropic area (tA) at the 10th percentile (blue) and the 90th percentile (red) of the range of observed tA values in each group of cells.} Each figure shows the model predictions for the different groups of cells classified on groups of 20 \% of their anthropic area percentage in 1985 (mostly anthropic to mostly natural from left to right). Grey shaded area shows the 80\% Highest Posterior Density Interval (HPDI).}
\label{fig:all_groups_tT_corrected}
\end{figure}

\begin{figure}
\centering
\includegraphics[width=\textwidth]{_chapter_01/all_groups_tAA_brmsTrends_corrected.pdf}
\caption[]{\textbf{Model predictions of trends in burned area from model M4 using trends in anthropic area (tA) as the predictor after fixing the value of trends in temperature (tT) at the 10th percentile (blue) and the 90th percentile (red) of the range of observed tT values in each group of cells.} Each figure shows the model predictions for the different groups of cells classified on groups of 20 \% of their anthropic area percentage in 1985 (mostly anthropic to mostly natural from left to right). Grey shaded area shows the 80\% Highest Posterior Density Interval (HPDI).}
\label{fig:all_groups_tAA_corrected}
\end{figure} 



\begin{figure}
\centering
\includegraphics[width=\textwidth]{_chapter_01/boxplot_natpercBA-clim_NCC.png}
\caption[]{\textbf{Distribution of natural area burned occurred in the period 1985-2020 for interval of each climatic variable}. \textbf{a} Temperature, \textbf{b}, VPD, \textbf{c} total annual precipitation, \textbf{d} dry-season precipitation, and \textbf{e} previous year total precipitation. Boxplots represent the distribution of the data points (cell-year pairs) falling in each interval.}
\label{fig:rq01_boxplots}
\end{figure}

\begin{figure}
\centering
\includegraphics[width=\textwidth]{_chapter_01/hot_cells_NCC_latlon.png}
\caption[]{\textbf{Figures presented to show why a small cluster of cells has been excluded from models M2 and M3.} \textbf{a} Burned natural area against annual average temperature for all the data points, we can observe how the cells experiencing the hottest years do not seem to follow the same behaviour as the rest. \textbf{b}, \textbf{c} Cells that experienced temperatures equal to or above 28.5 $^\text{o}$C five times or more in the period 1985-2020, and those that experienced MCWD values more extreme than -800 mm ten times or more respectively.}
\label{fig:rq01_hotcells}
\end{figure}

\clearpage

\newpage

\subsubsection{Supplementary Tables}
\begin{table}[h]
\centering
\caption[]{R-squared coefficient adapted for Bayesian inference \citep{gelman_r-squared_2019} calculated for model M1 (see Main) with its 95\% confidence interval, where anthropic area percentage is the causal predictor of interest.}
\begin{tabular}{c}
\begin{tabular}[c]{@{}c@{}}\textbf{Anthropic area}\\ \textbf{percentage}\\ \textbf{(Model M1)}\end{tabular} \\
\midrule
0.544 {[}0.542, 0.547{]}   \\
\bottomrule
\end{tabular}
\label{tab:ch1_01}
\end{table}

\begin{table}[h]
\centering
\caption[]{R-squared coefficients adapted for Bayesian inference \citep{gelman_r-squared_2019} calculated for models M2 (with temperature and Vapour Pressure Deficit, VPD, as the causal predictors of interest) and M3 (with dry season and previous year precipitation as causal predictors of interest) for each anthropic area decile (see Methods), with their 95\% confidence interval.}
\centerline{\begin{tabular}{lcccc}
\multicolumn{1}{l}{} & \begin{tabular}[c]{@{}c@{}}\textbf{Temperature}\\ \textbf{(Model M2)}\end{tabular}   & \begin{tabular}[c]{@{}c@{}}\textbf{VPD}\\ \textbf{(Model M2)}\end{tabular}   & \begin{tabular}[c]{@{}c@{}c@{}}\textbf{Dry season}\\ \textbf{precipitation} \\ \textbf{(Model M3)}\end{tabular} & \begin{tabular}[c]{@{}c@{}c@{}}\textbf{Previous year}\\ \textbf{precipitation} \\ \textbf{(Model M3)}\end{tabular} \\
% \cmidrule(rl){2-5}
\midrule
\textbf{{[}0, 10) \%    }     & 0.43 {[}0.41, 0.46{]} & 0.54 {[}0.52, 0.56{]} & 0.35 {[}0.32, 0.37{]}    & 0.43 {[}0.40, 0.45{]}       \\
\textbf{{[}10, 20) \%   }     & 0.39 {[}0.36, 0.43{]} & 0.47 {[}0.43, 0.50{]} & 0.29 {[}0.26, 0.33{]}    & 0.43 {[}0.39, 0.47{]}       \\
\textbf{{[}20, 30) \%   }     & 0.31 {[}0.27, 0.36{]} & 0.47 {[}0.43, 0.51{]} & 0.26 {[}0.22, 0.30{]}    & 0.34 {[}0.30, 0.38{]}       \\
\textbf{{[}30, 40) \%   }     & 0.24 {[}0.21, 0.27{]} & 0.34 {[}0.30, 0.38{]} & 0.23 {[}0.20, 0.27{]}    & 0.29 {[}0.26, 0.33{]}       \\
\textbf{{[}40, 50) \%   }     & 0.24 {[}0.21, 0.27{]} & 0.29 {[}0.26, 0.33{]} & 0.27 {[}0.23, 0.30{]}    & 0.27 {[}0.24, 0.30{]}       \\
\textbf{{[}50, 60) \%   }     & 0.29 {[}0.26, 0.33{]} & 0.27 {[}0.24, 0.31{]} & 0.29 {[}0.26, 0.33{]}    & 0.32 {[}0.28, 0.36{]}       \\
\textbf{{[}60, 70) \%   }     & 0.24 {[}0.21, 0.27{]} & 0.28 {[}0.24, 0.32{]} & 0.32 {[}0.28, 0.36{]}    & 0.28 {[}0.25, 0.32{]}       \\
\textbf{{[}70, 80) \%   }     & 0.29 {[}0.26, 0.32{]} & 0.29 {[}0.26, 0.33{]} & 0.34 {[}0.30, 0.37{]}    & 0.35 {[}0.32, 0.38{]}       \\
\textbf{{[}80, 90) \%   }     & 0.38 {[}0.34, 0.41{]} & 0.33 {[}0.29, 0.37{]} & 0.38 {[}0.34, 0.42{]}    & 0.39 {[}0.35, 0.42{]}       \\
\textbf{{[}90, 100{]} \%}     & 0.28 {[}0.17, 0.40{]} & 0.11 {[}0.05, 0.17{]} & 0.12 {[}0.07, 0.19{]}    & 0.15 {[}0.08, 0.23{]}      \\
\bottomrule
\end{tabular}}
\label{tab:ch1_02}
\end{table}

\begin{table}[h]
\centering
\caption[]{R-squared coefficients adapted for Bayesian inference \citep{gelman_r-squared_2019} calculated for model M5 - where trends in temperature (tT) and trends in anthropic area (tA) are the causal predictors of interest - for each 20\% interval of anthropic area group (see Methods), with their 95\% confidence interval.}
\begin{tabular}{lc}
\multicolumn{1}{l}{} & \begin{tabular}[c]{@{}c@{}}\textbf{Trends in anthropic area}\\ \textbf{percentage and in temperature}\\ \textbf{(Model M5)}\end{tabular} \\
\midrule
\textbf{{[}0, 20{]} \%}          & 0.40 {[}0.37, 0.44{]}                                                                             \\
\textbf{(20, 40{]} \%}           & 0.17 {[}0.09, 0.26{]}                                                                             \\
\textbf{(40, 60{]} \%}           & 0.08 {[}0.01, 0.18{]}                                                                             \\
\textbf{(60, 80{]} \%}           & 0.11 {[}0.01, 0.25{]}                                                                             \\
\textbf{(80, 100{]} \%}          & 0.32 {[}0.05, 0.45{]}        \\
\bottomrule
\end{tabular}
\label{tab:ch1_03}
\end{table}

\clearpage

\begingroup
\setstretch{1.2}
% \setlength\bibitemsep{0pt}
\printbibliography[heading=subbibnumbered, title = {Supplementary Information References}]
\endgroup

\end{refsection}