Monitoring the impacts of severe drought on southern California Chaparral species using hyperspectral and thermal infrared imagery

: Airborne hyperspectral and thermal infrared imagery acquired in 2013 and 2014, the second and third years of a severe drought in California, were used to assess drought impacts on dominant plant species. A relative green vegetation fraction (RGVF) calculated from 2013–2014 Airborne Visible Infrared Imaging Spectrometer (AVIRIS) data using linear spectral unmixing revealed seasonal and multi-year changes relative to a pre-drought 2011 reference AVIRIS image. Deeply rooted tree species and tree species found in mesic areas showed the least change in RGVF. Coastal sage scrub species demonstrated the highest seasonal variability, as well as a longer-term decline in RGVF. Ceanothus species were apparently least well-adapted to long-term drought among chaparral species, showing persistent declines in RGVF over 2013 and 2014. Declining RGVF was associated with higher


Introduction
California has been subjected to a multi-year, extreme drought, with tree ring records indicating that the drought may be the most severe in the last millennium [1]. Both warmer temperatures and greatly reduced wet season precipitation over multiple years [2] have increased stress on vegetation. Long-term, severe drought can cause canopy dieback and mortality in woody plants, resulting from loss of hydraulic conductance, depletion of non-structural carbohydrate reserves, and secondary stress by insect attack [3][4][5][6][7]. Plant senescence, canopy dieback, and mortality result in a reduction of leaf area and fraction of the ground surface covered by green vegetation. Drought impacts will vary by species, however, based on factors such as leaf structure, root depth, and cavitation resistance [8][9][10][11][12][13][14].
Hyperspectral data have capabilities for estimating vegetation biophysical and biochemical parameters [15][16][17][18] and mapping vegetation species [19][20][21][22][23][24]. These complementary capabilities give hyperspectral data strong potential for monitoring drought impacts on vegetation [25,26], but a paucity of hyperspectral time series datasets has prevented the development of drought monitoring applications. Thermal infrared time series data also have considerable potential for monitoring drought due to drought impacts on evapotranspiration and canopy temperature [27][28][29][30]. NASA has proposed the Hyperspectral Infrared Imager (HyspIRI) satellite mission, which would combine a hyperspectral imaging spectrometer covering the visible, near infrared, and shortwave infrared spectral regions (VSWIR) with a multispectral thermal infrared (TIR) sensor [31]. HyspIRI would provide new opportunities for global monitoring of terrestrial ecosystems [32], potentially including seasonal monitoring of drought stress on vegetation.
The California drought has coincided with data collection for NASA's HyspIRI Preparatory Campaign, a campaign designed to acquire data using VWSIR-and TIR-like airborne instruments. Data from the Airborne Visible Infrared Imaging Spectrometer (AVIRIS) [33] and the MODIS-ASTER Airborne Simulator (MASTER) [34] were acquired from an airborne platform over key ecosystems in California in 2013 and 2014, providing a valuable dataset for quantifying the impacts of the persistent drought on vegetation. This effort uses these newly available data to identify the initial impacts of a long-term drought on dominant plant species in the Santa Barbara region of southern California. Seasonal changes in fractional vegetation cover calculated from AVIRIS data and land surface temperature (LST) retrieved from MASTER data reveal differences in how individual species are responding to long-term drought.

Study Area
The study area includes the Santa Barbara Coast and portions of the Santa Ynez Mountains and Santa Ynez Valley (Figure 1). This area spans an elevation range from sea level to 1300 m, with coastal sage scrub and grassland communities dominating at lower coastal elevations, oak forest, and chaparral shrublands at higher elevation, riparian forest in mountain drainages, and oak savanna, coastal sage scrub, and grassland in drier inland areas. Precipitation and soil moisture data collected from a station in the study area [35] demonstrate reduced monthly precipitation and decreasing soil volumetric moisture content (VMC) since the last winter with above average precipitation, 2010-2011 ( Figure 2).

Image Data
Data from AVIRIS and MASTER were acquired over the study area on 19 July 2011, providing a pre-drought baseline. New data acquisitions occurred in April, June, and November/December 2013, and April, June, and August 2014. All data acquisitions occurred on a single date except for November/December 2013 (two flights separated by nine days) and June 2014 (two flights separated by two days). AVIRIS data cover a spectral range of approximately 350-2500 nm with 224 10 nm spectral bands [33]. An AVIRIS surface reflectance product at a uniform 18 m spatial resolution was provided by the Jet Propulsion Laboratory (JPL). Apparent surface reflectance was retrieved using ATREM as described by [36]. Surface reflectance images for each date were registered to an 18 m spatial resolution orthoimage basemap to improve co-registration between dates. Registration spatial error was required to be less than one pixel for all areas containing reference polygons. MASTER has 50 bands covering the visible through TIR [34]. A land surface temperature (LST) product retrieved from the TIR bands at a mean 35 m spatial resolution was provided by JPL. Temperature-emissivity separation required for retrieving LST used a water vapor scaling method described by [37]. Root mean square error of LST retrieval has been estimated at 1.2 K for ASTER data by [38] and at 0.7 K for MASTER data [39]. LST data were registered to the same 18 m orthoimage basemap and were resampled to 18 m spatial resolution using nearest neighbor resampling.

Ground Reference Data
Patches dominated by 13 vegetation species, including one case of two intermixed species, were used to construct a set of reference polygons covering 12 species classes (Table 1). Reference polygons were mapped using methods adapted from [40]. Homogeneous vegetation patches were manually delineated on high spatial resolution imagery, and a high-powered spotting scope was combined with in situ inspection to estimate species dominance within each patch. Reference polygons were required to be at least 75% dominated by a single species (or in one case by intermixed Artemisia californica and Salvia leucophylla) and to have a minimum approximate size of 40 m by 40 m. Reference polygons were mapped within the study area in 2003, 2009, and 2012 ( Figure 1), and were collected from outside of areas burned by recent fires in the Santa Ynez Mountains. Reference data collection is described in further detail in [21] and [41].

Analysis
Linear spectral unmixing (also known as spectral mixture analysis) models reflectance spectra as a linear mixture of endmembers representing pure cover types [42]. Unmixing assigns a fractional composition to each endmember, with green vegetation (GV), non-photosynthetic vegetation (NPV), soil, and shade representing typical endmember types used for vegetation analysis [43]. In a previous adaption of linear spectral unmixing, relative endmember fractions were modeled by unmixing time series data [44].   Changes in GV fraction over the AVIRIS time series were modeled using a relative linear spectral unmixing technique. The July 2011 AVIRIS image was selected as a reference date because it represents a high rainfall year for which shrub and tree canopy spectra would reflect peak growth conditions expressed by high near infrared/red contrast and minimum levels of exposed stems, branches, and soil. For each pixel in the study area, the July 2011 spectrum was used as a GV endmember representing a baseline green vegetation cover. An empirically-determined normalized difference vegetation index (NDVI) threshold of 0.3 was used to mask out non-vegetated pixels ( Figure 3). For each date in 2013 and 2014, each pixel was modeled as a linear combination of its July 2011 spectrum representing the GV endmember, a non-photosynthetic vegetation (NPV) endmember from a lab spectrum of bare branches, and a zero reflectance shade endmember. Singular value decomposition was used to calculate the fractions of GV, NPV, and shade endmembers for the 2013 and 2014 dates, providing a relative GV fraction (RGVF) based on change relative to the July 2011 spectrum. RGVF was shade-normalized to correct for solar zenith effects by dividing by the sum of the non-shade fractions. . An NDVI threshold of 0.3 was empirically determined using spectra extracted from species reference polygons (GV) and from areas with senesced grassland and soil cover (non-GV). Non-GV cover was identified based on appearance in high resolution orthoimagery and manual inspection of spectral shape in the July 2011 AVIRIS image. 95% of GV-type spectra were above the 0.3 NDVI threshold, while only 12% of non-GV-type spectra were above the same threshold.
Use of a relative fraction calculated on a per-pixel basis accommodates spectral variability that occurs by species and over space [45]. Thus RGVF represents change from a previous baseline for each pixel, rather than an absolute measure of vegetation cover. RGVF was used rather than a relative NPV or soil fraction due to the high spectral contrast between green vegetation and other sub-pixel components. This minimizes the impact of spectral confusion between NPV and soil on RGVF [43,46]. Since the GV endmember used to model each pixel is fixed, RGVF avoids problems with best fit endmembers potentially changing over time in a multiple endmember spectral mixture analysis-type GV non-GV model [20,47]. Use of a four-endmember model for RGVF calculation (including a soil endmember in addition to the GV endmember from the 2011 image plus NPV and shade endmembers) was attempted, but similarity between soil and NPV endmembers resulted in unrealistically large changes in soil and NPV fractions from date to date in response to changes in the spectral shape of modeled pixels. Since the NPV endmember best captured the dominant signal of canopy dieback and mortality, it was chosen over the soil endmember for the three endmember model used to calculate RGVF. Increased soil cover associated with canopy dieback and mortality will increase NPV fraction in the three endmember model due to the spectral similarity of soil and NPV. Severe drought impacts on vegetation canopies, representing die-back or mortality, would result in decreased RGVF compared to the 2011 baseline. The reference polygons were used to extract RGVF values for each species. Distributions of RGVF for each date were calculated. Decreased GV cover relative to July 2011 would result in a RGVF below one, while increased GV cover would result in a RGVF greater than one. RGVF was also compared to LST extracted from MASTER data using the same reference polygons.  QUDO demonstrated the largest decreases in RGVF from spring to summer and fall. QUDO is drought deciduous, but RGVF values larger than one in the April images indicate that the green-up and senescence of grasses surrounding and underneath QUDO canopies may be the primary contributor to seasonal trends. QUAG and PLRA, frequently found in riparian areas, had stable RGVF distributions across all dates.

Results and Discussion
Coastal sage scrub species (ARCA-SALE and ERFA) had the most variable RGVF distributions. ARCA-SALE exhibited RGVF values above one in April 2013. These species are drought deciduous, and higher GV fraction in spring 2013 relative to summer 2011 would be expected. ARCA-SALE and ERFA both had RGVF below one in April 2014, indicating that reduced winter precipitation in successive years may have reduced GV cover. ARCA-SALE displayed an unusual RGVF distribution centered below zero in November 2013. These extremely low relative fractions resulted from a total lack of chlorophyll absorption and red edge in ARCA-SALE spectra, indicating complete senescence.
Coastal sage scrub species had the highest LST values, likely associated with lower vegetation cover, lower evapotranspiration, and south-facing slopes typical of coastal sage scrub communities (Figure 7). QUAG, a broadleaf evergreen tree species, anchored the lower end of the LST scale. Roberts et al. [30] previously found that many species in this study area formed unique clusters in the space defined by GV fraction and LST. Species also formed clusters in the space defined by RGVF and LST (Figure 7). In April 2013, these clusters had a RGVF close to one and a wider range of LST values. Over the remainder of 2013 and all of 2014, the relationship between RGVF and LST skewed negative. Coastal sage scrub and chaparral species with large reductions in RGVF relative to July 2011 tended to have higher LST values, while tree species demonstrated smaller changes in RGVF and lower LST values.
Variability in RGVF corresponded to physiological and environmental characteristics. PLRA and UMCA are tree species found in the most mesic areas of the Santa Ynez Mountains, and had minimal changes in RGVF during 2013 and 2014. QUAG is very deeply rooted [12], and also demonstrated minimal change in RGVF. Among chaparral species, ADFA had less extreme declines in RGVF over 2013 and 2014 compared to CECU, CEME, and CESP. ADFA has been found to be more deeply rooted than Ceanothus species [48]. While Ceanothus species may have a higher cavitation resistance than ADFA [14], this did not prevent widespread dieback of Ceanothus species and corresponding declines in RGVF. Coastal sage scrub species are the most shallowly rooted [48] and have low cavitation resistance [13]. These species exhibited the largest seasonal swings in RGVF.
Deeply rooted species that were least affected by drought in 2013 and 2014 are likely to become more severely impacted as the drought continues. In situ data reveal that the continued drought has resulted in soil moisture depletion ( Figure 2). The percent of each month with soil volumetric moisture content (VMC) less than an empirically determined wilting point of 0.08 has declined as the drought has persisted, especially at 46 cm depth (Figure 2c). Winter precipitation has not percolated deep enough to substantially recharge soil moisture at depth. Continuation of the drought beyond 2015 is also likely to worsen drought effects beyond those seen in 2013 and 2014. The lagged response of shrub species to precipitation in February/March 2014 demonstrates that both timing and amount of precipitation are important for determining drought impacts.
This study did not differentiate between plant senescence, canopy dieback, and mortality. Senescence and dieback should produce seasonal cycles in RGVF, although long-term declines in vegetation cover may result from long-term drought stress. Mortality should result in more persistent reductions in RGVF, however, recruitment of herbaceous vegetation following mortality could introduce more seasonal variability. Hyperspectral VSWIR data have demonstrated the ability to map dominant plant species [19][20][21][22][23][24]30], which could be used to separate senescence and dieback (persistence of the same species) from mortality and replacement of one dominant species with another. Dominant species cover mapped using hyperspectral data could be used to extend these results beyond a limited number of field-assessed reference polygons. In doing so, variation in drought impacts on individual species may be found to vary by elevation, aspect, slope, and soil type. Accuracy of species mapping, especially in areas where no single species dominates, will need to be assessed to validate trends in RGVF and LST over larger areas. RGVF should be compared to ground measures such as leaf area index [49] and canopy cover to determine how canopy structural changes correspond to drought-associated decreases in RGVF. The HyspIRI Preparatory Campaign acquired AVIRIS and MASTER data over an area totaling approximately 130,000 km 2 . Due to the limited availability of pre-drought baseline imagery, our study area was limited to a much smaller 650 km 2 . A satellite mission could provide longer-term monitoring of drought conditions, leveraging the complementary abilities of VSWIR data to measure vegetation condition and TIR data to measure LST. Models combining VSWIR and TIR data may be useful for constraining reductions in evapotranspiration resulting from drought [50]. The capabilities of a combined hyperspectral VSWIR and TIR mission for mapping species and monitoring vegetation drought impacts could be used to investigate species response to drought along environmental gradients.  Table 1. Density is the kernel density estimate calculated using the "density" function in R statistical software (http://www.R-project.org/).  Table 1.

Conclusions
RGVF provides a method for measuring per-pixel change in GV fraction over time relative to a baseline image. Change in RGVF and LST were examined for dominant plant species in the second and third years of a continuing, severe drought. Species showed a range of responses to drought stress, corresponding with physiological and environmental characteristics. Deeply rooted tree species and tree species found in mesic areas showed the least change in RGVF relative to a pre-drought image. Coastal sage scrub species had the highest seasonal variability during 2013 and 2014, and both chaparral and coastal sage scrub species demonstrated long-term decline in RGVF. If the drought persists, cumulative, more severe impacts can be expected across a wider range of species. Hyperspectral VSWIR and TIR data hold promise for monitoring drought impacts. Improved knowledge of species-specific drought impacts may improve understanding of changes in ecosystem structure and function resulting from drought-induced mortality of select species [51][52][53].