Short-Term Recovery of the Aboveground Carbon Stock in Iberian Shrublands at the Extremes of an Environmental Gradient and as a Function of Burn Severity

Abstract

The degree to which burn severity influences the recovery of aboveground carbon density (ACD) of live pools in shrublands remains unclear. Multitemporal LiDAR data was used to evaluate ACD recovery three years after fire in shrubland ecosystems as a function of burn severity immediately after fire across an environmental and productivity gradient in the western Mediterranean Basin. Two large mixed-severity wildfires were assessed: an Atlantic site, dominated by resprouter shrubs and located at the most productive extreme of the gradient, and a Mediterranean site, dominated by obligate seeders and located at the less productive extreme. Initial assessment of burn severity was performed using the differenced Normalized Burn Ratio index computed from Landsat imagery. Thresholds for low and high burn severity categories were established using the Composite Burn Index (CBI). LiDAR canopy metrics were calibrated with field measurements of mean shrub height and cover at plot level in a post-fire situation. Pre-fire and post-fire ACD estimates, and their ratio (ACDr) to calculate carbon stock recovery, were computed from the predictions of LiDAR grid metrics at landscape level using shrubland allometric relationships. Overall, ACDr decreased both with high burn severity and low productivity, although the burn severity impact was not homogeneous within the gradient. In the Atlantic site, ACDr was similar under low and high burn severity, whereas it decreased with burn severity in the Mediterranean site. These results suggest that carbon cycling models could be biased by not accounting for both fire severity and species composition of shrublands under different environmental conditions.

1. Introduction

Shrubland ecosystems store about half of the global terrestrial carbon [1,2] and thus their potential feedback with the atmosphere, associated with changes in carbon stocks, has significant implications for global carbon cycling [3]. Much information exists about the estimation of carbon stocks in forest ecosystems in the Mediterranean Basin [3]. However, the function of the shrublands as a carbon sink has hardly been investigated [4], even though these ecosystems are widely distributed across the Mediterranean Basin [5], and, according to [6], they constitute an important carbon sink that is usually underestimated by carbon storage assessments. Although the biomass per unit area might be relatively small in shrubland ecosystems, land degradation due to wildfires represents a significant change in live aboveground carbon pools in this region [7]. Wildfire disturbance causes direct pyrogenic carbon release of the aboveground and belowground pools to the atmosphere through combustion and modifies the distribution of live and dead aboveground carbon stocks and their associated fluxes through mortality and regeneration [8,9,10,11,12,13,14]. In shrubland ecosystems, indirect carbon emissions resulting from wildfire-induced mortality and subsequent decomposition are not expected to be determinant for net ecosystem carbon balance since net primary productivity (NPP) would offset in the short to medium-term the decomposition carbon flux, which is relatively small to shape carbon exchange at that temporal scale [15]. Thus, post-fire mortality and regeneration in shrublands are the processes playing a key role in carbon transfer between aboveground pools.

Current shifts in the Mediterranean Basin fire regime [16] as a consequence of land use change, anthropogenic climate warming and lack of adaptive management strategies to global change [17,18,19], have promoted the expansion of shrublands prone to large, recurrent and severe wildfires [20,21]. Fire regime shifts in shrubland ecosystems also involve an increase in landscape heterogeneity [22], leading to mixed-severity wildfires that could entail unpredictable changes in ecosystem carbon cycling [8,23]. Burn severity, denoted as the total amount of biomass consumed [24,25], is one of the most crucial factors determining ecological effects of fire [26,27]. Monitoring the burn severity impact on carbon stock recovery in shrubland ecosystems is thus crucial for predicting climate impacts associated with changes in the ecosystem carbon balance [28,29], and also for developing adaptive management strategies to promote ecosystem resilience [11] and reflect changes in terrestrial carbon inventories [6].

Aboveground carbon density (ACD; tn C ha⁻¹) has been traditionally estimated using field data at plot level by correlating vegetation structural characteristics to aboveground biomass (AGB; t ha⁻¹) using allometric models and carbon content factors specific to each species or ecosystem [30]. Although field inventories are highly reliable in providing realistic ACD estimates, this approach is labor-intensive and time-consuming when applied over large burned areas [31]. In this sense, remote sensing techniques, coupled with field calibration plots, offer nowadays the most feasible alternative for scaling ACD regionally [9,32,33] and creating spatially explicit estimates [29]. Passive optical sensors, such as those on-board Sentinel-2, Landsat and MODIS platforms, have been widely used to monitor AGB and ACD together with field plot data in forest ecosystems e.g., [34,35,36,37,38]. However, the reflectance signal captured by passive sensors is mostly determined by the top-of-canopy structural properties, especially under canopy closure conditions [21]. Therefore, the AGB/ACD estimation will be limited to secondary correlations with canopy traits, such as shadowing or moisture content [39]. Conversely, active remote sensors, such as airborne Light Detection and Ranging (LiDAR) and C-band/L-band synthetic aperture radar (SAR), enable accurate three-dimensional quantification of the vegetation canopy structure at both plot and stand levels [32,40,41,42]. This technology enables establishing direct relationships between vegetation structure and AGB/ACD using allometric models at the plot level e.g., [31,41,42,43,44,45,46,47].

Several studies have evaluated the role of burn severity on ACD variability in forest ecosystems e.g., [8,10,48,49]. Nevertheless, there is a research gap in the evaluation of ACD recovery of aboveground live pools in shrublands as a function of burn severity, particularly using remote sensing techniques. This scientific knowledge will be essential to assess the need for mitigation strategies and support adequate management decisions in burned areas at large spatial scales, focused on stabilizing aboveground carbon stocks in these ecosystems [50]. The primary objective of this research was to evaluate the ACD recovery in live pools of shrub species on the short term after fire (less than five years; [51]) as a function of initial burn severity, by means of an area-based approach with airborne LiDAR data validated with field plot surveys, within the perimeter of two large mixed-severity wildfires affecting shrubland ecosystems across an environmental gradient in the western Mediterranean Basin. We have focused exclusively on aboveground carbon stock recovery of living shrub biomass because of the following reasons: (i) aboveground carbon stocks of live vegetation pools (i.e., NPP) are the driving force of the terrestrial carbon cycling, a major determinant of carbon sinks [52] and an essential governor of ecological processes [53]; (ii) the relationship between carbon pools in necromass and remote sensing data in shrublands cannot be directly estimated with ecological sense operationally, as opposed to field carbon inventories [54,55]; (iii) NPP would largely exceed the carbon flux due to dead biomass decomposition in shrublands, which would be much more extended in time than the fire return interval in Mediterranean shrublands [15].

Specifically, we seek to address the following research questions: (i) do the products derived from LiDAR data provide enough accuracy for estimating shrubland vegetation structural variables as a proxy of ACD at plot level in post-fire environments? (ii) is the ACD recovery conditioned by burn severity on the short-term after fire? and, (iii) does the ACD recovery pattern change across the gradient?

2. Material and Methods

The methods comprised four steps that were accomplished independently in each study site (Figure 1): (i) burn severity mapping immediately after fire through a threshold-based classification of spectral indices derived from Landsat satellite data; (ii) LiDAR data acquisition and processing for calculating pre and post-fire shrub canopy metrics at plot level; (iii) acquisition of field plot data for validating LiDAR plot metrics; (iv) estimation of carbon stock recovery three years after fire through ACD allometric relationships applied to LiDAR metrics at landscape level; and (v) data analysis.

2.1. Study Sites

Two wildfires were selected at the extremes of an Atlantic-Mediterranean climatic and productivity gradient on a west-east axis in the Iberian Peninsula (Figure 2).

The Atlantic site is located within the perimeter of a mixed-severity wildfire that burned 2523 ha in September 2013 in the northwestern Iberian Peninsula (Figure 2A). The site has a complex orography and elevation ranges between 0 and 628 m above sea level (ASL). The climate of the region is Atlantic, with a 50-year period mean annual temperature and precipitation of 13 °C and 1655 mm, respectively. Summers are temperate with no drought period. Soils are acidic and mostly classified as Umbrisols. The pre-fire landscape was a mosaic of Eucalyptus globulus Labill. and Pinus pinaster Ait. stands, gorse shrublands dominated by Ulex europaeus L. and Ulex gallii Planch., and heathlands dominated by Erica australis L., Erica umbellata Loefl., Pterospartum tridentatum (L.) Willk. and Halimium lasianthum subsp. alyssoides (Lam.) (Figure 3A). The post-fire landscape (three years after the wildfire) consisted of tree regeneration stands of Pinus pinaster in a seedling and sapling growth stage, resprouting stands of Eucalyptus globulus and shrublands dominated by the same pre-fire species. The area within the perimeter of the 2013 wildfire had partially burned in 2005, 2004, 2001, 2000, 1999 and 1995 [56].

The Mediterranean site is located within the perimeter of a mixed-severity wildfire that occurred in June 2012 in the eastern Iberian Peninsula (Figure 2B), which burned 29,752 ha of Pinus halepensis Mill. stands, as well as a large extension of garrigue (limestone shrubland with presence or dominance of Quercus coccifera L., Pistacia lentiscus L., Rosmarinus officinalis L., Cistus monspeliensis L. and Ulex parviflorus Pourr.) (Figure 3B). Terrain presents wide valleys and prominent crests, the elevation ranging between 114 and 995 m ASL. Climate is Mediterranean, with a mean annual temperature of 16 °C and a mean annual precipitation of 582 mm. Summers are hot and dry, with three months of drought. Soils are basic and classified as Haplic Calcisol and Calcari-Lithic Leptosol. The post-fire landscape was dominated by the above-mentioned pre-fire shrublands and Pinus halepensis regeneration stands. In addition to the 2012 wildfire, this site had been burned in 1978, 1991 and 1994 [56]. Within the fire perimeter, a study framework of 2917 ha representative of the mixed-severity conditions (Figure 2B) was selected to make the study surface comparable to that of the Atlantic site.

2.2. Burn Severity Mapping

Remote sensing data to estimate burn severity (Landsat Collection 1 Level-1 imagery) were retrieved from the USGS Earth Explorer data portal (http://earthexplorer.usgs.gov/ accessed on 11 May 2021). Landsat 8 Operational Land Imager (OLI) scenes for both 30th August 2013 at 11:21:33 UTC (pre-fire) and 15 September 2013 at 11:21:28 UTC (post-fire) were used in the Atlantic site (Path 205/Row 30). Landsat 7 Enhanced Thematic Mapper Plus (ETM+) scenes for both 23th August 2011 at 10:37:00 UTC (pre-fire) and 25 August 2012 at 10:38:50 UTC (post-fire) were used in the Mediterranean site (Path 199/Row 33). The optical bands of Landsat 8 OLI and Landsat 7 ETM+ scenes were atmospherically and topographically corrected to surface reflectance with the ATCOR algorithm [57], obtaining spectral products at 30 m of spatial resolution. Input data for ATCOR parameters (aerosol model, atmospheric model and visibility) were obtained from the MODIS water vapor product (MOD05) and meteorological data from the National Oceanic and Atmospheric Administration (NOAA) and the State Meteorology Agency of Spain (AEMET) [21]. The failure of Scan Line Corrector (SLC) occurred in Landsat 7 ETM+ imagery was corrected using a Delaunay interpolation method [26]. Finally, surface reflectance of each Landsat 7 ETM+ band was transformed to comparable Landsat 8 OLI surface reflectance with the functions provided by [58]. This operation minimizes significant reflectance differences between Landsat 7 ETM+ and Landsat 8 OLI, especially in the near-infrared (NIR) and the shortwave infrared (SWIR) regions due to the dissimilarities in the spectral response functions between both sensors [58].

Despite the influence of soil background signal and transferability issues that may hinder the performance of Normalized Burn Ratio (NBR)-based indices in burn severity initial assessments [21,59,60], threshold-based classification of the differenced NBR (dNBR; [61]) index is the most widely accepted approach for this purpose [62,63]. In fact, the dNBR is the primary spectral index within the Rapid Damage Assessment (RDA) module of European Forest Fire Information System (EFFIS) and the Monitoring Trends in Burn Severity (MTBS) project (together with RdNBR; [64]) in the United States. Therefore, the remote sensing-based estimation of the total amount of biomass consumed immediately after fire was conducted through the dNBR, using Landsat 7 ETM+ (Equations (1) and (3) and Landsat 8 OLI Equations (2) and (3) pre-fire and post-fire surface reflectance data from the NIR and SWIR regions.

$${\mathrm{NBR}}_{\mathrm{ETM}+}=\left(\mathrm{Band} 4-\mathrm{Band} 7\right) / \left(\mathrm{Band} 4+\mathrm{Band} 7\right)$$

(1)

$${\mathrm{NBR}}_{\mathrm{OLI}}=\left(\mathrm{Band} 5-\mathrm{Band} 7\right) / \left(\mathrm{Band} 5+\mathrm{Band} 7\right)$$

(2)

$$\mathrm{dNBR}=1000\left({\mathrm{NBR}}_{\mathrm{pre}}-{\mathrm{NBR}}_{\mathrm{post}}\right)-\mathrm{offset}$$

(3)

The offset term in Equation (3) is the average dNBR value from pixels in homogeneous and unchanged areas outside the fire perimeter [61,65] to make the index comparable between wildfires [61]. The dNBR thresholds were identified on the basis of field initial assessments (two months after fire in early fall) of the Composite Burn Index (CBI; [66]) in Atlantic and Mediterranean ecosystems [67,68] similar to those of the present study. Two field burn severity categories were established based on one of the CBI thresholds proposed by [69]: low (CBI ≤ 2.25) and high (CBI > 2.25). Using these CBI thresholds and the linear models proposed by [67] and [68], with coefficients of determination (R²) of 0.67 and 0.79, respectively, two dNBR burn severity categories (low and high) were established in the Atlantic and Mediterranean sites (Figure 2). Although we used external field burn severity data, the CBI approach provides a consistent and interpretable method for broad-scale comparisons of burn severity across time and regions with similar vegetation types and environmental conditions [70,71,72].

2.3. LiDAR Data Acquisition and Processing

LiDAR data were collected by the Spanish National Plan for Aerial Orthophotography (PNOA). Pre-fire data were collected for the Atlantic site on 5 February 2011 and for the Mediterranean site between 17 August 2009 and 3 September 2009, using a RIEGL LMS-Q680i sensor aboard a fixed-wing aircraft. Post-fire LiDAR data were collected for the Atlantic site between 1 and 19 August 2015, using a Leica ALS60 sensor, and for the Mediterranean site on 9 November 2015, using a RIEGL LMS-Q780 sensor. There were no fire events between pre-fire and post-fire LIDAR data acquisitions apart from the target wildfires. The discrete-return LiDAR sensors captured up to four returns per pulse. The mean point cloud density was 0.73 m⁻² (pulse spacing of 1.17 m) and 0.79 m⁻² (pulse spacing of 1.13 m) for the Atlantic and Mediterranean sites, respectively. The overall vertical accuracy (RMSE_Z) reported by the PNOA was lower than 0.2 m.

The point cloud classification into ground and non-ground returns was performed using the multiscale curvature classification algorithm [73] implemented in MCC-LIDAR 2.1 software. This algorithm has better performance compared to other open-source filtering algorithms in shrubland ecosystems and burned areas [74]. LAStools software (rapidlasso GmbH, Germany) was used to create a digital terrain model (DTM) with a grid size of 2 m, using a triangulated irregular network (TIN) interpolated from the filtered ground returns of the LiDAR point cloud. Then, LiDAR returns were normalized to heights above the ground using the underlying DTM and point clouds were clipped to the spatial extent of the field plots (30 m × 30 m). Several LiDAR metrics with ecological sense and closely related with canopy height and cover (input variables in the considered allometric equations for estimating shrub AGB/ACD) were computed from the height-normalized returns at plot level using US Forest Service’s FUSION software package version 3.80 [75], instead of the traditional approach of computing a wide battery of metrics and perform model selection procedures, which complicates the development of robust and generalizable models [76]. Additionally, canopy metrics were computed at landscape level across the study sites with a grid size of 30 m. A minimum height and cover threshold of 0.2 m was implemented to remove the influence of laying woody debris and other ground features such as rocks on the computed metrics. The metrics comprised: (i) 95th percentile of LiDAR 1st returns and (ii) average height of LiDAR 1st returns, both as a representative measure of canopy mean height in the plot [77,78,79]; and, (iii) LiDAR canopy cover, computed as the proportion of 1st returns above the cover threshold [80,81].

2.4. Field Plot Data and Relationship with LiDAR Plot-Level Metrics

Field data were collected in square plots of 30 m × 30 m within the wildfire perimeter of the Atlantic (40 plots) and Mediterranean (30 plots) sites. The plots were established in the summer of 2015, following a stratified random sampling design, using burn severity categories (low and high) as strata. Their center was located using a sub-meter accuracy GNSS receiver in post-processing mode. Each plot contained a cluster of four subplots of 2 m × 2m, whose center was located 6.5 m away from the center of the plot at azimuths of 45°, 135°, 225° and 315°. Within the subplots, we measured the mean height and cover as the vertical projected area occupied by the shrub stratum using a visual estimation method in steps of 5% [82,83,84]. The presence of dead standing biomass was negligible. Vegetation parameters quantified within the subplots were averaged at the plot level to be correlated with post-fire LiDAR data.

Univariate linear regression models were calibrated between mean shrub height measured in the field (dependent variable) and both 95th percentile 1st return heights and average height of 1st returns (independent variables) to identify the most contributing LiDAR plot-level metric through the coefficient of determination (R²) and the root-mean-squared error (RMSE). The same procedure was used to evaluate the relationship between field-estimated shrub cover and LiDAR cover metric at plot level.

2.5. Estimation of Carbon Stock Recovery and Data Analysis

Pre-fire areas dominated by trees were discarded from subsequent analyses since burn severity was assessed through passive optical remote sensing data and, therefore, the understory reflectance signal corresponding to the fire impact in the shrub stratum is more or less occluded depending on canopy openness [21,85,86]. Pre-fire LiDAR data were used to compute a canopy height model (CHM) with a grid size of 2 m generated from the TIN interpolation of the highest LiDAR returns and the subsequent rasterization onto a grid using LAStools software. The tree stand mask was computed from the grid cells with a height greater than 3 m, based on field expert knowledge. Additionally, grid cells lower than 0.2 m were masked to exclude other non-vegetation ground features. CHM thresholding has been recognized as a reliable approach for masking non-interest surfaces [47]. Mask accuracy was measured through a confusion matrix based on one thousand pixels randomly sampled within the Atlantic and Mediterranean sites, using as reference dataset several pre-fire PNOA orthophotographs at a spatial resolution of 0.5 m. True positive (TP), true negative (TN), false positive (FP) and false negative (FN), reflect accurate shrub pixel extraction, accurate rejection of non-shrub pixel, inaccurate shrub pixel extraction and inaccurate shrub pixel rejection [87]. From the confusion matrix, we computed the overall accuracy (ACC) and the Kappa index (κ), as well as commission (CE) and omission (OE) errors. The mask was expanded 15 m to discard regeneration stands close to the burned tree legacies in post-fire situation.

Pre-fire and post-fire ACD maps were computed within the non-masked areas from the LiDAR grid metrics at landscape level selected as proxies of mean shrub height and cover. We used the AGB allometric equations proposed by [88] for gorse shrublands (Equation (4)) and heathlands (Equation (5)) in the Atlantic site, and garrigue shrublands (Equation (6)) in the Mediterranean site, as well as the mean carbon content (%) for the dominant species of each formation [88].

$${\mathrm{ACD}}_{\mathrm{gorse}}={\displaystyle 0.4815}*\left[1.4902{*\mathrm{H}}_{\mathrm{m}}{}^{0.9889}*\arcsin \left(\sqrt{{\left(\mathrm{CC}/100\right)}^{1.0276}}\right)\right]$$

(4)

$${\mathrm{ACD}}_{\mathrm{heath}}={\displaystyle 0.5106}*\left[2.0484*{\left({\mathrm{H}}_{\mathrm{m}}*\arcsin \left(\sqrt{\left(\mathrm{CC}/100\right)}\right)\right)}^{0.8614}\right]$$

(5)

$${\mathrm{ACD}}_{\mathrm{garrigue}}={\displaystyle 0.4943}*\left[4.7723*{\left({\mathrm{H}}_{\mathrm{m}}{}^{0.5918}*\arcsin \left(\sqrt{\left(\mathrm{CC}/100\right)}\right)\right)}^{1.2271}\right]$$

(6)

where ACD is the aboveground carbon density in tons per hectare (t ha⁻¹), H_m is the mean shrub height at plot level in decimeters (dm) and CC is the percentage of shrub canopy cover at plot level. Bold numbers indicate the percentage of carbon content. There is a different temporal mismatch between pre and post-fire LiDAR data acquisition and the fire date in both sites. We applied an annual biomass growing equation for each shrub formation (see [88] for more details) to match the ACD to one year before the fire and three years after fire and make the data comparable. The map of aboveground carbon stock recovery of living shrub biomass was computed as the ratio of post-fire to pre-fire ACD estimates (ACD_r).

In each study site, ACD_r was sampled in one thousand pixels within the non-masked areas in the pre-fire situation, following a stratified random sampling design, using burn severity categories (low and high) as strata. A minimum distance of 30 m between points was ensured. A two-way ANOVA (2w-ANOVA) was performed to evaluate the effect of burn severity and environmental conditions (i.e., site, Atlantic or Mediterranean), and their corresponding interactions on ACD_r. Subsequently, a Tukey’s HSD test was implemented to decompose the interaction by means of pairwise comparison of mean ACD_r as a function of burn severity within each site. Statistical significance was determined at p < 0.05.

3. Results

Post-fire shrub height measured in the field at plot level was best correlated to the mean height of LiDAR 1st returns metric, both in the Atlantic (R² = 0.69; RMSE = 0.20 m) and Mediterranean (R² = 0.74; RMSE = 0.15 m) sites (Figure 4). The relationship was statistically significant at the 0.001 level. Conversely, the 95th percentile of LiDAR 1st returns provided the worst results in both sites (R² = 0.62–0.69; RMSE = 0.15–0.26 m). Therefore, the mean height of LiDAR 1st returns was selected as a proxy of mean shrub height in pre and post-fire ACD estimations. LiDAR canopy cover metric was correlated (p < 0.001) with field-estimated shrub canopy cover at plot level (Figure 4) under Atlantic (R² = 0.65; RMSE = 8.21%) and Mediterranean (R² = 0.78; RMSE = 7.41%) conditions. However, LiDAR metrics derived from the height distribution of the returns slightly underestimated mean shrub height and shrub canopy cover in both study sites (Figure 4).

The overall accuracy of the land cover mask in the pre-fire condition was 85.10% and 91.90% for the Atlantic and Mediterranean sites, respectively. The Kappa index was 0.70 for the Atlantic site and 0.83 for the Mediterranean site. Commission errors in the delineation of shrubland ecosystems were similar in both sites (between 7–9%). Conversely, the omission error was significantly higher in the Atlantic (17.67%) than in the Mediterranean site (11.50%), which indicates slight underestimation of the surface occupied by shrublands in the former site with respect to the reference data (

Table 1. Confusion matrix and accuracy metrics (ACC—overall accuracy; κ—Kappa index; CE—commission error; OE—omission error) of the pre-fire land cover mask (shrubland and non-shrubland) in the Atlantic and Mediterranean sites.

Atlantic Site		PNOA Orthophotograph
		Shrubland	Non-Shrubland
Mask	Shrubland	480	46
	Non-shrubland	103	371
ACC = 85.10%; κ = 0.70; CE = 8.75%; OE = 17.67%
Mediterranean Site		PNOA Orthophotograph
		Shrubland	Non-Shrubland
Mask	Shrubland	377	32
	Non-shrubland	49	542
ACC = 91.90%; κ = 0.83; CE = 7.82%; OE = 11.50%

).

ACD_r was significantly conditioned by both burn severity (F = 33.56; p < 0.001) (Figure 5A) and environmental conditions, i.e., site (F = 114.92; p < 0.001) (Figure 5B). Although burn severity had an overall negative impact on ACD_r, this effect was not homogeneous within the environmental gradient, as revealed from the significant interaction between burn severity and site in the 2w-ANOVA (F = 56.75; p < 0.001) (Figure 5C). In the Atlantic site, ACD_r was similar under low and high burn severity situations (p = 0.658), whereas it was significantly lower with increased burn severity in the Mediterranean site (p < 0.001). Finally, high burn severity reduced ACD_r to a greater extent in the Mediterranean site as compared to the Atlantic site (p < 0.001). The effects of low burn severity were similar in both sites (p = 0.084).

These results are consistent with the predicted ACD_r at landscape scale computed from pre and post-fire LiDAR canopy metrics with a grid size of 30 m (Figure 6). In the Atlantic site, spatial patterns of ACD_r were similar in both low and high severity burned areas, with recovery values predominantly higher than 70%. In addition, about 20% of the surface in the Atlantic site have accumulated aboveground carbon stocks in the post-fire situation on the short-term (Figure 6A). In contrast, 50% of the surface in the Mediterranean site burned at high severity featured an ACD_r lower than 20%, whereas such conditions only represented less than 10% of the areas burned at low severity in this site (Figure 6B).

4. Discussion

Quantifying aboveground biomass and carbon stocks in shrubland ecosystems through remote sensing techniques at large spatial scales is crucial for understanding wildfire impacts on carbon storage and cycling of terrestrial ecosystems [44,89,90,91], particularly in the context of changing fire regimes in Mediterranean ecosystems [17,21]. The results of this study demonstrate that the use of multi-temporal LiDAR data is a reliable approach for identifying the recovery patterns of aboveground carbon stocks in shrubland ecosystems as a function of burn severity under different environmental conditions.

4.1. LiDAR Data as a Proxy for Shrubland Carbon Stocks in Post-Fire Environments

Plot-level metrics computed from low pulse density LiDAR data in post-fire scenarios under Atlantic and Mediterranean environmental conditions successfully accounted for shrubland structural variables (i.e., shrub height and canopy cover) with a strong linkage with shrub biomass [90,92]. Despite the low point cloud density (between 0.7–0.8 points m⁻² for both sites), the use of a plot-stand approach to estimate shrubland structure minimizes this concern [44], since the point density translates in this case to a mean of 675 points per 900 m² (30 m × 30 m) plot. This is a reasonable number of points for computing shrubland plot metrics assuming a stable height distribution [45], as evidenced by previous studies in complex Mediterranean ecosystems [46].

The mean height of LiDAR 1st returns and LiDAR canopy cover (computed as the proportion of 1st returns above 0.2 m) metrics have been evidenced as strong proxies of field-measured vegetation structure at plot level, both in forest and shrubland ecosystems [44,45,46,47,92,93]. The overall accuracy and error in estimating field shrub height in the Atlantic and Mediterranean sites (0.69–0.74 and 0.15–0.20 m, respectively) agreed with the results of previous research in sagebrush and grassland communities in North America [89,94], and shrublands of the western Mediterranean Basin [95,96]. Similar to height, the LiDAR data estimation of mean shrub cover was significantly related to field reference cover, with errors of 7–9% for both study sites. This is consistent with [92], who used multi-scale LiDAR data for assessing shrub biomass across drylands ecosystems in the United States. The authors compared the airborne LiDAR canopy cover metric using a cover threshold of 0.2 m with the reference cover product derived from terrestrial LiDAR, yielding an error of 10%. In similar ecosystems, [97] also found that shrub cover can be accurately predicted with errors lower than 10%, the canopy cover metric being ranked as one the most important in the models using LiDAR data alone.

The estimation of ACD_r as a function of burn severity could be affected by the observed underestimation of mean shrub height and cover (also noticeable in the land cover classification mask in pre-fire condition). However, this error was systematically distributed across the study sites and burn severity scenarios (0.1–0.2 m for mean shrub height and around 10% for mean shrub cover). Underestimation of shrub canopy structure has been widely reported in the literature e.g., [89,91,92,96,98,99] and may occur because of (i) field measurement related errors and the vertical accuracy of the sensor [100]; (ii) low probability of the laser pulses hitting the top of the canopy depending on its architecture [99]; (iii) laser pulses penetrating the canopy [90]; and, (iv) misclassification of ground and non-ground returns by the ground filtering algorithm to generate a DTM in areas dominated by low and sparse vegetation, as well as in steep slopes [78]. In this sense, the slightly worse results in the Atlantic site could be related to misclassification of ground returns because of higher terrain complexity compared to the Mediterranean site.

4.2. Impact of Burn Severity and Environmental Conditions on Carbon Stocks Recovery

Regenerating vegetation with strong post-fire responses, such as shrub and herbaceous species, rapidly recovers from the loss of carbon stocks induced by fire [8,101]. However, this study revealed that short-term ACD_r in shrubland ecosystems is strongly conditioned by both burn severity and environmental conditions.

Previous research has demonstrated that environmental conditions play a key role in post-fire recovery trajectories in the first growing seasons e.g., [51,102]. According to [103], post-disturbance regeneration is expected to be very fast in productive environments with high resource availability, following the resource-productivity model [104]. The conditions of the Atlantic site, located at the most productive extreme of the climatic gradient in the Iberian Peninsula [56], will be translated into dramatic increases of NPP in the early post-fire periods, which is consistent with the higher overall shrub ACD_r (i.e., regardless of burn severity) compared to the Mediterranean site, located at the less productive extreme of the gradient [56]. In addition, shrubs with post-fire resprouting regeneration strategy are favored as resource availability increases [103,105]. The absence of significant ACD_r differences in the Atlantic site under low and high burn severity can also be explained by the above-mentioned resource-productivity model and the dominance of resprouter shrubs. Several authors evidenced that, to a certain extent, resprouter shrubs are less vulnerable to burn severity than post-fire obligate seeders e.g., [21,106,107]. In fact, the resprouting response might only be restricted whether the belowground bud bank is injured by extremely severe fires or the fire-free intervals are insufficient for recovering the reserves invested in resprouting [108,109]. Likewise, [110] indicated that promoting resprouter shrubs might increase the resilience to fire, ecosystem productivity and, therefore, terrestrial carbon sinks. In contrast, the ACD_r of living shrubs in the Mediterranean site was significantly lower as burn severity increases. This behavior may be associated with the post-fire seedling recruitment strategy of the dominant shrub species (obligate seeders) in this site [111]. Although obligate seeder shrubs are more tolerant to water deficit than resprouters and their germination is stimulated by heat [103,112], their post-fire recovery is slower [113,114] and the seed bank could be soundly hindered by severe fires [21], noticeably influencing short-term NPP trends [8].

Despite the lack of research about the burn severity influence on ACD_r in shrublands, several studies assessed wildfire effects on the carbon pools of understory shrub communities, among other compartments, in forest ecosystems. [8] found a strong post-fire shrub response in the understory of ponderosa pine (Pinus ponderosa Douglas ex C. Lawson) and mixed-conifer forests in Oregon, reaching or exceeding pre-fire NPP levels four or five years after the fire regardless of burn severity. In mixed-conifer forests in California, [48] also evidenced that wildfire had little effect on the recovery of shrub carbon stocks. However, under harsh post-fire physical environments in fire-prone open shrublands [114], with marginal microclimate amelioration and usually less soil nutrient availability compared to forest ecosystems [115,116], the influence of burn severity and ecosystem productivity interaction on aboveground carbon pools has turned out to be remarkable. This implies that the effects of fire disturbance in carbon cycling models and environmental policies at regional scales could be systematically biased if the related analyses are not focused on all landscape components [8], including open shrublands prone to mixed-severity wildfires.

4.3. Uncertainties and Research Implications

Despite the relevant findings reported in this study, the approach featured several uncertainties related to remote sensing data products, field measurements and allometric relationships. First, ACD proxies estimated from LiDAR data were only validated with field data at plot level in a post-fire situation, as in other related multitemporal LiDAR studies e.g., [45,49]. This should not pose a severe concern since the pre-fire situation would be even more favorable for LiDAR remote sensing technology in terms of increased vegetation height and less sparse vegetation cover [89] than the post-fire scenario. Second, allometric relationships were extensively validated by [88] using field plot data from the same shrublands occurring in our study sites. ACD relationships computed from LiDAR metrics in this study were not field validated because of the unavailability of shrub biomass data. Therefore, the propagation of error sources from LiDAR data (mainly underestimation of shrub height and cover) and allometric equations to aboveground carbon uncertainty [117] cannot be measured. Third, the lag between pre and post-fire LiDAR surveys and fire date in both sites was addressed by means of biomass growing equations for each specific shrubland ecosystem [88]. In the Atlantic site, the predicted post-fire ACD at landscape scale (year 2015) was matched to three years after fire (year 2016) to be comparable with the predicted stocks in the Mediterranean site. This approach does not completely reflect actual post-fire conditions since vegetation growth rates in fire-prone shrub ecosystems are influenced by the fire regime [118,119]. However, the results obtained in the Atlantic site would not be biased by this behavior since ACD_r was comparable under low and high burn severity scenarios and the burn severity trends will remain the same.

The reported uncertainties are common to any LiDAR project dealing with multitemporal data from several sites [120] and do not undermine the main trends observed. Future efforts are needed for quantifying the above-mentioned sources of propagating errors [76,121,122] of remotely sensed ACD in shrubland ecosystems either through Monte Carlo simulations e.g., [123,124], recommended by [125], or analytical approximations e.g., [122,126]. Meanwhile, the minimization of uncertainties related with carbon stock estimates in shrublands should be addressed by means of remotely sensed data of higher quality [127], such as LiDAR data with higher pulse density in spite of its direct implications on acquisition costs at large spatial scales [76]. In this sense, LiDAR sensors on-board unmanned aerial vehicles (UAV-LiDAR) allows for acquiring higher point density than traditional airborne LiDAR [128,129]. Simultaneous acquisition of UAV multispectral imagery (UAV-MS) at very high spatial resolution also allows for improving the classification of ground/non-ground returns of UAV-LiDAR data in areas with low vegetation and rough topography [95,130] and discriminating dead standing and lying biomass [131,132]. Therefore, UAV consumer-grade technology may enhance ACD estimation in shrubland ecosystems and should be considered in future research.

Despite the limitations outlined in this section, the proposed approach and the results evidenced in this study could provide reliable datasets to understand the sensitivity in specific ecosystems of operational products providing global carbon estimates [120]. Additionally, these datasets will improve the refinement of carbon budget inventories, particularly those based on fixed land cover [133], as well as the parametrization of carbon cycle models [134] based on the critical role of the fire regime in the changes of carbon stocks in open shrublands prone to mixed severity wildfires.

5. Conclusions

The quantification of carbon stocks through remote sensing techniques is crucial for understanding wildfire impacts at large spatial scales in shrubland ecosystems. The proposed remote sensing methodology, based on multitemporal airborne LiDAR surveys together with field plot data, has proved its applicability to quantify shrubland structural attributes as a proxy of live vegetation ACD computed from allometric relationships. This study showed for the first time that the ACD_r in open shrublands after wildfire is strongly conditioned, in the short-term, both by the interaction of burn severity and environmental conditions or ecosystem productivity. Burn severity had a major impact on ACD_r under Mediterranean conditions, corresponding to less productive environments dominated by obligate seeder shrubs, whereas it was not determinant in most productive areas under Atlantic climatic conditions dominated by resprouter shrubs. In addition, ACD_r was significantly higher in the Atlantic than in the Mediterranean site. Improvements of the proposed approach for future research include quantifying propagating errors by means of uncertainty analyses and the use of LiDAR data of higher quality to minimize remote sensing uncertainties related to carbon stock estimations in shrubland ecosystems. Our results can increase the understanding of the sensitivity of operational carbon monitoring products in shrubland ecosystems, and refine carbon inventories and the parametrization of carbon terrestrial models considering the fire-regime impact.