Monitoring Post-Fire Recovery of Chaparral and Conifer Species Using Field Surveys and Landsat Time Series

: Recovery trajectories derived from remote sensing data are widely used to monitor ecosystem recovery after disturbance events, but these trajectories are often retrieved without a precise understanding of the land cover within a scene. As a result, the sources of variability in post-disturbance recovery trajectories are poorly understood. In this study, we monitored the recovery of chaparral and conifer species following the 2007 Zaca Fire, which burned 97,270 ha in Santa Barbara County, California. We combined ﬁeld survey data with two time series remote sensing products: the relative delta normalized burn ratio (RdNBR) and green vegetation (GV) fractions derived from spectral mixture analysis. Recovery trajectories were retrieved for stands dominated by six different chaparral species. We also retrieved recovery trajectories for stands of mixed conifer forest. We found that the two remote sensing products were equally effective at mapping vegetation cover across the burn scar. The GV fractions (r(78) = 0.552, p < 0.001) and normalized burn ratio (r(78) = 0.555, p < 0.001) had nearly identical correlations with ground reference data of green vegetation cover. Recovery of the chaparral species was substantially affected by the 2011–2017 California drought. GV fractions for the chaparral species generally declined between 2011 and 2016. Physiological responses to ﬁre and drought were important sources of variability between the species. The conifer stands did not exhibit a drought signal that was directly correlated with annual precipitation, but the drought likely delayed the return to pre-ﬁre conditions. As of 2018, 545 of the 756 conifer stands had not recovered to their pre-ﬁre GV fractions. Spatial and temporal variation in species composition were important sources of spectral variability in the chaparral and conifer stands. The chaparral stands in particular had highly heterogeneous species composition. Dominant species accounted for between 30% and 53% of the land cover in the surveyed chaparral patches, so non-dominant land cover types strongly inﬂuenced remote sensing signals. Our study reveals that prolonged drought can delay or alter the post-ﬁre recovery of Mediterranean ecosystems. It is also the ﬁrst study to critically examine how ﬁne-scale variability in land cover affects time series remote sensing analyses.


Introduction
Wildland vegetation in Southern California frequently experiences intense crown fires that significantly alter ecosystem function by destroying aboveground biomass and altering soil properties [1]. Fires affect carbon cycling, nutrient cycling, soil water repellency, runoff, erosion, streamflow, and sediment export [2][3][4]. Burned vegetation is also the largest driver of carbon stock loss in California [5]. These effects generally diminish as vegetation regrows [2,4], although fire-induced changes in the composition of vegetation patches can lead to long-term changes in ecosystem structure and function [6].
Monitoring burn severity and post-fire vegetation regrowth are essential for detecting the ecological effects of wildfires, especially as large fires become more frequent in the western United States [7]. Field surveys are a common method of assessing ecological changes after fires, but they are expensive and time consuming, especially for continuous landscape-scale monitoring. Remote sensing is another method of monitoring land cover change, which provides broad spatial and temporal coverage. Landsat imagery covering the entire globe has been collected regularly for over 45 years and is now freely available for public use [8]. Recent studies have taken advantage of this data availability by developing time series data sets to monitor ecosystem recovery after wildfires [9][10][11]. The recovery trajectories derived from these data sets are an improvement over simple two-date change detection models because they capture the temporal variability that occurs after fires, including inter-and intra-annual trends [12].
Various remote sensing products have been used to develop post-fire recovery trajectories. Spectral indices are the most common method of assessing burn severity and post-fire vegetation regrowth using Landsat imagery. Many studies rely on the normalized burn ratio (NBR) and its variants for these purposes [13,14]. Spectral indices such as NBR primarily measure changes in greenness compared to pre-fire conditions. These indices are easy to calculate and are computationally efficient, however they produce unitless values that do not directly measure any biophysical property [15]. Spectral mixture analysis (SMA) is another method of assessing fire effects using remote sensing data [10,[16][17][18]. SMA estimates the fractional cover of different materials, such as green vegetation, within a pixel. While spectral indices such as NBR produce unitless measures of vegetation greenness, SMA can generate absolute estimates of green vegetation (GV) cover inside of a pixel. As a result, SMA fractions are physically meaningful values that can be compared with observations by field ecologists [15]. This has obvious advantages for studies that integrate remote sensing and field ecology data. Fire effects also vary at a scale that is finer than the spatial resolution of many spaceborne sensors [19]. Spectral indices that assign a single value to a pixel do not capture this subpixel variability, but SMA can decompose the remote sensing signal in order to determine a pixel's composition. As a result, SMA is more sensitive to the variability that occurs within vegetation patches after fires. That being said, SMA models require extensive calibration in order to produce accurate retrievals.
While remote sensing trajectories are often used to monitor ecosystem recovery after wildfires, previous studies have retrieved recovery trajectories without a precise understanding of the underlying land cover. Plant species exhibit different physiological responses to fires, so the types and number of species present in a patch affect its recovery behavior. The structure and composition of vegetation patches also change over time as vegetation regrows [20,21]. These factors affect the timing and magnitude of vegetation greenness signals in remote sensing data, but few studies have examined these sources of variability [17]. Likewise, few studies have examined whether individual species or functional types exhibit unique recovery trajectories after the same disturbance event [17,22,23]. This study fills that gap. We combine multi-temporal remote sensing imagery with field survey data to retrieve post-fire recovery trajectories for chaparral species and stands of mixed conifer forest. We also examine how fine-scale variability in species composition influences post-fire recovery trajectories. In doing so, we answer the following research questions: 2. How variable are post-fire recovery trajectories for different plant species and functional types?
What drives this variability? 3. How heterogeneous are the structure and composition of vegetation patches in chaparral and mixed conifer ecosystems? How does this heterogeneity affect post-fire recovery trajectories?

Study Area
This study examines vegetation regrowth following the 2007 Zaca Fire in Santa Barbara County, California ( Figure 1). The Zaca Fire burned 97,270 ha starting on 4 July 2007. It was fully contained on 2 September 2007. The majority of the burn scar is located in undeveloped wilderness areas inside of Los Padres National Forest. The elevation of the burn scar ranges from 254 m to 2077 m. Mean annual precipitation ranges from 324 mm to 1260 mm depending on elevation and topographic position, with 95% of the precipitation occurring between November and March [24]. The majority of vegetation cover inside the burn scar is chaparral shrubland, which is characterized by the presence of species such as Ceanothus spp., Arctostaphylos spp., and Adenostoma fasciculatum. Coastal sage scrub species such as Eriogonum fasciculatum occur in xeric locations. Oak species including Quercus berberidifolia, Quercus chrysolepis, and Quercus wislizenii are also prevalent across the landscape. Stands of conifers occur intermittently throughout the burn scar, especially at higher elevations. These stands include species such as Pinus coulteri, Pinus jeffreyi, Pseudotsuga macrocarpa, Abies concolor, and Calocedrus decurrens.
The dominant chaparral and conifer species are evergreen, but many herbaceous species senesce during the summer dry season. Botanical nomenclature follows Baldwin et al. [25]. Recovery from the Zaca Fire was affected by the 2011-2017 California drought. Annual precipitation in California is summarized using water years that last from 1 October until 30 September. Santa Barbara County received abnormally low rainfall during the water years of 2011-2015, when county-wide rainfall was 40 to 69 percent of average ( Figure 2) [26]. The water year of 2013-2014 was the worst drought year in the last 1200 years [27], and the period of drought from 2012 to 2015 was unprecedented in the historical record [28]. An atmospheric river provided substantial drought relief in January and February 2017 [29]. County-wide rainfall in the 2016-2017 water year was 136 percent of average [26]. Dry conditions returned in the next water year, and Santa Barbara County remained in some degree of drought condition until February 2019 [30].

Imagery
Landsat images from 2000 to 2018 (path 42, row 36) were used to monitor vegetation changes across the burn scar. Images from 2012 were omitted from the analysis because Landsat 5 failed in 2012 and Landsat 8 was not launched until 2013. Landsat 7 data were not used to fill the gap because its scan line corrector failed in 2003 [31]. The Level 2 surface reflectance products were acquired from the U.S. Geological Survey [32]. All images were captured in September, except for the 2015 image, which was captured in August because there was no cloud-free image available for September. September imagery was selected to maximize vegetation senescence at all elevations across the burn scar, which minimized variability in phenology caused by topography. As a result, GV cover is at its annual minimum in the September imagery. September image dates also coincide with the end of the California water year. The images from 2000 to 2006 were used to calibrate baseline metrics of vegetation cover. Based on fire history data, 98.3% of the burn scar had not burned in at least 40 years at the time of the Zaca Fire, so the pre-fire images represented mature vegetation [33]. The images from 2007 to 2018 were used to monitor vegetation regrowth after the fire. The measurements from 2013 to 2016 were affected by severe drought.
Imagery from the Airborne Visible/Infrared Imaging Spectrometer (AVIRIS) was used to build spectral libraries for SMA. AVIRIS is a 224 channel imaging spectrometer that was flown on a Lockheed ER-2 [34]. AVIRIS flightlines have been flown over the Zaca burn scar regularly since 2013 as part of the HyspIRI airborne campaign [35]. This study utilized five AVIRIS flightlines flown on 24 August 2015 and 17 June 2016. The five flightlines covered the majority of the burn scar. The Level 2 surface reflectance products for these images were acquired from NASA [36]. Additional pre-processing was performed on the 2015 flightlines in order to improve the reflectance retrievals and ensure consistent georegistration. See Meerdink et al. [37] for more information.

Relative Delta Normalized Burn Ratio
The normalized burn ratio (NBR) is commonly used to monitor land cover change after wildfires. NBR is derived from Landsat imagery using the formula where NIR is the near-infrared band and SWIR2 is the second shortwave infrared band [13]. The difference between pre-fire and post-fire NBR, known as delta NBR (dNBR), is often used as a metric of burn severity and vegetation regrowth [13]. NBR is sensitive to the amount of green vegetation in a pixel, so the largest dNBR values occur in pixels that had large amounts of green vegetation before a fire. To compensate for this bias, Miller and Thode [14] proposed the relative dNBR (RdNBR), which is calculated using the formula RdNBR normalizes dNBR by pre-fire vegetation cover, so it is a more reliable metric for comparing disturbance and recovery patterns across disparate sites. RdNBR values are correlated with burn severity. Values near zero indicate that a pixel is undistinguishable from pre-fire conditions [14]. Previous studies of forested ecosystems have reported that RdNBR values greater than 650 indicate high burn severity. The maximum observed RdNBR values in those studies were approximately 1300 [14,38]. In a study of a chaparral ecosystem, Hubbert et al. [39] reported that 99% of pixels had RdNBR values greater than 1000.
NBR was calculated for all of the Landsat images. RdNBR was then calculated for each year from 2007 to 2018. The median NBR value from 2000 to 2006 was used as NBR prefire for all RdNBR calculations. The NBR values from 2007 to 2018 were used as NBR postfire . This produced an RdNBR recovery trajectory for each pixel inside of the burn scar.

Multiple Endmember Spectral Mixture Analysis
SMA estimates the fractional cover of materials inside of a pixel [40]. SMA assumes that a pixels' spectrum is a linear combination of the spectra of the materials inside of the pixel, weighted by the fractional cover of each material [41]. SMA estimates the fractional cover of different materials by decomposing the observed spectrum using reference spectra of pure materials, known as endmembers. Endmembers can be extracted from images or collected in situ using a field spectrometer. Field spectra can be captured at the leaf scale, branch scale, or stand scale. Image spectra are typically stand-scale measurements [42]. The SMA decomposition is based on the equations where S is an observed spectral mixture for location i at wavelength λ. E is the set of N endmembers, k, which are weighted by the fractions f ki . ε is the residual between the modeled and measured spectra. The selected model for a pixel minimizes the root mean squared error (RMSE) between the modeled and measured spectra. Multiple endmember spectral mixture analysis (MESMA), a variant of SMA, allows the number and types of endmembers to vary by pixel. It also places additional constraints on acceptable models [43]. Because the number and types of endmembers vary by pixel, MESMA spectral libraries can contain multiple endmembers for each class and endmembers for classes that are not present across the entire scene. As a result, both the size of the spectral library and the number of potential models for each pixel can be very large. To address this, several techniques have been developed to select potential endmembers and choose the ones that form the best possible spectral library [44][45][46]. Endmember selection is an important step that determines the accuracy of the final product [47,48].

Spectral Library Development
For this study, MESMA GV fractions were used to monitor vegetation regrowth inside of the burn scar. Endmembers for the MESMA spectral library were extracted from the AVIRIS images and used to unmix the Landsat images. AVIRIS endmembers were used because the spectra were acquired at an appropriate spatial scale for ecological studies [42] and can be convolved to the spectral resolutions of broadband sensors such as Landsat. While endmembers can be extracted directly from Landsat imagery, the higher spectral resolution of AVIRIS data makes it possible to identify narrow absorption features, which may reveal the contents of the pixel in more detail [49]. The higher spatial resolution of AVIRIS also increases the likelihood of extracting spectrally pure endmembers.
The goal of the MESMA analysis was to quantify the fractional cover of GV, non-photosynthetic vegetation (NPV), and soil within each pixel, so potential endmembers were stratified across those classes. Classification and Assessment with Landsat of Visible Ecological Groupings (CALVEG) data were used as ground reference data during endmember extraction. CALVEG is a GIS database of California land cover that is maintained by the U.S. Forest Service [50]. We extracted endmembers from areas dominated by six CALVEG alliances that represent different plant species and functional types: 1. Ceanothus chaparral alliance 2. Chamise alliance 3. Scrub Oak alliance 4. Coast Live Oak alliance 5. Bigcone Douglas-fir alliance 6. Black Cottonwood alliance GV spectra were extracted from the June 2016 AVIRIS flight lines. Approximately 15 GV spectra were selected from each flight line for each alliance, totaling 480 potential GV endmembers. NPV spectra were extracted from the August 2015 flight lines in order to maximize vegetation senescence across the scene. We extracted 495 potential endmembers from areas dominated by annual grasses and forbs, which were senesced at the time the images were captured. Soil spectra were also extracted from the August 2015 flight lines in order to minimize spectral contamination from photosynthetic vegetation. We extracted 559 potential soil endmembers from areas covered by barren ground, agricultural soil, and streambed alluvium.
Iterative endmember selection (IES) was used to select the most representative endmembers from each subclass (i.e., the six GV alliances, the three soil types, and NPV) [44,51]. For each subclass, IES picks the pair of endmembers that has the highest kappa statistic for that subclass. It then iteratively adds and removes endmembers from the selection in order to increase the value of the kappa statistic. IES is completed when it is not possible to increase the value of the kappa statistic by adding or removing any more endmembers from the spectral library. IES is an automated process that selected 26 GV endmembers, 17 soil endmembers, and 3 NPV endmembers for inclusion in the spectral library.
Eight additional endmembers were manually added to the library. Three of the eight endmembers (two NPV endmembers and one soil endmember) were extracted from WestUSA, which is a reference spectral library that was created using a field spectrometer [42]. These endmembers were brighter than the image spectra and helped ensure that the spectral library captured the full range of illumination variability that occurred within the scene. Three ash endmembers were extracted from an AVIRIS image flown shortly after the 2004 Gaviota fire in Santa Barbara County, California. A fourth ash endmember was extracted directly from the September 2007 Landsat image. The burn scar was covered in ash when the 2007 Landsat image was collected, so the ash endmembers helped model parts of that image. Ash is short-lived on the landscape after fires [52], so the four ash endmembers were only included in the spectral library for the 2007 image. Finally, photometric shade with zero reflectance at all wavelengths was added to the spectral library [41]. This enabled MESMA to model observed spectra that are proportionally darker than mixtures of endmembers because of illumination variability. After unmixing, the fractional cover allocated to the shade endmember was proportionally reallocated to the other materials in each pixel in order to estimate subpixel land cover [41]. This process is known as shade normalization.
After preliminary endmember selection was completed, the endmembers were manually inspected to ensure that they were pure representations of their class. Several endmembers were removed from the library because they represented a combination of multiple land cover types, which caused the models to overestimate GV fractions. In order to guide this process, we compared the modeled GV fractions for 2018 with ground reference data of GV cover from the field survey transects [53,54]. The selected endmembers maximized the agreement between the two data sets. The creation of the ground reference data is described in Section 2.3. In total, 18 GV endmembers, 2 soil endmembers, and 1 NPV endmember were manually removed from the library to ensure that it only contained pure endmembers.
The final MESMA spectral library contained 29 endmembers, including 8 GV endmembers, 16 soil endmembers, 4 NPV endmembers, and photometric shade. The library for the 2007 image also contained the four ash endmembers. All plant species and all soil types were represented in the final library. The spectral libraries were convolved to the spectral resolutions of Thematic Mapper (TM), Enhanced Thematic Mapper Plus (ETM+), and Operational Land Imager (OLI) before unmixing.

Spectral Unmixing
MESMA was performed using VIPER Tools version 2.0 [55], an extension for ENVI (Harris Geospatial Solutions, Broomfield, CO, USA). All Landsat images were unmixed using the spectral library described in the previous section. The potential models for each pixel contained between two and four endmembers from different classes (i.e., GV, NPV, soil, and shade). The selected model minimized the RMSE between the modeled and measured spectra. Several constraints were placed on MESMA in order to ensure that the selected models were physically reasonable: endmember fractions were constrained between 0 and 1, shade fractions were constrained between 0 and 0.8, and the maximum allowable RMSE was 0.025 [56]. All fractions were shade normalized after unmixing. The shade-normalized GV fractions were extracted for each image date in order to create a GV trajectory for each pixel inside of the burn scar.

Assessment of Remote Sensing Products
The selected MESMA model maximized the agreement between the 2018 GV fractions and ground reference data of GV cover from the field survey transects. RMSE and Pearson's correlation coefficient were used to quantify model fit. The units of the GV fractions are the same as the units of the field observations (i.e., percent cover of vegetation scaled from 0-1), so the two data sets could be directly compared. NBR and RdNBR produce unitless values, so they cannot be directly compared with field observations. We assessed their ability to model changes in vegetation cover by calculating the correlation between the 2018 NBR values and the ground reference data. While RdNBR was ultimately used to quantify post-fire vegetation regrowth, RdNBR values do not represent land cover at one point in time; RdNBR is a multi-temporal metric that measures changes in NBR values between two image dates. Thus, it is appropriate to use NBR, but not RdNBR, to assess model fit. We also calculated the correlation between the MESMA and NBR products to quantify the strength of their agreement.

Vegetation Field Surveys
In order to characterize the land cover in the remote sensing images, we conducted vegetation field surveys in the summer of 2018. These data enabled us to analyze how fine-scale variability in land cover influences remote sensing signals of post-fire vegetation regrowth. Species cover data were collected from eighty-two 45 m transects located throughout the burn scar. The transect locations were selected based on six criteria: 1. in an area of high burn severity (2007 dNBR > 600) 2. homogenous vegetation composition 3. homogenous slope and aspect 4. minimum 40 m from a road to minimize the impacts of the road on surrounding vegetation (e.g., invasive grass) 5. maximum 150 m from an accessible road or trail 6. estimated slope less than 35 degrees Because species composition is an important control on post-fire recovery behavior, we also stratified the survey sites by dominant vegetation type using CALVEG data. Approximately half of the sites were dominated by mixed chaparral (e.g., Ceanothus spp., Arctostaphylos spp., and A. fasciculatum), a quarter of the sites were in oak woodlands, and a quarter of the sites were located at the interface of chaparral shrubland and mixed conifer forests. The sites located near the forests did not have conifer overstories and were primarily covered by conifer saplings, shrubs, or herbaceous species.
Species cover data were collected from the 45 m transects using a digital camera mounted on a 5 m tall boom (Wonder Pole, Salem, OR, USA). The camera was adjusted so that the lower border of each image aligned with the transect tape and the upper border of the image extended 5 m downslope of the transect tape. Photographs were taken every 5 m along the transect tape, totaling 10 photographs for each transect. A small number of survey sites had no slope. For those sites, photographs were taken on both sides of the transect tape at five evenly spaced points along the transect. Every plant species within each transect was identified and extensive notes were taken recording the locations of each species within the transect. Grass species were not identified individually and were assigned to a single taxonomic group.
To estimate the fractional cover of different species within each transect, points were sampled from each photograph using a regular sampling scheme [57,58]. In the laboratory, we overlaid a 7 × 6 point grid onto each photograph and labeled the species cover at each grid point. Non-biotic land cover types (e.g., rock and charcoal) were also labeled accordingly. In total, there were 42 grid points for each photograph and 420 grid points for each transect. The proportion of grid points assigned to each species was used as an estimate of its fractional cover within the transect.
GV, NPV, and soil fractional cover were estimated using the same method. These estimates were used to calibrate endmember selection for the spectral mixing models, as described in Section 2.2.3. Because these cover types vary at a finer spatial resolution than species cover, a 10 × 10 grid was overlaid onto each photograph, totaling 1000 grid points for each transect. The proportion of grid points assigned to each cover type was used as an estimate of its fractional cover within the transect.
The species cover data were analyzed using principal components analysis (PCA) to identify species assemblages across the transects. We only analyzed cover data for dominant species, defined as species that had a fractional cover greater than 25% in at least one transect [59]. PCA was performed on cover data for all of the transects, even if a species was not dominant in a particular transect.
We also performed two additional analyses to quantify fine-scale variability in land cover. We calculated the mean percentage of land cover that was identical between all pairs of transects dominated by each species. This enabled us to test whether stands dominated by the same species have consistent land cover across the landscape. We also calculated the percent cover of dominant species in each transect. Species exhibit different physiological responses to disturbance, so compositional heterogeneity within and across vegetation stands may inhibit our ability to attribute remote sensing signals to particular plant types.

Conifer Canopy Cover
The field survey protocol could not quantify the canopy cover of mature conifer trees, so conifer canopy cover was determined by manually interpreting high-resolution imagery. Conifer stand boundaries were extracted from CALVEG. We extracted all 756 stands within the burn scar dominated by the "Mixed Conifer-Pine" alliance. This alliance was selected because its stands are located in areas of known conifer mortality after the fire (Dr. Nicole Molinari, personal communication), and because it contains most of the conifer species found within the burn scar. The stands are primarily located in the vicinity of Big Pine Mountain in Los Padres National Forest.
The conifer stand boundaries were overlaid onto high-resolution imagery available through the National Agriculture Imagery Program and Google Earth (Google, Mountain View, CA, USA). We utilized images captured from 2002 to 2018. Percent canopy cover from before and after the fire was manually assessed and recorded for each stand. Canopy cover was classified into four classes: no trees (0% canopy cover), sparse (1-25% canopy cover), open (25-75% canopy cover), and dense (75-100% canopy cover). Thresholds in 25% increments were selected to minimize subjectivity and error during the manual interpretation process [60]. They are similar to thresholds used by Canada's National Forest Inventory [61] and previous studies [62].

Recovery Trajectories
RdNBR values were computed for each year from 2007 to 2018. GV fractions were computed for each year from 2000 to 2018. Values for the remote sensing products were extracted from the centroid of each 45 m transect. In order to determine whether different species exhibited distinct recovery trajectories, the transects were grouped by dominant species. Mean RdNBR and GV trajectories were then computed for each dominant species. Median RdNBR values and GV fractions were also extracted from the 756 stands of mixed conifer forest. The stands were then grouped based on pre-fire and post-fire canopy cover and the mean trajectories were computed for each group.

Land Cover
Species cover data for the chaparral shrubland were collected from 82 vegetation transects between June and August of 2018. Sixty-nine plant species and four non-biotic land cover types were observed within the transects. There were between 3 and 19 species in each transect, with a mean diversity of 8.4 species. In total, there were 14 dominant species, which had more than 25% cover in at least one transect.
PCA of the transect data for the dominant species revealed that the ground cover of one dominant species often had low or negative correlations with the ground cover of other dominant species (Figure 3). This was particularly true for the four most prevalent species: Q. berberidifolia, Ceanothus cuneatus, A. fasciculatum, and Arctostaphylos glandulosa. This suggests that these species are diagnostic of the different vegetation alliances that exist across the landscape. These chaparral alliances are well documented in literature, and their distributions respond to fine-scale variability in environmental conditions. A. fasciculatum dominates in many arid locations. Obligate-seeding Ceanothus spp. (e.g., C. cuneatus) and Arctostaphylos spp. dominate on ridge lines and arid south-facing slopes. Q. berberidifolia dominates on mesic north-facing slopes [63,64].
We also calculated the mean percentage of land cover that was identical between all pairs of transects dominated by each species (Table 1). This provided a metric of compositional similarity that revealed whether stands dominated by a given species had consistent land cover across the landscape. For the four most prevalent species, the compositional similarity ranged from 54.27 to 60.84%.
Nearly half of the land cover within those transects was different, even though the dominant species were the same. We also calculated the dominant species cover in each transect (Table 1). For the four most prevalent species, the mean dominant species cover ranged from 32.53% to 47.94%. The remainder of the land cover belonged to non-dominant land cover types. P. coulteri and Quercus agrifolia were excluded from these analyses because they only dominated one transect each.  Based on these findings, Q. berberidifolia, C. cuneatus, A. fasciculatum, A. glandulosa, Ceanothus leucodermis, and Cercocarpus betuloides were selected for remote sensing analysis. The first five species were selected because they were the most prevalent species across the landscape. C. betuloides was selected because previous studies reported that it had extremely low mortality during the 2011-2017 California drought [65]. The six selected species also represent different functional responses to fire. After wildfires, chaparral species either resprout from surviving underground plant tissue (obligate resprouters), recruit from seeds (obligate seeders), or both (facultative seeders) [66]. Q. berberidifolia and C. betuloides are obligate resprouters. C. cuneatus is an obligate seeder. A. fasciculatum, A. glandulosa, and C. leucodermis are facultative seeders [67][68][69].
In order to assess fire impacts on the mixed conifer forest, pre-fire and post-fire canopy cover were manually labelled for all 756 CALVEG stands dominated by the "Mixed Conifer-Pine" alliance ( Figure 4). Saplings could not be resolved in the imagery, so the classifications were based on the fractional cover of mature trees. In total, 254 stands (34.6%) changed canopy cover class as a result of the fire. 481 stands (65.4%) did not change canopy cover class. Sixty-two stands (8.4%) experienced complete conifer mortality. No stand had more canopy cover after the fire than it did before the fire. Twenty-one stands were classified as having no trees before the fire. These stands appear to have been misclassified by CALVEG and were excluded from further analysis.

Assessment of Remote Sensing Products
The GV fractions for 2018 achieved strong agreement with the ground reference data of vegetation cover from the transect photographs (RMSE = 0.096, n = 80; Figure 5). The 2018 GV fractions (r(78) = 0.552, p < 0.001) and 2018 NBR values (r(78) = 0.555, p < 0.001) had nearly identical correlations with the ground reference data. The two data sets also exhibited very strong linear relationships with each other. We calculated the correlation between the GV fractions and NBR values for each year from 2000 to 2018. The annual correlation coefficients ranged from 0.9 to 0.98. The mean annual correlation coefficient was 0.952. The four lowest correlation coefficients occurred in the four years after the fire.

Relative Differenced Normalized Burn Ratio
Mean 2007 RdNBR values for the six chaparral species ranged from 1275 to 1423. All six species recovered rapidly between 2007 and 2011, but recovery stalled between 2013 and 2016 ( Figure 6). These dates coincided with a period of extreme drought [70]. C. leucodermis generally had the highest mean RdNBR values throughout the time series, followed by A. glandulosa. C. betuloides had the lowest mean RdNBR values after 2011, suggesting that it was more resilient during drought than the other species. The trajectories for A. fasciculatum, C. cuneatus, and Q. berberidifolia had intermediate RdNBR values and were similar throughout the time series. The RdNBR values for the conifer stands were inversely correlated with the degree of canopy loss following the fire (Figure 7). In stands that experienced complete canopy loss (n = 62), the mean 2007 RdNBR was 1114. Partial canopy loss leading to sparse (n = 108) or open (n = 84) canopies resulted in RdNBR values of 701 and 351, respectively. For stands that did not experience canopy loss (n = 481), the mean 2007 RdNBR value was 302. The degree of canopy loss also influenced the shapes of the recovery trajectories. In stands that experienced complete canopy loss (i.e., no mature trees survived), the RdNBR values decreased rapidly during the first four years of recovery and then more slowly after that. Other stands had more linear trajectories. Average 2018 RdNBR values were relatively similar across all combinations of pre-fire and post-fire canopy cover. The mean values for the different combinations only ranged from 150 to 272.

Green Vegetation Fractions
The six chaparral species had similar GV trajectories before the fire ( Figure 8). The greenness of all six species fluctuated because of inter-annual variability in precipitation. There were substantial declines in vegetation greenness in 2002 and 2004, which were years that had below average rainfall in Santa Barbara County. C. leucodermis had the highest average GV fractions from 2001 to 2004.
The fire destroyed practically all green vegetation in the chaparral shrubland: mean 2007 GV fractions for the six species ranged from 0 to 0.01. GV fractions increased rapidly in the four years after the fire, but the recovery of C. leucodermis lagged behind that of the other species. GV fractions generally declined between 2011 and 2016, presumably as a result of the drought. C. betuloides had the highest GV fractions during the drought years, while GV fractions for the other species were substantially lower. C. leucodermis and A. glandulosa exhibited the least change during the drought, but the post-fire GV fractions for both species were much lower than their pre-fire values. By September 2017, GV fractions for four of the six species had recovered to the values observed in the dry pre-fire years of 2002 and 2004. In certain cases, GV fractions recovered to the values observed in wetter years such as 2003, but this pattern is not consistent across the data set. C. leucodermis and A. glandulosa did not recover to pre-fire values.  In the conifer stands, GV fractions declined between 2006 and 2007 as a result of the fire and increased between 2007 and 2018 as vegetation regrew (Figure 9). Conifer stands that experienced partial canopy loss exhibited linear recovery trajectories, while stands that experienced complete canopy loss exhibited logarithmic recovery trajectories. In stands that had small amounts of canopy cover after the fire, the post-fire GV signal must have been strongly influenced by shrubs and herbaceous species in the understory.
The maximum possible canopy loss occurred in conifer stands that had dense canopies before the fire and no trees after the fire. In those stands (n = 26), 2007 GV fractions were 92% less than their pre-fire median and 2018 GV fractions were 24% less than their pre-fire median. In stands that did not change canopy cover class (n = 481), 2007 GV fractions were 20% less than their pre-fire median and 2018 GV fractions were 12% less than their pre-fire median. 2018 GV fractions were less than pre-fire median GV fractions in 545 of the 756 stands that we analyzed.

Discussion
Field survey data are essential for understanding and evaluating remote sensing signals of disturbed ecosystems [71]. Likewise, remote sensing data can help interpolate field observations across an entire landscape. In this study, we combined field survey data with remote sensing imagery to retrieve post-fire recovery trajectories for different plant species and functional types. We also used the field survey data to identify sources of variability in the recovery trajectories that are not apparent from the imagery data alone.

Remote Sensing Products
The RdNBR and GV trajectories provide different insights into the post-fire recovery of burned vegetation. RdNBR measures the normalized difference between pre-fire and post-fire NBR values, so it represents the change in vegetation cover compared to pre-fire conditions. GV fractions are an absolute measure of vegetation cover on a single image date. We found that the NBR values and GV fractions had nearly identical correlations with the estimates of GV cover from the transect photographs. This finding suggests that they are equally effective at mapping post-fire vegetation cover. However, NBR and RdNBR produce unitless values, while GV fractions are physically meaningful values that can be compared with observations by field ecologists [15]. Unitless spectral indices can be combined with field survey data to model vegetation cover and other ecosystem traits [72], but it is not clear when such models can be transferred from site to site without additional ground reference data [73]. The NBR correlation with vegetation cover that we observed is similar to the correlation reported by Morresi et al. [74], but additional work is needed to validate the consistency of these relationships. Few studies have assessed the ability of NBR and other spectral indices to map fractional cover of vegetation.

Chaparral Species
The RdNBR and GV trajectories for the chaparral species revealed rapid vegetation regrowth from 2007 to 2011. The trajectories for the species were initially very similar, but diverged over time as vegetation regrew. The exception to this trend was C. leucodermis, whose recovery lagged one year behind that of the other species. Five of the six species exhibited substantial declines in GV fractions during the 2011-2017 California drought. By September 2017, GV fractions for four of the six species had recovered to the values observes in the dry years of 2002 and 2004, but the GV fractions did not consistently recover to the values observed in wetter years. We propose three factors that may account for the observed variability between species: post-fire ecological succession, the physiological responses of species to fire and drought, and inter-annual variability in precipitation.
The initial similarity between the trajectories may be attributable to a succession of ephemeral herbaceous species that germinate after fires. Immediately after a fire, the pre-fire dominant shrubs only comprise a small fraction of the land cover. It takes at least four to five years for shrubs to regain their pre-fire cover, and a succession of herbaceous species dominates vegetation patches while the shrubs regrow [21,75]. These ephemeral species likely influenced the remote sensing estimates of vegetation cover, and they may have obscured differences between the RdNBR and GV trajectories for the chaparral species. The herbaceous plants decline in cover as the shrubs regain dominance [21,75], so their influence on remote sensing signals diminishes over time. The trajectories began to diverge around 2010, and they exhibited distinct trends in subsequent years.
The physiological responses of species to fire and drought were also important drivers of variability in the recovery trajectories. The composition of patches dominated by obligate resprouters (e.g., Q. berberidifolia and C. betuloides) changes little after fires because resprouting plants emerge from established root systems and are extremely resilient. Resprouting plants also benefit from having established root systems because it decreases their vulnerability to water stress during post-fire recovery [76,77]. Obligate and facultative seeders (e.g., C. cuneatus, A. fasciculatum, A. glandulosa, and C. leucodermis) recruit from dormant seed banks and usually experience substantial population increases after fires. As a result, patches dominated by seeding species change in both structure and species composition for several years after the disturbance [21]. Unlike resprouting plants, seedlings must establish new root systems after a fire, which makes them more vulnerable to water stress until they reach maturity [76][77][78].
The physiological differences between the species were not expressed in the recovery trajectories until the onset of the drought in 2011, when the trajectories began to exhibit sensitivity to inter-annual variability in precipitation. C. cuneatus was the most sensitive to precipitation variability throughout the time series. Its GV fractions declined substantially during the drought years of 2011-2017 and increased rapidly during the wet growing seasons of 2004-2005 and 2016-2017. Its RdNBR values also increased during the drought and decreased rapidly between 2016 and 2017. It is likely that C. cuneatus was the most sensitive to drought because it is an obligate seeder, and seedlings have less extensive root systems than resprouting plants. C. betuloides, an obligate resprouter, was the most resilient species during the drought. It had the highest GV fractions and lowest RdNBR values after 2011. C. betuloides benefits from an established root system and is highly resistant to drought-induced cavitation, which enables it to thrive under dry conditions [79]. Our finding is also consistent with that of Venturas et al. [65], who reported that C. betuloides did not experience any mortality during the drought. A. glandulosa and C. leucodermis exhibited the least change during the drought, but their post-fire GV fractions were substantially lower than their pre-fire GV fractions. As a result, they generally had the highest RdNBR values throughout the time series. A. glandulosa is well adapted to arid microclimates, which may explain its relative stability during the drought [80]. The explanation for the drought response of C. leucodermis is less clear. C. leucodermis exhibited the slowest recovery, so its drought response might have been affected by the timing of its regrowth relative to the onset of the drought.
In 2018, the RdNBR values for all six species remained above zero, indicating that they did not recover to median pre-fire conditions 11 years after the fire. The GV fractions for four species recovered to the values observed in the dry years of 2002 and 2004, but C. leucodermis and A. glandulosa did not recover to pre-fire GV fractions. Several studies have used Landsat time series to analyze the post-fire recovery of chaparral ecosystems, but no studies have examined the impact of a prolonged drought on recovery trajectories [9,22,81,82]. The previous studies reported full recovery in 7-15 years, which may suggest that the 2011-2017 California drought delayed the recovery to pre-fire conditions within the Zaca burn scar. For example, Peterson and Stow [81] reported that GV fractions stabilized approximately seven years after a fire. Hope et al. [22] reported that the normalized difference vegetation index (NDVI) returned to pre-fire values in 10 years. Storey et al. [9] reported that the enhanced vegetation index, soil adjusted vegetation index (SAVI), and modified SAVI returned to pre-fire values in four years. They also reported that NBR and NDVI returned to pre-fire values in 9 and 12 years, respectively. McMichael et al. [82] reported that leaf area index stabilized approximately 15 years after a fire. That being said, there are important methodological differences between the previous studies and our study. The previous studies did not stratify their samples based on burn severity, while our study only examined high burn severity sites. Thus, burn severity may also account for some of the difference in the reported recovery times.

Mixed Conifer Forest
We also examined the 756 stands of mixed conifer forest within the Zaca burn scar. There was significant variability in stand density before and after the fire, so we grouped the stands based on their pre-fire and post-fire canopy cover. This enabled us to retrieve RdNBR and GV trajectories for varying degrees of canopy loss and vegetation type conversion. Two distinct patterns emerged among the recovery trajectories. First, as the extent of canopy loss increased, the 2007 RdNBR values increased and the magnitude of change between the 2006 and 2007 GV fractions also increased. Second, as post-fire canopy cover decreased, the shapes of the trajectories transitioned from being linear (high post-fire canopy cover) to logarithmic (low post-fire canopy cover).
The first pattern highlights a difference between the recovery trajectories for chaparral and conifer species. While chaparral stands are comprised of low-stature vegetation that is prone to intense crown fires, conifer forests have pronounced vertical structure that protects foliage from being scorched by flames. Low severity fires generally burn foliage that is close to the ground. As a fire's severity increases, it burns foliage higher in the canopy [83]. Foliage above the maximum scorch height will still be green, so conifer stands with surviving trees can have non-zero GV fractions immediately after a fire. The 2007 GV fractions were approximately proportional to post-fire canopy cover, and the RdNBR values exhibited similar trends. Some stands that did not experience substantial canopy loss had higher than expected 2007 GV fractions, presumably due to surviving understory vegetation.
The second pattern can be attributed to the influence of understory vegetation on remote sensing signals after the fire. In stands with closed canopies, the recovery trajectories represent changes in the spectral properties of the tree crowns. In areas without closed canopies, the remote sensing signal is strongly influenced by understory cover [84]. Stands with dense canopy cover after the fire exhibited linear increases in GV fractions and little change in RdNBR values. Presumably, these trajectories represent the regeneration of scorched foliage within the tree crowns. Stands with no trees after the fire exhibited logarithmic trajectories that are qualitatively similar to the trajectories of the chaparral shrub species. These trajectories likely represent shrubs and other understory vegetation because the tree crowns were completely destroyed by fire. Stands with intermediate amounts of canopy cover after the fire had recovery trajectories that were neither perfectly linear nor logarithmic. These findings suggest that the conifer and shrub species have distinct recovery patterns, which are caused by their canopy architecture and physiological responses to fire.
Unlike the chaparral species, the conifer stands did not exhibit a drought signal that was directly correlated with annual precipitation. Both the RdNBR and GV trajectories indicated continuous recovery throughout the 2011-2017 drought. Nonetheless, 545 of the 756 conifer stands did not return to their pre-fire values within 11 years, indicating that there was long-term change as a result of the fire or drought or both. The most obvious source of long-term change is the fire-induced mortality of mature conifer trees. Prolonged drought can also reduce growth, cause crown dieback, and inhibit seedling establishment in conifer species [85][86][87]. Based on previous studies of fire effects, the drought likely delayed the return to pre-fire levels of greenness in the conifer stands. Fernandez-Manso et al. [10] retrieved GV trajectories for field survey plots in a Mediterranean forest. Plots with low burn severity recovered in seven years and plots with moderate burn severity recovered in 13 years. Moressi et al. [74] performed a similar analysis using spectral indices. They found that modeled recovery times varied between 8 and 12 years, depending on the forest type, burn severity, and spectral index. Almost none of the stands in our analysis recovered to pre-fire levels of greenness within 11 years, which suggests that fire effects alone do not account for the observed recovery times.

Species Composition
We also used the field survey data to quantify fine-scale variability in species composition within and across vegetation stands. Many studies rely on thematic land cover maps to analyze fire behavior and effects, but the biophysical attributes of vegetation patches vary continuously in both space and time [88,89]. Because the spatial variability occurs at a scale that is finer than the spatial resolution of most sensors, reflectance signals for individual pixels are influenced by numerous plant species and functional types. We observed between 3 and 19 species in each field survey transect, with a mean diversity of 8.4 species. The mean dominant species cover ranged from 30.16 to 52.94%, depending on the species. While classification schemes based on dominant species are simple and easy to understand, non-dominant species inevitably introduce spectral variability into remote sensing signals. If linear mixing is assumed, this variability is approximately proportional to the non-dominant species' ground cover within a pixel. Dominant species accounted for less than half of the land cover in most transects, so non-dominant land cover types were as influential on the remote sensing signals as the dominant species. These findings suggest that the recovery trajectories may represent vegetation alliances dominated by particular species, rather than the dominant species themselves. We also analyzed whether stands dominated by a given species had consistent land cover across the landscape. On average, between 51.39% and 82.69% of the land cover was identical across transects dominated by the same species. This suggests that vegetation alliances are somewhat stable across the landscape, but there is still substantial within-class variability in patch composition.
Temporal variability in species composition also affects remote sensing trajectories. A succession of ephemeral species dominates vegetation patches after fires, as discussed in Section 4.2.1. Even after a vegetation patch has completed its initial recovery, it may not return to the state that was present before the fire. The post-fire distributions of seeding species are a function of fire intensity and the locations of dormant seed banks at the time of the fire [21]. Secondary disturbances such as drought can also affect the plant species and life forms that regenerate within burn scars [6]. Without pre-fire species cover data, it is difficult to determine where and to what extent the long-term species composition changed as a result of the fire. This introduces variability into remote sensing metrics that rely on pre-fire baseline values, such as dNBR and RdNBR. It also affects the interpretation of remote sensing trajectories that include pre-fire years. Few studies have considered temporal variability in species composition as a source of variability in remote sensing trajectories [17]. These factors merit further consideration in literature because remote sensing trajectories are often used to infer the drivers and effects of ecological change [90]. Misattributing remote sensing signals to different vegetation types could lead to erroneous conclusions about ecosystem dynamics.

Methodological Considerations
This study used field survey data to summarize and interpret remote sensing signals of post-fire vegetation regrowth. Certain limitations should be considered when interpreting the results. For example, the field survey data were collected at a different spatial scale than the Landsat spectra. The field survey data were collected from 45 m × 5 m transects, while Landsat spectra were extracted from 30 m × 30 m pixels. Documenting all of the plants along a 45 m transect is substantially more efficient than documenting all of the plants in a 30 m × 30 m area. This is especially true in chaparral ecosystems, where shrubs often form impenetrable canopies. The transects effectively subsampled the Landsat pixels, and we assumed that the subsamples were representative of the Landsat pixels as a whole. To the extent possible, we selected survey sites that had homogenous slope, aspect, and vegetation composition, which bolstered the validity of that assumption. For safety reasons, we also selected survey sites that were located near accessible roads and did not have steep slopes. Large areas of the burn scar are not located near roads, so certain vegetation alliances may be overor underrepresented in the field survey data. The field survey data were also used to calibrate the MESMA model. It is possible that the MESMA model was overfit to a particular sample of ground reference data. However, the MESMA GV fractions exhibited a very strong linear relationship with the NBR values, which were derived from an entirely different modeling framework. This suggests that the field survey data were adequate to calibrate the MESMA model reliably. There may still be some error in the MESMA GV fractions, as they are modeled and not measured values.

Conclusions
Integrating field survey data into remote sensing analyses makes it possible to characterize the remote sensing signals for land cover types that would otherwise be difficult to differentiate. It also makes it possible to identify sources of variability in remote sensing data that are not apparent from imagery alone. In this study, we used time series imagery and field survey data to retrieve post-fire recovery trajectories for chaparral species and stands of mixed conifer forest. In doing so, we answered three research questions.
First, we examined whether NBR and GV fractions were effective for mapping vegetation cover across a burn scar. The selected MESMA model maximized the agreement between the modeled GV fractions and ground reference data of GV cover from the field survey transects. The two data sets achieved strong agreement, indicating that the GV fractions were a reliable estimate of GV cover across the burn scar. The NBR values and GV fractions had very similar correlations with the ground reference data, indicating that they were equally effective at mapping vegetation cover.
Second, we examined sources of variability in post-fire recovery trajectories for chaparral species and stands of mixed conifer forest. For the chaparral species, precipitation variability and species-level physiology were the most obvious sources of variability throughout the time series. The chaparral species exhibited rapid regrowth between 2007 and 2011, but GV fractions for five of the six species declined substantially during the drought. In 2018, RdNBR values for the six species remained above zero, indicating that they had not recovered to pre-fire conditions. For the conifer stands, the shapes of the recovery trajectories were determined by post-fire canopy cover. Stands with high post-fire canopy cover had linear recovery trajectories, while stands with low post-fire canopy cover had logarithmic recovery trajectories. The conifer stands did not exhibit a drought signal that was directly correlated with annual precipitation, but in a majority of stands the GV fractions did not recover to their pre-fire values.
Finally, we examined variability in the structure and composition of chaparral and conifer patches. In the surveyed chaparral patches, dominant species only covered 30 to 53% of the land area. Even within patches dominated by the same species, there was significant variability in species composition across the landscape. Non-dominant species inevitably affected the remote sensing signals, and these effects may not be consistent from patch to patch. In the conifer stands, canopy dieback exposed understory land cover, including shrubs, grasses, and soil. In stands with substantial amounts of canopy loss, the recovery trajectories were strongly influenced by understory vegetation. For both chaparral and conifer stands, these effects vary over time because a succession of ephemeral species populate vegetation patches as the long-term dominant species recover.
Our study expands on previous literature by exploring the opportunities and limitations of using Landsat time series to monitor ecosystem recovery after disturbance. This technique has gained popularity since the Landsat archive became freely available. Few studies have used species cover data to retrieve post-disturbance recovery trajectories for different plant species and functional types. The trajectories we retrieved were sensitive to changing environmental conditions and some physiological differences between the species. The pre-fire and post-fire values were also useful for assessing whether the ecosystem recovered to pre-fire conditions. However, our findings also suggest that there may be fundamental limitations in our ability to infer the traits of individual species using Landsat imagery. Species cover often varies at a scale that is finer than the spatial resolution of spaceborne sensors, so researchers must be cognizant of the impact of non-dominant land cover types on remote sensing signals. These considerations will become increasingly important as hyperspectral spaceborne sensors become operational and researchers attempt more precise retrievals of ecosystem properties.
Our study also reveals that prolonged droughts can slow or even reverse the post-fire recovery of Mediterranean ecosystems, which can lead to long-term changes in ecosystem structure and function [6]. Both fire activity and drought conditions have increased under anthropogenic climate change [7,91], so it is increasingly likely that land managers will be forced to address multiple overlapping disturbances in sensitive ecosystems. Without targeted management and restoration efforts, vulnerable species may disappear from the landscape if they cannot maintain reproducing populations under changing environmental conditions [86].