Fostering Carbon Credits to Finance Wildﬁre Risk Reduction Forest Management in Mediterranean Landscapes

: Despite the need for preserving the carbon pools in ﬁre-prone southern European landscapes, emission reductions from wildﬁre risk mitigation are still poorly understood. In this study, we estimated expected carbon emissions and carbon credits from fuel management projects ongoing in Catalonia (Spain). The planning areas encompass about 1000 km 2 and represent diverse ﬁre regimes and Mediterranean forest ecosystems. We ﬁrst modeled the burn probability assuming extreme weather conditions and historical ﬁre ignition patterns. Stand-level wildﬁre exposure was then coupled with fuel consumption estimates to assess expected carbon emissions. Finally, we estimated treatment cost-efﬁciency and carbon credits for each fuel management plan. Landscape-scale average emissions ranged between 0.003 and 0.070 T CO 2 year − 1 ha − 1 . Fuel treatments in high emission hotspots attained reductions beyond 0.06 T CO 2 year − 1 per treated ha. Thus, implementing carbon credits could potentially ﬁnance up to 14% of the treatment implementation costs in high emission areas. We discuss how stand conditions, ﬁre regimes, and treatment costs determine the treatment cost-efﬁciency and long-term carbon-sink capacity. Our work may serve as a preliminary step for developing a carbon-credit market and subsidizing wildﬁre risk management programs in low-revenue Mediterranean forest systems prone to extreme wildﬁres. pole-stage (LP), high pole-stage (HP), and timber stage (T).


Introduction
Carbon emissions from fires in European countries, and Mediterranean areas, in particular, represent about 4.03 Tg C per year [1,2]. However, a few extreme fires account for the bulk of burned areas and carbon dioxide emissions to the atmosphere [3,4]. These are extreme escaped fires that easily overwhelm firefighting capacity, and the containment is restricted to strategic locations in backing and flanking fire spread areas [5,6]. This is the reason why recent studies advocate for a comprehensive long-term solution to better coexist with fire, and forest fuel management is emerging as a fundamental strategy complementing fire suppression and ignition prevention in populated areas [7][8][9].
Fuel treatments mitigate wildfire intensity on treated locations and reduce wildfire spread and overall wildfire likelihood at the landscape scale [10][11][12]. Commonly imple-Land 2021, 10, 1104 2 of 23 mented prescriptions in Mediterranean forests are low pruning, thinning from below, prescribed fires, mechanical mastication, and targeted grazing [13][14][15]. On the other hand, the most frequent treatments at the landscape scale consist of fuel break networks dividing the landscape into large planning areas, systematic buffer clearings within the wildlandurban interface, and intensively managed scattered treatment units known as strategic management points (SMPs) [16][17][18]. Wildfire managers design SMPs based on expert criteria, heading-fire pathway analysis, and fire-weather synoptic conditions [19]. The implementation of SMPs is gaining extensive adoption in many fire-prone Mediterranean areas, such as Catalonia, where fuel treatments have become the primary risk reduction strategy [20].
Primary fuel treatment objectives in Mediterranean areas are protecting human communities, timber harvesting, restoring fire-adapted ecosystems, reducing extreme wildfire potential, and facilitating safe and effective fire suppression [21,22]. Moreover, some studies have shown that forest management increased carbon sequestration [23,24] and developed management prescriptions under different wildfire hazard scenarios based on stand structure conditions [25,26]. However, protecting carbon stocks from extreme fires was often considered a secondary objective in southern European regions, and the reduction in emissions remains unknown in large-scale ongoing fuel reduction programs. Previous works conducted in the western US implemented quantitative assessment methods that coupled fire simulation modeling with stand-level carbon loss functions to estimate the effect of fuel treatments in reducing carbon emissions from wildfires [27][28][29]. These studies found that thinning weight plus prescribed fire treatment intensity and the frequency and severity of future wildfires determine the benefit of increasing carbon pools in managed forests.
Multifunctional Mediterranean forests currently provide meager economic revenues from timber harvesting. Over the last 50 years, limited management plus the prevalent fire exclusion policy favored fast fuel buildup in open woodlands and dense forest regeneration in clear areas [30,31]. Moreover, large-scale afforestation campaigns planted extensive areas over marginal lands with pines to facilitate forest development on poor soils [32]. As a result, the fires rapidly evolved from fuel-limited short-events to weather-driven catastrophic events [33]. In central Catalonia (northeastern Spain), these new forest structures fostered stand-replacing high severity wildfires associated with increased mortality rates on the non-serotinous sub-Mediterranean conifers forests [34,35]. Similarly, extreme fires that burned broadleaved forests triggered significant changes in stand structures (i.e., conversion to coppice forest or scrublands), hampering woodlands' management and economic viability of many products, such as cork oak production and extraction [36].
The European Union assumed the commitment to reduce emissions by at least 40% by 2030, and forest ecosystems were identified as core natural carbon pools with the potential capacity to store and compensate carbon emissions considerably (COM/2019/640). As a result, the European Union green deal articulated the action plan towards a green and climate-neutral sustainable economy in the EU by 2050. It will mobilize €100 billion over 2021 to 2027 to assist the most affected regions. Among many other actions, promoting a bio-economy in local markets is one of the main strategies for overcoming the fossil fuel-dependent economy [37,38]. Likewise, encouraging the assumption of voluntary carbon credits by large corporations to compensate for their emission may represent an additional funding source to reverse the landscape degradation occurring in many EU rural areas [39,40]. Wildfire managers, in turn, require a quantitative assessment to evaluate carbon sink capabilities and request financial compensation for fuel management programs in low-revenue, poorly managed fire-prone areas.
This study explores the opportunity for carbon-credit-oriented forest management in diverse fire-prone Mediterranean areas across Catalonia (northeastern Spain), where large historical wildfires caused substantial losses in the last decades [36,41]. Providing standlevel and spatially explicit quantitative results is essential to identify high emission hotspots and redirect management efforts at sufficient intensities. Specifically, we addressed the following research questions: (i) What is the expected wildfire carbon emission under Land 2021, 10, 1104 3 of 23 the current conditions? (ii) What is the overall reduction in emissions per treated area for planned fuel reductions? and (iii) where could carbon credits help subsidize ongoing risk-reduction fuel management programs? The methods presented in this study may help the small local landowners increase their revenues for the services provided to the community and vindicate large-scale fuel reduction programs in low timber revenue Mediterranean landscapes and elsewhere. Our main goal is to advance the ecosystem service carbon markets that incentivize farmers and rural communities to generate fireresilient landscapes adapted to extreme wildfires. We ultimately attempt to develop and implement a sound framework for assessing emission reductions and provide the technical blueprint required by large corporations and EU authorities to prove the emission reductions by fuel treatment programs.

Study Area
We conducted this study in different planning (n = 6) areas of Catalonia (northeastern Spain) (Figure 1). The extent of the planning areas corresponded to the existing landscape units of Catalonia [42]. The climate is predominantly Mediterranean, with increasing rainfall on pre-littoral mountain ranges, milder winters closer to the coastline to the east, and a transition to sub-Mediterranean and mountain climate in the Pyrenees. Catalonia is one of the most significant fire-prone regions in the Mediterranean Basin and encompasses various physiographic gradients and fire regimes ( Figure 1B). On average, 650 fires burn 11.5 thousand ha per year, from which 2% of large fires (>100 ha) account for more than 88% of the burned area (1983 to 2015). A few extreme events (i.e., >1000 ha, catastrophic wildfire episodes of 1986,1994,1998,2003, and 2012) make up the bulk (>65%) of the affected areas [43]. Most fire ignitions (>90%) are caused by humans and show apparent occurrence patterns clustering close to anthropogenic features, such as communication corridors and urban areas [44]. See vegetation description and wildfire history in Appendix A for further details.

Forest Fuel Treatments
Treatment spatial allocation corresponded to the strategic management points (SMP) provided by the Bombers GRAF team (Group of Support to Forest Actions) of the Firefighters of the Catalan Government [22]. These typically locate on watershed divides or ravine bottoms that intersect with primary flow paths to restrict the fire potential and

Landscape Data
We generated the landscape file assembling topography, surface fuel, and canopy metric grids at a 40 m resolution as required for fire spread simulation modeling using ArcFuels [45]. The modeling domain encompassed the extent of the planning areas plus a 10 km buffer to account for the fires incoming from the neighboring regions and predict realistic wildfire likelihood estimates. Adjacent planning areas (i.e., 2 and 3, and 5 and 6) were merged in the same fire modeling domain. Topographic data, including elevation, aspect, and slope grids, were derived from a 25 m resolution digital terrain model (ign.es). We assigned surface fire behavior models [46] to the 1:5000 scale land-cover vector map of Catalonia [47] to generate the surface fuel model grid. For that purpose, we considered vegetation characteristics, such as species composition, tree cover, thickness, and shrub, plus herbaceous fuel thickness and heights gathered from the 4th National Forest map and the 2012 habitat map of Catalonia [48]. The canopy metrics consisted of canopy height, canopy cover, canopy base height, and canopy bulk density and were obtained from the LiDAR-derived 20 m resolution biophysical variable grids of Catalonia [49].

Forest Fuel Treatments
Treatment spatial allocation corresponded to the strategic management points (SMP) provided by the Bombers GRAF team (Group of Support to Forest Actions) of the Firefighters of the Catalan Government [22]. These typically locate on watershed divides or ravine bottoms that intersect with primary flow paths to restrict the fire potential and create better opportunities for an effective fire response using a tactical fire [21,50]. Likewise, the design accounts for the existing operational constraints, including the prescribed fire implementation and the accessibility to the site. In addition to treatment units in the planning areas, we also considered the SMPs in neighboring areas within a 10 km buffer ( Figure 2).    We assumed the most widely implemented treatment prescriptions, including a tree thinning followed by a prescribed fire in conifer forests or mechanical mastication in broadleaved forests. Specifically, we considered a heavy-weight thinning from below, where ladder fuels (i.e., suppressed and dominated trees) were cleared to a 60% residual basal area target [26]. The canopy metrics were reduced accordingly to a canopy cover of 60%, a canopy bulk density of 0.09 kg m −3 , and a canopy base height of 3.6 m [51]. We assigned timber litter TL 1 and TL 2 surface fuel types for conifer and broadleaved treated forests [46]. The treatment cost was provided by local forest managers [52] (Table 1). We annualized the costs assuming a treatment duration of 8 years for surface fuel reductions and 14 years for the thinning, as described in previous works conducted in the study area [53,54]. The stand-level treatment cost varied between 100 and 228.75 € ha −1 per year. Table 1. Summary table with treatment cost data for fuel reductions in SMPs. The overstory is first reduced to a 60% canopy cover implementing a thinning from below, followed by prescribed fire (in conifer forest) or mechanical mastication (in broadleaved forests). We assumed average cost conditions for ongoing fuel reduction plans.

Wildfire Occurrence
We trained a machine learning Random Forest model [55] using fire ignition location data ( Figure 1B) retrieved from the EGIF Spanish fire database [43] to predict fire ignition probability from a set of geospatial data grids. These layers included proxies for accessibility (distance to roads, distance to forest tracks, and distance to trails), agricultural activities (distance to croplands), human pressure on wildlands (distance to the wildland-urban interface), and potential sparks from power-lines (distance to power lines) [56]. The model was calibrated using a k-fold cross-validation (k = 4) procedure using a sample of 10,835 fires from 1998 to 2015. At each step or fold, we split the sample into a training sample with 75% of fires to fit the model and a 25% validation sample to evaluate the model, calculating the area under the receiver-operating (ROC) curve (AUC) [57]. The sub-models were then combined into the final model used to predict the probability of fire occurrence, presenting an average AUC > 0.73. We generated a 40 m resolution ignition probability grid, where values ranged between 0 and 1 ( Figure 3). This raster grid was used as input in fire simulation modeling to locate required fire ignition patterns on the different fire modeling domains.

Fire-Weather Scenarios
The fire-weather conditions present substantial region-wide differences across Catalonia [58]. Therefore, we considered multiple fire-weather macro-areas or pyromes in the planning areas to accurately replicate the changing gradients on historical burn patterns ( Figure 1B). The delimitation of the fire-weather macro-areas was based on climatic and physiographical conditions and historical large fire footprints [59]. The wildfire season was the annual period concentrating 90% of the burned area from large fires (>100 ha; see Appendix A). We considered all the large fire ignitions within a 10 km buffer to compute the wildfire season's duration, and the analysis was conducted together for neighboring planning areas [43].
We then characterized the fire-weather conditions occurring during the wildfire season from hourly temperature records, rainfall, wind speed and direction, relative humidity, and solar radiation from different automatic stations using the Fire Family Plus v 4.2 software [60]. For each macro-area, we identified a representative weather station with a long data record (>15 years). Since containment efforts are very effective under mild weather conditions [61], we assumed 97th percentile extreme weather conditions in wind speed for most frequent wind directions and ERC-G fuel moisture content to generate the fire modeling weather scenarios. Specifically, we developed a set of probabilistic fireweather scenarios in terms of fire spread duration, wind speed, wind direction, and fuel moisture content combinations ( Table 2).  Figure 1B) and geospatial features using Random Forest. The numbers refer to the planning area codes (Appendix A; Figure 1).

Fire-Weather Scenarios
The fire-weather conditions present substantial region-wide differences across Catalonia [58]. Therefore, we considered multiple fire-weather macro-areas or pyromes in the planning areas to accurately replicate the changing gradients on historical burn patterns ( Figure 1B). The delimitation of the fire-weather macro-areas was based on climatic and physiographical conditions and historical large fire footprints [59]. The wildfire season was the annual period concentrating 90% of the burned area from large fires (>100 ha; see Appendix A). We considered all the large fire ignitions within a 10 km buffer to compute the wildfire season's duration, and the analysis was conducted together for neighboring planning areas [43].
We then characterized the fire-weather conditions occurring during the wildfire season from hourly temperature records, rainfall, wind speed and direction, relative humidity, and solar radiation from different automatic stations using the Fire Family Plus v 4.2 software [60]. For each macro-area, we identified a representative weather station with a long data record (>15 years). Since containment efforts are very effective under mild weather conditions [61], we assumed 97th percentile extreme weather conditions in wind speed for most frequent wind directions and ERC-G fuel moisture content to generate the fire modeling weather scenarios. Specifically, we developed a set of probabilistic fire-weather scenarios in terms of fire spread duration, wind speed, wind direction, and fuel moisture content combinations (Table 2).  Figure 1B) and geospatial features using Random Forest. The numbers refer to the planning area codes (Appendix A; Figure 1).

Fire Spread Simulation Modeling
We used the minimum travel time (MTT) algorithm [62] as implemented in FCon-stMTT to model fire spread [59]. The algorithm calculates a two-dimensional fire growth at a resolution set by the user by minimizing fire travel time from the cell corners based on Rothermel's fire spread model [63]. The MTT has been widely used in previous studies assessing wildfire transmission and fuel treatment effects in complex terrains worldwide [12,64,65]. The fire ignitions were first distributed within the modeling domain according to the ignition probability grid ( Figure 3). Then each fire was independently modeled considering subarea-level 97th percentile extreme fire-weather conditions for the wildfire season (Table 2).
To calibrate the surface fire spread model (see Appendix B), we replicated the actual large fire size (>100 ha) distribution as well as average fire size in every planning area separately, except for adjacent units ( Figure 1A). Multiple day events were decomposed into daily blow-up progressions and constant fire-weather conditions during the fire's duration [66]. We set the fire spread durations that better replicated the historical fire size distribution under extreme weather conditions. The fire modeling domains were saturated with thousands of fires that burned each pixel more than 20 times on average and a total area equivalent to 10,000 years of iterations. Fire suppression efforts were omitted due to their limited containment capability during extreme fire events in the study area [67]. As a result, we obtained the pixel-level annual burn probability as [68]: where aBP is the annual burn probability determined at a 40 m resolution as the n number of times a given xy pixel burned divided by the Y total number of modeled wildfire years or iterations. Table 2. Fire-weather scenarios for extreme weather conditions (97th percentile) during wildfire season for the different planning areas. While the neighboring planning areas ( Figure 1) were merged in a single fire modeling domain or landscape file, we considered multiple fire weather macro-areas. We used wildfire-season automatic weather station data to generate fire-weather scenarios. We set the fire spread duration probabilities that replicated historic fire distributions under constant extreme fire-weather conditions. See modeled and predicted fire size distributions in Appendix B.  (12) 12 (28) 12 (34) 10 (20) 10 (6) --Mountain peaks (Cadí-Nord) --12 (8) 12 (9) 14 (20) 14 (34) 12 (29) (7) 37 (35) Plain of Lleida (Tárrega) -19 (6) 20 (24) 19 (17) 19 (17) 18 (36) --Inner mountains (St. Salvador Guardiola) -10 (5) 15 (32) 10 (41) 10 (13) 10 (9) --Pre-coastal depression (Font-Rubí) --14 (7) 19 (58) 19 (23) 18 (7) 19 (5) -Coastal belt (Torredembarra) -20 (21) 15 (20) 16 (36) 18 (23) ---

Expected Carbon Emissions
We estimated pixel-level expected emissions from wildfires combining modeled wildfire likelihood estimates with conditional fire effects from surface fires as [69]: where eE CO 2 is the expected carbon dioxide gas emission (T ha −1 yr −1 ), aBP is the annual burn probability (%), cFC is the conditional fuel consumption (T ha −1 ), and EF CO 2 is the emission factor (T T −1 , dry basis). Emissions are, therefore, expectations estimated as the burn probability times consequences [70]. We used EF for stand-level dominant tree species, ranging between 1415 and 1879 g kg −1 [71]. The conditional fuel consumption (cFC) for surface fuels was first estimated by fuel category, including litter, duff, 1 h, 10 h, 100 h, and 1000 h using available models from the literature [72], and then summed to a 40-m resolution grid. Next, we assumed extreme fire-weather conditions (Table 1) and stand-level data on the forest floor and above ground dead fuel biomass dry weight per ha for the dominant forest types in the study area to calculate the fuel consumption (see Appendix C). Finally, the reduction in carbon dioxide gas emissions was estimated as the difference between the current conditions and the managed scenario or treated landscape in SMPs (Figure 2). We focused on assessing the reductions due to aBP differences and did not consider emissions from treatments due to the lack of reliable data for ongoing works.

Cost-Efficiency Analysis
First, we estimated the pixel-level difference, as CO 2 T ha −1 per year, between the current conditions and the managed scenario. Then, the CO 2 emission results in managed scenarios were determined after implementing fuel reductions in SMPs (Figure 2), with the treatment cost and prescriptions described in previous sections. Specifically, we modeled the reduction in burn probability for the managed scenarios to assess the treatment effects across the landscape. These results were summarized by the planning area and treatment unit (i.e., SMP). We then estimated the potential revenue from carbon credits considering a €13.18 T −1 reference market price for prevented emissions [73]. Finally, we determined the potential contribution of the carbon credits to financing the treatment cost, as a percentage, in each planning area. See the methodological flowchart on Appendix D.

Expected Carbon Emissions
Expected carbon emissions from wildfires varied largely within and among the different planning areas ( Figure 4). The average results ranged between 0.003 and 0.060 CO 2 T ha −1 per year at the landscape scale. Primarily, burn probability was the main causative factor explaining high expected emissions, despite the existing differences in stand-level conditional fuel consumption. We found the highest contrast between the neighboring forest lands of the planning areas 2 and 3, where similar stand structures (e.g., timber-stage conifer forests) located at less than 15 km of distance showed differences greater than 100 times in expected emissions. The strong north-south climatic gradient explains this difference because southern portions show drier fire-weather conditions ( Table 2) and a higher wildfire occurrence ( Figure 3). The open plains of planning area 4 and high-fuel load unburned black pine mature forests of central areas in planning area 3 presented the highest emission values from among all forest types. Most black pine mature stands of planning area 4 (i.e., mature timber-stage black pine forests) showing remarkably high conditional and expected emissions (>0.2 CO 2 T ha −1 yr −1 ) were located in unburned patches or islands of the latest extreme fire episodes ( Figure 1B). Conversely, recently burned areas in planning area 4, including low pole-stage Aleppo pine and mixed forests, presented a lower potential emission (<0.02 CO 2 T ha −1 yr −1 ). On the other hand, incoming fires arrived from all around in planning areas 5 and 6 ( Figure 4), where high emission areas tend to concentrate in the central portions. These areas were primarily dense and mature Aleppo pine stands exposed to long-distance spreading extreme wildfire events.
As expected, we found the lowest values in planning area 2 (Figure 4), where more than 80% of the mountainous areas in the north presented very shallow emission values (<0.001 CO 2 T ha −1 yr −1 ). Except for some rare events occurring in conifer plantations, the wildfires in high-elevation areas (>1500 m.a.s.l) of the Pyrenees hardly ever burn more than 100 ha of forest lands. Winter fires associated with dry and gusty foehn winds can also burn reduced forest patches in mountainous alpine areas, but these are isolated low-frequency episodes with a limited potential compared to Mediterranean fires [74]. Nonetheless, southern wind-driven fires occasionally exposed south-facing lower slopes of the neighboring east-west orientation mountain range in planning area 2. We also Land 2021, 10, 1104 9 of 23 found this topographic pattern in planning area 1, but the transition was smoother because the central valley and dominant southern winds occurring during the wildfire season presented the same orientation. In planning area 4 (Figure 4), a large portion to the west showed shallow values (<0.001 CO 2 T ha −1 yr −1 ) mainly because of the lower wildfire occurrence and a northern wind direction channeling the fires through the central valley.
than 80% of the mountainous areas in the north presented very shallow emission values (<0.001 CO2 T ha −1 yr −1 ). Except for some rare events occurring in conifer plantations, the wildfires in high-elevation areas (>1500 m.a.s.l) of the Pyrenees hardly ever burn more than 100 ha of forest lands. Winter fires associated with dry and gusty foehn winds can also burn reduced forest patches in mountainous alpine areas, but these are isolated lowfrequency episodes with a limited potential compared to Mediterranean fires [74]. Nonetheless, southern wind-driven fires occasionally exposed south-facing lower slopes of the neighboring east-west orientation mountain range in planning area 2. We also found this topographic pattern in planning area 1, but the transition was smoother because the central valley and dominant southern winds occurring during the wildfire season presented the same orientation. In planning area 4 (Figure 4), a large portion to the west showed shallow values (<0.001 CO2 T ha −1 yr −1 ) mainly because of the lower wildfire occurrence and a northern wind direction channeling the fires through the central valley.

Carbon Emission Reduction in Managed Scenarios
The landscape-scale emissions in managed scenarios were reduced by between 49 (planning area 6) and 444 CO 2 T per year (planning area 3), which resulted in a reduction between 11 (planning area 6) and 35% (planning area 2) compared to the non-managed conditions (Table 3). Interestingly, a high percentage reduction was not associated with a high decrease in CO 2 T emission per year (e.g., see planning area 2). In planning area 2, most wildfire activity concentrated in the southern portion, where we concentrated many of the treatments, thus explaining the high reduction when presented as a percentage. On the other hand, the highest reduction per planning area extent corresponded to the planning area 4 and the highest reduction per treated area to the planning area 6. The higher overall performance in those units is explained by the high wildfire activity and the higher chance of SMPs encountering a fire in the future. However, we note that the treatment intensity (i.e., % treated area), spatial patterns (i.e., treatment unit size and clustering degree), and topographic position (e.g., water divides vs. valley bottoms) varied among and within the different planning areas. Treatment design is a significant factor affecting emission reductions that should be considered while interpreting the results [10,11].
Emission reduction maps showed very complex patterns and significant differences within the planning areas ( Figure 5). Not surprisingly, the treated stands were the highest emission reduction areas, but the performance varied widely within and among the planning areas. As expected, the most notable emission reductions concentrated within a 5 km buffer area of influence around treated SMPs. Overall, the effect reduction decreased with the increase in the distance to the SMPs. Nonetheless, this influence was broader in highly fire-affected regions where the effect resulted very substantially (>0.01 CO 2 T yr −1 ) far away from treated units ( Figure 5). Planning area 3 is a good example where massive catastrophic fires can potentially spread as much as 10 km ( Figure 1B). Table 3. Summary table with landscape-scale expected carbon dioxide emissions from wildfires and emission reductions after treatments. The emission reductions represent landscape-scale results accounting for the fuel treatment effects on treated SMPs (Figure 2) and neighboring lands ( Figure 5).

Planning Areas (Code)
Area ( (planning area 6) and 444 CO2 T per year (planning area 3), which resulted in a reduction between 11 (planning area 6) and 35% (planning area 2) compared to the non-managed conditions (Table 3). Interestingly, a high percentage reduction was not associated with a high decrease in CO2 T emission per year (e.g., see planning area 2). In planning area 2, most wildfire activity concentrated in the southern portion, where we concentrated many of the treatments, thus explaining the high reduction when presented as a percentage. On the other hand, the highest reduction per planning area extent corresponded to the planning area 4 and the highest reduction per treated area to the planning area 6. The higher overall performance in those units is explained by the high wildfire activity and the higher chance of SMPs encountering a fire in the future. However, we note that the treatment intensity (i.e., % treated area), spatial patterns (i.e., treatment unit size and clustering degree), and topographic position (e.g., water divides vs. valley bottoms) varied among and within the different planning areas. Treatment design is a significant factor affecting emission reductions that should be considered while interpreting the results [10,11]. Emission reduction maps showed very complex patterns and significant differences within the planning areas ( Figure 5). Not surprisingly, the treated stands were the highest emission reduction areas, but the performance varied widely within and among the planning areas. As expected, the most notable emission reductions concentrated within a 5 km buffer area of influence around treated SMPs. Overall, the effect reduction decreased with the increase in the distance to the SMPs. Nonetheless, this influence was broader in highly fire-affected regions where the effect resulted very substantially (>0.01 CO2 T yr −1 ) far away from treated units ( Figure 5). Planning area 3 is a good example where massive catastrophic fires can potentially spread as much as 10 km ( Figure 1B). The shadow effect in emission reductions obtained in the study areas was directly associated with the lower burn probability of the managed scenario ( Figure 6). Specifically, the footprint effect expanded on the opposite side to the large fire arrival border. For instance, the SMPs disrupting southern wind-driven wildfires exhibited a reduction in the treatment's northern side. Previous studies also found this reduced burn probability shades close to treatments due to smaller fire footprints obtained in modeled managed scenarios [69,75]. Nonetheless, very high aBP reductions in large, treated patches did not necessarily produce a significant shade effect in the adjacent lands if the polygons represent a fire sink area. We found this fire-sink effect in the large SMPs of central planning area 5, where the aBP footprint reduction in neighboring regions ( Figure 6) was minimal compared to the treated areas and emission reduction ( Figure 5). These were the treatment blocks border. For instance, the SMPs disrupting southern wind-driven wildfires exhibited a reduction in the treatment's northern side. Previous studies also found this reduced burn probability shades close to treatments due to smaller fire footprints obtained in modeled managed scenarios [69,75]. Nonetheless, very high aBP reductions in large, treated patches did not necessarily produce a significant shade effect in the adjacent lands if the polygons represent a fire sink area. We found this fire-sink effect in the large SMPs of central planning area 5, where the aBP footprint reduction in neighboring regions ( Figure  6) was minimal compared to the treated areas and emission reduction ( Figure 5). These were the treatment blocks with a low wildfire occurrence encircled by fire source communication infrastructure and urban development areas (Figure 3). The treatment unit level emission reductions showed variable results within and among the different planning areas ( Figure 7A). They revealed that SMPs were not necessarily allocated in the high-emission hot spots (Figure 4). The bulk of the SMPs attained reduction values between 0.01 and 1 T CO 2 per year on average, and only less than 5% of the treated areas was above 0.08 T CO 2 ha −1 per year. The top SMPs in the planning areas 1, 4, and 3 showed reductions greater than 1 T CO 2 per year, but this high performance shifted drastically when we analyzed the reduction per treated area. In planning area 4, for instance, the reduction per treated area varied more than five times between the treatment units with a similar total reduction. Other planning areas, such as the 1 and 2, presented a narrower difference of two times while addressing the performance per treated areas. The treatments implemented in planning area 2 were the worst among the different treatment units and rarely surpassed 0.01 T CO 2 ha −1 per year.
The widest variation in performance was obtained in planning area 5 ( Figure 7B). Here, the difference between the first and third quartiles was close to 0.05 T CO 2 per year and reflected how the average value alone is not enough to describe the performance results within planning areas. These results presented in box plots could assess the potential performance improvement for the best set of SMP treatment units after excluding those located in rare remote areas with a reduced potential to restrict fire spread and carbon emission (<0.01 CO 2 T ha −1 yr −1 ). Likewise, we found that some planning areas, such as 4 and 5 in particular, contained a significant number of observations above the third quartile ( Figure 7B) implemented in very high-emission areas (>0.05 CO 2 T ha −1 yr −1 ) where the project-level average reduction disguised a great allocation. These units were the high-performance strategic management points or fuel break, or treatment barriers implemented in the central part of the planning area 4 (Figure 3). planning areas 1, 4, and 3 showed reductions greater than 1 T CO2 per year, but this high performance shifted drastically when we analyzed the reduction per treated area. In planning area 4, for instance, the reduction per treated area varied more than five times between the treatment units with a similar total reduction. Other planning areas, such as the 1 and 2, presented a narrower difference of two times while addressing the performance per treated areas. The treatments implemented in planning area 2 were the worst among the different treatment units and rarely surpassed 0.01 T CO2 ha −1 per year. Figure 7. Summary statistics for the strategic management points. The performance is shown by treatment unit (A) and planning area (B) level. We obtained treatment unit-level wildfire occurrence and wildfire likelihood (C). The conditional fuel consumption is presented by planning area (D). Treatment implementation (Figure 2) operational constraints in terms of the average slope and distance to rode are presented by treatment unit (E). Finally, we show the burn probability reduction in treated SMP ( Figure 6) by planning area (F). In box plots, the boxes indicate the first/third quartiles, the whiskers indicate 10th and 90th percentiles, the black line is the median, and the dots indicate values below the 10th percentile or above the 90th percentile. The colors refer to planning areas, as shown in boxplots.
Although an increasing wildfire occurrence (Figure 3) correlated with a higher overall wildfire likelihood, high ignition probability values did not necessarily connote a high burn probability ( Figure 7C). For instance, many SMPs in planning area 4 with ignition probability values of about 60% sowed a very variable burn probability, even up to three times. These complex patterns were also found in other Mediterranean areas [76,77]. Wildfire burn probability was a significant factor explaining expected emissions because conditional fuel consumption ( Figure 7D) alone would have provided a misleading result.
While planning areas 1 and 4 presented similar fuel consumption results ( Figure 7D), expected emissions were much higher in planning area 4 ( Figure 7B). Most treatment units were located close to roads (<250 m) on gentle slopes (<30%) ( Figure 7E). Overall, burn probability reductions were higher than 80%, except in planning area 6, where the size of the treatment units was tiny, and large fires easily surpassed the SMPs ( Figure 7F).

Financing Fuel Treatment Cost with Carbon Credits
Forest management provided a significant gross benefit or revenue from preventing wildfire carbon emissions in some planning areas (Table 4). Specifically, the annual income was between 645 and 5853 € for the carbon emission reductions estimated in this study (Figure 5), assuming a €13.18 per CO 2 T market price [73]. While the lowest attained revenue normalized per planning area extent was 0.37 € ha −1 in planning area 2, the highest value resulted in 4.70 € per ha in planning area 4. Despite the low total revenue of the planning area 6 (643 € yr −1 ), the attainment per treated area at the landscape scale was above the average (2.35 € ha −1 ). Since the treatment allocation showed highly variable stand-level results (Figure 7), we calculated the percentage treated area in high emission hot spots (>0.01 CO 2 T ha −1 yr −1 ). This was an alternative performance metric to determine the SMP area with a significant contribution in reducing emissions. In contrast with the good overall allocation of treatments in planning areas 5 and 6 (>90%), in some planning areas, less than half of the treated area would be considered attractive for carbon credit investment. By contrast, a significant revenue emerged when we computed most carbon credits to a much lower treated area (e.g., 3245 € yr −1 to just 29% of the SMP area in planning area 4). Table 4. Revenue from a potential carbon credit market in the different planning areas. We assumed a fuel treatment effective duration between 8 to 14 years. The wildfire managers implementing fuel reduction programs in the study area provided the treatment cost ( Table 2). The SMP area over high emission hotspots (>0.01 CO 2 T ha −1 yr −1 ) was computed as a performance metric for the fuel reduction programs. −1 )   1  628  64  164  26  602  74,972  2  1980  0  1099  129  1851  278,880  3  3775  74  1133  711  3064  471,821  4  691  29  24  310  382  76,698  5  3171  96  1657  50  3121  436,342  6  274  93  104  1  The total treatment cost by planning areas was between 34,770 and 471,821 € per year. The areas requiring a thinning to a greater extent presented the highest cost. These SMPs were predominant in planning areas 2 and 5, where thinning was implemented in more than 50% of the area, and the average cost was above 135 € ha −1 per year. On the other hand, the most economical treatments were located in planning area 4 (110 € ha −1 yr −1 ) despite requiring the more expensive mechanical mastication instead of prescribed fire in 45% of the area. The carbon credit revenue's potential contribution in financing the fuel treatment cost was between 0.3 in planning area 2 and 4.2% in planning area 4. Considering only the SMP area allocated in carbon emission hotspots (Table 4), the contribution increased substantially. It would potentially cover up to 14.5% of the cost in the best case (i.e., planning area 4). Conversely, the potential revenue from carbon credits in planning areas 2 and 5 was below 1% in any case. While we expected poor results in planning area 2, the contribution in planning area 5 was exceedingly low for the observed high fire activity (i.e., annually burned area). We presume that this was due to the design of exceeding large treatment blocks in planning area 5 (Figure 2), where much narrower barriers or fuel breaks perpendicular to dominant winds, such as those implemented in planning area 4, would have obtained the same effect for less than half of the cost.

Discussion and Conclusions
Despite the European Union's accession to the Paris Agreement [78] and commitment to drastically cut CO 2 emissions, severe wildfires continue to burn vast areas and pose a significant threat to the Mediterranean forests [4,79,80]. As a result, forest management works developing wildfire risk reduction plans have substantially increased in recent years [9,81,82], but the emission reduction effects are still largely unknown. However, sustainable forest management support for preserving carbon pools in fire-prone forest ecosystems lacks specific financing lines. To our knowledge, this work is the first study conducted in Mediterranean areas estimating prevented CO 2 emissions in ongoing risk reduction plans through modeling. Precisely, we assessed the emission reduction effect of strategic management points (SMPs) and revenue from carbon credits in a wide range of fireprone Mediterranean landscapes. This study provides a valuable baseline for developing a carbon credit market intended to economically compensate small forest landowners for preserving fire resilient cultural landscapes in fire-prone southern European regions.
The fuel consumption quantification from large fires is essential to calculate potential carbon credits in Mediterranean landscapes. Large fires caused substantial losses in the past decades in Catalonia and will likely continue to burn vast portions due to limited management and young forest expansion in marginal agricultural lands [30,41]. Therefore, landscape planning efforts derived from stand-level conditional effects would largely ignore expected CO 2 emissions from large fires ignited elsewhere. In other words, the treatment units or stands present near-impossible odds of getting burned by a fire ignited within the same forest. Indeed, humans ignite most fires close to developed sites and communication infrastructure, and wildfires then hit forest lands after traveling long distances [36,83]. On the other hand, CO 2 emission estimates at regional scales overlook the complex patterns of the burned areas [84,85]. In contrast to most previous works in the Mediterranean region, we implemented stochastic fire simulation to replicate historic fire footprint distributions across the different planning areas and model the annual burn probability at high resolution (40 m). We accounted for the potential sources of uncertainty and variability in model inputs (i.e., ignition probabilities, extreme fire-weather conditions, and fuel loads in forest types) by modeling thousands of iterations (10,000 years) and summarizing the wildfire likelihood pixel-level results to an overall annual burn probability.
We identified the areas where ongoing fuel reduction programs are most likely to generate long-term carbon benefits. The results demonstrated how fuel management programs implemented in fire-prone areas reduced expected carbon emissions on treated stands and the neighboring forest lands. However, many SMPs located in remote regions would hardly encounter a wildfire and presented poor performance results. Thus, we strongly suggest prioritizing SMP treatments in high emission spots and postponing or excluding the fuel reductions in remote areas. Focusing on high emission areas would reduce the fuel treatment program cost drastically. In addition, our high-resolution expected carbon emission maps may help forest managers determine which suitable sites meet these conditions (Figure 4). The SMPs that we tested in this study were designed based on the expert criteria [20,22]. These SMPs were specifically designed to reduce large-fire potential, increase firefighting contention capacity, and did not explicitly address CO 2 emission reduc-tions [21]. Some previous studies estimated encounter rates between fire perimeters and treated areas to evaluate the fuel treatment effectiveness [86], but we decided to implement a high-resolution quantitative assessment based on the burn probability to measure CO 2 emission reductions [69]. Likewise, the burn probability has been widely used in previous works to prioritize and test fuel treatment effects at the landscape scale and select the most convenient design [87][88][89].
We conducted this analysis as a preliminary step to estimate spatially explicit results for surface fires. Other carbon emissions sources, including emissions during prescribed fires, crown fires, and the use of machinery, were not considered in our analysis. Previous studies assessing the prescribed fire combustion emissions reported a wide range of values between 2 and 10 CO 2 T ha −1 depending on the forest systems and treatment intensities [90,91]. Conducting prescribed fires in early spring would represent the year's proper timing for a light burn, preventing duff and heavy log consumption due to a higher fuel moisture content [92]. Considering an eight-year fuel treatment rotation interval, these light burns would represent about 0.375 T CO 2 ha −1 per year. Concerning crown fires, these only consume the thin branches and leaves (<6 mm) and modeling these emissions in high resolution at large scales would result in an exceedingly complex calculation. Nonetheless, high-resolution LiDAR and remote sensing-derived estimates would allow for computing the crown fire biomass in the following studies [93,94]. And lastly, we understand that chainsaws and mastication equipment also represent a CO 2 emission. Still, the manually implemented prescribed fire was the dominant treatment type (88% of the area vs. mechanical mastication in 12% of the treated area). We dismissed the use of machinery to extract thinning due to the problematic access to SMPs and the lack of commercial interest. In all, frequent but light surface burns in high emission spots, implemented in spring under high moisture content conditions for the soil, would deliver the most increased net carbon sequestration at the landscape scale [51].
We can conclude from these results that few landowners would be suitable to receive economic compensation for treating SMPs. Furthermore, the potential financial revenue from carbon credits was only substantial in extremely high CO 2 emission areas, such as the central fuel break barrier of the planning area LU4 (Figure 2). At best, the revenue from carbon credits could fund up to 14% of the treatment cost if we focus on SMPs with expected emissions above 0.01 CO 2 T ha −1 per year and exclude all the rest. Nonetheless, we need to emphasize that treatment costs will decrease by a factor of three after the first intervention because the following maintenance treatments are much easier to implement [52,54]. Moreover, we noted that results require careful consideration due to the wide range of treated intensities (i.e., % treated areas) and SMP size and shapes, which may affect the modeling outcomes. Likewise, treated SMP in large fire areas interact closely, and individualizing the contribution to the total emission reduction is highly challenging. Indeed, each treatment would require a separate modeling analysis to precisely determine reduced emission at the landscape scale.
We believe that protecting carbon stocks on forest systems may become a significant management objective in future projects as emission neutrality and the promotion of a bio-based economy gain momentum in the European Union. Nonetheless, future efforts should be oriented toward determining the tipping points between carbon benefits and losses from a broader range of forest management actions, which should provide better and more accurate information to compute realistic compensations [95]. However, optimal solutions may compete with other existing objectives (e.g., reducing carbon emissions vs. reducing wildfire transmission to communities), but the trade-off analysis allows assessing treatment co-location opportunities (i.e., forest stands where treatments can meet multiple goals) on vast landscapes [96,97]. Ultimately, treatment prioritization works would provide a valuable set of treatment solutions that would require implementing a risk assessment framework to evaluate cost-efficiency.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Dominant Vegetations Types and Historic Fire Activity
Irrigated agricultural lands, mosaics of low shrublands, and herbaceous xerophytic vegetation cover the central depression below 450 m. The higher altitude and coarse reliefs to the north confine cultivated plots to valley bottoms, with forested areas dominated by Mediterranean oaks and low shrublands on slopes (Lavandula angustifolia Mill., Rosmarinus officinalis L. and Quercus coccifera L.). These shrublands and forests are gradually replaced by tall-shrubland species (Buxus sempervirens L. and Juniperus communis L.), mid-mountain oak (Quercus pubescens Willd.), and conifer species (Pinus nigra Arn. and Pinus sylvestris L.) first on north-facing slopes, and then on higher elevations (Pinus uncinata Ram.). Mosaics of rocky outcrops, low shrublands (Genista balansae Boiss.), and pastures cover the high mountain tops above 2400 m. The Mediterranean maquis (Pistacia lentiscus L. and Arbutus unedo L.) appear in combinations with densely regenerated young Aleppo on the pre-littoral mountain ranges pine cohorts (Pinus halepensis Mill.). Silicicolous shrublands (Cistus ssp. and Erica ssp.) are found in coastal lowlands, sometimes with stone pine (Pinus pinea L.). The cork oak (Quercus suber L.) forests are limited to the northeastern lowlands. Table A1. Summary table with main features of the planning areas in the study and observed wildfire activity from 1983 to 2015. See planning area locations in Figure 1. We considered 100 ha as the large fire (LF) threshold. The wildfire activity data was summarized collectively for adjacent planning areas. The wildfire season was determined from the cumulative burned area ( Figure A1).

Appendix C. Dead Fuel Loads in Dominant Forest Types
Required dry biomass data for assessing the FC was obtained from field sampling campaigns conducted on the dominant forest typologies (Table 1; Figure A2) and complemented with fuel loading data from standard fuel models in herbaceous type models and live woody components [46]. We used different methods to estimate dead fuel loadings for the different fuel fractions. The field sampling combined litter and duff extractions on 30 cm radius circular plots (n = 16 plots per forest structure type) and 1, 10, and 100 h fuel extractions on 2 m square plots (n = 3 per forest structure type). All these fuel samples were then oven-dried and added to estimate the dry weight biomass per ha. To calculate the amount of 1000 h dead biomass, we conducted a set of 30 m sampling strips (n = 3 per forest structure type) where we measured the diameter of the fallen logs on the center and the length within a 1 m width. Then, we used species-specific allometric equations [98] to estimate the total amount of biomass from the volume of the logs, i.e., 1000 h fuel category. and 100 h fuel extractions on 2 m square plots (n = 3 per forest structure type). All these fuel samples were then oven-dried and added to estimate the dry weight biomass per ha. To calculate the amount of 1000 h dead biomass, we conducted a set of 30 m sampling strips (n = 3 per forest structure type) where we measured the diameter of the fallen logs on the center and the length within a 1 m width. Then, we used species-specific allometric equations [98] to estimate the total amount of biomass from the volume of the logs, i.e., 1000 h fuel category.  Table 1. Regular stand ages included low polestage (LP), high pole-stage (HP), and timber stage (T).

Appendix D. Methodological Flowchart
The cost-efficiency process was conducted in four main steps ( Figure A4). We first assessed expected carbon emissions as the annual burn probability times stand-level conditional carbon emissions for current conditions. Then, we implemented the landscape fuel reduction treatments in SMPs. Next, we predicted the annual burn probability reduction for the treated LCP to assess the expected carbon emission for the managed scenario. Finally, we calculated the carbon emission reduction for every Euro invested in treatments. The carbon credits of reduced carbon emissions were considered as revenue in managed scenarios.  Table 1. Regular stand ages included low pole-stage (LP), high pole-stage (HP), and timber stage (T).

Appendix D. Methodological Flowchart
The cost-efficiency process was conducted in four main steps ( Figure A4). We first assessed expected carbon emissions as the annual burn probability times stand-level conditional carbon emissions for current conditions. Then, we implemented the landscape fuel reduction treatments in SMPs. Next, we predicted the annual burn probability reduction for the treated LCP to assess the expected carbon emission for the managed scenario. Finally, we calculated the carbon emission reduction for every Euro invested in treatments. The carbon credits of reduced carbon emissions were considered as revenue in managed scenarios.