Quantifying biomass consumption and carbon release from the California Rim fire by integrating airborne LiDAR and Landsat OLI data

Quantifying biomass consumption and carbon release from the California Rim fire by integrating airborne LiDAR and Landsat OLI data Mariano Garcia 1,2 , Sassan Saatchi 1 , Angeles Casas 3 , Alexander Koltunov 3 , Susan Ustin 3 , Carlos Ramirez 4 , Jorge Garcia-Gutierrez 5 , and Heiko Balzter 2,6 1 Jet Propulsion Laboratory, California Institute of Technology, Pasadena, California, USA, 2 Centre for Landscape and Climate Research, University of Leicester, Leicester, UK, 3 Center for Spatial Technologies and Remote Sensing, University of California, Davis, California, USA, 4 Region 5 Remote Sensing Lab, McClellan, USDA Forest Service, Vallejo, California, USA, 5 Department of Computer Science, University of Seville, Seville, Spain, 6 National Centre for Earth Observation, University of Leicester, Leicester, UK Abstract Quantifying biomass consumption and carbon release is critical to understanding the role of fires in the carbon cycle and air quality. We present a methodology to estimate the biomass consumed and the carbon released by the California Rim fire by integrating postfire airborne LiDAR and multitemporal Landsat Operational Land Imager (OLI) imagery. First, a support vector regression (SVR) model was trained to estimate the aboveground biomass (AGB) from LiDAR-derived metrics over the unburned area. The selected model estimated AGB with an R 2 of 0.82 and RMSE of 59.98 Mg/ha. Second, LiDAR-based biomass estimates were extrapolated to the entire area before and after the fire, using Landsat OLI reflectance bands, Normalized Difference Infrared Index, and the elevation derived from LiDAR data. The extrapolation was performed using SVR models that resulted in R 2 of 0.73 and 0.79 and RMSE of 87.18 (Mg/ha) and 75.43 (Mg/ha) for the postfire and prefire images, respectively. After removing bias from the AGB extrapolations using a linear relationship between estimated and observed values, we estimated the biomass consumption from postfire LiDAR and prefire Landsat maps to be 6.58 ± 0.03 Tg (10 12 g), which translate into 12.06 ± 0.06 Tg CO2 e released to the atmosphere, equivalent to the annual emissions of 2.57 million cars. 1. Introduction Forest fires are a natural process playing a critical role in the structure and functioning of many terrestrial ecosystems that are adapted to fires. They are a key element of forest composition and succession, resulting in stand thinning, removing understory vegetation, and regulating patterns of carbon (C) accumulation by promoting vertical stratification of the forest canopy [Chuvieco, 2008; French et al., 2003]. Whereas low to moderate severity fires may increase nutrients and promote vegetation renovation, high burn severity fires reduce organic matter and deteriorate soil structure and porosity [Certini, 2005; Neary et al., 1999]. Forest fires also impact the C cycle through the direct release of C into the atmosphere during the biomass burning [Andreae, 1991; French et al., 2003; van der Werf et al., 2006]. Therefore, quantifying the impact of fires on ecosystem structure and C emissions from biomass burning are important elements of evaluating the interactions of climate and ecosystems. Climate condition is a key driver of fire occurrence. In areas where climate change leads to warmer and drier future conditions, the frequency and the intensity of fires are expected to increase, causing significantly larger impacts on ecosystem function and its role on the global carbon cycle [Liu et al., 2013; Stocks et al., 1998]. In addition, areas with increased fuel load, due to the loss of traditional woodcutting activities and fire exclusion policies, may adversely favor the occurrence of extreme fires with large socioeconomic and ecological consequences. These so-called megafires [Attiwill and Binkley, 2013; San-Miguel-Ayanz et al., 2013] release to the atmosphere huge amounts of C accumulated during decades of tree growth. The expected increase in the intensity and frequency of megafires [Attiwill and Binkley, 2013] requires improved accounting of the C released by these fires and better characterization of their role in the climate system. Quantification of the C emissions from biomass burning can be achieved using remote sensing by either direct top-down measurement of trace gases released during the fire [Arellano et al., 2004; Kaufman et al., 1992] or by GARCIA ET AL. C RELEASED BY THE CALIFORNIA RIM FIRE 340 PUBLICATION S Journal of Geophysical Research: Biogeosciences RESEARCH ARTICLE 10.1002/2015JG003315 Key Points: •Quantification of biomass consumption and carbon release by a megafire •Integration of airborne LiDAR and Landsat OLI data to estimate biomass consumption •Uncertainty of estimates calculated and provided for each biomass consumption estimation Supporting Information: •Supporting Information S1 •Figure S1 •Figure S2 •Figure S3 •Figure S4 Correspondence to: M. Garcia, [email protected]; [email protected] Citation: Garcia, M., S. Saatchi, A. Casas, A. Koltunov, S. Ustin, C. Ramirez, J. Garcia-Gutierrez, and H. Balzter (2017), Quantifying biomass consumption and carbon release from the California Rim fire by integrating airborne LiDAR and Landsat OLI data, J. Geophys. Res. Biogeosci.,122, 340–353, doi:10.1002/ 2015JG003315. Received 22 DEC 2015 Accepted 19 JAN 2017 Accepted article online 26 JAN 2017 Published online 18 FEB 2017 ©2017. The Authors. This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. an indirect bottom-up modeling approach in which burned area observations are combined with ecosystem biomass and fuel data [Seiler and Crutzen, 1980] Cr¼ABfcβ;(1) where Ais the area burned (in hectare), Bis the biomass density (Mg ha 1 ), f c is the fraction of the biomass that is carbon, and βis the combustion completeness (CC) or burning efficiency, i.e., the fraction of biomass consumed during burning. Uncertainties in the C release derived from the model by Seiler and Crutzen [1980] result from the uncertainty in each of its parameters. van der Werf et al. [2006] pointed out Aas the most uncertain parameter at a global scale; however, uncertainties in biomass density and combustion completeness significantly impact emission estimates as well. Fuel load, i.e., the amount of biomass available for burning, is usually estimated based on biome-averaged values, which have large uncertainties due to the spatial variability in biomass within the biome [van der Werf et al., 2006; van Leeuwen et al., 2014] and may not represent well the biomass density of the affected area. More recently, biogeochemical models that account for the effect of herbivory and fuel wood collection have been used to estimate fuel load [van der Werf et al., 2006, 2010]. Nevertheless, the accuracy of the fuel load estimates is subject to the difficulty of parameterizing the model appropriately. Thus, development of more accurate methods for estimating aboveground biomass and the fraction of the biomass consumed during the fire is crucial to improve and validate fire emission estimates, as well as to correct estimates of fuel loads yielded by biogeochemical models used at regional and global scales. Remote sensing provides a robust way of deriving fuel load. In particular, airborne and satellite LiDAR measurements have allowed development of novel techniques to accurately quantify the vertical and horizontal structures and the aboveground biomass over a wide variety of forest ecosystems [García et al., 2012; Hudak et al., 2012; Næsset et al., 2013]. As most LiDAR measurements are currently acquired from airborne platforms, the data are often limited in terms of spatial and temporal coverage. This hinders the analysis of prefire and postfire ecosystem status, in most cases due to the lack of data that would enable the characterization of aboveground biomass before the fire. On the other hand, multispectral optical and Radar sensors onboard satellites provide wide-area coverage and frequent observations, which makes these sensors particularly well suited for monitoring fuel load dynamics. As a result, integration of airborne LiDAR with optical or Radar satellite observations provides a convenient alternative to overcome the shortcomings of LiDAR data availability, allowing for more accurate representation of the fuel load and the biomass dynamics, thus improving fire emission estimates. Here we explore both the spatial and temporal extrapolation of the LiDAR data with multispectral Landsat OLI (Operational Land Imager) data over large areas to quantify the biomass before and after the fire. We also quantify the burned fuel consumption and the amount of carbon released during the Rim megafire in California. The LiDAR data were acquired after the fire event covering the entire burn and a 2 km buffer of unburned vegetation. We developed a methodology to address several specific objectives as follows: (1) to estimate aboveground biomass (AGB) from LiDAR data, (2) to quantify prefire and postfire biomass over the entire burned area, (3) to quantify the biomass consumed and the C released during the megafire, and (4) to assess estimation uncertainty. 2. Materials and Methods 2.1. Study Area The study site comprises the footprint of the Rim fire, in the Sierra Nevada Mountains, California, USA. This megafire started on 17 August 2013 and was not contained until 24 October 2013, burning more than 104,000 ha of the Stanislaus National Forest and Yosemite National Park (Figure 1). The area is topographically rough, with elevations ranging from 60 to 2400 m and slopes of up to 90%. The fire burned through a mosaic of vegetation types, which included low-elevation grasslands, chaparral, and foothill-oak woodland savanna habitat; mixed conifer-broadleaf forests dominated by pines in the lower montane zone; and mixed conifer forests in higher elevation areas dominated by firs. A more detailed description of the study site can be found in Casas et al. [2016]. Journal of Geophysical Research: Biogeosciences 10.1002/2015JG003315 GARCIA ET AL. C RELEASED BY THE CALIFORNIA RIM FIRE 341 According to the burn severity map developed by the U.S. Forest Service Rapid Assessment of Vegetation Condition after Wildfire processing system [Miller and Quayle, 2015], the area experienced different levels of burn severity, with 35.8% in high burn severity, 25.4% in moderate severity, 30.9% in low severity, and 7.9% remained unburned. 2.2. Data Sets 2.2.1. Reference Aboveground Biomass Data Sets As part of the study, 71 circular plots of 0.09 ha (equivalent to the Landsat pixel resolution) were collected across a 2 km buffer zone around the fire perimeter. The plots were located following a stratified random sampling scheme using a Landsat-based prefire vegetation map provided by the U.S. Forest Service. The strata were defined by vegetation type (softwood (conifers), hardwoods (broadleaf trees), and mixed forests) and diameter classes (12.7–25.2 cm, 25.2–50.6 cm, 50.6–76 cm, and >76 cm). Although the randomness of plot locations was slightly perturbed by accessibility constraints, the plots represented well the range of available biomass and were used to calibrate and validate the models based on the LiDAR data. For each tree with diameter at breast height (DBH) greater than 10 cm, the species was recorded and the DBH tallied. The center of each plot was positioned using a differential GPS, with a horizontal accuracy after Figure 1. Study area comprising the footprint of the Rim fire, in the Sierra Nevada Mountains, California, USA. Journal of Geophysical Research: Biogeosciences 10.1002/2015JG003315 GARCIA ET AL. C RELEASED BY THE CALIFORNIA RIM FIRE 342 postprocessing better than 0.5 m. The spatial error in two plots exceeded 2 m, and these were discarded from further analysis. Additionally, fractional cover (FC) was estimated from the LiDAR data, and those plots having a FC less than 10% were also removed from the analysis. Thus, the final data set consisted in 65 field plots. For each tree the total aboveground biomass (AGB) was computed using the National Biomass Estimator Library (NBEL) developed by the Forest Management Service Center (FMSC). The NBEL synthesizes published biomass equations for the U.S. and also unpublished allometric models developed by FMSC that are stored along with their associated metadata in a SQLite database [Wang, 2014]. The models used to estimate AGB were functions of DBH and tree height. The height was estimated from the measured DBH using the species level equations compiled by Keyser and Dixon [2012]. Table 1 summarizes the structural attributes of the sampled plots. To extrapolate AGB values using Landsat OLI data, we used the LiDAR-derived AGB as a reference. More than 500 pixels were randomly selected based on the histogram for the LiDAR-derived AGB to ensure that the samples covered the full range of AGB values in the study area. Our sampling was limited to unburned and lowseverity areas. In this way, the effect of fire-induced structural changes was eliminated from the sample. 2.2.2. LiDAR Data and Processing LiDAR data were collected on November 2013 by the National Center for Airborne Laser Mapping using an Optech Gemini Airborne Laser Terrain Mapper instrument that recorded up to four returns per pulse. These data covered the extent of the fire plus a 2 km buffer. The site was flown at a mean elevation of 2200 m above sea level with a maximum scan angle of ± 14° and a nominal 50% overlap between flight lines. The average point density was approximately 20 points m 2 . The vendor provided the point cloud data set and a 1 m digital elevation model (DEM) that was used to normalize the height of each return. We computed 16 variables related to canopy structure derived from the vertical structure of the forest measured by the LiDAR data, including height percentiles (H 25 ,H 50 ,H 75 ,H 90 , and H 99 ), the mean, standard deviation, kurtosis, skewness, and coefficient of variation for the height of the returns above 2 m. Biomass distribution within the canopy was estimated as follows: H 99 –H 50 ,H 99 –H 25 ,H 90 –H 50 , and H 90 –H 25 , as well as the canopy depth (maximum minus minimum vegetation height). Finally, the ratio of hits above 2 m to all returns was computed to represent canopy fractional cover [Hopkinson and Chasmer, 2009; Morsdorf et al., 2006] (Table S1 in the supporting information). These metrics have been shown to have strong relationships with the aboveground biomass [Hudak et al., 2012; Næsset et al., 2013]. In addition, we derived 19 metrics from the intensity of LiDAR returns after normalization to a standard range [García et al., 2010]. These metrics included intensity percentiles (I 25 ,I 50 ,I 75 ,I 90 , and I 99 ), the mean, kurtosis, range, skewness, coefficient of variation, and standard deviation, as well as accumulated intensity at H 25 , H 50 ,H 75 ,H 90 , and H 99 . The canopy FC was also estimated as the ratio of the canopy energy to the total energy [García et al., 2010; Hopkinson and Chasmer, 2009], as further detailed in Table S1 in the supporting information. We did not apply a correction factor to the intensity of the ground energy to account for the differences in reflectance between canopy and ground at the wavelength used by the LiDAR system, as suggested by other authors [Lefsky et al., 1999b; Morsdorf et al., 2006], because this factor is site dependent and it was not available for our study site. In addition, we computed the canopy reflection sum [Means et al., 1999] by summing the intensities of the canopy returns within the plot [Hall et al., 2006] and the density-weighted canopy reflection sum [García et al., 2010], which accounts for variation in point density throughout the study site that results from topography or scan angle differences. Following the approach described by Muss et al. [2011], we also constructed pseudowaveforms for each plot and derived six additional metrics related to the vegetation spatial distribution. These metrics included the Table 1. Properties of the Plots Measured in the Field a DBH (cm) Height (m) AGB (Mg ha 1 ) Stem Density (tree ha 1 ) Maximum 142.5 53.6 645.4 1455.6 Minimum 10.0 4.9 24.21 22.2 Average 30.4 16.7 195.8 358.8 Standard deviation 19.9 8.5 143.1 263.7 a The plot size is 0.09 ha (Landsat OLI pixel size on the ground). Journal of Geophysical Research: Biogeosciences 10.1002/2015JG003315 GARCIA ET AL. C RELEASED BY THE CALIFORNIA RIM FIRE 343 height of the median energy and the height/median ratio, which are related to AGB [Drake et al., 2002; Hyde et al., 2005], the canopy height profile (CHP) describing canopy structure [Lefsky et al., 1999a], the mean canopy height, and the quadratic mean canopy height. In addition, we included the CHP variation coefficient and the area under the canopy waveform (AUCW) representing vertical heterogeneity and the amount of canopy material, respectively [Bouvier et al., 2015; Muss et al., 2011]. 2.2.3. Landsat OLI Data and Processing We selected two orthorectified Landsat OLI reflectance images (path/row: 043/034) to estimate the biomass before and after the fire. The image dates were chosen based on the following two criteria: (1) minimize time gaps between the fire event and Landsat and LiDAR data acquisitions and (2) cloud and snow cover over the burned area must be zero. As a result, we chose the prefire image acquired on 30 July 2013 and the postfire image acquired on 18 October 2013. Although the postfire image was collected 6 days before the fire was fully contained, this was the only image that satisfied the above criteria. The next available image was acquired 1 year after the fire, and by then, cleaning and salvage logging operations had introduced structural changes that were undesirable for our purpose. The imagery, downloaded from the U.S. Geological Survey Earth Explorer web site (http://earthexplorer.usgs. gov/; accessed: 21 August 2015), had been processed based on the Landsat Ecosystem Disturbance Adaptive Processing System (LEDAPS) atmospheric correction [Masek et al., 2006]. To ensure an appropriate coregistration between the different data sets used, the Landsat images were geocorrected by collecting ground control points (n>25) using an intensity LiDAR image as reference and applying a linear transformation model, yielding an RMSE <1 pixel. Estimation of AGB and its dynamics from multispectral data relies on establishing relations between biomass and vegetation spectral response [Foody and Mathur, 2004; Meng et al., 2007; Pflugmacher et al., 2014], which is usually performed using spectral indices and/or transformations. We selected Landsat OLI reflectance bands 2 through 7 and derived four vegetation indices: the normalized difference vegetation index (NDVI) [Tucker, 1979], Normalized Difference Infrared Index (NDII) [Hunt and Rock, 1989], Enhanced Vegetation Index [Huete et al., 2002], and Visible Atmospherically Resistant Index [Gitelson et al., 2002]. We also used brightness, greenness, and wetness components of the tasseled cap transformation [Kauth and Thomas, 1976] using the OLI-specific coefficients [Baig et al., 2014]. Tasseled cap components, particularly wetness, have been used to assess forest structure and forest change [Cohen et al., 1995; Hansen et al., 2001; Pascual et al., 2010]. From the brightness and greenness component we computed the tasseled cap angle and tasseled cap distance metrics [Powell et al., 2010] that are related to vegetation cover and vegetation composition and structure, respectively [Pflugmacher et al., 2014]. In Table S2 in the supporting information, we provide additional details. Finally, for each band we calculated a set of isotropic texture metrics from the gray level co-occurrence matrix for window sizes varying from 3×3 to 9×9 pixels. These metrics included the following: homogeneity, contrast, standard deviation, entropy, angular second moment, and correlation. They were previously found effective for biomass estimation in different forest ecosystems [Kelsey and Neff, 2014; Lu, 2005]. 2.3. Ancillary Data The DEM provided with the LiDAR data was resampled to 30 m using the average of all pixels included within the 30 m cell. Subsequently, the slope and aspect were computed for each pixel. A Landsat-based prefire vegetation map was provided by the U.S. Forest Service and reclassified into three vegetation type groups as follows: softwood, hardwoods, and mixed forest. These auxiliary layers were used in subsequent statistical analyses. A burn severity map derived from Landsat data, which was sensitive to different burn severity levels [Key and Benson, 2006], was also provided by the U.S. Forest Service, separating unburned areas from low, medium, and high severity burn levels. 2.4. AGB Estimation In order to model AGB, we used a least squares support vector regression (SVR) approach (see supporting information for more details), using the LS-SVMlab toolbox developed by De Brabanter et al. [2011] and implemented in MATLAB [MATLAB, R2010a]. Nonparametric regression techniques are attractive because they do not make explicit assumptions about data distributions and thus can model complex relationships between dependent and independent variables [Zhao et al., 2011]. We used a Gaussian radial basis Journal of Geophysical Research: Biogeosciences 10.1002/2015JG003315 GARCIA ET AL. C RELEASED BY THE CALIFORNIA RIM FIRE 344 function kernel, controlled by its bandwidth (h). The second parameter determining the estimation accuracy is the penalty parameter (γ), which represents a trade-off between the model complexity and error on the training data. These two parameters were obtained by a grid search with tenfold cross validation. A different support vector regression (SVR) model was developed for each of the remote sensing data sets we used. Given the large number of metrics derived from each data set, we used three methods for feature selection to reduce the processing time, increase model generalization performance, and help with the interpretation (see supporting information for more details) [Weston et al., 2000]. The model based on the metrics (features) derived from the LiDAR data was calibrated and validated using the field measurements described in section 2.2.1. However, to develop and validate the models based on the Landsat data (prefire and postfire biomass models), we used the LiDAR-based AGB estimates, which allowed us to construct a much larger reference data set. Thus, for each date, we created a data set of 514 samples, with the features derived from the Landsat OLI image and the ancillary data. In both cases, approximately 70% of the samples were used for model training/calibration and the remaining samples were used for independent validation. An initial evaluation of the LiDAR-based postfire AGB estimates showed the presence of high AGB values within the high burn severity class. These areas had high H 50 values, while the amount of canopy material (AUCW) was very low. Due to the strong correlation between H 50 and AGB, this resulted in an overestimation of AGB in high burn severity areas. Therefore, we modeled the relationship between H 50 and AUCW and derived a correction factor that was applied to the H 50 metric to account for the changes in the distribution of canopy returns as a result of the fire. Subsequently, we recomputed AGB based on the weighted H 50 (see supporting information for more details). Although the overestimation for H 50 was most evident in high burn severity areas, we applied the AGB correction over the entire burn, including all severity levels. This was motivated by the fact that burn severity maps derived from Landsat data may not accurately capture the actual damage caused by the fire to understory and midstory vegetation in low and moderate severity areas [Miller and Quayle, 2015], especially if canopy cover is high. Hereinafter, the corrected AGB estimates are referred to as the “corrected LiDAR-based AGB”or the AGB obtained with corrected LiDAR-based data. 2.5. Biomass Consumption and Carbon Release Estimation The amount of biomass consumed by the Rim fire was computed as the difference between the prefire AGB estimated from Landsat and the postfire AGB estimates derived from LiDAR and Landsat, as discussed above in section 2.4. Furthermore, De Santis et al. [2010] related the changes in different vegetation indices to combustion completeness, i.e.,the fraction of biomass consumed by the fire.Therefore,inadditiontothecorrected LiDAR-based postfire AGB, we also computed the difference between the prefire and postfire Landsat AGB estimates. In this way, we could evaluate the extent to which optical data alone can capture changes in biomass. Finally, after computing the consumed biomass we estimated the released C by multiplying the biomass by a C fraction coefficient 0.5, as it is commonly done [Marklund and Schoene, 2006]. We reported the consumed biomass for the entire study area and also per burn severity levels. 2.6. Uncertainty Analysis The performance of the SVR models was evaluated using the coefficient of determination (R 2 ), the adjusted coefficient of determination (R 2 adj), the root-mean-square error (RMSE), and the relative root-mean-square error (relRMSE) (Table S3 in the supporting information). The AGB extrapolation to the entire study area was based on a two-step regression approach as follows: (i) first, field measurements were used to model AGB from LiDAR data and then (ii) the LiDAR-based AGB so obtained was used to train the Landsat AGB model. The errors in these steps, denoted by RMSE 1 and RMSE 2 , propagated to the final estimates. To calculate the overall error (variance) of the two-step AGB estimator (i.e., σ 2 modeling ), we assumed the error of the individual steps to be independent, which leads to σ 2 modeling = RMSE 12 + RMSE 22 . Other sources of error contributing to the overall uncertainty, such as remote sensing measurement error, sensor noise, field measurement error, and other sources (see supporting information), were not included in the uncertainty analysis and considered negligible compared to the error due to modeling. A comprehensive analysis of the error propagation and accounting for all error sources requires additional data not available for this study. The Journal of Geophysical Research: Biogeosciences 10.1002/2015JG003315 GARCIA ET AL. C RELEASED BY THE CALIFORNIA RIM FIRE 345 ground-based AGB estimation error was considered negligible, as there was no uncertainty associated with the allometric equations used in NBEL. We assessed the errors in the AGB and the biomass consumption estimates at a burn severity class level and combined these class-wide errors into the overall, site-wide quantities. This analysis accounted for spatial autocorrelation of the AGB error at the pixel scale using the techniques described by McRoberts [2006] and Weisbin et al. [2014] (details are provided in the supporting information). As the biomass consumption estimate is the difference between two correlated AGB estimates, its variance was estimated following McRoberts et al. [2014]. 3. Results 3.1. AGB Estimates Using LiDAR Data Table 2 shows the results for the SVR model based on the metrics selected with one of the feature selection methods, the expert knowledge method. This model yielded the best results over the validation data and did not include correlated metrics as happened with the evolutionary algorithm. Therefore, this model was used to generate the LiDAR-based AGB map for the study area. The results of the models built using the other two feature selection methods, the stepwise regression and the evolutionary algorithm, are given in Table S4 in the supporting information. Figure 2 shows the scatterplot of the LiDAR-based versus field-based AGB estimates. The application of the SVR model using the H 50 corrected based on the AUCW (section 2.4) showed significant differencescomparedtotheoriginalimage.AftercorrectingH 50 ,thenewAGBspatialdistributionagreedbetter with the burn severity map (as illustrated in Figure3), with an ~80% decrease in AGB in high burn severity areas and an ~50% reduction in uncertainty (Table 3). Between the two LiDAR-basedmethods of AGB estimation, the AGB based on the corrected LiDAR data reflects more accurately the variation in the postfire AGB over the landscape (Figure 3). Table 3 also presents the mean AGB estimates and their uncertainties for different burn severity levels, vegetation types, and AGB estimation methods. 3.2. AGB Estimates Using Landsat OLI Data Our first attempt to model AGB from Landsat data was based on a model that included texture that is the entropy of band 3 in a 9 × 9 window. This model yielded an R 2 of 0.73 and 0.61 and a relative RMSE of 46% and 56% for the calibration and the validation data, respectively. Nevertheless, when the AGB map was produced, it exhibited a clear pattern related to the entropy band rather than to the actual AGB distribution (Figure S4 in Table 2. Accuracy of the SVR Model of AGB Based On Metrics Selected by Expert Knowledge a Feature Selection Selected Variables R 2 R 2 Adj RMSE (Mg ha 1 ) relRMSE (%) Expert knowledge AUCW H 50 0.82 0.81 59.98 30.63 0.82 0.79 67.18 34.31 0.81 0.80 62.17 31.76 a For each accuracy measure, the first line corresponds to the calibration data set, the second line corresponds to the validation data set, and the third line (in bold) corresponds to the accuracy for the combined data sets. Figure 2. LiDAR-based AGB estimates compared to field-based AGB estimates. The solid line represents the 1:1 line. Journal of Geophysical Research: Biogeosciences 10.1002/2015JG003315 GARCIA ET AL. C RELEASED BY THE CALIFORNIA RIM FIRE 346 the supporting information). Therefore, we performed a second analysis that excluded texture measures to avoid this artifact. The final features selected to model AGB from Landsat were the OLI reflective bands 2 to 6, the NDII, and the elevation. Table 4 summarizes the performance of the prefire and postfire models. Figure 4 shows the scatterplots of the Landsat-based versus LiDAR-based AGB for each model. The combined uncertainty of the Landsat AGB estimates was computed using the RMSE of the models developed at each step (section 2.6) using the complete reference data set, i.e.,the 65 plots for the LiDAR model and the 514 pixels for the Landsat model. This resulted in an overall relative RMSE of 55.3% and 49.9% for the postfire and prefire AGB, respectively. We also computed the uncertainties for the estimates of the mean AGB over areas with different burn severity levels and vegetation types. The uncertainty of the aggregated estimates was much smaller due to the large number of pixels and relatively weak spatial correlation (Table 3). However, the overall uncertainty of the postfire AGB estimates obtained with Landsat was double the uncertainty of corrected LiDAR-based AGB, particularly in high severity regions where Landsat significantly overestimated LiDAR AGB by more than 150%. 3.3. Estimation of the Carbon Released by the Rim Fire Table 5 presents the estimates for the biomass loss and C released to the atmosphere during the Rim fire using the three different postfire AGB estimates derived from the LiDAR and the Landsat data. In all cases, the AGB estimates were corrected for bias based on a linear regression between the observed and the estimated values (Figures 2 and 4). This process improves the dilution bias existing in the results but may introduce some random errors due to adjustments and changes introduced on all pixels, even with high estimation accuracy [Xu et al., 2016]. The estimates of the consumed biomass changed from 2.75 ± 0.21 Tg obtained using the original LiDAR-derived AGB to 6.58 ± 0.03 Tg derived using the corrected LiDAR data. Applying a 50% biomass to carbon conversion factor, we obtained the estimates of the released C that Figure 3. Spatial distribution of the postfire AGB estimated from LiDAR data: (a) burn severity map (added for comparison), (b) AGB using the original H 50 metric, and (c) AGB using a corrected H 50 metric (also referred to as the “corrected LiDAR AGB”). Table 3. AGB Mean ± Std Error (in Units of Mg ha 1 ) Within the Strata Defined by Burn Severity and Vegetation Type for Estimates Obtained Using LiDAR and Landsat OLI Data Strata LiDAR-Based AGB Corrected LiDAR- Based AGB Landsat-Based Prefire AGB Landsat-Based Postfire AGB Unburned Coniferous 216.77 ± 0.64 216.77 ± 0.64 214.07 ± 1.49 209.39 ± 2.24 Deciduous 74.81 ± 0.44 74.81 ± 0.44 81.64 ± 0.78 95.32 ± 1.55 Mixed 147.56 ± 0.87 147.56 ± 0.87 115.74 ± 0.57 118.01 ± 0.70 Low severity Coniferous 250.08 ± 0.45 183.81 ± 0.32 235.34 ± 0.79 232.69 ± 0.57 Deciduous 101.42 ± 0.80 94.48 ± 0.42 96.77 ± 2.40 102.47 ± 0.93 Mixed 193.71 ± 0.62 146.36 ± 0.54 162.81 ± 1.04 154.96 ± 0.96 Moderate severity Coniferous 190.29 ± 0.54 90.65 ± 0.22 184.56 ± 1.09 122.58 ± 0.36 Deciduous 76.74 ± 0.45 65.38 ± 0.28 82.86 ± 1.05 77.25 ± 0.62 Mixed 144.28 ± 0.61 88.4 ± 0.35 136.58 ± 2.63 90.04 ± 0.57 High severity Coniferous 151.05 ± 0.32 34.97 ± 0.14 172.13 ± 0.45 92.66 ± 0.28 Deciduous 64.34 ± 0.36 30.60 ± 0.22 81.50 ± 1.84 92.53 ± 0.89 Mixed 132.23 ± 0.44 33.89 ± 0.20 133.71 ± 0.55 90.46 ± 0.45 Journal of Geophysical Research: Biogeosciences 10.1002/2015JG003315 GARCIA ET AL. C RELEASED BY THE CALIFORNIA RIM FIRE 347 ranged from 1.37 ±0.11 TgC to 3.29 ±0.02 TgC. For the Landsat-based approach the consumed biomass estimate was 3.93 ± 0.17 Tg or 1.87 ±0.08 TgC. Figure 5 shows the spatial distribution of AGB difference between the Landsat prefire and the Landsat (Figure 5b) and LiDAR (Figure 5c) postfire estimates, representing the biomass consumption across the Rim fire domain. The difference between the LiDAR and Landsat data used for postfire AGB estimates can be readily observed in the significant underestimation by Landsat of the fire impact in terms of the biomass loss, particularly in the high severity areas. The AGB loss image derived from corrected LiDAR data shows the best visual agreement with the burn severity map (Figure 5a). Nonetheless, it can be seen that the match between the biomass difference map and the burn severity map is not perfect. One reason for this is that the burn severity map is only an approximate qualitative representation of the fire damage. Besides that, the burn severity map reflects a relative change in the overall vegetation spectral signal rather than an absolute measure of the biomass loss. Therefore, different burn severity levels could represent similar biomass losses depending on the prefire biomass levels and the extent and severity of forest disturbance, and the actual biomass loss can vary substantially within a burn severity class. Furthermore, we quantified the biomass loss for each burn severity class and estimated the proportion of biomass consumed by the fire, i.e., the combustion completeness factor, for each severity level (Table 6). It can be seen from Table 6 that for all burn severity levels, both, Landsat and the uncorrected LiDAR data underestimated the CC values provided by De Santis et al. [2010] (last row in Table 6), while the corrected LiDAR data consistently overestimated those values. Table 4. Features Selected From Landsat OLI for the SVR AGB Model (Left Most Column) and the Model Accuracy for Each Date a Variables Selected R 2 R 2 adj RMSE (Mg ha 1 ) relRMSE (%) Prefire B2-B6, NDII, elevation 0.79 0.79 75.43 36.02 0.72 0.71 87.82 41.94 0.76 0.76 80.65 38.52 Postfire 0.73 0.72 87.18 41.63 0.6 0.58 105.01 50.15 0.67 0.67 94.77 45.26 a For each accuracy measure and for each date, the first line corresponds to the calibration data set, the second line corresponds to the validation data set, and the third line (in bold) corresponds to the combined data set. Figure 4. Scatterplots of LiDAR-derived versus Landsat-derived AGB for the (a) prefire and (b) postfire images. Journal of Geophysical Research: Biogeosciences 10.1002/2015JG003315 GARCIA ET AL. C RELEASED BY THE CALIFORNIA RIM FIRE 348