Estimating Individual Tree Mortality in the Sierra Nevada Using Lidar and Multispectral Reflectance Data

Widespread tree mortality events occur during periods of severe drought in temperate conifer forests and are expected to become more frequent in many areas due to climate change. Improved mapping of individual tree mortality is needed to identify risk factors and design effective conservation strategies. In this study, we used National Ecological Observatory Network (NEON) lidar and multispectral reflectance airborne observations to map individual tree mortality over a 160 km2 area during and after the 2012–2016 drought for two sites in California's Sierra National Forest. We used NEON lidar to derive tree locations and crown perimeters and multispectral data to map tree mortality for more than 1 million trees. We found that 25.4% of the trees in our study area died between 2013 and 2017, with considerably higher mortality at the lower‐elevation Soaproot Saddle site. Between 2017 and 2019, an additional 2.0%–2.8% of the trees died each year. Two wildfires in 2020 and 2021 increased tree mortality within burned area perimeters by 49%–89% between 2019 and 2021. Consistent with previous work, we found that tree mortality risk increased as a function of tree height. Tree mortality was positively associated with distance from rivers, trees per hectare, and decreasing slope at the lower elevation site. In contrast, increasing slope was positively associated with tree mortality at the higher elevation site. Our approach and dataset provide a means to study the combined effects of drought and wildfire on tree mortality and may improve projections of forest resilience under a changing climate.

During the 2012-2016 California drought, evapotranspiration in the Sierra Nevada exceeded precipitation in mid-elevation forests by about 1,500 mm in the 4 years from 2012 to 2015 (Goulden & Bales, 2019). Canopy water content losses were observed throughout California's forests (Asner et al., 2016). Previous studies have found that topographic features related to hydrology play a key role in mortality risk during drought conditions in the Sierra Nevada. For example, lower elevation trees on shallower slopes are at higher risk of dying during drought (Paz-Kagan et al., 2017;Stovall et al., 2019). In addition, trees surrounded by granite outcrops, farther from rivers, and on southwestern slopes have been shown to have a higher risk of tree mortality (Paz-Kagan et al., 2017). Stovall et al. (2019) used 2013 lidar from the National Ecological Observatory Network (NEON) and imagery between 2009 and 2016 from the National Agricultural Imagery Program (NAIP) to map and classify 1.8 million individual trees as dead or alive over two NEON sites in the Southern Sierra Nevada . Large trees (greater than 30 m in height) were two times as likely to die as smaller trees (trees between 5 and 15 m in height) . In contrast, a nearby study in Sequoia National Park over the same drought period showed that trees belonging to the genus Pinus were the only group to show an increasing trend in mortality with height, whereas most species showed an increase in resilience with height (Stephenson & Das, 2020). Since the western pine beetle prefers large Pinus trees, the authors reasoned that the pine beetles were responsible for the trend in increasing mortality with height for Pinus trees (Stephenson & Das, 2020).
To test hypotheses about drivers of tree mortality during drought conditions, we need accurate tree mortality maps that rely on robust and reproducible classification algorithms. In this study, we use an updated set of lidar and hyperspectral reflectance observations available from NEON for the years 2013, 2017, 2018, 2019, and 2021 to create maps of tree mortality during and after the 2012-2016 California drought. Our approach makes use of these new data resources to build on the work of Stovall et al. (2019) within an open-source framework. We use all available lidar to identify the location of treetops and then identify the perimeters of individual tree crowns based on the lidar data for each year. Tree mortality is based on co-registered hyperspectral reflectance data from the same flights which are more closely aligned with the lidar data than the NAIP imagery used in past work. NEON data providers also correct their hyperspectral reflectance data using a surface bidirectional reflectance distribution function. As a result, vegetation indices such as relative greenness (green surface reflectance divided by the sum of red, green, and blue surface reflectance) and normalized difference vegetation index (NDVI) derived from their hyperspectral datasets are potentially more robust with respect to variations in solar zenith angle, collection angle, and surface topography. In contrast, for NAIP imagery, little information about the time of day and angle of data collection is publicly available to correct for these variations.
In this study, we compare our revised method for identifying individual tree mortality following the 2012-2016 California drought to the method developed by Stovall et al. (2019) over a similar region and time interval, quantify tree mortality in subsequent years (2017)(2018)(2019)(2020)(2021) in response to drought legacy effects and wildfire, and evaluate environmental geospatial controls on tree mortality. The method we develop may be useful for tracking tree mortality in other forests where aircraft lidar and spectral data are available. Insights from the regional maps of tree mortality over time may help inform forest stewardship, models simulating forest ecosystem processes, and environmental policy.

Study Site
We selected NEON's Soaproot Saddle and Lower Teakettle sites in the Sierra National Forest in the Sierra Nevada which experienced severe drought between 2012 and 2016 ( Figure S1 in the Supporting Information S1). The spatial extent of these sites is shown in Figure 1. The elevations range from 1,000 to 1,400 m at Soaproot Saddle which has a mean precipitation of about 900 mm per year (Krauss et al., 2018). Because the site is located along the rain-snow transition, winter snow melts quickly (Krauss et al., 2018). Dominant tree species at Soaproot Saddle include Ponderosa pine (Pinus ponderosa), incense cedar (Calocedrus decurrens), canyon live oak (Quercus chrysolepis), and California black oak (Quercus kelloggii) with high mortality observed among Ponderosa pine from bark beetle infestations during the drought (Krauss et al., 2018). Compared with Soaproot Saddle, Lower Teakettle has a higher elevation and wetter climate with respect to both precipitation and snow accumulation which provides more moisture during the dry season. The elevation for Lower Teakettle ranges from 2,000 to 10.1029/2022JG007234 3 of 18 2,800 m with a mean precipitation of 1,220 mm per year (Krauss et al., 2018). In addition, snow accumulates at Lower Teakettle with a 30-year average maximum depth of 1,140 mm (Krauss et al., 2018). In contrast to Soaproot Saddle, the forest of Lower Teakettle is characterized by red fir (Abies magnifica), white fir (Abies concolor), Jeffrey pine (Pinus jeffreyi), and lodgepole pine (Pinus contorta) (Krauss et al., 2018).
The yellow polygons in Figure 1 show the overlap of high-resolution lidar from NEON's airborne observation platform for the years 2013, 2018, 2019(NEON, 2021. NEON also has co-registered, high-resolution (1 m by 1 m) hyperspectral data (NEON, 2022b) corrected for atmospheric effects and topography (Karpowics & Kampe, 2022) over the same spatial extent and years except for some missing data in 2013. These rich datasets provide the opportunity to make detailed tree mortality maps over time.

Lidar and Canopy Height Models
We downloaded lidar point cloud tiles collected from NEON's airborne observation platform (NEON, 2021). The lidar resolution for 2013 is 3 points/m 2 and increases to 6 points/m 2 in 2017-2019 and 18 points/m 2 in 2021. Each lidar point is described by an x-, y-, and z-coordinate. The empirically-derived horizontal error is 3-5 cm, while the vertical error is estimated to be 5-10 cm (Krause & Goulden, 2022). To reduce noise, we filtered out any altimeter values more than five standard deviations from the mean. We also removed altimeter values greater than 120 m in height after normalizing the lidar to account for the surface topography.
In addition, we downloaded NEON's canopy height model raster data (NEON, 2022a) to help create vegetation masks to aid in the granite-detection algorithm described in Section 2.4 and to compute the canopy cover fraction. A recent study found significant effects of canopy cover on microclimate using plots ranging from 5 to 1,300 m 2 in area (Zellweger et al., 2020). We chose to use the 2013 canopy cover fraction within a 20-m radius of each Survey. The background raster image shows the NASA Shuttle Radar Topography Mission elevation from low elevation (dark green) to high elevation (light purple). We shaded the background from no shade (northern aspect) to light shade (southern aspect) to illustrate the topography of the region. The center location of the zoomed-in image is shown with a black point in the California inset inside the Sierra Nevada ecoregion (black hashed polygon). tree's location (an area of about 1,260 m 2 ) to investigate tree mortality related to initial microclimate conditions. Because we find a polygon describing the crown perimeter of each tree for each year of data (Section 2.2.5), we exclude the individual tree's 2013 crown perimeter from the 20-m buffer when calculating the canopy cover fraction.
We also downloaded NEON's digital terrain model raster data (NEON, 2023) to derive slope and aspect to assess potential tree mortality drivers. Both the canopy height and digital terrain model data have a resolution of 1 m by 1 m. To reduce potential noise in the high-resolution digital terrain model, we calculated the median slope and aspect within a 10-m radius of each tree. The algorithm we used to find each tree location is described in Section 2.2.5.

Hyperspectral Reflectance Data
We downloaded mosaic orthorectified hyperspectral reflectance data from NEON's airborne observation platform for the years 2013, 2018, 2019(NEON, 2022b. The surface radiance collected by the airborne observation platform is transformed to surface reflectance by removing atmospheric effects from haze, cloud shadow, waterbodies, water vapor and topography and applying a bidirectional reflectance distribution function (Karpowics & Kampe, 2022). This dataset includes 426 spectral bands and is distributed at a spatial resolution of 1 m by 1 m. The wavelength interval for band lengths between 340 and 1,349 nm is 0.40 nm, while the interval from 1,351 to 1,449 nm is 0.98 nm. We approximated the Landsat 8 spectral bands by taking the mean surface reflectance for blue (452-512 nm), green (533-590 nm), red (636-673 nm), near infrared (851-879 nm), and shortwave infrared 1 (1,566-1,651 nm) to create spectral indices and natural-color (red-green-blue or RGB) images for training labels. We also selected two water absorption bands, water band 1 (w 1 ) from 1,440 through 1,460 nm, and water band 2 (w 2 ) from 1,935 through 1,955 nm, because the water in leaves is a strong absorber of radiation in these wavelengths (Carter, 1991). Reflectance in the water bands should be lower where there is higher leaf water content. All wavelength ranges listed here are inclusive of the boundaries.
From the spectral bands, we created five vegetation indices including relative greenness, NDVI, normalized difference moisture index (NDMI), and two additional indices like NDMI created with the water absorption bands instead of the shortwave infrared 1 band. These indices are described in more detail in Table S1 in the Supporting Information S1. Because the relative greenness yielded the highest accuracy in the algorithm described in Section 2.3 (Table 1), we focus on relative greenness for this study. Diurnal variations in canopy moisture during the air campaigns may have contributed to lower performance of the moisture indices.

National Hydrography Dataset
We downloaded the National Hydrography Dataset Plus High-Resolution for region 1803 which encompasses our study site from the United States Geological Survey (USGS) Earth Explorer website (U.S. Geological Survey, 2019). While we show named rivers within this dataset in Figure 1, we used all rivers in the dataset for our analysis of the relationship between distance to rivers and tree mortality.

Landsat Surface Reflectance
We retrieved Landsat Collection 2 Level 2 Tier 1 surface reflectance data from the USGS EarthExplorer website. Path 42, row 34 completely covers our region of interest, so we used this tile in our analysis. To compute NDMI, we collected the near infrared band (NIR; Band 4 for Landsat 5 and Band 5 for Landsat 8) and shortwave infrared 1 band (SWIR1; Band 5 for Landsat 5 and Band 6 for Landsat 8) and applied the formula (1) NDMI ranges from −1 to 1 where −1 represents dry conditions (low near infrared surface reflectance relative to shortwave infrared 1 surface reflectance) and 1 represents moist conditions. The change in dry season NDMI has been found to correspond to tree mortality (Goulden & Bales, 2019), so we compared our estimates of tree mortality against dry season NDMI. To obtain dry season NDMI, we selected the four Landsat tiles temporally nearest 15 September for a given year, filtered each tile for clouds, and took the median of the surface reflectance pixels for the surface reflectance bands before computing NDMI. We aggregated 2013 and 2017 tree mortality to the Landsat resolution and took the difference to compare to the corresponding change in NDMI from Landsat between late summer 2013 and 2017. Because some hyperspectral data is missing in 2013, this analysis does not span the entire region of interest shown in Figure 1.

Tree Polygon Algorithm
Using a 96-m-by-96-m square moving window with a 24-m buffer, we combined all the lidar point clouds from 2013, 2017, 2018, 2019, and 2021 for the Soaproot Saddle and Lower Teakettle sites into one combined point cloud of x, y, and z locations. We combined the point clouds so that we would detect any trees that were standing at any point between 2013 and 2021. If a tree fell or was logged between 2013 and 2017, for example, we would have the treetop location from the combined point cloud and use the corresponding change in crown area from 1 year to the next to classify it as dead. To filter noise in the lidar point clouds, we dropped any z-values greater or less than 5 standard deviations from the mean. Following lidR package documentation, we first rasterized the canopy height to 1-m resolution and smoothed the canopy height raster using a 3 × 3 median filter. Elsewhere in our methods, we used the NEON-processed canopy height models; however, for this portion of the methods, we generated the canopy height model corresponding to our moving window as we processed the lidar point cloud. We used the find_trees function (R package: lidR) to identify the treetops of the combined point cloud using the local maxima filter lmf (R package: lidR).
The lmf function uses its own moving window on the smoothed canopy height raster to identify local maxima and classify them as treetops. We wanted the length of the local maxima filter window to grow with the height of the trees as suggested by the lidR package documentation, so we adjusted the window size according to the following function: where w is the width of the window and h is the height of the canopy, both in meters. We decided on the parameters for this function through trial and error on a subset of trees to reduce the number of trees that are counted more than once when the window is too small and that are missed when the window is too large.
After finding the locations of the trees from the whole point cloud, we went through each year of lidar data (2013, 2017, 2018, 2019, and 2021) to assign lidar points to each tree for each year. To do this, we used the function segment_trees with the dalponte2016 algorithm (R package: lidR) which has three adjustable parameters (Dalponte & Coomes, 2016). We selected a minimum threshold of 3 m and adjusted the rules for adding returns to a given tree so that returns need to be at least 65% of the height and 75% of the mean collection of returns assigned to a given tree. Because there are partial trees around the edges of the 96-m-by-96-m window, we filtered out any trees with treetop locations in the 24-m buffer to prevent partial trees or duplications in the segmentation process. We used the concaveman function to draw polygons around the x-and y-coordinates of each tree for each year (R package: concaveman) and adjusted the height of each tree using the canopy height model for each lidar year. The resulting attributes for each tree are the crown perimeter polygons and height for each year surveyed by NEON air campaigns. All analyses were completed using R, version 4.1.2.

Training Data
We created an evenly spaced training grid. First, we drew a bounding box around all the NEON lidar data files for all years (2013, 2017, 2018, 2019, and 2021). Then, we selected a reference Landsat 8 image for path 42, row 34, which completely covers our region of interest. We transformed the raster image to the same coordinate reference system by 300 m in both the x-and y-directions, resulting in 120 points in the x-direction and 72 points in the y-direction for a total of 8,640 grid points. Each grid had a 30-m by 30-m square perimeter coinciding with Landsat data pixels.
We manually labeled trees inside these training squares using 2017 RGB imagery that we derived from the NEON hyperspectral data following our flowchart shown in Figure S2 in the Supporting Information S1. First, for any tree that had a crown area greater than 1 square meter that drops below 1 square meter in 2017, 2018, and 2019, we labeled the tree dead for the corresponding year. Next, if a quarter or more of the tree crown polygon was taken up by road, ground, or part of another tree, we labeled it as uncertain. If neither of those criteria were met, we considered evidence that the tree was alive or dead in 2017. If at least one third of the crown area of the tree was red or gray, we labeled it as dead if the same was true in 2018 and 2019. If the tree appeared dead by crown area color in 2017 (i.e., if at least one third of the crown area was red or gray) but more than two-thirds of the crown area in 2018 or 2019 was green, we labeled it as alive. Alternatively, if the crown area in 2017 was more than two-thirds green, we labeled it as alive. If we could not see that the tree was alive in 2017, we labeled it as live if at least two thirds of the crown area was green in 2018 or 2019 under the assumption that a tree alive in 2018 or 2019 must be alive in 2017.
We also recorded three levels of certainty (0, 0.5, and 1) where 0 refers to no certainty, 0.5 is less certain, and 1 is relatively certain. While our analysis considered 0.5 and 1 acceptable levels of certainty, future analyses could filter out the trees with certainty values of 0.5. Additional detail for the assignment of labels and certainty values is included in Figure S2 in the Supporting Information S1. We visually inspected 9,025 trees and labeled 8,897 as dead or alive using our protocol. The life status of 128 trees was unclear, and these trees were not included in fitting parameters for the tree mortality algorithm. We labeled 6,377 trees (71.7%) as alive and 2,520 trees (28.3%) as dead. When we filtered out trees with no corresponding hyperspectral data, we had a total of 8,873 labeled trees for our algorithm.

Wildfire Perimeters
We downloaded the National United States Forest Service Final Fire Perimeters to quantify the number of trees which died within wildfire perimeters in our region of interest (USDA Forest Service National Forest System Lands GIS and Fire personnel, 2023). The 2020 Creek Fire and 2021 Blue Fire both burned through the Soaproot Saddle site, affecting tree mortality between 2019 and 2021.

Tree Mortality Classification Algorithm
Our algorithm for generating tree crown perimeter polygons and classifying the trees and live or dead is shown in Figure 2. We explored five vegetation indices (relative greenness, NDVI, NDMI and two additional normalized difference moisture indices derived from two water absorption bands) described in Table S1 in the Supporting Information S1 to classify trees as dead or alive. Previously, Stovall et al. (2019) masked NAIP imagery where the canopy heights were less than 5 m to isolate tree canopies and limit tree height to 5 m and greater. We followed the 5-m minimum height for our tree dataset, but we masked the hyperspectral data corresponding to canopy height model values less than 4 m so that we could still collect some lower canopy returns within the crown perimeters, particularly for trees near 5 m in height. In addition, we found that snow confounded tree mortality estimates in high elevation regions in early summer, so we also masked pixels where the mean luminosity (mean of the red, green, and blue bands) was greater than 0.2. For reference, we set the minimum reflectance to 0 and the maximum to 0.2 when viewing the RGB imagery in the Quantum Geographic Information System (QGIS) software (which is free and publicly available) for labeling trees for training data. Reflectance values greater than 0.2 for the three visible light bands were bright surfaces such as snow or granite.
We partitioned our 8,873 labeled trees into a training, validation, and test dataset by randomly sampling 60% of the labels (5,331 trees) for the training set, 20% of the labels for the validation dataset (1,777 trees), and held out the remaining 20% of the labels (1,777 trees) for the testing dataset. We used both a validation dataset and testing dataset, so that we could use the validation dataset to select the best spectral index while setting the test dataset aside for the final accuracy assessment. We set a seed of zero for the random number generator within R to ensure reproducibility for the training, validation, and testing split.
Our model for classifying dead and live trees has two parameters that use an individual vegetation index derived from NEON's hyperspectral reflectance data. For all five indices we considered, increasing values of the index toward 1 indicate healthier vegetation (e.g., more green or moist), so we used the same algorithm to explore all five indices. The first parameter is a vegetation index threshold below which we label a pixel within a tree's crown perimeter dead. For example, healthy tree canopies typically have high NDVI values for which near infrared surface reflectance is high relative to red surface reflectance. The second parameter follows the methods of Stovall et al. (2019) and is a threshold of dead pixels within a tree crown perimeter over which the tree is classified as dead. For example, Stovall et al. (2019) classified a tree as dead if more than 37.5% of its pixels within the crown area were classified as dead with respect to falling below an NDVI threshold.
We considered vegetation index thresholds from 0 to 1 with an increment of 0.01 and crown area fraction thresholds of 0-1 with an increment of 0.01. We computed the accuracy on the training dataset for every combination of these two parameters. The grid search for the vegetation indices with the highest training dataset accuracy (NDVI and relative greenness) are shown in Figure 3. Here, accuracy is defined to be the number of correct classifications of a tree as dead or alive over the total number of classifications. The dark blue regions of the parameter space maps indicate an accuracy resulting from guessing that most of the trees are dead, while the green regions result from guessing that most of the trees are alive. When we classify trees with at least 37% of pixels under a relative greenness of 0.37 as dead, we achieve an accuracy of 93.3% on the training dataset.

Application, Validation, and Testing of Method
Because the number of labels is much smaller than the total number of trees within the domain, we analyzed the accuracy of the parameters in our validation dataset in relation to the number of labeled training data values used to fit the parameters to see if we labeled enough trees. We set a seed in R and randomly sorted the training labels. Then we selected the first n labels to use as training data where n begins at 10 and increases by increments of 10. We ran the whole algorithm apart from evaluating the performance on the testing dataset for each quantity of cumulative labels. We repeated this process for 30 randomly-generated sequences of the training data to capture the range in accuracy generated by the order in which the training data are selected. The outcome from this analysis is shown in Figure 4. We found that the accuracies on the training and validation dataset converge after including just a couple thousand trees for training data.
The optimal parameters and accuracy analysis of the indices explored in this study are arranged in order of training dataset accuracy in Table 1. We consider a true positive to be an instance where a dead tree is correctly classified as a dead tree. Sensitivity refers to proportion of true positives which are the total number of accurate classifications of dead trees over the total number of dead trees. Specificity refers to the number of correctly identified live trees over the total number of live trees. Of the five indices explored, relative greenness yielded the highest accuracy on the training and validation datasets and also had the highest sensitivity and specificity on the validation dataset.
The confusion matrices for the training, validation, and testing datasets are shown for relative greenness in Table 2. The accuracy of our relative greenness model was 92.1% on the dataset held out for testing. A sample image of trees labeled by our algorithm is shown in Figure S3 in the Supporting Information S1. Shadows visible in the RGB image ( Figure S3a in the Supporting Information S1) appear to be removed by the normalization in our relative greenness metric ( Figure S3b in the Supporting Information S1).

Granite-Detection Algorithm
We hypothesized that higher granite cover fraction near a tree would be associated with negative health outcomes for trees based on soil quality as in a previous study (Paz-Kagan et al., 2017) and due to higher relative sensible heating and less latent cooling among trees surrounded by granite outcrops compared with those surrounded by vegetation. To explore this hypothesis, we developed a granite-detection algorithm. We classified all pixels with a canopy height less than 0.5 m and a luminosity greater than 0.11 as granite pixels. Lastly, we rasterized National Hydrography Dataset waterbody polygons and United States Forest Service roads polygons to filter out any pixels that are associated with waterbodies or roads. To investigate the role of granite outcrops in tree mortality, we calculated the fraction of ground pixels within 20 m of a tree that are classified as granite.

Tree Matching Algorithm
We downloaded the Stovall (2019) dataset from figshare and matched the trees to trees in our dataset to enable a comparison of our method with the previous work it builds upon. For each tree in our dataset, we found the distance from the x-and y-coordinates of the 2013 treetop location and each tree in the Stovall dataset. If the closest tree in the Stovall dataset fell inside the 2013 tree crown perimeter for our tree, we considered it a match. If multiple trees met this criterion, we selected the one closest in height to our tree. We matched 714,656 trees out of 1,011,634 trees located in our region of interest with overlapping NEON spectral data. Our comparison is based on this set of 714,656 trees.
We used the matched dataset to compare the distribution of the mean NDVI within the crown area of trees in 2017 and over time for each of the following categories: (a) trees where both studies label the tree as alive, (b) trees where both studies label the tree as dead, and (c) trees we labeled as live in 2017 but Stovall labeled as dead in 2016. We use the first two categories as reference datasets to determine whether the trees in the third category are more likely to be alive or dead in 2017.

Tree Mortality During and After the 2012-2016 Drought
We estimate that 324,301 of a total 1,011,577 trees (32.1%) were dead in the study domain by 2017 with 25.4% of the total trees dying between 2013 and 2017 ( Figure 5). Tree mortality was considerably higher at Soaproot Saddle than Lower Teakettle (Figures 5a and 5b for cumulative numbers and Figure 6 for spatial maps). Of the 275,103 trees we mapped at Soaproot Saddle, we classified 49.6% as dead by 2017 (Figure 5a). In contrast, of the total 736,474 trees we mapped at Lower Teakettle, only 25.5% were dead by 2017 (Figure 5b). The differences between the two sites were further amplified when considering the change in mortality between 2013 and 2017. During this interval, tree mortality increased by 45.8% at Soaproot Saddle and by 17.4% at Lower Teakettle.
Maps of the temporal evolution of tree mortality are shown for Soaproot Saddle in Figure S4 in the Supporting Information S1 and for Lower Teakettle in Figure S5 in the Supporting Information S1. We found very little additional tree mortality at Soaproot Saddle between 2017 and 2019; however, additional mortality was visible by 2021 when absolute tree mortality reached a cumulative level in non-fire affected areas of 55.9% (Figures 5 and 6). In Lower Teakettle, visible increases in tree mortality occurred each year between 2017 and 2019, with little additional change in the final interval between 2019 and 2021. There are some areas that appear misclassified in 2019 because trees labeled dead in 2019 are labeled live later in 2021. We explore some key areas of inconsistent classification in our 2019 tree mortality estimates in Figure S6 in the Supporting Information S1. One region in Lower Teakettle ( Figure S6a in the Supporting Information S1) appears to have high levels of snow cover and low luminosity, while the second region ( Figure S6b in the Supporting Information S1) has trees that were brown in 2019 but green in 2021.
The Creek Fire in 2020 and the Blue Fire in 2021 contributed to significant additional tree mortality at Soaproot Saddle. While most of the Creek Fire affected areas of the Sierra National Forest to the northwest of our region of interest, part of the wildfire extent intersected our study domain, affecting 76,644 trees ( Figure S4d in the Supporting Information S1). Within the Creek Fire perimeter intersecting our study domain, 45.2% of the trees were dead in 2019 prior to the start of the fire. After the fire in 2021, 85.5% of the trees within the fire perimeter were dead (an increase of 89%). We found that tree mortality within the 2021 Blue Fire perimeter increased from 64.3% to 95.6% between 2019 and 2021 (an increase of 49%) among the 11,068 trees within the Blue Fire's perimeter. We show a zoomed in image of the 2021 Blue Fire burn area in Figure S7 in the Supporting Information S1.
We plotted the mean relative greenness and mean NDVI density distributions for our estimates of live and dead trees in Figure 7. The relative distributions for mean relative greenness reach a maximum at 0.42 for live trees and 0.33 for dead trees. The mode for mean NDVI is 0.76 for live trees and 0.37 for dead trees.

Biophysical Drivers of Tree Mortality
Next, we consider biophysical drivers of tree mortality within our study region. We explored tree height, canopy cover fraction within 20 m (excluding an individual tree's canopy within its crown perimeter), distance from the nearest river, trees per hectare, the mean distance of the 10 nearest trees, slope, aspect, and the fraction of granite within 20 m to investigate the role of tree height and canopy cover, proximity to rivers, tree population density, and topography. We used the tree locations derived from the combined point cloud to compute the tree population density, mean distance of the 10 nearest trees, distance from rivers, and fraction of granite within 30 m. For tree height and canopy cover, we used the tree heights and crown perimeters from 2013. Figure 8 shows our analysis of each of these feature variables for Soaproot Saddle and Lower Teakettle individually, while the analyses for the combined dataset are shown in Figure S8 in the Supporting Information S1. The slopes, y-intercepts, and mortality ranges for each feature variable are shown in Table S2 in the Supporting Information S1.
To assess the correlations among feature variables and our target variable, we computed the correlations between each pair of feature variables and the change in mean relative greenness between 2013 and 2017 for each tree ( Figure S9 in the Supporting Information S1). We chose change in mean relative greenness within the crown perimeter of the tree to have a continuous variable for the tree's health status rather than a binary 0 or 1 for dead or alive when computing the correlations among variables. We found that tree height has the strongest correlation with change in mean relative greenness at the site level ( Figure S9 in the Supporting Information S1). The likelihood of mortality increased by 8.4% for every 10-m increase in height at Lower Teakettle and by 10% for every 10-m increase in height at Soaproot Saddle (Table S2 in the Supporting Information S1). The drier, lower elevation site (Soaproot Saddle) also had stronger positive relationships between tree mortality and trees per hectare, distance to rivers, and canopy cover fraction than did Lower Teakettle (Table S2 in the Supporting Figure 6. The spatial pattern of tree mortality in 2017 averaged on a 30-m-by-30-m grid for (a) Soaproot Saddle and (b) Lower Teakettle. Light yellow pixels indicate 100% tree mortality for the trees in that pixel, while dark blue indicates 0% tree mortality. The aspect is converted from 0 to 360° to a scaled value reflecting North (1) to South (0) and overlayed in white to gray to show the underlying topography.
Information S1). Whereas tree mortality fraction decreases with increasing slope at Soaproot Saddle, tree mortality fraction increases with increasing slope at Lower Teakettle.

Comparison With Previous Work
To make a fair comparison with Stovall et al. (2019), we only compared trees from Stovall et al. (2019) that intersected a crown perimeter of one of our trees. If multiple tree top locations intersected our crown perimeter, we chose the one with the closest tree height to ours. Both studies only consider trees above 5 m in height. The number of trees we found in the overlapping domain of the two studies, the fraction of those trees we were able to match, and the corresponding mortality fraction of each group is shown in Table S3 in the Supporting Information S1. While the tree top locations relative to the lidar point clouds in the two studies are similar in a high-resolution example ( Figure S10 in the Supporting Information S1), Stovall et al. (2019) found more trees overall. The spatial patterns and differences in tree top locations at Soaproot Saddle and Lower Teakettle are shown for representative regions in Figure S11 in the Supporting Information S1.
Because different years were studied by Stovall et al. (2019) compared with this study, the closest comparison of tree mortality after the drought is the cumulative tree mortality in 2016 from Stovall et al. (2019) compared to our 2017 mortality. While Stovall et al. (2019) found a cumulative mortality of 41% for 2016 compared with our 32% from 2017, we can compare only the region of overlap. Using the values from Table S3 in the Supporting Information S1, we find that our cumulative tree mortality fraction in 2017 in the overlapping domain of the two studies is 33% versus 40% for Stovall et al. (2019).
To consider spatial differences in mortality between the two studies, we rasterized the Stovall et al. (2019) dataset to the Landsat grid. Figure 9 shows the mortality fraction from our study, that of Stovall et al. (2019), and the difference between the two studies. In the difference maps (Figures 9c  and 9f), brown areas represent pixels where our tree mortality estimates are higher than those of Stovall et al. (2019), while green areas represent pixels where our estimates of tree mortality are lower.  Stovall et al. (2019), the trees that they labeled as live that we labeled as dead are not necessarily incompatible. These trees may have been alive in 2016 and then died by 2017. However, 125,785 of the trees that Stovall et al. (2019) labeled as dead in 2016 were labeled as alive in 2017 in this study. In other words, 45.4% of the matched trees labeled as dead in the Stovall et al. (2019) study are labeled as alive in our study.
To better understand the differences between the two approaches, we compared the mean NDVI within the crown perimeters of trees we labeled as live, but Stovall et al. (2019) labeled as dead ( Figure 10). The trees where both studies agree are well-separated by mean NDVI. Overall, the distribution of the mismatched trees shown in Figure 10a more closely matches the distribution of the live trees shown in Figure 10b. Because we used the combined point cloud to collect our tree locations, we can track each individual tree over time. We plotted the mean NDVI within the crown perimeters between 2013 and 2021 and the change in the crown area over the same time frame (Figure 11). The mean relative greenness for the trees we labeled as live and Stovall et al. labeled dead more closely resembles the trees both studies labeled as live than those labeled dead. In Figure 11b, we see higher variance in crown area. However, the trees that both studies agree are dead by 2017 show a decline in crown area, with a crown area of 0 m 2 falling within one standard deviation of the mean. The other two categories (where both studies agree the trees are live in blue, and where only our study labels trees as live in gray) decline more slowly. We would still expect to see some decline from 2013 to 2021, since this group of trees includes trees that die between 2017 and 2021.
We compared the tree mortality between 2013 and 2017 rasterized to the Landsat grid to the change in late summer NDMI over the same years. We used the change between 2013 and 2017 to capture the large tree mortality event at the end of the 2012-2016 drought and due to the years of data availability of the NEON datasets (2013, 2017, 2018, 2019, and 2021). The result is shown in Figure S12 in the Supporting Information S1. We found that there is a negative linear association with an R 2 value of 0.37. This indicates that higher tree mortality is associated with a decrease in late summer NDMI consistent with work by Goulden and Bales (2019).

Discussion
We mapped the spatial patterns of tree mortality under drought conditions in the Sierra Nevada. We loosely followed the methodology charted by Stovall et al. (2019) but made several key changes with respect to the datasets used and the ways in which trees were segmented and labeled. We found evidence that fewer trees may have died by 2016 than previously estimated, although there were broadly similar patterns with higher levels of mortality at the lower elevation Soaproot Saddle. Our lower estimate of cumulative tree mortality (32% vs. 41%) could be due to differences in the total number of trees analyzed (we analyzed half the number of trees) and proportion of  low-to-high elevation regions that we considered. For example, trees from the lower elevation site (Soaproot Saddle) were more likely to die, so considering fewer trees from that region will lead to a lower tree mortality estimate. However, within the region where our study overlaps with previous work, our method-which had 93% accuracy on our training dataset-disagreed with 45.4% of the labels for trees labeled dead by Stovall et al. (2019). Comparison of NDVI and crown area time series for this set of mismatches revealed a pattern that was more consistent with the characteristics of live trees identified by both studies.
Differences in spatial patterns from previous research efforts may be due to differences in the preprocessing of the spectral reflectance data. Spectral reflectance data are sensitive to the time of day that the returns are collected. While NEON hyperspectral data is corrected for atmospheric effects and topography (Karpowics & Kampe, 2022), there is not enough information provided with the NAIP imagery to apply this correction. Previous work correcting NAIP imagery using Landsat data revealed significant differences in reflectance indices such as NDVI before and after correction (Zhang et al., 2019). The differences in reflectance in the NAIP imagery at different times of day and different angles of reflection in the highly variable terrain in the Sierra Nevada mountains may account for many of the differences between the tree mortality map presented in this study and previous results by Stovall et al. (2019).
Consistent with previous work, we found that tree mortality risk increases with height. We found a steep increase in mortality risk with height (about 10% per 10 m of height) at the lower elevation site (Soaproot Saddle) up to a tree height of 30 m. Trees taller than 30 m maintained a 70% likelihood of mortality after the drought. At the higher elevation site (Lower Teakettle), we found a more moderate increase in mortality risk with height (about 8.4% per 10 m of height). There was much more tree mortality at the lower elevation Soaproot Saddle site, consistent with previous work on drought-associated tree mortality (Byer & Jin, 2017;Stovall et al., 2019).
After the drought, 2.0% of the trees died between 2017 and 2018, and 2.8% died between 2018 and 2019. Legacy effects from the 2012-2016 drought may have kept mortality levels elevated in subsequent years. Most of these additional increases occurred at Lower Teakettle. When we mapped additional tree mortality for each year of data, we found some inconsistencies, particularly between 2019 and 2021 where trees that we labeled as dead in 2019 were later labeled as live in 2021. On closer examination in QGIS, one region is in an area where there is both visible snow and low luminosity in the surface reflectance, potentially due to a late evening or early morning data collection flight. The low luminosity poses a challenge to remov ing snow-filled pixels by our methods which impose a luminosity threshold. Other areas in the eastern regions of Lower Teakettle appear brown in 2019 but look green and healthy in spectral imagery from July 2021. Because our algorithm only looks at the spectral signature of the trees for a given year, we would not expect these trees to be accurately classified without additional processing.
Because biophysical drivers of tree mortality vary significantly by species (Stephenson & Das, 2020), improvements in species classification models from high-resolution spectral and lidar data may help to create accurate maps of tree mortality risk. Future directions for this work include using the individual tree mortality dataset to build tree mortality models which incorporate non-linear relationships and interactions among biophysical features. Raster data provided with this study may be helpful to scientists working to improve Landsat-based estimates of tree mortality. In addition, these estimates of the impact of drought and subsequent wildfire may assist Figure 10. Comparison of mean crown area normalized difference vegetation index (NDVI) distributions for trees labeled differently and the same between our study and that of Stovall et al. (2019). The x-axis shows mean NDVI of individual trees, while the y-axis shows the relative frequency distributions. (a) We show the mean NDVI relative frequency distributions for trees that were identified as live in this study but dead in Stovall et al. (2019). The boundaries for the middle 99% of the mean NDVI values of these trees are shown by the black vertical bars that carry down into the next panel. (b) The mean NDVI relative frequency distributions of trees that both our study and Stovall et al. (2019) labeled as live are shown in blue, while the relative frequency distributions for trees we both labeled dead are shown in light yellow.

10.1029/2022JG007234
16 of 18 land managers in forest conservation efforts and may provide helpful parameters for dynamic vegetation models which include these effects.

Conclusion
Our study provides an open-source methodological framework for estimating tree mortality from high-resolution lidar and spectral data that can be applied to other forests using hand-labeled training data and parameters specific to the trees in new study sites. We estimate that 49.6% of trees were dead by the end of the 2012-2016 California drought at Soaproot Saddle and 25.5% were dead at Lower Teakettle. While tree mortality at Soaproot Saddle appears to saturate in 2017, we found that trees at higher elevations continue to die at a rate of about 2% per year after the drought subsided. Moreover, we found that subsequent wildfire at Soaproot Saddle increased tree mortality fraction by 49%-89% within burned areas. The datasets provided here may help to constrain dynamic vegetation models to improve our understanding of forest disturbance now and in the future.