Estimating the Trade-Offs between Wildfires and Carbon Stocks across Landscape Types to Inform Nature-Based Solutions in Mediterranean Regions

Abstract

Climate and land-use changes have been contributing to the increase in the occurrence of extreme wildfires, shifting fire regimes and driving desertification, particularly in Mediterranean-climate regions. However, few studies have researched the influence of land use/cover on fire regimes and carbon storage at the broad national scale. To address this gap, we used spatially explicit data from annual burned areas in mainland Portugal to build a typology of fire regimes based on the accumulated burned area and its temporal concentration (Gini Index) between 1984 and 2019. This typology was then combined with carbon stock data and different landscapes to explore relationships between landscape types and two important ecosystem services: wildfire reduction and carbon stock. Multivariate analyses were performed on these data and the results revealed a strong relationship between landscapes dominated by maritime pine and eucalypt plantations and highly hazardous fire regimes, which in turn hold the highest carbon stocks. Shrubland and mixed landscapes were associated with low carbon stocks and less hazardous fire regimes. Specialized agricultural landscapes, as well as mixed native forests and mixed agroforestry landscapes, were the least associated with wildfires. In the case of agricultural landscapes, however, this good wildfire performance is achieved at the cost of the poorest carbon stock, whereas native forests and agroforestry landscapes strike the best trade-off between carbon stock and fire regime. Our findings support how nature-based solutions promoting wildfire mitigation and carbon stock ecosystem services may prevent and revert land degradation harming Mediterranean regions.

1. Introduction

Climate and land-use changes have made Mediterranean-climate regions increasingly prone to extreme fire events and aridity, leading to carbon emissions and land degradation [1,2,3,4,5,6,7]. On the one side, terrestrial ecosystems are a global greenhouse gas (GHG) sink but, in contrast, agriculture and forestry are major sources of anthropogenic emissions, making land-use change the key driver of every terrestrial environmental syndrome including wildfires, desertification, and climate change [8]. Land use may paradoxically prevent or contribute to wildfires, desertification, and climate change, depending on the quantity and longevity of either the carbon that is emitted into the atmosphere or stored in ecosystems [2,6,9,10]. Most Mediterranean countries have experienced a forest transition, which resulted in forest expansion during the 20th century [11]. Still, in some regions, such as mainland Portugal, the effect is being reversed as severe wildfires are mostly affecting the northern regions of the country, while desertification is growing in the south, leading to large carbon emissions [2,12,13,14,15,16].

In Portugal, as in other Mediterranean countries, wildfires are becoming more frequent, larger, and more hazardous in the last few decades, due to land-use changes and warmer and dryer climates [2]. In fact, its mainland territory is mostly affected by two types of fires, the very frequent fires, and the extraordinary wildfires that also tend to be extreme in size and severity [17]. Extraordinary wildfires become extremely large, intense, and above suppression capacity, causing significant ecological and socio-economic impacts, including relevant human fatalities and carbon emissions [7,15,18]. Frequent fire impacts include forest regeneration shortages, even for fire adapted species, and transformation of forests into shrublands, since fire recurrence might be shorter than trees’ reproductive maturity [13,19]. Besides direct ecological effects, frequent wildfires greatly contribute to lack of forestry practices and rural abandonment [20,21].

Landscape specialization, in Portugal and other Mediterranean countries, is characterized by agricultural intensification on fertile lands and farmland retreat on marginal lands, where forests have regenerated or were installed, often with fast growing species and originating in large forest patches, increasing fuel continuity and fire risks [6,10,22]. Nowadays, agricultural and forest plantations have created large-scale specialized landscapes, while the rest of the territory is abandoned where no viable socioecological solution is possible [23]. Native forests became only residual, and mixed landscapes with agroforestry systems or non-timber forestry systems—including silvopastoralism, cork, nuts and fruits—are mostly confined to marginal areas [24]. Although intensive production systems provide a significant part of food and timber to society, specialized landscapes produce negative impacts on the environment and society [25]. Without regulation services like fire protection (in plantation-dominated landscapes) or carbon storage (in farmland-dominated landscapes), these specialized landscapes are prone to disturbances like extreme wildfires, floods, droughts, ecosystem degradation, biodiversity loss, desertification, and land degradation, which are all amplified by climate change [3,26,27].

Land use is not the only driver involved in wildfire regimes and carbon storage, as there are other well-known drivers also affecting both processes (e.g., topography, climate, human activities). However, in Mediterranean-climate regions, land use/land cover (LULC) is the major driver that must be taken into account to address fire protection and carbon storage from a management perspective [3,10,17,22,28,29,30,31,32,33].

Previous research has pointed out the advantages of using landscape approaches to address climate and fire issues [31,32,34,35], with nature-based solutions being recommended to explore trade-offs and synergies across multiple ecosystem services [5,32,35,36,37,38]. According to most studies in Mediterranean regions, trade-offs between fire regulation and climate services are expected since increasing carbon storage implies forest and shrubland expansion, which may lead to increased fire risks [5,32,35,37,39]. However, the influence of specific LULC on fire regimes has not yet been sufficiently explored, mostly because of the limited geographic scope of such studies. The role of LULC types on fire regimes is currently under dispute, with some studies pointing to a strong association with extreme wildfire events [3,40,41,42,43], while others assign LULC types a minor influence on fire severity and size [44,45].

A wide region in central Portugal was recently studied focusing on the landscape level (represented by parishes), with the authors reporting a clear association between LULC and fire regimes, but even in these cases, there was a limited number and range of landscape types [22,30,46]. In another recent study for mainland Portugal assessing 12 fire regime drivers, the role of LULC was variable, depending on fire regimes, but very influential in determining extremely hazardous regimes [17].

This paper aimed to compare different landscape types at the national level regarding the fire regime and carbon storage levels that characterize each type. Therefore, the objective of the paper is not to explore the dynamics between fire regime and carbon stock, but rather to analyze how each of them change across landscape types. This is expected to help identify which landscape types strike the best trade-off between more favourable fire regimes and higher carbon stock. As the associations between landscape composition, diversity and configuration, on the one hand, and fire regime, on the other, have been previously explored [6,10], our goal here is to explore how these trade-offs vary across increasing levels of landscape specialization. We then aim to apply LULC landscape shares as modelling variables to improve our understanding of landscape behaviour for both fire protection and carbon stock. This framework involves a three-step approach: (1) it starts with an extensive characterization of fire regimes and carbon stock at the landscape level to find the landscape association with both ecosystem services; (2) secondly, following previous results, we build a predictive model associating LULC with the most hazardous fire regime, coupled with a carbon stock model; and finally, (3) we use the models to simulate scenarios of LULC change to promote wildfire mitigation and carbon storage. The results are used to explore potential nature-based solutions to improve fire regimes integrating carbon stock to find trade-offs at the landscape level [5,32,35,37]. We assumed carbon stock to be a proxy of climate-resilience ecosystem services [39], considering that maximizing carbon storage in the biosphere is a main objective of climate change mitigation [38]. For the carbon stock setup, we attributed to each LULC class a figure for the average of three carbon pools, namely the above- and below-ground biomass, and litter. Average carbon stock has the advantage of reflecting the net carbon sequestration balance, including all the anthropogenic and natural disturbances related to LULC (e.g., fire, droughts, harvesting) [9,36].

To our knowledge, the simultaneous assessment of fire regimes and carbon stock at such a wide regional (national) level is innovative and valuable in nature-based solutions for wildfire mitigation research, and due to its extensive geographical scope, it may contribute to uncovering the most relevant landscape associations with both fire regimes and carbon stock.

Our results are expected to bring new insights into landscape planning policy-making and contribute to reverting the reinforcing feedback loops of climate change, wildfires, and desertification, particularly in Mediterranean regions.

2. Materials and Methods

2.1. Study Area

Mainland Portugal is located in the SW of the Iberian Peninsula, occupying a surface area of c.a. 89,000 km², and its highest altitudinal point is 1993 m above sea level. The climate in Portugal is predominantly Mediterranean, characterized by dry and warm summers with north–south variations from the Atlantic to Mediterranean, respectively, from NW to SE, which splits the country in half (Figure S1 in supplementary material). Climate aridity in the south and east parts of the country has been increasing and is projected to be the main climate trend in the future [16]. Land use is mostly forest areas (forests, shrubs, and unproductive lands), occupying about 69.4% of mainland Portugal. Agriculture uses 23% of the total area and has been decreasing in the last few decades, dropping 13% in the period 1995–2015 [24]. However, the last agriculture census has recorded an expansion of the utilized agricultural area between 2009 and 2019 (8.4%), most of which is due to the expansion of permanent crops and grasslands [47]. According to the last forest inventory held in 2015 [24], maritime pine (Pinus pinaster) and eucalypt plantations (mainly Eucalyptus globulus) represent 48% of the woodland area, while evergreen broadleaf forests (mostly Quercus suber and Quercus rotundifolia) occupy 33%. The remaining woodlands represent only 18%. Shrublands and grasslands occupy nearly half of the forest areas (46%), but they seem to have declined for the first time since 1995 (−2.3%) [24].

2.2. Landscape Spatial Data

The spatial baseline data used in this study consisted of official sources, namely the 2018 Land Use/Land Cover map produced by the National Territory Directorate (https://www.dgterritorio.gov.pt/Carta-de-Uso-e-Ocupacao-do-Solo-para-2018, accessed on 5 May 2021) and the Sixth National Forest Inventory—2015 [24] from the National Forest Services. To simulate a landscape scope, a regular grid of hexagons was created, with a constant area of 3091 ha, which corresponds to the average size of the lowest-level administrative unit in Portugal (freguesia) for which some statistical data are available. The hexagonal shape was chosen based on several topological and compactness advantages associated with representing landscapes [48]. These hexagonal polygons were used as the analysis unit for all subsequent calculations. To avoid topologic inconsistencies and minimize the border effect, all hexagons less than 2600 ha within the study area (terrestrial Portugal Mainland) were discarded from further analyses.

• Landscape composition

The characterization of landscape composition was based on the 2018 LULC map, which was used in a raster format with a 100 × 100 m resolution pixel and simplified to consider only five major LULC classes: farmland (FAR), agroforestry (AGF), forest plantations (PLA), native forests (FOR), and shrubland (SHR). FAR included all farmland classes (i.e., croplands and grasslands); AGF includes all open woodlands (hereafter agroforestry); FOR indicates native forests (i.e., broadleaf, evergreen broadleaf, conifers and mixed) hereafter forests, and PLA stands for industrial forest plantations (i.e., Eucalyptus spp., Pinus pinaster and invasive species), hereafter plantations. This last LULC class includes eucalypt (mainly E. globulus) and maritime pine plantations, since they are intensively managed for timber production because they have shared the same spaces over the studied 36-year period, and eucalypt has been expanding since 1970, mostly at the expense of maritime pine [44]. Typical stands of both species are even-aged, understorey is intensively controlled, and clear cuts include the harvesting mode, with the main differences being the autochthonous character of maritime pine against the exotic eucalypt, which is coppiced for at least three rotations of 10–13-year intervals. Invasive species (mainly Acacia spp.) were aggregated in this class because their spatial distribution and total proportion, below 1% of the forest areas [24], does not contribute any relevant results. Other LULC classes (water bodies, urban, coastal, and rocky areas) as well as hexagons with a significant area (>20%) in these classes were not included in the analysis.

• Landscape diversity

Landscape composition provided the data to compute landscape diversity variables, keeping the pixel resolution and analysis units the same. The characterization of landscape diversity was focused on the number and proportions of land-use classes, corresponding to alpha diversity at the landscape/hexagon level [49]. Landscape diversity was summarized using three indicators: landscape richness (LR—the number of land-cover classes), the Shannon Diversity Index (SDI), and Shannon Evenness Index (SHEI) [50,51].

• Landscape configuration

The way these LULC categories are spatially distributed in the landscape was used to compute landscape configuration. Landscape configuration was characterized using five landscape metrics: Largest Patch Index (LPI), Edge Density (ED), Number of Patches (NP), Cohesion, and Contiguity [52]. All configuration metrics were computed using fuel presence as positive value (plantations, forests, agroforestry and shrublands), and null values for fuel absence (croplands and grasslands). A principal components analysis (PCA) following the eigenvalue >1 Kaiser criterion was performed to summarize the landscape configuration variables. We used the Kaiser–Meyer–Olkin (KMO) measure to assess the sampling adequacy of PCA, and Bartlett sphericity test for its significance.

• Landscape typology

A hierarchical cluster analysis (HCA) using the composition, diversity, and configuration variables was performed to build a landscape typology. Ward’s method [53] and the Squared Euclidean distance metric were used to produce a range of solutions from 4 to 11 clusters. The best solution in terms of cluster numbers was selected through expert judgement to gauge the sufficient variability of wildfire regimes and carbon stock patterns in the country.

2.3. Fire Regime

In this study, the fire regime followed the sensu stricto concept from [54], including only temporal and spatial distribution to determine the fire characteristics. The fire regime characteristics in our modelling setup included the cumulative percentage of burned areas (CPBA), computed over a 36-year period (1984–2019), and a temporal concentration index of the total accumulated burned area at the landscape level, assessed through the Gini Concentration Index (GCI), following recommendations from previous studies [22,30]. Both variables were used to create a two-axis system to establish meaningful fire regime categories at the landscape level, [55] using the same method used by [46] for Central Portugal [30,46,56]. The CPBA was compiled from the Portuguese Forest Services geographic database (https://www.icnf.pt/florestas/gfr/gfrgestaoinformacao/dfciinformacaocartografica, accessed on 20 February 2021) in order to compute the summation of all the burnt areas in each hexagon across the 36-year period. The burnt area measure is a simple metric that is accurate and available for the considered period for the study area. The minimum burned area (BA) criteria included 5 ha from 1984 to 2004 and 1 ha since 2005, and these were considered consistent for the present study. The CPBA represents the cumulative percentage of hexagon areas over the 36-year period, with values from 0% to 550.8% (e.g., 100% means that the area burnt during the 36-year period is equal to the hexagon’s total potential burning area). GCI was used as a concentration index, taken as a proxy of the forest fire intensity or large fire events [56]. The GCI values ranged from 0 to 1, representing the minimum to maximum temporal concentration of cumulative burning areas in each hexagon. A GI equal to 1 means that all the burnt area occurred in a single year.

By integrating GCI with CPBA, it was possible to compare hexagons with the same burnt areas but different concentrations over time, defining fire regimes based on the values of both indicators. Four quadrants were first obtained using the arithmetic averages of CPBA (52.4%) and GCI (0.84), and then converted to 3 categorical classes representing fire regimes in terms of both the cumulated burned area and concentration of fires. The first fire regime (FR1) was designated as small burned areas, corresponding to areas with minimal to no fires, where the CPBA is below 52.4%. A second fire regime (FR2) identifies hexagons with significant burned areas but with low temporal concentrations and less hazards (CPBA > 52.4%; GCI < 0.84), and a third fire regime (FR3) was assigned to hexagons with both significant CPBA and high GCI, which was therefore identified as the most hazardous fire regime.

2.4. Carbon Stock

Carbon stock was used as a proxy of climate mitigation ecosystem services [37,39]. The carbon information was used to estimate carbon stock for all LULC classes, namely farmland (6 classes), agroforestry and forests (11 classes) and semi-natural (grasslands and shrublands). Three carbon pools were considered, including the above- and below-ground biomass, and litter. Because of the lack of robust information, soil pool was not considered in our study, although above- and below-ground carbon represent the main sources for soil organic carbon [57,58]. Forest carbon data were sourced from the 6th National Forest Inventory [24], which measured 7627 field forest plots with inventory direct measures, from stratification data using 2015 aerial photographic land cover, from 356.359 photo-points, and included estimates of the volume and biomass affected by the wildfire events from 2016 to 2018 [24]. For this study, we used national carbon values averages for each species and expanded the field plot values to photo points. For agriculture and semi-natural areas, carbon information was taken from the 2018 NIR—National Inventory Report 2020 [15] (

Table 1. Above- and below-ground carbon values (Ton C ha⁻¹) and sources.

Species	Ton C ha−1	Source
Chestnut tree	86.76	IFN 1
Acacia	54.54
Oaks (deciduous)	45.8
Other conifers	41.06
Other broadleaves	33.5
Umbrella pine	32.39
Maritime pine	31.32
Eucalypt	25.76
Cork oak	25.53
Holm oak	22.97
Carob tree	16.99
Shrubs	13.72	NIR, 2020 [15]
Olive orchards	10.72
Other permanent cultures	9.94
Vineyards	3.67
Grasslands	1.47
Annual agriculture crops	0.62

¹ National average of carbon measurements from national inventory 2015 calculated including biomass wildfire events from 2016 to 2018 [24].

).

2.5. Comparing the Relative Performance of Each Landscape Type for Fire Regimes and Carbon Stock

To compare the performance of different landscape types for fire regimes, we tested the null hypothesis that landscape type and fire regimes are independent events, whose probability can be estimated by multiplying the probabilities of both events (Pearson’s chi-square statistics at 95% significance). Cramér’s V test was used to assess the degree of association between landscape types and fire regimes. For this purpose, we performed crosstabulation analysis, using the SPSS software (IBM Corp. version 28). Adjusted standard residuals were also used to test whether the observed frequencies were above/below those that would be expected by chance under the independent event hypothesis (p < 0.05).

To compare the performance of different landscape types regarding the carbon stock, we simply compared the averages of carbon stock across landscape types and performed a one-way Anova to test the null hypothesis that the observed differences were due to chance (<0.01). The squared ETA indicator [59] was also computed to assess the proportion of the variation in the computed carbon stock for all hexagons that is due to the differences across landscape types (as opposed to differences within types).

2.6. Simulation of the Probability of Occurrence of Hazardous Fire Regimes and Carbon Stock Levels Based on the Proportion of Forest Plantations and Its Conversion to Other LULC Classes to Obtain Corresponding Trade-Off Curves between the Probability of Hazardous Fire Regimes and Carbon Stock

After assessing fire regimes across landscape types, and according to our focus on extreme wildfire events, we refined the analysis of the strongest associated LULC class with the most hazardous fire regime (FR3), namely the forest plantations. This assessment allowed us to first understand fire regime relation, and then the effects on carbon stock. After that, a conversion of forest plantation to the other LULC classes (farmland, agroforestry, native forests and shrubland) allowed us to simulate the improvement of both ecosystem services at the landscape level. This was achieved after a simple mathematical transformation, according to [60], where a binary logit model was used to obtain the probability (p) of forest plantations associated with the most hazardous FR3, as follows:

$$p=\frac{1}{1+e^{-({\beta }_0+{\beta }_1{X}_1+\dots +{\beta }_n{X}_n)}}$$

(1)

where β_i are the logit function parameters, i = {0,…,n}, and x_i (i = {1,…,n}) are the corresponding independent variables related to LULC classes as percentages.

A linear regression assessed the forest plantation’s carbon stock, with successful tests for the tolerance and variance inflation factor, which rejected multicollinearity between the independent variables. Integration of the two models (binary logistic and linear regression) was performed after a simple complementary conversion to one of FR3 to represent the probability of avoidance of most hazardous FR3. Leaving the forest plantation class out of the function, as the reference class, the same methodology as multiple logistic binary and multiple linear regressions assessed the effects of plantation conversion into other LULC classes on the fire regime and carbon stock, respectively.

3. Results

3.1. Landscape Characteristics

The hexagonal grid of landscape units led to a total of 2595 valid analysis units, which were used in subsequent analyses. The landscape composition results, computed at the level of each of these hexagonal units, are presented in five maps representing the proportion (0–1) of each of the main LULC classes (Figure 1).

The landscape diversity results are represented in three maps, corresponding to the landscape richness, Shannon Diversity Index and Shannon Evenness Index (Figure 2).

The PCA performed on the configuration metrics variables yielded two PCAs with eigen values above 1, explaining 85.2% of the total cumulative variance. Both PCs were retained and used to describe the landscape configuration in subsequent statistical analysis (

Table 2. Loadings and percentage of variance explained for the 2 principal components extracted to describe landscape configuration. The total percentage of variance explained was 85.2%.

Variable	PC1	PC2
Cohesion	−0.880	0
Mean contiguity	0.121	0.990
Largest Patch Index	−0.905	0
Edge density	0.958	−0.106
Number of patches	0.856	0
Explained variance (%)	64.843	20.436

).

Sampling adequacy was confirmed (KMO = 0.754), and a significant correlation within the variables was found (Bartlett chi-squared = 9810.9; <0.0001). Also, all the communalities extracted above 0.5 confirmed the adequacy of the procedure. The first component explained 65% of the variance, while the second explained only 20%, hereafter represented by regression factors 1 and 2, which were renamed edge density (ED) and mean contiguity, respectively. The spatial distribution of these two new composite landscape configuration variables was assessed in two maps, depicting the scores of each analysis unit (hexagon) on both edge density (ED) and the mean contiguity axes (Figure 3).

The six-cluster solution was selected from the results of the hierarchical cluster analysis, based on the inspection of the resulting dendrogram and expert judgment, with each cluster representing a landscape type (

Table 3. Six landscape types with composition, diversity and configuration information. Maximum variable values are shown in bold.

Landscape Types	N	Composition					Diversity			Configuration
		Farmland	Agroforestry	Plantations	Native Forests	Shrubland	Shannon	Shannon Evenness	Richness	Edge Density	Mean Contiguity
1-Specialized agricultural landscapes	524	0.718	0.095	0.052	0.085	0.050	0.776	0.522	4.416	1.231	0.267
2-Mixed agroforestry landscapes	153	0.322	0.552	0.008	0.114	0.004	0.924	0.650	4.235	−0.445	−0.094
3-Mixed landscape with shrubland, farmland and native forests	609	0.290	0.009	0.162	0.183	0.356	1.183	0.817	4.291	−0.266	−0.283
4- Mixed native forest landscapes	598	0.250	0.155	0.123	0.416	0.056	1.164	0.733	4.906	−0.564	0.286
5-Mixed landscapes with plantations and farmland	459	0.287	0.003	0.494	0.108	0.109	1.065	0.738	4.283	−0.177	0.039
6-Specialized plantations landscapes	252	0.125	0.001	0.778	0.048	0.048	0.671	0.481	4.032	−0.830	−0.541
Total	2595

, Figure 4 and Figure S2). In fact, this was the solution with a better compromise to represent the major landscape types with sufficient segregation of different LULC categories, and to reflect the relevant variability of wildfire regimes and carbon stock.

Cluster 1 grouped the landscapes dominated by farmland (72%), while the remaining area was occupied by agroforestry (10%) and native forests (9%), with very low values of both the Shannon Diversity (0.78) and Shannon Evenness (0.52) indexes and medium landscape richness (4.4), having been designated as “specialized agricultural landscapes”. It is worth noting that this assessment does not include diversity within agricultural crops, but only between farmland and other LULCs. These landscapes were mostly found in the central and southern parts of the country (Figure 4) and include intensive irrigated areas as well as dryland low-intensity farming, being associated with the highest average edge density, reflecting high patch perimeter to area ratios and patch numbers (woody patches).

In cluster 2, agroforestry represents 55% of the land, farmland represents 32%, and native forests represents 11%. These landscapes presented low richness and diversity (LR = 4.2; SDI = 0.92; SHEI = 0.65), and intermediate ED, although it must be highlighted that agroforestry often includes two different vegetation strata, and thus it may be vertically diversified if not horizontally diversified. It was, therefore, designated as “mixed agroforestry landscapes”.

Cluster 3 corresponded to a mixed landscape, combining shrubland (36%), farmland (29%), native forests (18%), and forest plantations (16%). It is a cluster with very high diversity (LR = 4.3; SDI = 1.18; SHEI = 0.82), and medium ED because it only contrasts the farmland composition with all other instances of LULC. This cluster was therefore named “mixed landscape with shrubland, farmland and native forest”.

In cluster 4, native forests covered 42% of the area, agroforestry covered 16%, and farmland covered 25%, presenting the highest richness and high diversity (LR = 4.9; SDI = 1.16; SHEI = 0.73), and an intermediate negative ED, revealing a low exposure to the farmland border. It was therefore named “mixed native forest landscapes”. This landscape cluster corresponded mostly to evergreen oak woodlands (Q. suber and Q. rotundifolia) and umbrella pine (Pinus pinea), but also to the rest of woodland species, such as deciduous oaks (Q. robur, Q. pyrenaica, Q. faginea, Q. canariensis), chestnut, other pines, and riparian species.

In cluster 5, forest plantations represented 49% of the area, although farmlands represent 29% and native forests and shrubs about 11% each. The cluster demonstrates low diversity, richness and edge density, and therefore was designated “mixed landscapes with plantations and farmland”.

Cluster 6 was the most specialized landscape type, with 78% of the area covered with forest plantations, while ca. 13% are farmlands, and all the other LULC instances are residual. It showed the lowest scores for diversity, richness, and edge density (SDI = 0.67; SHEI = 0.48; LR = 4.03; ED = −0.83), and therefore was named “specialized plantations”.

3.2. Fire Regime

Burnt areas from the years 1984 to 2019 showed a divided country where northern half landscapes experienced burning much more than the southern landscapes. SW mountains were the only exception in the southern half (Figure 5a).

FR1, here designated as small burned areas, dominated most of the southern half of the country, but also the coastal areas in the centre and north, and the extreme NE. FR2, meaning frequent but time-dispersed fire events and thereby related to less hazardous effects, occurs mostly in the northern parts of the country. FR3 was present in the centre and SW mountains and it was related to significant burned areas that are highly concentrated in time, which were associated with the most hazardous effects on the environment and societies (Figure 5b).

3.3. Associations of Fire Regimes with Landscape Types

According to our cross tabulation results, we rejected the null hypothesis that landscape types and fire regimes were independent events, with high confidence (X²(10) = 1332.769; p < 0.001). In addition, a strong association was found between landscape types and fire regimes (Cramér’s V = 0.52). The adjusted standardized residuals have shown that, for all cells, the magnitude of the deviation between the observed value and the value that would be expected based on the independent event hypothesis was such that we can also reject this hypothesis for each cell; the signal of the residual told us whether the observed value was above (+) or below (−) the expected value (

Table 4. Comparison of observed frequency of the three fire regimes and average carbon stock across landscape types.

Landscape Types		Fire Regime (FR)			Ton C ha−1
		1	2	3
1-Specialized agricultural landscape	Proportion of FRs	97.1%	1.5%	1.3%	8.3
	Adjusted Residuals	+	−	−
2-Mixed agroforestry landscape	Proportion of FRs	100%	0%	0%	13.4
	Adjusted Residuals	+	−	−
3-Mixed landscape with shrubland, farmland and native forest	Proportion of FRs	31.4%	46.5%	22.2%	14.6
	Adjusted Residuals	−	+	+
4- Mixed native forest landscapes	Proportion of FRs	90.0%	1.7%	8.4%	16.2
	Adjusted Residuals	+	−	−
5-Mixed landscapes with plantations and farmland	Proportion of FRs	46.2%	24.0%	29.8%	17.8
	Adjusted Residuals	−	+	+
6-Specialized plantation landscapes	Proportion of FRs	28.2%	8.3%	63.5%	22.8
	Adjusted Residuals	−	−	+
TOTAL	Proportion of FRs	64.5%	16.6%	18.8%

Adjusted residual values are all significant (p < 0.05).

).

Our results showed that specialized agricultural landscapes and mixed agroforestry landscapes were mostly associated with FR1, which represented 97% and 100% of the hexagons in these two landscape types, respectively. The positive signs of the adjusted standardized residuals showed that FR1 was more frequent in these landscape types than would be expected if fire regimes and landscape types were independent events. Mixed landscapes with shrubland, farmland and native forests were mostly associated with FR2, with 46% of hexagons in this landscape type demonstrating this fire regime, and a positive sign of the residual; FR1 only represents 31% of the cases and occurs less than expected (negative residual). The most hazardous FR3 represented 22.2% of the cases in this type but occurred there more often than expected based on random values (positive residual).

Mixed native forest landscapes were mostly associated with FR1, which occurred in 90% of the cases and exhibited a positive residual; the other two fire regimes occur less often in this landscape type than would be expected by chance and represent a small frequency of cases. In mixed landscapes with plantations and farmland, FR1 represented 46% of the cases, but was less frequent here than expected (negative residual); the other two fire regimes were reported more frequently than expected if we assumed fire regimes and landscape types were independent of each other (positive residuals) and represent 30% (FR3) and 24% (FR2) of the cases.

Specialized plantation landscapes were mostly associated with FR3 (64% of the cases), which were reported more frequently than expected by chance. The other two fire regimes had lower frequencies in this landscape type and represented less than what could be expected if fire regimes and landscape types were independent of each other. Comparing the two previous landscape types suggests that increasing the weight of forest plantations in the landscape clearly increases the probability of occurrence of FR3.

3.4. Carbon Stock

The computed values of landscape-level carbon stocks varied between 0.6 Ton C ha⁻¹ and 62.3 Ton C ha⁻¹ (Figure 6). Carbon peaked in the NE region of the country, corresponding to a mixed landscape with chestnut, oak and maritime pine stands, followed by a coastal centre region up to 31.7 Ton C ha⁻¹ dominated by maritime pine, and the SW Atlantic coastal triangle, mainly composed of mixed pine stands (P. pinea and P. pinaster) and eucalypt up to 29 Ton C ha⁻¹. Comparing carbon stock averages across landscape types revealed that all the landscape types had statistically different average carbon stocks (p-value < 0.001). Almost half of the variance in carbon stocks across hexagons (ETA² = 0.45) was explained by the landscape type (Table 4).

The lowest carbon stock average was found in specialized agricultural landscapes (8.3 Ton C ha⁻¹) and the maximum was found in specialized plantations (22.8 Ton C ha⁻¹), with the medium value reported in mixed agroforestry landscapes (13.4 Ton C ha⁻¹), mixed landscapes with shrubland, farmland and native forests (14.6 Ton C ha⁻¹), mixed native forest landscapes (16.2 Ton C ha⁻¹) and mixed landscapes with plantations and farmland (17.8 Ton C ha⁻¹). Our integrated carbon stock and fire regime results showed a negative trade-off across specialized forest plantation landscapes (Table 4).

3.5. Associations between Forest Plantations, Hazardous Fire Regimes and Carbon Stock

Our binary logistic regression results have shown that there was a positive and significant relation between the most hazardous FR3 and the proportion of forest plantations (PLA) in the landscape (

Table 5. Estimated coefficients of the binary logistic regression.

Variable	Beta	SE	95% CI		EXP (B)	p
			LL	UL
(Intercept)	−2.88	0.098			0.056	0.000
PLA %	0.042	0.002	1.038	1.047	1.043	0.000
N	2766
Nagelkerke’s R2	0.268
Prediction accuracy	0.846
Chi-square sign	0.000

). The model showed an 85% accuracy in predicting the occurrence of FR3, based on the proportion of forest plantations alone. Nagelkerke’s R² indicated that 27% of the variance was explained by the independent variable.

Based on the estimated parameters of the logistic model, the probability of the fire regime FR3 occurring in a hexagon (dependent on the proportion of forest plantations, PLA%) was given by:

$$p_{FR3}=\frac{1}{{1+e}^{-(-2.88+0.042\ast PLA\%)}}$$

(2)

A graphical representation of the fitted values for plantation proportions up to 100%, showed a positive association with FR3, rising in line with the specialization of the landscape in forest plantations (Figure 7).

The results of the linear regression that predicts carbon stock from the proportion of forest plantations in the landscape revealed a positive and significant relationship between both variables:

$$\begin{gathered} C\left(\mathrm{T}\mathrm{o}\mathrm{n}.{\mathrm{h}\mathrm{a}}^{-1}\right)=11.368+0.142\ast PLA\% \\ \left(\mathrm{F}\right(1.2764)=1567.596;p<0.001;{\mathrm{R}}^2=0.362 \end{gathered}$$

(3)

Using functions (2) and (3), we simulated the joint change in the probability of avoiding (not having) FR3 and the expected value (average), as we increase the proportion of forest plantations in the landscape. This simulation resulted in a trade-off curve that shows how increasing the proportion of forest plantations in the landscape will increase the carbon stock, but at the expense of the probability of having a fire regime more benign than FR3 (Figure 8).

The results of the multiple logistic binary regression revealed a significant relationship between FR3 and LULC (chi-squared = 778.083; p < 0.001) with a percentage of correct predictions of FR3 occurrence of 87%, and Nagelkerke’s R² of about 40% (Table S1 in supplementary material). The graphical representation of the fitted values confirms that increasing plantation conversion to any of the other LULCs improves the probability of avoiding FR3 (Figure 9a). The results also showed that the improvement of the fire regime of the landscapes dominated by forest plantations depends on replacing forest plantations by other LULC classes, and that the efficiency of such improvement depends on the LULC that is replacing the forest plantations. In fact, the cost of a similar reduction in the probability of avoiding FR3, in terms of the loss of the forest plantation area, is smaller when we convert it to farmland rather than agroforestry, and is smaller with a conversion to agroforestry than to native forest. So, in a landscape with 100% forest plantations, reducing the probability of FR3 to 20% (cf. solid line for 80% avoidance of risk of FR3 in Figure 9a) would mean losing almost 40% of the forest plantation area if we are converting forest plantations into farmland, 50% if we convert the area to agroforestry, and almost 80% if we convert the area to native forests. The same reasoning would apply to the more ambitious target of reducing the probability of FR3 to 10% (avoidance of risk of 90%, dashed line in Figure 9a). Assuming that forest plantations have been chosen by the owners because they are more profitable in that context, this area cost will translate into a financial loss. Note that the curves in Figure 9a (and also Figure 9b and Figure 10) do not take into account the stochasticity of the dependent variables, which are regression-based predictions, and thus integrate an error term. This could be accounted for by representing confidence intervals for those predictions. However, with so many curves, it would be too cumbersome, and thus we opted for this simpler representation to improve the communication of the main messages.

However, replacing the forest plantation area with other LULC classes also leads to a reduced carbon stock. A significant multiple regression was estimated with carbon stock as the dependent variable and the proportions of different LULC classes as the independent variables (F (4.2761) = 1251.644; p < 0.001; R² = 0.645) (Table S2). This model was used to predict the reduction in carbon stock as we replaced forest plantations with each of the other LULC classes, starting with a landscape with a 100% forest plantation proportion (see the four lines in Figure 9b for each of the other LULC classes). Here, we can see that the cost of converting forest plantations into other LULC classes, in terms of the carbon stock loss, also depends on the alternative LULC we are considering. In this case, the highest cost is when we use farmland as an alternative and the lowest is when we use native forest as the alternative.

Using the relationships represented in Figure 9a,b and keeping the forest plantations as the reference class (now omitted), we simulated the joint change in the probability of avoiding FR3 and the expected average carbon stock value, as we increased the proportion of forest plantations converted to each of the other LULCs at a specific time (Figure 10). These simulations reveal the marginal carbon cost of increased fire security. There is a particular trade-off curve for each LULC class used to replace forest plantations, which underlines the dependence of such a marginal cost on the alternative LULC class. The optimum for each alternative LULC class occurs when the slope of the corresponding trade-off curve (marginal cost) is equal to the ratio between the per-unit benefits of fire risk avoidance and carbon stock. The choice between these four optima (one for each alternative LULC class) will depend on the market (food, timber, etc.) and other non-market (biodiversity, water, etc.) values associated with each LULC class.

4. Discussion

4.1. Methodological Prospects

Our results suggest that robust models representing a landscape’s fire protection and carbon storage performance can be built with basic fire and carbon information, which are usually available in Mediterranean-climate regions. This approach allowed the decoupling of the two studied ecosystem services to optimize models for each component, and then integrate them to explore trade-offs in respect to different LULC classes to address climate change mitigation in the land system and inform nature-based solutions [61].

Our findings on fire regimes confirmed the advantages of integrating the total burned area (CPBA) with the temporal concentration of fires (GCI), and the use of averages in both indicators, following the work of [46] and their results that fully agree with studies with different scales and time periods, including [22,56] and [17], who designed a complete and detailed fire regime map for mainland Portugal, accepted by the Portuguese fire authority Agency for Integrated Rural Fire Management [17]. Considering that both ecosystem services are priorities for the society and the environment, our focus on FR3 comes from the recent 2017 extreme wildfires that had extensive impacts, including over 100 casualties [62].

Regarding landscape types, our cross-tabulation results have shown that the FR is highly dependent on landscape composition, diversity, and configuration, in agreement with most previous studies [3,17,22,42,43,46]. Our carbon stock results, including above- and below-ground biomass and litter, confirmed their robustness as a proxy for climate research, where soil carbon information is not available or accurate enough. Although the innovative approach of joining forest and agriculture biomass carbon does not allow direct comparisons, our findings confirmed their dual relevance in sourcing soil organic carbon, in line with desertification and aridification studies [57,58,63,64], and also in feeding GHG emissions [61].

The estimated logistic and linear models proved their effectiveness in exploring the relationships between forest plantations and the most hazardous fire regime and carbon stock, confirming positive and significant results with each one, and a strong negative trade-off between them. Similar but multiple models were also correctly used to assess the joint effect on FR3 and carbon stock of forest plantation conversion to the other LULC classes.

4.2. Fire Regime and Carbon Stock Associations across Landscapes

Our results showed that landscape types can be divided in three groups according to fire regimes; the most hazardous fire regime is strongly associated with forest plantations, suggesting the effects of its extensive distribution with very low diversity and a fire-prone configuration; the less hazardous fire regime is associated with mixed landscapes with shrubland, farmland and native forests; a third group including specialized agricultural landscapes, mixed agroforestry landscapes and mixed native forest landscapes is associated with small burned areas.

Integrating fire regimes and carbon stock showed that the proportion of forest plantations in the landscape have a marked inverse relationship with carbon stock and FR3. But that was not the case when comparing agriculture with agroforestry landscapes and the latter showed a better trade-off between the two indicators than the former (higher carbon stock and lower FR3 odds). Comparing forest plantation landscapes with native forests also revealed a better trade-off for native forests than for plantations at the landscape level.

Our findings on forest plantations agree with most of the previous research, namely [17,19,22,40,41,43,46,65,66] for maritime pine evaluation. Assessing only eucalypt [44] found no effect on the burned area within three decades of its expansion but attributed the prevalence of mega-fires (>10 Kha) to regional forest specificities, including the extension, continuity, and homogeneity of eucalypt stands. A partial agreement is found with [30] who did not find the eucalypt contribution to increase the GCI, although it was justified by its high post-fire regeneration rate. However, comparison with our study should be carried out with caution because these authors used a different grouping of forest classes [30]. Specialized landscapes of forest plantations reach the highest average carbon stock of 22.8 Ton C ha⁻¹ due to their significant presence, although the individual values for each of the species (maritime pine and eucalypt) are only intermediate compared with other forest species. On the contrary, the country’s carbon peaks are found in hexagons dominated by native forests, but they are too scarce to constitute a landscape typology. We recall that forest plantation landscapes group eucalypt and maritime pine stands in pure and mixed stands with both species, where management efforts can range from intensive to virtually absent across the entire territory [44]. For these reasons, maritime pines and eucalypt have often been considered together, including in more recent studies [17].

Mixed landscapes with shrubland, farmland and native forests were mostly associated with less hazardous fire regimes, in agreement with other authors [17], who found this very frequent fire regime associated with traditional fire use on pastures and shrublands. Although serving as a mosaic builder and preventing forest encroachment, this fire regime may also start wildfires if uncontrolled [17], as these land covers are strongly associated with wildfire initiation, especially in peri-urban or farmland territories [44].

We found a strong association of specialized agricultural landscapes with small burned areas, in full agreement with previous research in the Iberian Peninsula, which attributed farmland a strong role in decreasing fire risk, burned areas and fire hazards [17,22,30,41,67]. However, as specialized agricultural landscapes are mostly occupied by farmlands, they hold the minimum carbon stock, the lowest diversity and the poorest richness, which may contribute to amplifying desertification processes in most of the semiarid southern half of Portugal and the NE region [16,63].

Mixed agroforestry landscapes have a similar fire regime as agricultural landscapes, followed by mixed native forest landscapes, with the latter demonstrating a higher carbon stock, in full agreement with previous literature [30,40,41,65].

4.3. Implications of Carbon Stock and Fire Protection Trade-Offs for Nature-Based Solutions

The outcomes of our research suggest that nature-based solutions should be based on changing LULC towards a better-balanced distribution of closed and open landscapes, decreasing their specialization and thus contributing to less hazardous fire regimes and increasing carbon stock to prevent land degradation processes.

Specifically, in forest plantation specialized landscapes, our findings show that converting forest plantations to farmland could have a successful and immediate effect on fire prevention but with a strong negative trade-off with carbon stock, which would be better-off if conversion to agroforestry or to native forests instead is performed. However, this native forest solution is not so efficient at improving fire regimes when forest plantations largely dominate landscapes, as a larger converted area is required for such an improvement.

Carbon stock could also benefit from a conversion of eucalypt plantations into maritime pine stands, as changing maritime pine management to close-to-nature sylviculture could promote longer carbon storage forest products [68], with potential additional benefits for biodiversity [69,70,71], although this should be further researched.

Our results on farmland performance in wildfire prevention agree with other authors [22] who found a minimum threshold of ca. 40–50% of farmland in the Central Portugal region to keep fire regime at low-hazard levels, and also with previous research who considered spatial compartmentation vital to minimizing wildfire impacts, not mere fuel brakes, and further conversion to native forests for which species have much lower flammability [2,3,17,22,29,32,34,35,72].

Concerning specialized agriculture landscapes, our results suggest that increasing forestland proportions in these areas would improve diversity and provide carbon inputs to soils [57], as biomass limitation is a land degradation indicator present in all Mediterranean and semi-arid regions [13,63]. Native forests are the best LULC to include because besides having the highest unitary carbon stock, they also increase landscape diversity and total biodiversity [49]. These carbon stock outcomes on specialized agricultural landscapes allow relevant nature-based solution planning, considering the improvement of most fragile soils demands carbon sourcing from vegetal residues like composts, biochar, or manures [57]. Improving carbon stock could also benefit from changing occupations within farmland LULC, where the categories exhibit a wide carbon stock range (orchards have a 17-fold higher carbon stock than annual crops). However, a holistic solution to prevent land degradation would need field and farm-level approaches to effectively couple carbon with biogeochemical cycles within vegetation, organic matter, and microbial biomass [73], which is beyond the scope of the present study.

Consideration of the shrubland class does not provide any relevant fire regime or carbon stock effect, but promotes landscape diversity and biodiversity [49] and our results advise that for mixed landscapes with shrubland, farmland and native forests, the promotion of natural succession, and fuel connection interruption segregating shrublands from forestlands to prevent small fires to expand are necessary. Regulating traditional fire use and complementing it with prescribed burning seems a broad recommendation from specific studies [74].

Regarding the mixed agroforestry landscapes and native forests, our findings suggest increasing the forest proportion would improve carbon stock without a significant decrease in fire regulation indicators, keeping all the advantages of the mosaic diversity and configuration [75]. At the stand level, native forests could also benefit CS by increasing density and deciduous hardwood presence in their compositions.

Native forests are, according to our results, the best land-use solution for fire regimes, carbon stock and diversity, but it is, however, time-consuming to regain market value, as most tradable products take longer to be achieved, compared to all alternative land uses. Considering most land is privately owned, priority should be given to policies promoting native forests, regarding these trade-offs.

Together, the effects of landscape specialization are increasing GHG emissions by both extreme fire events and soil fertility decay, implying a net territorial decarbonization that already affects half of the country, and tend to form positive feedback loops with climate change and desertification processes [2,3,16,26,63]. The lack of landscape diversity across specialized landscapes does not help the recovery from extreme disturbances and ecosystem degradation, as resilience is based on biodiversity [76,77].

In addition to wildfire prevention, climate change may imply broad bioclimatic shifts in forests involving productivity reductions for forest plantations and the development of evergreen oaks, carob trees, and umbrella pine in Atlantic areas where they would substitute forest plantations [78]. Public perception [44] and official reports from the previous 2017–2018 wildfires in Portugal are in general agreement with our findings [66,79,80]. Forestry and climate policies have already been approved to stop eucalypt expansion (Law 2017/2017, 17 August), and decreasing plantation areas within the Portuguese Carbon Neutrality Roadmap [81], but these are still very modest steps towards solving the magnitude of the problem. To reverse rural abandonment, another recent policy measure in Portugal was the creation of integrated areas for landscape management AIGP (Decree-Law 28-A/2020, 26 June), which formalized in only 2 years, with 70 AIGPs already in operation, covering a total area of 140 Kha (https://www.dgterritorio.gov.pt/paisagem/ptp/aigp accessed on 18 January 2023).

4.4. Limitations and Uncertainties

The scope of this research included the provision of new evidence on the performance of different landscape types regarding wildfire avoidance and carbon stocks. Considering the landscape level for all calculations, this may result in biased values for fire regimes if the burned areas are divided by several polygons, or if it concentrates only on one and recommending groups of polygons should be considered rather than isolated polygons.

In terms of statistical limitations, the values of Nagelkerke’s R² for the two logistic regressions revealed that part of the observed fire regime variance was dependent on factors besides LULC, which were not considered here. In fact, these accuracy values were expected since similar studies demonstrated that LULCs were predominant in most hazardous fire regimes modelled by our regressions, in agreement with [3,10,17,28,29,30].

Further attention on the fire regimes should be given to potential uncertainties due to the temporal mismatch between the fire information over a 36-year period (1984–2019), and the fire regime characterization because the LULC data are only from the year 2018. Minimization of this question includes the species that change the most within plantation LULC classes, with eucalypt mostly replacing maritime pine [44], and the fact that past fire history shapes the current landscape [21,46].

Our carbon research was focused on the terrestrial phase of the carbon cycle, referring to the current average (assumed close to the steady state) biomass of each LULC, so it does not include the forestry carbon footprint or forest-based generated products. Short-lived forest products include bioenergy, pulp and paper, while building, furniture and non-timber products are associated with longer carbon lifetimes. Life cycle assessments should be carried out to provide robust information on complete carbon cycles, which is out of the scope of this paper.

5. Conclusions

We showed that land-use changes improving fire regulations and carbon storage may contribute to reverting the Mediterranean land degradation trend. The three main findings seem relevant when implementing nature-based solutions at the landscape level to effectively address both ecosystem services.

Firstly, the requirement of assessing several interrelated ecosystem services to find integrated solutions for environmental problems, especially when they are prone to complex trade-offs and feedback, was noted. In such cases, pragmatic information may require individual analysis (e.g., integrating the Gini Index and total burned area), keeping land use as the independent variable, given that it drives both fire regimes and carbon stock.

Secondly, the key role of landscape type showed that highly specialized landscapes degrade at least one of the studied ecosystem services, in contrast with mixed landscapes, which seem to strike a better balance between ecosystem services. Specialized maritime pine and eucalypt plantation landscapes store higher amounts of carbon but are strongly associated with the most hazardous fire regime, transforming carbon in atmospheric emissions and significantly damaging people and ecosystems. Agricultural-dominated landscapes may not burn as much, but unless carbon is increasingly stored in agroecosystems, they will not help prevent desertification.

Finally, LULC trade-off analysis confirmed the need to drastically reduce forest plantation proportions in landscapes to improve its fire profile, keeping at least 50% of landscapes as farmlands (crops and grasslands) and agroforestry. From this point, conversion to native forests seems the long-term solution, because they reach the best trade-off between the carbon stock and fire regime, followed by agroforestry, which is also associated with minimal fires and high biodiversity. The major opportunity cost for forest owners of conversion to native forests may be the loss of market value of the systems’ output, which may require public support. Shrublands have low to intermedial carbon stock, but high biodiversity and are mostly associated with frequent but less hazardous fire regimes. These five major LULC classes are needed to conveniently monitor improvement for both ecosystem services at the landscape level, including farmlands, forest plantations, native forests, agroforestry and shrublands.

In place of individual fire or desertification measures, all these findings should integrate nature-based solutions to promote a resilient landscape to cope with climate change and effectively revert extreme wildfire and land degradation occurrences.