Prescribed burning mitigates the severity of subsequent wildfires in Mediterranean shrublands

Background Prescribed burning (PB) is becoming relevant in fuel reduction and thus fire hazard abatement in fire‑prone ecosystems of southern Europe. Yet, empirical evidence on the effectiveness of this practice to mitigate wildfire severity in Mediterranean shrublands is non‑existent, despite being the focus of PB efforts in this region. Here, we intended to quantify the protective effect of PB treatment units (2005–2021) to subsequent wildfire severity in shrublands across mainland Portugal, as well as the relative contribution and complex interactions between drivers of wildfire severity in PB‑treated areas and untreated neighboring counterparts through Random Forest regression. We leveraged cloud‑computing remote sensing data processing in Google Earth Engine to estimate fire severity (PB and wildfire) as the Relativized Burn Ratio (RBR) using Landsat data catalog. Results PB treatment was particularly effective at mitigating wildfire severity at the first PB‑wildfire encounter in shrublands, with a mean reduction of around 24% in RBR units. Fuel age (i.e., time since prescribed burning) in PB‑wildfire intersection areas overwhelmed to a large extent the effect of fire weather, burning probability, and PB severity. The mitigating effect of PB on wildfire severity persisted for a fuel age of around 5 years. However, this effect decreased with increasingly adverse fire weather conditions, such that variation in wildfire severity was somewhat insensitive to fuel age under extreme fire weather. Similarly, the lowest wildfire severity experienced in sites with high burning probability, along with the interaction effect observed between burning probability and fuel age, suggest that repeated PB treatments may be useful in controlling fuel accumulation and mitigating wildfire severity. The relative contribution of fire weather in explaining wildfire severity was exceedingly high in untreated areas, doubling that of the other variables in the model in the absence of PB treatment variables. Conclusions Our results suggest that the implementation of PB treatments at intervals of less than 5 years is of para‑ mount importance to control fuel build‑up and fire hazard under extreme fire weather in productive Mediterranean shrublands. Further research on this topic is warranted in other shrublands worldwide, namely in Mediterranean‑type climate regions.


Results
Los tratamientos de PB fueron efectivos para mitigar los efectos de severidad del fuego en las primeras interacciones de PB con áreas incendiadas no tratadas, con una reducción media de alrededor del 24% en unidades de RBR.La edad del combustible (i.e.el tiempo transcurrido desde la PB) en la intersección PB-incendio superó en gran medida los efectos de la meteorología del incendio, la probabilidad de quema, y la severidad de la PB.El efecto mitigador de la PB en la severidad de los incendios persistió en el combustible por unos 5 años luego de la PB.Desde luego, este efecto decreció con el incremento de condiciones adversas en la meteorología del incendio, tal como que la variación en la severidad del incendio fue no sensible a la edad del combustible bajo condiciones meteorológicas extremas.Similarmente, las severidades más bajas experimentadas en sitios con alta probabilidad de quema, junto con los efectos de las interacciones observadas entre probabilidad de quema y la edad del combustible, sugiere que la repetición de los tratamientos de PB puede ser útil para controlar la acumulación de combustible y mitigar la severidad de los incendios.La contribución relativa de las condiciones meteorológicas en la severidad del fuego fue notablemente alta en áreas no tratadas, duplicando el efecto de otras variables del modelo en ausencia de las variables de las PB.

Background
Wildfires are a natural and recurrent disturbance in the terrestrial ecosystems of Mediterranean-type climate regions (Pausas et al. 2008), and have shaped landscape pyrodiversity and adaptive plant traits in these regions for thousands of years (Keeley et al. 2012;Jones and Tingley 2022).However, strong fire regime shifts have been evidenced in Mediterranean ecosystems worldwide (Pausas and Fernández-Muñoz 2012;Fernandes et al. 2014;Lestienne et al. 2022).In southern Europe, land use and land cover changes in recent decades, including rural abandonment, and a deficit of adaptive fuel management have resulted in highly flammable landscapes (Moreira et al. 2011;Pausas and Keeley 2021), coupled with fire suppression policies and a consequent increase in large wildfires potential (Fernandes et al. 2020;Moreira et al. 2020).Concurrently, the effects of anthropogenic climate change, including increased heat waves and prolonged droughts (Giorgi and Lionello 2008;Molina et al. 2020), have led to fuel dryness conditions conducive to extreme fire events (Ruffault et al. 2018).As a result, large areas are expected to burn severely owing to the connection between large wildfire growth and extreme fire behavior (Parks et al. 2014a;Harvey et al. 2016).These events can have unprecedented impacts on ecosystem multifunctionality (Lasslop et al. 2019;Fernández-Guisuraga et al. 2023a), typically imply rapid changes in plant community structure and dynamics (Fernández-Guisuraga et al. 2019;Nolè et al. 2022), and are a significant threat to human lives and assets (Wunder et al. 2021), particularly in the wildland-urban interface (Bowman et al. 2017).To further complicate matters, the frequency of extreme wildfires is expected to increase in the future, and mitigation of their impacts should become a priority for land managers and environmental agencies (Duane et al. 2021;Jain et al. 2022).
Fire suppression has been the main focus of fire management policies in southern Europe since the second half of the 20th century, but it can exacerbate fuel accumulation and landscape scale fuel connectivity, precluding firefighting efforts by promoting extreme fire behavior under adverse fire weather (Moreira et al. 2020).In this context, the only feasible solution to constrain wildfire size and severity is to reduce fuels such that potential fire-growth rate and energy release are reduced (Fernandes 2015;Prichard et al. 2021).Prescribed burning (PB) offers a range of ecological benefits, suggesting higher cost-effectiveness than other fuel treatment options (Hesseln 2000;Hunter and Taylor 2022), and is used in a broad variety of fire-adapted vegetation types worldwide (e.g., Burrows and McCaw 2013;Ryan et al. 2013;Dems et al. 2021;Kupfer et al. 2022).In single-layered vegetation types, i.e., grasslands and shrublands, PB treatment is stand-replacing and fireline intensity tends to be high.In forests, PB is usually performed as a low-intensity understory fire to reduce crowning likelihood in subsequent wildfire encounters (Fernandes et al., 2015;Westlind and Kerns 2021).
The degree to which fuels are consumed and then reaccumulate is commonly used to assess the immediate and longer-term effectiveness of PB and will be reflected in modified fire behavior during wildfire reburns, decreasing wildfire extent passively or actively, i.e., through facilitated wildfire suppression operations (Fernandes and Botelho 2003;Espinosa et al. 2019).However, the effect of fuel reduction on fire-spread rate, and hence on burned area, may be negligible under extreme fire weather (Reinhardt et al. 2008;Cruz et al. 2022).Consistent with the need to consider fire-caused damage rather than wildfire extent (Moreira et al. 2020), PB effectiveness is more adequately captured by assessing changes in wildfire severity within PB-wildfire intersection areas (Reinhardt et al. 2008;Fernandes 2015), either in-situ or through remote sensing techniques (e.g., Arkle et al. 2012;Prichard and Kennedy 2014;Espinosa et al. 2019;Collins et al. 2023).
PB in southern Europe primarily targets shrubland and is characterized by low-treatment effort at the landscape level, consisting of small-sized (typically < 30 ha) dispersed treatments (Fernandes et al. 2013(Fernandes et al. , 2022)).Previous research on the effectiveness of PB in the Mediterranean basin is extremely scarce when examined from the angle of wildfire reburns as compared to other regions worldwide, such as Australia or North America (Fernandes 2018).Davim et al. (2021) assessed the likelihood of a fire intersecting a PB treatment in Portugal, and the unburned proportion of the treatment, finding that survival of PB treatments to wildfire is typically low.Moving to the landscape level beyond the limitations of study cases (Moreira da Silva 1997;Molina-Terrén et al. 2006) and simulation studies (Alcasena et al. 2018;Benali et al. 2021), the empirical analysis of Davim et al. (2022) reported a modest return-for-effort for the effect of PB on subsequent wildfire extent in Portugal.Fernandes (2018) indicates a number of European studies that have assessed how PB affects modeled fire behavior or severity, but the ability of current models to reliably depict how fuel treatments affect fire characteristics is debatable (Fernandes 2009;Cruz et al. 2014Cruz et al. , 2022)).Empirical evidence in this regard is quite scarce and is restricted to treateduntreated comparisons of observed fire behavior during experiments in a northern Portugal maritime pine (Pinus pinaster) stand, 2 to 13 years after prescribed burning (Fernandes 2009).Espinosa et al. (2019) examined wildfire severity in prescribed-burnt maritime pine stands in Portugal and Italy by sampling a limited number of sites, and concluded that PB reduces crown scorch height for 2-6 years after treatment; spatial variability due to variation in other factors influencing fire behavior and in stand height was manifest.Field studies on PB effectiveness in decreasing wildfire severity are often constrained by small sample sizes that may not encompass most of the variability in environmental conditions, as opposed to when relying on landscape-scale assessments based on remote sensing techniques (Arkle et al. 2012;Collins et al. 2023).
As the support to PB by policies and land managers increases in southern Europe, it is crucial to understand its role in mitigating the ecological impacts of extreme wildfires to inform adaptive management strategies, namely in the shrublands and grasslands that are the focus of PB efforts in the region (Fernandes et al. 2022).And as the incidence of extreme wildfire events increases in the context of global change, so will the need to scale up fuel treatment efforts and the frequency of wildfiretreatment intersections (Buma et al. 2020).Here, we intend to quantify the effect of PB in mitigating wildfire severity in shrublands across mainland Portugal by leveraging remote sensing data processing in geospatial cloudcomputing platforms.Specifically, we aim to answer the following questions: (i) does wildfire severity decrease within PB treatments compared to untreated neighboring vegetation?; (ii) does the contribution of fire weather to wildfire severity differ in treated versus untreated areas?; and (iii) do fuel-related attributes of PB (i.e., fuel reduction effectiveness and fuel recovery) interact with site productivity and topography in shaping wildfire severity?

Study area and data sources
The study area comprises mainland Portugal, an area of about 89,100 km 2 in the Iberian Peninsula (western Mediterranean Basin; Fig. 1).Elevation in the country ranges from sea level to around 2000 m.Topography is rough, particularly in the mountain ranges of the north and central regions, predominantly occupied by shrublands.The climate in mainland Portugal is mostly humid Mediterranean, with warm and dry summers, and cool and wet winters.Mean annual temperature and precipitation range between 7 and 18 °C and 400-2800 mm, respectively (Mora and Vieira 2020).The spatial database of PB units from 2005 to 2021 (Table 1; Fig. 1) was provided by the Portuguese Forest Service (ICNF) and derived from official treatment reports.PB is concentrated in central and northern Portugal (but also present in the southern region; < 1% of PB units), mainly in shrublands.PB units located in other fuel types (i.e., grasslands and forests; 11.47%) were discarded through photointerpretation of mainland Portugal orthophotographs (acquisition period 2004-2021; spatial resolution of 25-50 cm).The Portuguese atlas of wildfire perimeters for the period 1990-2022 (Fig. 1), also provided by the ICNF, was generated by supervised image classification and subsequent photointerpretation of Landsat 4-5 Thematic Mapper (TM), Landsat 7 Enhanced Thematic Mapper Plus (ETM+), and Landsat-8 Operational Land Imager (OLI) surface reflectance products (Oliveira et al. 2012).The minimum remotely-sensed mapping unit was 1-5 ha, but the fire atlas also included smaller fires from ground GPS records, namely during the PB period (2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018)(2019)(2020)(2021).
We obtained the PB-wildfire intersection areas by overlapping the spatial database of PB units and the atlas of wildfire perimeters.We discarded PB-wildfire intersections that occurred after the first PB-wildfire encounter.We retained for subsequent analyses intersections larger than 10 ha to be consistent with the spatial resolution of the remote sensing products and the sampling methodology used in this study.A 200-m buffer around the intersected PB areas was established to evaluate PB effectiveness in reducing wildfire severity and differences between wildfire severity drivers in treated areas and in their untreated neighboring counterparts.This limited buffer size was selected to minimize differences in fuel types and topography with the PB areas, which was visually assessed by means of mainland Portugal orthophotographs for the period 2004-2021, and by the European Digital Elevation Model (EU-DEM), version 1.1 with a 25-m grid size.There were no significant differences between the size of the PB-wildfire intersections and the buffers (p-value > 0.35 of paired sample t-test).
Wildfire severity was estimated using the Relativized Burn Ratio (RBR; Parks et al. 2014b) computed from Landsat TM, ETM +, or OLI Level 2 (Collection 2, Tier 1) surface reflectance products in Google Earth Engine (GEE; Gorelick et al. 2017) (https:// devel opers.google.com/ earth-engine/ datas ets/ catal og/) at a spatial resolution of 30 m.Despite the recent availability of multispectral satellite products at higher spatial resolution (e.g., Sentinel-2 since 2015), we have chosen Landsat to maintain spectral consistency throughout the study period (2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018)(2019)(2020)(2021).We chose a relativized metric of wildfire severity to ensure more accurate estimates than metrics of absolute change, e.g., the differenced normalized burn ratio (dNBR; Key 2006), in areas with low vegetation cover and in heterogeneous landscapes encompassing multiple vegetation types (Miller and Thode 2007).Contrary to other relativized metrics such as the relative dNBR (RdNBR; Miller et al. 2009), the RBR is mathematically more stable (Parks et al. 2014b).However, we also computed the dNBR as a wildfire severity metric.We used the code of Parks et al. (2018) as a basis in GEE, but adapted to initial assessments of wildfire severity, to compute the RBR and the dNBR indices from pre and post-fire, cloudfree satellite scenes as close as possible to the wildfire date.
We computed fire severity for the PB treatment (PB severity), fuel age, burn probability (BP), fire weather, and topographic texture as putative drivers of wildfire severity in PB treatments.The same variables, except PB severity and fuel age, were considered as potential drivers of wildfire severity in the surrounding untreated vegetation.Fuel age and severity of previous wildfires were not explicitly considered in buffer areas since fire intensity and severity are expected to be exceedingly higher in wildfires than in PB treatments (Arkle et al. 2012) and thus the leverage in fuel reduction.Also, fuel age is the same in intersection and buffer areas prior to the implementation of the PB treatment, and we discarded intersections after the first encounter as stated above.
PB severity for each treatment unit as a proxy for fuel reduction was calculated in GEE following the same procedure as for wildfire severity.Time since prescribed burning (TSPB), or fuel age, was the time (years) elapsed between PB and the wildfire interception.We calculated BP with a grid size of 30 m within each wildfire scar as a proxy for cumulative losses in site productivity in reburning areas (Davim et al. 2021) using the ratio of the number of times burned to the time elapsed from 1990 to the wildfire year, i.e., the inverse of the fire return interval (Fernandes et al. 2012).Fire weather description was based on the Canadian Forest Fire Weather Index System ( Van Wagner 1987).We retrieved the Fire Weather Index (FWI), which rates potential fireline intensity, for each wildfire from the ERA-5 FWI reanalysis dataset available from the Copernicus Data Store as 0.25° gridded data (Vitolo et al. 2020).Finally, the terrain ruggedness index (TRI; Riley et al. 1999) and topographic position index (TPI; Guisan et al. 1999) were calculated as proxies for topographic texture (Fernandes et al. 2016).TRI and TPI were computed from the EU-DEM version 1.1.

Data analyses
Wildfire severity and their drivers were extracted for each PB-wildfire intersection area and its untreated neighboring area using a regular grid of points separated at least 100 m.We fitted a Mann-Whitney test due to the violation of parametric test assumptions to evaluate PB effectiveness in reducing wildfire severity by determining the significance of the differences in the mean wildfire severity between treated and untreated areas.Statistical significance was determined at the 0.05 level.The dataset showed a Moran's I value for the response variables lower than 0.037, corroborating the lack of significant spatial autocorrelation (Moran's I < 0.1; Diniz-Filho et al. 2012).Random Forest regression (RFR; Breiman 2001) is an ensemble learning algorithm that was used to disentangle the relative contribution of the wildfire severity drivers in PB treatments.Therefore, the dependent variable was wildfire severity and the putative independent variables were the severity of the PB treatment, TSPB, BP, the FWI, and topographic texture.We opted for the RFR algorithm because it makes no assumptions regarding the distribution of dependent variables, provides relative variable contribution estimates, it can efficiently handle overfitting issues, outliers and spatial autocorrelation, generates unbiased error estimates, and can reveal non-linearities between the dependent and independent variables, as well as complex interactions among the independent variables (Cutler et al. 2007;Gigović et al. 2019;Fernández-Guisuraga et al. 2023b).The Boruta feature selection technique (Kursa and Rudnicki 2010) is a wrapper algorithm around RFR based on a permutation test using a holdout approach for variable importance measures and was used to determine which predictors are essential and at what magnitude (Hornero et al. 2021).The algorithm was implemented to ensure low redundancy among model predictors, while improving model interpretation, robustness, and predictive performance (Speiser et al. 2019).In this sense, Boruta has been widely implemented as a feature selection algorithm in previous fire ecology research (e.g., Vijayakumar et al. 2016;Eskandari et al. 2020), including prescribed burning (Martín-Pinto et al. 2023).Boruta provides Z-scores as important measures of predictors marked as unconfirmed, tentative, and confirmed against shadow variables (Kursa and Rudnicki 2010).We thus retained the predictors marked as confirmed by the Boruta algorithm in the prediction subset of the RFR model.The ntree RFR hyperparameter was set to 2000 to ensure stable predictions (Probst and Boulesteix, 2018), whereas the optimum mtry hyperparameter value was found by tuning through 10-fold crossvalidation.The internal out-of-bag error rate was used to assess the variance explained by the RFR model, and the percent increase in mean square error (%IncMSE) to quantify the importance of the predictors.The relationship between each predictor and wildfire severity in the model was examined through partial dependence plots.We inspected two-way interactions between the predictors through the H-statistic.The H-statistic measures the degree to which the interaction between the predictors contributes to the variability of the dependent variable, and ranges between 0 (no interaction strength) and 1 (the entire variability of the dependent variable is determined by a specific interaction).The strongest interactions were represented through three-dimensional partial dependence plots.
We followed the same procedure (Boruta algorithm and RFR) to unravel the relative contribution of the wildfire severity drivers in untreated buffer.We also fitted an RFR model calibrated solely with a binary treatment variable (PB, no treatment) and FWI to investigate if fire weather modulates the mitigating effect of PB on wildfire severity.

Results
PB treatments in mainland Portugal in the 2005-2021 period were frequently intersected by subsequent wildfires, 22.44% (excluding subsequent encounters) out of a total number of 3,124 PB units.78 of the PB-wildfire overlaps were larger than 10 ha (mean = 28.95ha; median = 22.65 ha; range = 10.27-165.72 ha; interquartile range = 15.59-30.82ha) and were thus considered in subsequent analyses.Summary statistics and distribution of the study variables in PB units and in the untreated buffers are displayed in Table 2 and Figure SM1 of the Supplementary Material.The PB-wildfire intersections covered a substantial range in TSPB, from 1 to 11 years, but 75% of the observations pertained to fuel ages lower than 6 years at the encounter time.PB severity was substantially lower than wildfire severity (t-value = 23.13;p-value < 0.001).FWI variation was wide both in the PB treatments and in the untreated buffers, ranging from 1 to 72, and was distributed continuously with no gaps up to FWI > 60.
Wildfire severity as per the RBR index was significantly lower (U = 113,471; p-value < 0.001) in PB treatments than in their untreated buffers (Fig. 2), with a mean 24% reduction.Untreated areas featured a high number of very low wildfire severity outliers.Results for the commonly-used dNBR index followed the same patterns than the RBR index (Fig. 2).

Wildfire severity in PB-wildfire intersections
The Boruta feature selection algorithm wrapped around RFR indicated that TSPB was the single most important driver of wildfire severity (median Z-score = 26.79) in PB-wildfire intersections, followed by FWI (16.21),BP (12.15) and PB severity (6.51) (Fig. 3).Topographic variables were deemed unimportant.
The RFR model calibrated from TSPB, FWI, BP and PB severity accounted for 57% of wildfire severity variability (RMSE = 112.79 in RBR units) in PB-wildfire intersection areas (Figure SM2 of the Supplementary Material).RFR variable importance (Figure SM3 of the Supplementary Material) followed the same pattern than the Boruta wrapper algorithm.Non-linear effects and potential thresholds are evident in the partial dependence plots for each wildfire severity predictor (Fig. 4).The mitigating   effect of PB on wildfire severity persisted up to a fuel age of around 5 years, but a slight increase in wildfire severity continued to be observed as TSPB increases above that threshold.Wildfire severity increased sharply with FWI, leveling off at FWI above 30.A very slight increase in wildfire severity up to BP values of 0.1 was evidenced, but above this value the relationship became inverse.The variability of wildfire severity with PB severity was the opposite of that observed for BP.Wildfire severity varied in response to the interaction between RFR predictors.Indeed, the interaction strength as measured by the H-statistic (Table SM1 of the Supplementary Material) was remarkably high for all four predictors, especially for TSPB and FWI (H-statistic ≈ 0.45), indicating that nearly 50% of variation in the predicted wildfire severity outcome depends on TSPB and FWI interactions, among them and with the other two variables.The mitigating effect of TSPB on wildfire severity decreased with increasing FWI (see an example in Figure SM4 of the Supplementary Material for PBwildfire intersections with low fuel ages and different FWI), such that variation in wildfire severity was somewhat insensitive to fuel age under extreme fire weather in 3-year-old fuels and beyond (Fig. 5).While wildfire severity was substantially decreased by high BP and both low TSPB and FWI, the effect of BP was undermined under high fuel age and extreme fire weather.Similarly, intermediate PB severity constrained wildfire severity more intensely when TSPB and FWI values were low.Finally, the combination of intermediate PB severity and low BP produced a slight decrease in wildfire severity.
Wildfire severity in the untreated buffers FWI was the most important driver of wildfire severity as represented by the Boruta feature selection algorithm within the untreated buffers surrounding PB-wildfire intersections (Fig. 6).The median Z-score for the FWI (40.17) was twice as large as that of BP (13.22) and TRI (12.70) variables.
The RFR model calibrated from FWI, BP and TRI for the untreated buffers accounted for lower wildfire severity variability (45%; RMSE = 121.83 in RBR units) (Figure SM5 of the Supplementary Material) than the RFR model in PB-wildfire intersections.RFR variable importance for the untreated neighborhood (Figure SM6 of the Supplementary Material) were similar to those provided by the Boruta wrapper algorithm.The non-linear effects and FWI and BP thresholds evidenced in the partial dependence plots for the untreated buffers (Fig. 7) were consistent with the findings for PB-wildfire intersections.Intermediate TRI values were related to the highest wildfire severity outcome, whereas a sharp decrease in wildfire severity was evidenced with TRI > 125.
As in PB-wildfire intersections, the interaction strength was remarkably high for the RFR predictors in the untreated buffer (Table SM2 of the Fig. 5 Three-dimensional partial dependence plots for the two-way interactions in prescribed burning (PB)-wildfire intersections using the Random Forests regression (RFR) algorithm Supplementary Material), particularly for FWI and BP (H-statistic > 0.45).Large increments in wildfire severity were predicted when high FWI and extremely low BP coincided with low TRI.The interaction between FWI and BP followed the same pattern but was less pronounced than in PB-wildfire intersection areas (Fig. 8).

Does fire weather modulate the mitigating effect of PB on wildfire severity?
The RFR algorithm calibrated from the binary treatment variable and FWI explained 38% of wildfire severity variability.The interaction strength between both variables was markedly high (H-statistic = 0.472).The mitigating effect of PB treatment on wildfire severity is strong even under severe fire weather and prevents the occurrence of high wildfire severity peaks (Fig. 9).In contrast, there are hardly any differences in the wildfire severity outcome under less elevated fire weather.

Discussion
PB treatments have local to sub-regional relevance in southern Europe to reduce fuels and thus fire hazard in fire-prone ecosystems (Fernandes et al. 2013(Fernandes et al. , 2022)).However, non-anecdotal evidence on the effectiveness of PB in mitigating wildfire severity in ecosystems other than maritime pine forest (Fernandes 2009;Espinosa et al. 2019) was absent up to now.This study is thus the first attempt to shed light into PB effectiveness in the region from this perspective, which is increasingly perceived as a more meaningful management target than decreasing wildfire extent on its own (Reinhardt et al. 2008;Moreira et al. 2020).
PB in shrublands was particularly effective at mitigating wildfire severity, with a mean reduction in wildfire severity of around 27% and 24% in dNBR and RBR units, respectively.Collins et al. (2023) conducted the first landscape-scale assessment involving empirical fire severity data in relation to PB effectiveness in shrublands.
The authors also found that treated areas in temperate shrublands in southeastern Australia experienced lower wildfire severity than those areas that were not previously treated.However, their results are not directly comparable with those of our study since Collins et al. (2023) analyzed categorized wildfire severity data that quantify the degree of foliage scorch and consumption.Our results are also consistent with reports of prior landscape-scale assessments based on remote sensing techniques in other vegetation types worldwide, mainly forests and woodlands.For instance, PB treatments in mixed-conifer forests of Washington, USA, reduced wildfire severity by around 25% in terms of RdNBR in the Tripod Complex fires (Prichard and Kennedy 2014).Similar results were reported in dry conifer forests in southern Idaho, where wildfire severity estimated using the dNBR index in three large PB treatments was about 20% lower than in untreated counterparts (Arkle et al. 2012).Hislop et al. (2020) also found a mean wildfire severity difference of Fig. 8 Three-dimensional partial dependence plots for the two-way interactions in untreated buffers in the vicinity of PB treatments using the Random Forests regression (RFR) algorithm Fig. 9 Three-dimensional partial dependence plots for the two-way interaction between binary treatment variable (prescribed burning (PB), no treatment (NT)) and fire weather index (FWI) using the Random Forests regression (RFR) algorithm around 20% in terms of dNBR between PB treatments and untreated counterparts in southeastern Australia's dry eucalypt woodlands and forests within the extent of the 2019 and 2020 Black Summer extreme fires.The congruence in PB effectiveness among our results and those of previous studies in forests and woodlands supports PB as an effective means of mitigating wildfire severity also in open Mediterranean vegetation, which precisely tend to reach high fuel consumption at less adverse fire weather conditions compared to forests (Keeley and Fotheringham 2001;Fernández-Guisuraga et al. 2021;Landesmann et al. 2021).Considering the simple structure of shrubland, i.e., fire spread is essentially supported by the elevated fuels (Anderson et al. 2015), the extent of the difference in wildfire severity between treated and untreated areas was higher than anticipated.This is particularly relevant since shrublands may be as important as forests regarding carbon sequestration in the Mediterranean Basin (Pasalodos-Tato et al. 2015), and are expected to be more vulnerable to postfire runoff and erosion than forests (Pierson and Williams 2016).
The high number of very low wildfire severity outliers in untreated buffers, manifest in the relativized metric of wildfire severity (i.e., the RBR), could be associated to the fire-shadow effect (Finney 2001;Loehle 2004), i.e., areas of lower than expected wildfire spread likelihood adjacent to treated areas, particularly on the treatment side opposed to the direction of fire spread.It is less likely that this phenomenon might be associated to non-vegetated and thus unchanged pixels at a spatial resolution of 30 m, due to the small buffer size chosen around PB-wildfire intersections.
Fuel age (i.e., TSPB) was the dominant driver of wildfire severity in PB-wildfire intersections, overwhelming to a large extent the effect of fire weather, which is consistent with the work previously mentioned and fire-size limitation in pyrodiverse shrub-dominated landscapes (Fernandes et al. 2016).However, the mitigation longevity of PB in shrublands barely extends beyond 5 years after treatment, presumably because of the fast fuel recovery enabled by the Mediterranean humid climate prevalent in northern and central Portugal.Generic total fuel loading equations estimate 12.9 Mg ha −1 for 5-yearold shrubland (Rosa et al. 2011), which corresponds to 43% of the modeled fireline intensity at steady-state (asymptotic) fuel accumulation (Fernandes et al. 2014).For a TSPB range of 0-8 years, Collins et al. (2023) also found that PB effectively decreased canopy defoliation for up to 3-5 years in temperate shrublands in southeastern Australia.Likewise, the effectiveness of PB treatments in conifer (e.g., Malone et al. 2011;Espinosa et al. 2019) and dry sclerophyll forests (e.g., Hislop et al. 2020;Leavesley et al. 2020) is greatly reduced for TSPB > 5 years because of litter and understory fuel build-up.Fire weather also related non-linearly with the severity of PB treatments reburned by wildfire, with a transition to high severity observed at FWI > 30.This threshold is lower than observed in prescribed-not burning at extremely high intensities under moderate fire weather conditions (Anderson and Anderson 2010).Short TSPB effectively mitigated wildfire severity, but this effect was overwhelmed by extreme fire weather for TSPB > 2 years as indicated by the strong interaction between the two drivers in PB-wildfire intersections, consistent with Collins et al. (2023) findings.Price and Bradstock (2012) reported that wildfire severity was also sensitive to fuel age at moderate fire weather conditions in a broad range of forest types in southeastern Australia.Regardless, the short-lived effect of PB under extreme fire weather conditions implies that wildfire severity mitigation will be less likely if such conditions become more frequent with climate change (Fernandes et al. 2013).This will be particularly relevant in more productive environments, as shown by the non-linear relationship between wildfire severity and BP as a proxy for loss in site productivity in recurrently burned areas.Specifically, the highest wildfire severity in treated areas is experienced in shrublands with BP < 0.15 and thus relatively low cumulative losses in site productivity and high fuel loading (Rosa et al. 2011).
Repeated PB treatments, with lower severity than wildfires, as evidenced here and elsewhere (Fernandes 2015), may be useful in controlling fuel accumulation and mitigating wildfire severity, along with the interaction effect observed between BP and fuel age.The relationship between wildfire severity and PB severity followed the opposite pattern of that observed for BP.PB tends to be more severe in shrubland with greater fuel load (Baeza et al. 2002;Knapp and Keeley 2006) and thus higher site productivity.Fuel buildup after PB is thus expected to be fast and reach higher levels, consistent with subsequent higher wildfire severity.This may account for the apparent interaction between PB severity and BP.The low relative contribution of topographic texture in explaining wildfire severity in PB treatments may be related to the overwhelming influence of the other drivers, but also because of the narrow TRI range, between level terrain and intermediately rugged, a feature of PB in Portugal (Davim et al. 2022).
The wildfire severity variability accounted for by the RBR model in the untreated buffer contiguous to PB units was moderately lesser than in the PB treatments.Jointly with the unbalanced, strong contribution of fire weather to wildfire severity in untreated areas (whereas it was overwhelmed by TSPB in PB-wildfire overlaps), demonstrates the robust modulating effect of PB on wildfire severity, particularly in productive landscapes (Prichard and Kennedy 2014).It could have made sense to consider fuel build-up in the untreated buffer, and therefore include time since the last wildfire as a variable that would improve model performance.However, fire severity and fuel consumption are remarkably higher in wildfires in comparison with PB (Fernandes 2015), and thus their post-disturbance fuel dynamics are expected to differ.Non-linear relationships and thresholds between wildfire severity and the FWI and BP were consistent with treated areas, but topographic texture, namely TRI, became important in the model probably due to the absence of the overwhelming contribution of PB treatment variables and higher topographic variability around PB units.A pronounced decrease in wildfire severity outcome was evident at higher TRI values, corresponding to possible fire behavior moderation by disruptions in fuel continuity by physical obstacles or changes in fuel moisture and wind speed (Guyette et al. 2002;Holden et al. 2009;Fernandes et al. 2016).On the contrary, flat or undulating terrain tends to experience higher wildfire severity (Chafer et al. 2004) as evidenced by the model outcome for lower TRI values, particularly under severe fire weather and high site productivity.
Variables that may be influencing wildfire severity (both in PB-wildfire intersections and the untreated buffers) have not been considered in this study, for example, fire spread patterns in relation to wind and slope directions, or wildfire suppression efforts (Prichard and Kennedy 2014;Viedma et al. 2020;Davim et al. 2021).These variables are not available at this scale of analysis and in an after the fact study as ours and would require direct observation during the wildfires.Future research could benefit from the recent availability of fine-scale products as proxies for fuel moisture such as the Landsat Collection 2 (C2) Level-3 Provisional Actual Evapotranspiration (ETa) Science Product (Senay et al., 2023) or the ETa product from the ECOsystem and Spaceborne Thermal Radiometer Experiment on Space Station (ECOSTRESS; Fisher et al. 2020).Although the effect of climatic conditions on fuel build-up after PB treatment has not been considered in this study and may be relevant variables given their high variability across mainland Portugal (Mora and Vieira 2020), less than 1% of the PB units are located in the less productive areas of the southern regions.Future efforts should also be targeted at assessing the effectiveness of PB treatments to mitigate soil fire severity in Mediterranean shrublands at regional scales.This is particularly relevant since spectral indices, such as those used here, are weakly sensitive to the variability of fire severity on the forest floor and substrate due to the background reflectance signal occlusion by the upper vegetation strata (Koetz et al. 2004).In this sense, future research could leverage the use of synthetic aperture radar due to its sensitivity to fire effects on soils (Fernández-Guisuraga et al. 2022).

Conclusions
We provided the first sound evidence of a mitigating effect of PB on wildfire severity in shrublands in Mediterranean-type climate regions.The mean reduction in wildfire severity at first encounter within PB treatment units indicates that fuel reduction treatments in shrubland are as effective as in woodlands and forests worldwide, which has substantial implications in Mediterranean regions where shrublands are the dominant land cover.The extent of the difference in wildfire severity between treated and untreated areas was higher than anticipated considering the less complex structure of shrublands in comparison to other vegetation types.
The longevity of the PB treatment effect evidenced here suggests that fuel-age mosaics resulting from repeated treatments and with a relevant fraction of < 6-yearold patches will succeed in moderating wildfire severity across the landscape.However, under more extreme fire weather, particularly in higher-productivity environments, it is likely that fire severity mitigation will be attained through avoided wildfire extent -through effective fire suppression operations enabled by fuel treatments -rather than by fire severity mitigation within treatments.Further research in this topic is warranted in other shrublands worldwide, namely in Mediterranean-type climate regions.Variable importance (%IncMSE) of the wildfire severity (RBR) descriptors in untreated areas of the vicinity of PB treatment units using the Random Forests regression (RFR) algorithm.Table SM1.Interaction strength of RFR predictors in prescribed burning (PB)-wildfire intersection areas.Table SM2.Interaction strength of RFR predictors in untreated areas of the vicinity of PB treatment units.

Fig. 2
Fig. 2 Variation in wildfire severity between prescribed burning (PB) treatments and their untreated (NT) buffers.Mann-Whitney test results are displayed in the upper area of each plot

Fig. 3
Fig. 3 Variable importance to wildfire severity assessed by the Boruta feature selection algorithm wrapped around Random Forest regression (RFR) against the shadow variables (shadowMin, shadowMean and shadowMax) in prescribed burning (PB)-wildfire intersections.The boxes in green represent variables marked as confirmed by the Boruta algorithm, in yellow as tentative, and in red as unconfirmed

Fig. 6
Fig. 6 Variable importance to wildfire severity assessed by Boruta feature selection algorithm wrapped around Random Forest regression (RFR) against the shadow variables (shadowMin, shadowMean and shadowMax) in the untreated buffers surrounding PB treatments.The boxes in green represent variables marked as confirmed by the Boruta algorithm, in yellow as tentative, and in red as unconfirmed

Table 1
Information on prescribed burning (PB) treatment units in shrublands across mainland Portugal for the period 2005-2021

Table 2
Summary descriptive statistics for the descriptors and putative determinants of wildfire severity RBR) in prescribed burning (PB)-wildfire intersection areas using the Random Forests regression (RFR) algorithm.FigureSM3.Variable importance (%IncMSE) of the wildfire severity (RBR) descriptors in prescribed burning (PB)-wildfire intersection areas using the Random Forests regression (RFR) algorithm.Figure SM4.Example of PB-wildfire intersections and untreated buffers with low fuel ages (2-years) and different FWI.Figure SM5.Relationship between observed and predicted wildfire severity (RBR) in untreated areas of the vicinity of PB treatment units using the Random Forests regression (RFR) algorithm.Figure SM6.