Detection of drought-induced blue oak mortality in the Sierra Nevada Mountains, California

. From 2011 to 2016, California experienced a millennial-intensity drought, generating high levels of tree mortality. Remote sensing has been used to monitor the long-term impacts of drought; how-ever, discriminating dead from live trees in arid and semiarid deciduous woodlands is challenging. The goals of this study were to assess and map the spatial patterns of drought-induced tree mortality in a blue oak ( Quercus douglasii ) woodland, a highly drought-tolerant species forming savannas along the lower foothills surrounding California ’ s Central Valley. Airborne hyperspectral imagery was used to identify the most important wavelength regions predicting drought-induced blue oak mortality. The best metric to predict canopy stress was a normalized ratio using the spectral bands 937.53 and 1100.08 nm with a correla-tion with tree mortality of R 2 = 0.83. The image prediction of mortality for nine ﬁ eld plots found that 16 of 98 trees died (17.9%) during the drought. We further evaluated tree mortality in 82 independent plots, and we found the greatest image predictive accuracy for tree mortality between 1% and 10%. When applied at the landscape level, the regression-based index found mortality ranged from 1% to more than 51% in oak stands with an average mortality of 10% over the entire study region. In addition, tree mortality at landscape level showed higher tree mortality in blue oaks on south-facing aspects probably because of higher insolation rates and in sites with low potential for accumulated drainage.


INTRODUCTION
Vegetation health is closely linked to the hydrological cycle, and intensified water deficits can cause ecosystems to flip from carbon sink to carbon source, as vegetation dieback changes the ratios of assimilation and decomposition (Paw et al. 2004, Ma et al. 2007, Huang et al. 2010, Hicke et al. 2012, Earles et al. 2014, Fisher et al. 2017). If changes in the water cycle are longlasting (from seasonal to inter-annual), they can alter biogeochemical cycles (Guardiola-Claramonte et al. 2011, Michaelian et al. 2011, and episodes of drought and heat can induce forest mortality (e.g., Allen et al. 2010).
Climate change predictions indicate increased risk of droughts with hotter temperatures and increased variability of precipitation (Allen et al. 2015). These changes are expected to increase forest and woodland mortality, especially in semiarid regions, because conditions are at or near tipping points for changes in ecosystem composition and function (van Mantgem et al. 2009, Allen et al. 2010, Weed et al. 2013, Trenberth et al. 2014. Extended drought impacts to vegetation begin with changes in physiological processes related to stomatal closure and reductions in transpiration and assimilation (Burkett et al. 2005, Hicke et al. 2012, Earles et al. 2014, Brodribb et al. 2020, and if persistent, can be followed by changes in structure, for example, reduced leaf size, loss of leaf area, stem dieback (Fensham and Holman 1999), mortality, and sudden collapse (Breshears et al. 2005, Matusick et al. 2013. These changes can cause ecosystem-scale shifts in species distributions (Mueller et al. 2005, Mueller et al. 2005, Engelbrecht et al. 2007, Bates et al. 2008, Willson et al. 2008) and altered plant structure and functions (Guardiola-Claramonte et al. 2011) ranging from declining leaf area index and plant density to species replacement and conversion to a different functional type (e.g., forest to grassland; Aguiar et al. 1996, Jacobsen and Pratt 2018, McLaughlin et al. 2020. Interest in extreme drought events (e.g., McDowell and Sevanto 2010, Allen et al. 2015), in particular hot droughts (AghaKouchak et al. 2015, Diffenbaugh et al. 2015, Teuling 2018, has focused on improving our understanding of forest traits associated with resilience and tolerance, on conditions that lead to large scale tree mortality (Allen andBreshears 1998, van Mantgem et al. 2009), and on ecosystem changing events, such as occurred in the southwestern U.S. drought in 2002-2003 that covered 12,000 km 2 (Breshears et al. 2005).
Plants in Mediterranean climate ecosystems are adapted to survive hot and dry summers when low soil moisture causes slow growth or induces partial or full dormancy. Thus, these ecosystems may provide insights into plant mechanisms of drought tolerance and endurance. California recorded a drought of millennial impact from December 2010 to March 2016 ( Fig. 1; Griffin and Anchukaitis 2014, Jezdimirovic and Hanak 2016, Williams et al. 2020, and such conditions are predicted to become more common in California by midcentury (Ullrich et al. 2018). The unprecedented combination of low precipitation and high air temperatures, with 2014 and 2015 being the hottest years in the state's history (Fettig et al. 2019b), generated an extreme level of tree mortality in conifer forests throughout the state, estimated at 129 million trees between 2010 and 2017 from the combined effects of drought, extreme temperature, and associated insect kill such as bark beetles (USFS 2017a, Fettig et al. 2019b). In this study, we focused on the use of hyperspectral imagery to track the responses to drought of blue oak (Quercus douglasii), an endemic tree that forms the extensive woodlands and savannas in California, covering about 12,000 km 2 (Sandiford 1997, Keeley 2002. Unlike conifers that are attacked by bark beetles and other insects when weakened by drought, blue oak are not subject to tree killing insect attacks and mortality would be a direct response to drought and extreme temperatures.
Remotely sensed monitoring of drought impacts and mortality Most regional-scale remote sensing studies of drought impacts have used spectral indices for leaf pigments or water content from multispectral data such as Landsat or MODIS (Fraser and Latifovic 2005, Anderson et al. 2010, Meigs et al. 2011). Fraser and Latifovic (2005 used NDVI and the ratios NIR/Red and NIR/SWIR to estimate tree mortality in a coniferous forest in Quebec, Canada. Anderson et al. (2010) used a combination of vegetation indices and climatological data to explore the relationship between remote sensing metrics and tree morality in a tropical forest. They found positive relationships with EVI and a negative relationship with NDWI. Byer and Jin (2017) assessed tree mortality in Sierra Nevada forests using time series of NDVI, NDWI, and EVI from MODIS data. However, these indices miss important plant physical and physiological processes in other wavelength regions (Blackburn 1998, Asner and Martin 2009, Schaepman et al. 2009). Airborne imaging spectroscopy provides contiguous narrow-band spectral information across the solar spectrum, often at high spatial resolution (Asner 1998, Ustin et al. 2004. Asner et al. (2016) used airborne imaging spectroscopy, LiDAR data, and satellite models trained on the spectroscopy data to estimate the water lost from California's forest ecosystems over the drought years between 2011 and 2015 and to detect synoptic spatial patterns in the decline of canopy water content. Numerous studies have demonstrated the capabilities of imaging spectroscopy data to evaluate vegetation injury (Zarco-Tejada et al. 2002, Coops et al. 2003, White et al. 2007) but few studies have discriminated dead and live vegetation (Meddens et al. 2011).
The spectra of live vegetation are well understood, and photon absorptions in the visible part of the spectrum and water in the near-infrared (NIR) and short-wave infrared (SWIR) region of the spectrum are well characterized. Dead stems and leaves scatter photons across the spectrum lacking strong absorptions over most of it (Ustin and Jacquemoud 2020). Because these signals are not distinctive, they are difficult to identify. Our research hypothesis is that, despite this, the difference between the spectral signals that discriminate dead (no green foliage) and live trees (with green foliage) will exceed canopy changes caused by other factors, such as those described by Macomber and Woodcock (1994), including high spectral variability in tree populations and understory components, spectral distortions caused by topographic position, and/or atmospheric effects (water vapor, aerosols, dust). Our assumption is that temporal changes observed in an arid woodland's spectral responses during the drought represent changes in structural canopy traits that can be correlated with the percentage of tree mortality.
The goals of this study were to (1) determine how well imaging spectroscopy data can be used to detect tree mortality in a sparse blue oak savanna ecosystem by identifying the most important wavelength regions that predict drought-induced blue oak mortality and (2) assess the spatial patterns of drought-induced tree mortality in a blue oak savanna in California.

Study site
The study site is the U.S. Forest Service's San Joaquin Experimental Range (SJER; 37°5 0 N, 119°43 0 W) located in the western foothills of the southern Sierra Nevada Mountains, California, USA. The site covers an area of 4994 ha, having elevations between 210 and 520 m above sea level. Summers are hot and dry with mean temperatures between 24°and 27°C. Winters are cool and wet with mean temperatures ranging from 4°to 10°C. Annual precipitation is around 486 mm and occurs between October/November and April/May (USFS 2017b (Griffin and Critchfield 1972). It is one of California's most xeric-tolerant oak species. The trees remain physiologically active at −4 MPa in the Coast Range (Griffin 1973, Weitz 2018) and at −6 MPa predawn water potentials in the Sierra Nevada Mts. (Baker et al. 1981, Xu andBaldocchi 2003). Blue oaks typically exhibit a slow decline in stomatal conductance and evapotranspiration and maintain traits related to photosynthesis, such as Vcmax, quantum yields, light compensation points, total chlorophyll and nitrogen, maximum assimilation rates, and lower dark respiration rates into late summer (Callaway and Mahall 1996, Matzner et al. 2003, Xu and Baldocchi 2003, Kueppers et al. 2005. Additionally, blue oak has physiognomic traits that enhance tolerance such as phenological timing, plasticity in leaf size, high leaf mass area, drought deciduousness, and heterogeneous root distributions with both deep and shallow roots that allow it to occupy a wide range of soil conditions (Baker et al. 1981, Knops and Koenig 1994, Callaway and Mahall 1996, Kueppers et al. 2005. These traits are considered to have a genetic basis and include acclimation from antecedent weather conditions (Garcia et al. 2010, Vincente-Serrano et al. 2012, Gouveia et al. 2017. Blue oak can survive several years of low rainfall but if drought is prolonged, it will eventually experience canopy damage (dieback) or death. Blue oak is one of the longest-lived oak species with tree ring ages of 500 years common and some exceptional ages up to 600 and possibly 700 years (Stahle et al. 2013). Insects and grazing are mostly focused on seeds and seedlings and with little impact on mature trees. Therefore, wildfire, with an estimated 25.2 years of return interval (McClaren and Bartolome 1989, Sandiford et al. 2012, was not a primary cause of tree mortality, at least before fire was largely suppressed.

Data
Imaging spectroscopy data.-NASA's Airborne Visible/Infrared Imaging Spectrometer (AVIRIS), onboard NASA's ER-2 aircraft, was collected by the Jet Propulsion Laboratory (JPL) during the NASA HyspIRI summer flight campaigns in 12 June 2013 between 10:42 and 10:57 Pacific Standard time (PST), 3 June 2014 between 12:16 and 12:31 PST, and 2 June 2015 between 10:41 and 10:58 PST, and for this study, the subtraction between 2013 and 2015 was used to evaluate the impact over the drought period. The larger study area (blue box on Fig. 2) was part of a NASA research program in preparation for the Hyperspectral Infrared Imager (HyspIRI) mission (National Research Council [NRC] 2007, Lee et al. 2015.
AVIRIS collects raw data that is calibrated by the Jet Propulsion Laboratory into level 1B surface radiance measurements in 224 spectral bands in the solar spectrum between 370 and 2500 nm with a nominal bandwidth of 10 nm and a spatial resolution of 18 m (Green et al. 1998). Several bad bands were identified by JPL and removed, leaving 193 bands that were retained in the Atmosphere Removal Algorithm (ATREM) model to predict Level 2 surface reflectance (Boardman 1999, Thompson et al. 2015. The bad bands might have been noisy because they occur in atmospheric water bands or are wavelengths of low detector sensitivity and/or in bands that measured low photon flux passing through the atmosphere.
The AVIRIS image was radiometrically calibrated and geographically rectified by JPL. A topographic correction was applied to the reflectance values for slope and aspect.
We secondarily applied the model proposed by Soenen et al. (2005), which is based on the sun-canopy-sensor (SCS) correction developed by Gu and Gillespie (1998). This correction considers that the orientations of trees do not grow perpendicular to the inclined terrain. The formulation is as follows (Eq. 1): where ρ n is the normalized reflectance, ρ is the uncorrected reflectance, α is the slope of the terrain, θ is the solar zenith angle, i is the illumination angle, and C is an empirical parameter v www.esajournals.org introduced by Soenen et al. (2005) to moderate the overcorrection of the SCS correction at large incidence angles. The parameter C (Eq. 2) is a function of the slope (b) and the intercept (a) of the regression line derived from the relationship between the reflectance and the cosine of the illumination angle.
Although JPL did a preliminary georectification of the data, the images showed residual misalignment that required additional correction. The images were co-registered using the automated image registration technique developed by Koltunov et al. (2012). This method combines band-wise compensation for radiometric differences between images (Koltunov et al. 2008) with an iterative gradient-based alignment method by Irani (2002), under the affine image motion model. The image from June 2015 was coregistered using June 2013 as the base image.
Auxiliary data.-A vegetation map called CAL-VEG (USDA 2017) with 30 m pixel resolution was used to identify areas dominated by blue oak and to remove the pixels dominated by grass, rock, or bare soil.
Multispectral imagery with 1 m resolution acquired in 28 June 2012 around 9 am and 0.6 m in 1 July 2016 acquired around midday from the National Agriculture Imagery Program (NAIP) was used to verify that the timing of observed mortality occurred during the 2013-2016 drought by careful evaluation of the two images (i.e., 2012 and 2016, before and after drought). The main limitation of this analysis was that the last AVIRIS image was recorded in summer 2015 and there was additional tree mortality in the fall of 2015 and spring of 2016. Thus, NAIP 2016 imagery may detect more tree mortality than AVIRIS 2015.
Field data.-Three field campaigns were used in this research, carried out in summer 2014, summer 2015, and spring 2017. The 2014 and 2015 field data collections were concurrent with AVIRIS overflights and were used to train the model, while the 2017 field data were used for validation (Appendix S1: Table S1).
Data collected during the field campaign of June 2014 included tree structure and composition in 21 circular plots of 20 m radius and three plots of 18 m radius. This inventory was repeated in June 2015 for nine plots (blue circles Fig. 2) dominated by blue oak that had no tree mortality during the 2014 campaign. The inventory in both field campaigns included the following data: ground cover fraction, tree species and location, diameter at breast height (DBH), canopy height, canopy base height, crown width, and health status, including whether each tree was alive or not (Table 1). For trees bigger than 10 cm of DBH, all measurements were recorded, while for small tress, only the health status and species were recorded. The plots can be grouped into sets of three. One set of plots (1 in Fig. 2) is situated in a very-dry, pure blue oak patch; set 2 is located in a patch dominated by blue oak with some gray pine present; and set 3 is situated in a pure blue oak patch with higher moisture content (Fig. 2). In total, 97 trees were measured. In addition, LAI measurements using digital hemispherical photographs were taken during these two field campaigns. LAI measurements were processed using Hemiview software.
The topography of each plot was characterized using a 4 × 4 pixel window to estimate the dominant aspect per plot. In a first step, we analyzed all individual pixels within the window to determine whether there was a dominant aspect. In a second step, the average aspect was calculated for the pixels within the dominant aspect. In most plots, between 80% and 90% of the pixels were within the dominant aspect. The set for plot 1 with 50% blue oak mortality was located on a steep slope; slope has been correlated with mortality in pine forests during this drought (Fettig et al. 2019a).
In spring 2017, we collected additional data on tree mortality for validation (Fig. 2). Trees were categorized into two DBH classes, below and above 10 cm, and in one of four health states, totally healthy with 0% dead, less than 50% of the canopy dead, more than 50% of the canopy dead, and canopy completely dead. One hundred 20 m radius plots were randomly selected along nineteen transects (one was removed because it was not dominated by blue oak) within pure blue oak patches, identified from the CALVEG map. In the field, however, we found that some plots were not blue oak dominated and these plots were removed, yielding 82 blue oak dominated plots. We identified tree mortality within each of these plots. In order to verify that tree mortality observed in the 2017 field collections occurred during our study period, we compared the NAIP imagery to the 82 plots measurements for June 2012 and 2016, collected at 1 and 0.6 m pixel resolution, respectively (https://earthexplorer.usgs. gov/). The color infrared NAIP images show healthy green vegetation in shades of bright red, while less healthy trees have smaller canopies and a darker, less red color. We overlayed the 82 field plot locations on both NAIP images to confirm that the 2017 recorded tree mortality occurred during the drought (Appendix S1: Fig. S1). Since the study area is an open woodland predominantly of blue oak and we accurately recorded the position of the trees (which are large compared to pixel size), we could identify each tree in the NAIP imagery. We visually compared the images 2012 and 2016 to identify any tree that died during the drought period. Spring 2017 data were not used to train the mortality model because it might overestimate tree mortality that was not present in AVIRIS 2015. We assume most of the tree mortality happened in winter 2015-2016 in our study region, but since we cannot be certain, we did not use these data for model training.

Analytical methods
The AVIRIS data acquired in June 2013 and June 2015 were analyzed for the most useful metrics, which were selected from a sensitivity analysis of 32 of spectral metrics (i.e., spectral indices, absorption-based features, and physically derived metrics) applied to the 2014 and 2015 training data set of nine plots (Table 1). The metrics selected are related to chlorophyll and water content, which previous studies found were related to tree mortality from drought (Byer andJin 2017, Sapes et al. 2019).
Optimized narrow-band normalized difference indices were then developed and compared for all possible combinations of two AVIRIS bands (n = 17,205) to determine whether the detection accuracy could be improved beyond that predicted by the 32 initial metrics. The best metric was then used to develop several models (i.e., lineal, polynomial, exponential regression) for mapping blue oak mortality, and the model with the highest accuracy was applied to the entire study site to predict percent tree mortality at the landscape level. The training-testing data for the model came from 2014 and 2015 field campaigns (Table 1) and were validated by the 82 blue oak dominated plots from the 2017 field campaign (Fig. 3).
Optical metrics computation.-Thirty-two optical metrics tested were grouped into (1) spectral indices (n = 16), (2) absorption-based features (n = 15), (3) a physically derived metric (n = 1), and (4) calculation of optimized band combinations. Sixteen spectral indices, both simple ratio and normalized difference ratio indices, were evaluated in this analysis (Appendix S1: Table S2). Fifteen absorption-based features (Appendix S1: Table S3) were calculated using the continuum removal approach (Clark and Roush 1984) to calculate local wavelength minima for the absorption maximum and the maxima at the edges of the absorption feature as adapted by Pu et al. (2003), using code from Clark and Roberts (2012). Location of absorptions features is shown in Appendix S1: Fig. S2. The physically derived equivalent water thickness (EWT) measures mm thickness of water per pixel (mm/m 2 ). It was computed using a Beer-Lambert light-extinction model on the 970 nm water absorption feature (Roberts et al. 1997). Optimized normalized difference indices were calculated by computing all possible combinations of two AVIRIS bands to identify band combinations with the strongest relationships to tree mortality and to identify the spectral region(s) most closely related to percent tree mortality. Narrow-band normalized v www.esajournals.org difference indices were computed based on untransformed spectral bands where ρ λa and ρ λb are the reflectance at bands λ a and λ b, respectively (Eq. 3). The indices were selected based on the maximum R 2 value between the index and the percent tree mortality within each plot. Only optimized indices that improved the R 2 of the previous 32 spectral metrics were considered.
Assessment of percent tree mortality.-In the first step, we computed the mortality percent by subtracting the optical metrics described above in years 2013 minus 2015 for each pixel in the AVIRIS images, using shapefiles of the 20 m radius plots measured during fieldwork. The study site is dominated by open areas where understory that might influence the tree's spectral reflectance was minimal. Most pixels at the study site are dominated by annual herbaceous grasses that were dry and dormant during the two flight periods; therefore, their contribution to change should be negligible. Given control for other likely confounding factors (the flights were approximately the same time of day, cloud free, low atmospheric water vapor, no/low smoke, etc.), we assumed that the temporal spectral change observed is due primarily to changes in the tree canopy.
We performed a sensitivity analysis on the correlations among all optical metrics and the field-measured percent of tree mortality by plot. Fig. 3. Methods flowchart: the rectangular box labeled tree mortality/field plot represents input from the original nine plots and the 82 plots (Fig. 2). The NAIP data used for validation were done at the step labeled Field Tree Mortality Survey NAIP imagery.
The best-fit model was applied to each correlation (e.g., lineal, polynomial, exponential regression). The models with the highest accuracy were compared to the plot data and applied to the entire study site to evaluate percent tree mortality at the landscape level. We applied a mask to remove non-tree pixels based on the CALVEG map. The 82 plots of 20 m radius collected in spring 2017 were used to validate the model at the landscape scale by estimating the prediction error and the standard error (Eqs. 4 and 5). The model was applied to a landscape of similar cover type using the CALVEG map to select areas dominated by blue oak forest.
where σ p and σ est represent the prediction and the standard error, respectively; Y R and Y P, the real and the predicted percentage of blue oak mortality, respectively; and N, the number of observations. Finally, we evaluated the relationship between topography and tree mortality by using aspect, the topographic position index (TPI), incoming solar radiation (ISR) for June, and the cumulative drainage model (CDM). TPI was generated using the Land Facet Corridor Designer Application ArcGIS extension. It was calculated using circular neighborhoods where neighborhood radius was based on pixel units. We used an area of 2 × 2 pixels thus generating neighborhood of 36 × 36 m. Each cell value represents the difference in elevation in meters between the cell studied and the average elevation of the neighborhood. The raster value is in TPI units. ISR was generated using Area Solar Radiation tool in ArcGIS Pro 2.5. Monthly incoming solar radiation for June was calculated at 7-day intervals, 0.5-hour intervals, default 1000 sky size (the upward looking viewshed from the cell of the sky map, and sun map raster), and standard overcast sky for a diffuse model. It calculated the insolation based on models developed by Rich et al. (1994) from the hemispherical viewshed algorithm and after developed by Rich (2000, 2002). The raster units are in watt hours per square meter (WH/m 2 ). CDM was generated in ArcGIS using the hydrology toolset of ArcToolbox to compute the potential uphill drainage area for each pixel. The flow accumulation tool calculates how many cells flow into each downslope cell of the output and it is used to represent the relative amount of water that can flow into each cell, thus the relative moisture availability of the pixel.
All image analyses were done with ENVI 4.7 and statistics (i.e., correlations, errors, regressions) with R.

Relationship between selected optical metrics and percent tree mortality
Among the 32 metrics used in the linear sensitivity analysis, only four (i.e., NDWI [R 2 = 0.76], RWI [R 2 = 0.74], MSI [R 2 = 0.68], and NDII [R 2 = 0.42]) had R 2 higher than 0.40. Wtr1AbWvl also had a high R 2 but the statistical relationship was not predictive (Appendix S1: Table S4 and Fig. S3). Two are simple ratios (RWI and MSI), and two are normalized indices (NDWI and DNII). The highest predictions were NDWI and RWI. All four spectral indices are related to canopy water content. The accuracy of the relationships using NDWI, RWI, and MSI was improved by a second-degree polynomial model, while the linear model was best for NDII. The polynomial metric with the highest accuracy (R 2 = 0.84) was for RWI, followed by NDWI (R 2 = 0.80) and MSI (R 2 = 0.77).

Results from the index optimization analysis
Among the 17,205 possible combinations of the 193 AVIRIS bands, the 100 best band indices have correlations of R 2 between 0.77 and 0.83 (Fig. 4a). The best correlations were for NIR bands around 752 and 956 nm (λ a ) and around 849 and 1262 nm (λ b ). The strongest relationship (i.e., R 2 = 0.83, Fig. 4b) was when bands at 937 nm and 1100 nm were combined (Eq. 6).
where λ937 and λ1100 represent reflectance value at 937 and 1100 nm.
The optimized narrow-band index analysis provided additional spectral regions related to percent tree mortality (Fig. 4a). The highest correlations were found with a combination of bands from 723-966 nm and 946-1472 nm (Fig. 4b), both of which relate to water content. Slightly lower correlations were found with a combination of bands from the visible region (433-608 nm) and SWIR region (1452-1800 nm) which are related to pigment content and water content respectively. Additional correlations were found with a combination of bands in the NIR (723-936 nm) and SWIR (2376-2445 nm) regions, again primarily related to water content.
The best optimized index (ND (λ 937.08_λ 1100.53) was fit to several models, and the best results were found using a second-degree polynomial function and a linear function providing an R 2 of 0.84 and 0.83, respectively. The optimized index had the highest accuracy for estimating tree mortality in relation to the previously described optical metrics. Since the v www.esajournals.org difference in accuracy among these models is negligible, we selected the simpler linear model to extrapolate results to the entire study area. The model formula is provided in the following equation (Eq. 7): Y ¼ 2018:1 Â X À 10:73 (7) where Y represents percent blue oak mortality and X the difference between the selected optimized index in 2013 and 2015.
Tree percent mortality assessment at regional scale Once we identified the best performing spectral metric (Eq. 6) for detecting declines in canopy moisture and tree mortality, it was applied to the 4994 ha (154,128 pixels) for which imaging spectroscopy data were available. From a total of 26,100 ha of pure blue oak woodlands, 562 ha of mixed blue oak and gray pine woodlands, and 1922 ha of non-tree dominated areas, we calculated that 1% of the area had tree high mortality between 51 and 100%, 8% had mortality between 26 and 50%, 30% had mortality between 11 and 25%, 33% had mortality between 1 and 10%, and in 28% no mortality was observed (Fig. 5).
The spatial pattern of tree mortality (Fig. 5b) is related to both the topography and canopy water loss ( Fig. 5a and c, respectively). Changes in water content are due to a combination of reduced water content per leaf area and a reduction in leaf area per pixel. Leaf water content, measured on oak leaves from this site in 2013 and 2014, did not show a significant change in water content on a leaf area (g/cm 2 ) or leaf mass (g/g) basis, nor did leaf mass area (dry weight, g cm 2 ) change between these years (S. L. Ustin, M. Huesca, K. L. Roth, et al. unpublished data). Field measurements of LAI from fisheye photographs measured in the same plots in 2014 and 2015 (n = 54 per year, 6 plots × 3 subplots × 3 photographs) show a large decrease in LAI (LAI = 0.75 in 2014 and 0.49 in 2015), suggesting that much of the decline in NDWI measured between 2013 and 2015 is due to lower LAI/pixel than to lower water content per leaf. As blue oak are a drought-tolerant deciduous species, a decline in LAI is an expected mechanism to reduce stress once soil moisture becomes limiting. Although the drought started in 2013, it was not until after June 2014 that leaves manifested significant reduction in size and in water content and canopies had significant declines in LAI.
Canopy water measured by NDWI declined everywhere, but the loss was greater in some areas than others. The condition of these sites is correlated with areas of high tree mortality, as the decline in water content is expected to precede tree mortality. In relation to the topography, there is more tree mortality on the southeasternto southwestern-facing slopes, compared to northern-facing slopes.

Map accuracy assessment
Accuracy of blue oak mortality estimates for the study site (Fig. 5b) was performed using the 82 field validation plots (distributed across all 19 transects) from 2017 by comparing mortality at plot locations with the 2012 and 2016 NAIP data (Fig. 2). The NAIP imagery was used to confirm whether the mortality in our field data occurred during the drought or occurred pre-2012. We deleted tree mortality for a small number of trees that died before the drought. With our adjusted true tree mortality numbers, we compared the true tree mortality to the mortality predicted by the model for each plot. In general, transects close to the training plots are expected to have higher accuracy. The difference between actual mortality and estimated mortality is less than 3% in transects 1a, 2b, 3a, 3b, 4a, and 10a of the validation plots. The largest difference between field-measured and imagery-predicted values was found in transect 5b with 21.64%, and the smallest was found in transect 2b with −0.06% (Appendix S1: Table S5). Average error for the 19 transects was a 2.08% underestimate, based on the field verified mortality.
Most plots were located on SE, S, or SW aspects, and none of the plots fell into our highest level of predicted mortality (Appendix S1: Table S5). On average, the field verified mortality was higher (mean = 14.55%) than the imagery-predicted mortality (mean = 12.47%) with under-prediction of 2.08%. The aspect for each measurement ranged between southeast and southwest, and no transects were located on slopes with northern aspects. The CalVeg map (USDA 2017) was used to identify locations for pure valley oak transects, and we did not sample in gray pine or mixed blue oak and gray pine stands from the CalVeg map. Later, we found this decision resulted in no pure blue oak pixels located on north-facing aspect orientations. The small sample size of the data in Table 1 and Appendix S1: Table S5 precluded further interpretation. The confusion matrix (Table 2) shows that the highest errors are found in the 0% and 26-50% mortality classes and the most accurate class was in 1-10% mortality. Plots with 0% mortality were wrongly classified more than 50% of the time into the 11-25% mortality category. Plots with 26-50% mortality were most often wrongly over-classified into the 1-10% mortality class. The best performing class (1-10%) had an even spread of errors across the other classes except for the 51-100% class, for Note: The values show the percentage of plots in each group classified at different levels of tree mortality (i.e., 0, from 1% to 10%, from 11% to 25%, from 26% to 50%, and 51% and above tree mortality). which they were never classified. The source of the classification errors in this table was likely from the discretization of continuous variables into classes but considering that the largest errors were for the 0% and 26-50% categories, and not in adjacent categories, it supports further model refinement is needed.

Mortality relationship to topography
We evaluated the relationship between mortality and south-facing slopes (Fig. 6a), which shows the relationship between aspect in degrees vs. predicted tree mortality. It clearly shows heightened mortality in the south-facing orientations relative to more northern-facing aspects. The frequency distribution of pixels reinforces our results that mortality was low across most aspects and within the range of 20%, consistent with Table 1, Appendix S1: Table S5, and our image data (Fig. 5b). Above 40% tree mortality, differences in aspect emerge with the hotter and drier south-facing slopes having the highest mortality. The low number of data points in the northern directions (337.5°to 22.5°) is due to masking out all non-oak dominated pixels (grasslands, meadows, shrublands, and mixed oak-pine woodlands). Clearly stands of pure blue oak on northern-facing slopes are uncommon, whereas on the south-facing aspects, only blue oak (with grass understory) is present. Pixels with higher solar irradiance (Fig. 6b) are associated with higher tree mortality, and higher solar irradiation is at its highest in the southwest orientation. TPI shows that mortality is associated with flatter topographic areas, which is a consequence of the low slopes throughout the region. Therefore, there were no significant differences using the TPI index (Fig. 6c). Pixels with higher potential for flow accumulation, from the cumulative drainage model analysis (Fig. 6d), are associated with pixels having lower tree mortality (less than 20%), while pixels having very low potential for flow accumulation have high mortality. The south-facing aspects combined with the low potential for water drainage all show that the driest areas experienced the highest oak mortality. If tree mortality that is independent of flow accumulation represents the random mortality rate, the pixels with higher mortality must represent the enhanced mortality caused by the drought.

DISCUSSION
California's 2012-2016 drought represents the types of increasingly severe impacts predicted for drought-prone ecosystems as climate change progresses (Allen et al. 2010, Ullrich et al. 2018. The extent of blue oak mortality in our study areas was surprising, given the droughthardiness of this species. The cause of the mortality appears to be a direct impact from a combination of water deficit, such as more negative water potentials, and indirect effects from loss of leaf area through the production of smaller and fewer leaves and early senescence, that possibly contributed to carbon starvation (McDowell et al. 2008, McDowell andSevanto 2010). These observations are consistent with Weitz (2018) who reported physiological and morphological changes in blue oak throughout the drought and McLaughlin et al. (2020) who correlated blue oak mortality with precipitation, low soil moisture, and lithology from a site in California's Inner Coast Range, at a comparable latitude to our site. We observed unusually small leaves and widespread leaf senescence at our site in mid-June 2015 (S. L. Ustin, M. Huesca, K. L. Roth, et al. unpublished data). This pattern of early season senescence followed by fall and winter mortality contrasts with mortality in the conifer forests found at higher elevations in the region, where massive die-off occurred from bark beetles and other insects that attacked drought-weakened trees (Fettig et al. 2019a, Restaino et al. 2019). In the blue oak populations, there was no outbreak of insect predation.
Our results show temporal, spatial, and spectral changes in blue oak woodlands occurred over the multiyear drought. Furthermore, we found the NAIP imagery confirmed the tree mortality that we documented using the AVIRIS imagery, and that mortality ranged between 1% and 60% of the blue oak canopy cover, as also confirmed by our comparison of the 82 plots in 2017 with the 2012 and 2016 NAIP imagery. We recognize that the randomly selected field sites missed the areas of highest mortality as shown in this analysis of the image data, which was a consequence of the plots only covering about 1% of the total area. However, our results are roughly equivalent to Das et al. (2019), who report 23% standing dead Q. douglasii in Sequoia National v www.esajournals.org Park (about 75 km south), but their methods could not confirm that all dead trees died in the drought.
Our research hypothesis overall, that changes in the spectral reflectance of arid and open woodlands (with low cover fractions of green plant understory) that discriminate dead and live trees will exceed canopy changes caused by other factors, is confirmed. Changes in reflectance of AVIRIS hyperspectral data are correlated with differences between live and dead blue oak trees that allow mortality to be detected and quantified. Although the model we used is empirical and not expected to provide accurate predictions under other conditions, it shows that a spatially explicit map of tree mortality at the landscape scale can be created and that the method could be extended across the entire range of blue oak.

Model performance
This study illustrates the utility of using imaging spectroscopy data for mapping percent of blue oak mortality at landscape levels (Fig. 5b). The model developed using training data from the inventory plots produced a map of tree mortality with high accuracy. The validation performed with 82 independent plots was predictive (R 2 = 0.61) for the plots containing between 5% and 50% tree mortality (Table 1 and Appendix S1: Table S5). The highest prediction accuracy was found in the validation plots with percent tree mortality greater than 0% but less than 10%. We did not identify high tree mortality >60% from our field surveys, likely because our randomly selected training data did not include areas with the higher ranges of mortality. Our fieldwork in 2014, 2015, and 2017 did not identify trees that died between 2015 and 2016, and the data in Table 2 indicate that the image data underestimated mortality compared to the field observed mortality by nearly 15%. The 82 plots measured for validation in summer 2017 data showed that there are a few areas (20 m radius) in the region with more than 50% tree mortality, in agreement with the inventory transect data and NAIP results.
The assessment between tree mortality and topography shows that blue oaks are surprisingly scarce in north-facing aspects. This may be an indication that they are less competitive in cooler and wetter north-facing slopes. In addition, our finding that pixels with a higher potential for flow accumulation are associated with lower tree mortality is consistent with the modeling study by Brown et al. (2018) who reported an increased probability of mortality in blue oaks with the combination of high spring, summer, and fall temperatures from 2013-2014, low precipitation (<363 mm) in 2015, and a loss of water storage 30 cm below baseline.
Our results could be improved with a more complete and comprehensive set of training data that covered the full range of mortality variability and if the AVIRIS data acquisition had continued into 2016 and 2017. The model would also be improved if it had been possible to identify the extent and pattern of canopy mortality at the individual tree scale in AVIRIS data. Our pixel size was 18 × 18 m, and it was difficult to tell whether detected mortality was comprised of multiple dead trees in a pixel or a single dead crown. Higher spatial resolution data may have helped resolve these differences. Thirdly, acquiring measurements of additional environmental variables that play a role in forest senescence and mortality processes would improve model accuracy, such as soil moisture and soil drainage patterns at the landscape scale (and shown for a subset of data in Fig 6a and 6b).

Climate change
Climate change predictions show increased risk of drying in arid ecosystems (Diffenbaugh et al. 2015, Thorne et al. 2015, Ullrich et al. 2018). These changes are expected to increase tree mortality, because conditions are already at tipping points, becoming more arid with drier conditions (van Mantgem et al. 2009, Allen et al. 2010, Weed et al. 2013, Trenberth et al. 2014. Increased temperatures during drought periods exaggerate protective responses and induce drought responses more quickly (McDowell andSevanto 2010, Sala andMencuccini 2014). A positive feedback is set up that results in more rapid ecosystem conversion to less resilient and more arid systems. We found that hyperspectral imaging (imaging spectroscopy) can provide a basis for monitoring change in this drought adapted ecosystem and that it could provide spatially and temporally resolved data with which to populate models for more accurate predictions at the local and watershed scales that forest managers need.
Consistent with our results, others have reported that blue oak mortality is concentrated in the driest, southern part of the Sierra Nevada Range (Brown et al. 2018, Jacobsen and Pratt 2018, Das et al. 2019, Fettig et al. 2019a) and southernmost interior Coast Ranges (Brown et al. 2018, McLaughlin et al. 2020. We show that blue oaks on south-facing aspects are more susceptible to drought stress, while the potential for soil moisture to drain into a pixel is correlated with less drought stress. We noted steeper terrain was also correlated with trees having greater drought stress. Greater ability to predict specific locations and conditions of where blue oak mortality is predicted to occur would provide forest managers essential information to adopt better management practices for vulnerable forest sites for this endemic forest type.

ACKNOWLEDGMENTS
This work was funded by the U.S. Geological Survey Southwest Climate Science Center, Grant ID GS15AP00173. We thank the field crews who collected the field data in the different years used in this publication, including Maria del Mar Alsina, Pia van Benthem, Mariano Garcia, Paul Haverkamp, Mui Lay, Keely Roth, and Michael Whiting. We also thank the students in NASA's Student Airborne Research Program who participated in fieldwork with us in 2013, 2014, and 2015. We thank Mr. Woody Turner, program manager for the HyspIRI planning grant # NNX12AP87G to S.U. and to Jack Kaye, NASA Associate Director for Research, Earth Science Division, for supporting the AVIRIS overflights.