Airborne and Spaceborne Lidar Reveal Trends and Patterns of Functional Diversity in a Semi-Arid Ecosystem

Assessing functional diversity and its abiotic controls at continuous spatial scales are crucial to understanding changes in ecosystem processes and services. Semi-arid ecosystems cover large portions of the global terrestrial surface and provide carbon cycling, habitat, and biodiversity, among other important ecosystem processes and services. Yet, the spatial trends and patterns of functional diversity in semi-arid ecosystems and their abiotic controls are unclear. The objectives of this study are two-fold. We evaluated the spatial pattern of functional diversity as estimated from small footprint airborne lidar (ALS) with respect to abiotic controls and fire in a semi-arid ecosystem. Secondly, we used our results to understand the capabilities of large footprint spaceborne lidar (GEDI) for future applications to semi-arid ecosystems. Overall, our findings revealed that functional diversity in this ecosystem is mainly governed by elevation, soil, and water availability. In burned areas, the ALS data show a trend of functional recovery with time since fire. With 16 months of data (April 2019-August 2020), GEDI predicted functional traits showed a moderate correlation (r = 41–61%) with the ALS predicted traits except for the plant area index (PAI) (r = 11%) of low height vegetation (<5 m). We found that the number of GEDI footprints relative to the size of the fire-disturbed areas (=< 2 km2) limited the ability to estimate the full effects of fire disturbance. However, the consistency of diversity trends between ALS and GEDI across our study area demonstrates GEDI’s potential of capturing functional diversity in similar semi-arid ecosystems. The capability of spaceborne lidar to map trends and patterns of functional diversity in this semi-arid ecosystem demonstrates its exciting potential to identify critical biophysical and ecological shifts. Furthermore, opportunities to fuse GEDI with complementary spaceborne data such as ICESat-2 or the upcoming NASA-ISRO Synthetic Aperture Radar (NISAR), and fine scale airborne data will allow us to fill gaps across space and time. For the first time, we have the potential to monitor carbon cycle dynamics, habitats and biodiversity across the globe in semi-arid ecosystems at fine vertical scales.


INTRODUCTION
Understanding the drivers of ecosystem processes and services at regional and global scales provide pivotal knowledge to assess ecosystem responses under changing conditions (Díaz et al., 2007;Isbell et al., 2015). Functional traits, including morphological, physiological and phenological characteristics, depict the growth, reproduction, and survival of individual species and individual plants (Bu et al., 2019). Further, variations in functional traits are shown to be linked to environmental conditions and their fluctuations over time (Schneider et al., 2017). For example, (Ahrens et al., 2020) demonstrated that functional traits such as leaf size and specific leaf area show variation with changes in temperature and precipitation. In addition to climate factors, other influences including human intervention and pathogens can alter functional traits and their spatial distribution (Oliva et al., 2020).
The distribution of functional traits within and between species of an ecosystem represent the demography, resilience, and response strategies to disturbance (Díaz et al., 2004;Poorter and Markesteijn, 2008;Serbin et al., 2019). Therefore, a wealth of research uses functional trait-based approaches to assess ecosystem processes and services including productivity, nutrient cycling and biodiversity (Hooper et al., 2002;Bardgett and van der Putten, 2014;Violle et al., 2014;Wieczynski et al., 2019). For example, direct observations of functional traits are widely utilized at various spatial and temporal scales, mostly in forested ecosystems, to elucidate overall structure, function, and diversity (Funk et al., 2017;Wieczynski et al., 2019). In addition, functional traits have been utilized in assessing community assembly processes across a variety of traits (e.g., Pakeman and Stockan, 2014;Medeiros et al., 2019).
Selection of representative functional traits and scaling those traits are critically important as ecosystem to global scale processes are a function of combined traits of co-occurring species and their abundance (Funk et al., 2017). However, large-scale functional diversity measurements are limited due to the lack of spatially continuous data sets (Jetz et al., 2016), where functional diversity is derived from the distribution of functional trait values at a given spatial grain and extent (Mouchet et al., 2010). Availability of remotely sensed data at fine to coarse spatial and temporal scales facilitates upscaling traits from local to regional and global scales relevant for ecosystem service management (Abelleira Martinez et al., 2016;Braun et al., 2018;Schimel, et al., 2013), and functional leaf traits are now available across biomes in North America (Wang et al., 2020). Among many, the traits that represent canopy architecture (morphological functional traits) show direct relationships between carbon storage (Rödig et al., 2019), habitat distribution and quality (Bae et al., 2019), and biodiversity (Bagaram et al., 2018), and are widely used to characterize regional to global scale ecosystem processes. Importantly, lidar remote sensing provides opportunities to accurately estimate the morphological functional traits in order to map morphological diversity at local and regional scales (Schneider, et al., 2017;Zheng et al., 2021). With the launch of the Global Ecosystem Dynamics Investigation (GEDI) mission, we have new opportunities to map functional traits and biodiversity at global scales in forest ecosystems (Marselis et al., 2019;Rödig et al., 2019;Schimel and Schneider, 2019;Dubayah et al., 2020;Schneider et al., 2020). Nevertheless, the applicability of GEDI in estimating functional traits and diversity in semi-arid ecosystems with short and sparse vegetation is still unknown.
Semi-arid ecosystems cover approximately 40% of the terrestrial landscape and show large dynamicity in ecosystem structure and function (Conti and Díaz, 2013). Hence, semi-arid ecosystems greatly impact global carbon dynamics, productivity, and habitat quality (Poulter et al., 2015). Structure-functioning relationships in frequently disturbed semi-arid ecosystems are unclear, largely due to gaps in spatially continuous data and the weak response of sparse and short height vegetation in optical remote sensing (Kulawardhana et al., 2017;Stavros et al., 2017). In particular, understanding semi-arid ecosystem responses to global change is challenging due to the complex and dynamic interactions among multiple ecosystem functions. Assessing the spatial patterns of functional diversity and its abiotic controls are critically important towards unraveling this complexity (Schlesinger et al., 1990;D´Odorico et al., 2013). Especially, understanding how community assembly is controlled by the balance of abiotic drivers is important in predicting the response of ecosystems to environmental change (Pakeman and Stockan, 2014). Further, functional diversity estimates provide the input data to help constrain model accuracies of ecosystem processes at landscape scales in different regions across the globe (Stavros et al., 2017;Braghiere et al., 2019;Dashti et al., 2021).
In this study, our objectives were to 1) evaluate the spatial pattern of functional diversity estimated from small footprint airborne lidar (ALS) with respect to abiotic controls and fire in a semi-arid ecosystem; and 2) assess GEDI's potential of capturing the semi-arid ecosystem diversity at regional scales using GEDI data from 16 months (April 2019-August 2020) over the same study area. We selected a range of abiotic controls including topography, distance to water, topographic wetness index, and soil. Finally, we explored how disturbance from fire regulates the functional diversity in this semi-arid ecosystem.

Study Area
The study was carried out in the Reynolds Creek Experimental Watershed (RCEW) in southwest Idaho (Figure 1). The study area is characterized by a range of topography (1,100-2,200 m) and vegetation communities. While many varieties of grass, forbs, and shrubs dominate the low elevations, trees of aspen (Populus tremuloides), juniper (Juniperus occidentalis), and Douglas fir (Pseudotsuga menziesii) mark the high elevations. Low stature sagebrush (Artemisia spp.), bitterbrush (Purshia tridentata), and grass of varying densities and cover are found throughout RCEW. In addition, riparian vegetation with cottonwood (Populus trichocarpa) and willow (Alnus, Ceanothus, and Salix species) are found within valleys, and mostly along streams across the watershed. The study area has experienced prescribed and natural fires and supports grazing. As a consequence, invasion of cheatgrass (Bromus tectorum) in native shrub areas and juniper encroachment have occurred in this study area during the last few decades.

Field Data
Reference field data were collected at 43 plots at 10 × 10 m ( Figure 1). The plots were selected using a stratified random sampling approach considering the complicated terrain conditions. We established 5 transects at 1, 3, 5, 7, and 9 m in each plot. We collected canopy heights of all the shrubs, plant area index (PAI) and images at 2 m intervals along each transect totaling 20 measurements per plot using a ceptometer (AccuPAR LP-80, Decagon Devices Inc., Pullman, WA, USA) and a camera (Nikon COOLPIX AW120) respectively. The ceptometer data provides estimates of PAI following (Norman and Welles, 1983). Plot scale PAI was estimated by averaging the 20 measurements (Glenn et al., 2017). The collected images were analyzed using the "SamplePoint" freeware program to estimate the species abundance presented within each plot (v1.59, Booth et al., 2006). Each photo approximately covered 2 m 2 on the ground. We placed 100 equally spaced grid points within each photo and identified the material under each grid point as vegetation species, dead wood, litter or as bare ground which were later used to estimate the percent cover of each material in those plots. Field estimated functional traits were used to assess the accuracy of ALS estimated functional traits.

Abiotic Factors
We selected six abiotic factors identified as potential factors driving the vegetation dynamics in the study area from previous research (e.g. Seyfried et al., 2018;Seyfried et al., 2001). They are elevation, slope, aspect, topographic wetness index (TWI), soil, and distance to water (DTW) (Supplementary Figure S1). The topographic variables of altitude, slope, aspect, and topographic wetness index (TWI) were estimated using the small footprint airborne lidar derived 3 m digital elevation model. The slope data were categorized into 10°groups between 0-90°. The aspect data were categorized into two major directions as north (+/− 90°-180°from south) and south (0°-+/− 90°from south). We used only north-south aspects as the study area shows vegetation growth dynamics primarily on these two aspects (Seyfried et al., 2001). We did not use the flat aspect category as there were not many pixels in this category due to the complex topography of the study area. Stream networks and soil types were retrieved from the Reynolds Creek Critical Zone Observatory database (GIS Server, Reynolds Creek Critical Zone Observatory, 2015). Euclidian distance from perennial and from all streams to each pixel were estimated and used as the distance to water (DTW). We did not use climate data such as precipitation or temperature as the local Frontiers in Remote Sensing | www.frontiersin.org November 2021 | Volume 2 | Article 743320 distribution of these two variables within the study area is primarily controlled by elevation and aspect (Seyfried et al., 2018).

Airborne Lidar Data
Small-footprint waveform lidar data were acquired August 2014 using the NASA Airborne Snow Observatory (ASO) (Riegl, 2017) (RIEGL Laser Measurement Systems GmbH, Horn, Austria) dual laser scanner (wavelength 1,064 nm, outgoing pulse width 3 ns). The study area was scanned at a pulse repetition rate of 400 kHz per laser and the backscattered signal was sampled at 1 ns interval. All the lidar waveforms in each flight line were Gaussian decomposed following the workflow of (Ilangakoon et al., 2018). We estimated spatial coordinates, incident angle, pulse width, amplitude, and scattering cross section of all echoes in each waveform. The resulting average point density across the watershed was 10-14 pts/m 2 with positional accuracies of 0.14 m in vertical and 0.11 m in horizontal directions. The vertical and horizontal positional accuracies were assessed using TerraScan software.

ALS Functional Traits
The point cloud resulting from the waveform Gaussian decomposition was used to estimate the ALS based plant functional traits. Finding the bare ground is critical yet challenging as the low-stature vegetation in the study area tends to widen the waveform without a separate vegetation pulse. We used a pulse width threshold (3.2 ns) to separate bare ground from vegetation based on (Ilangakoon et al., 2018). Then, all lidar returns were height normalized using the ground plane and were aggregated into 10 × 10 m spatial pixels. A spatial resolution similar to the field plot size was used to accommodate both trees and shrubs while also allowing us to identify small variations of functional diversity in shrubdominated areas. Canopy height is the maximum height from all above ground returns from each 10 × 10 m pixel area. PAI was estimated as, where, P(θ) and G(θ) are the gap fraction at incident angle θ and projection function, respectively (Zheng et al., 2016). Ω is the clumping index and L is the true leaf area index (Norman and Campbell, 1989). The gap fraction was estimated as, where σ g is the scattering cross sections from ground returns and σ all is the scattering cross sections from all returns at a given angle. Vegetation sparseness of the study area resulted in very similar spectral signatures for both the bare ground and the canopy (Dashti et al., 2019). We did not correct the PAI for canopy and ground reflectance. Estimates of the scattering cross section from small footprint waveform lidar are explained in Wagner (2010) and Ilangakoon et al. (2018). We predicted the foliage height diversity (FHD) at plot resolution (10 × 10 m). We used plot scale LAI profiles developed only from the above ground point cloud. The above ground lidar points were fragmented into 20 cm layers. The 20 cm layer thickness was used given that it is approximately 2 times the point cloud vertical uncertainty. For each layer, the ratio between numbers of points in each layer to the total number of points was calculated. The total foliage height diversity was then estimated using the following equation (MacArthur and MacArthur 1961).
where, p i is the fraction of foliage in ith canopy layer. Once the three traits were estimated at 10 × 10 m pixel resolution for the whole study area, the trait values were normalized to vary between 0 and 1 following (Schneider et al., 2017).

Functional Diversity
We used three morphological traits (CH, PAI, and FHD) to estimate measures of functional diversity for the study area. We selected CH, PAI, and FHD as they can be directly estimated from both discrete return and waveform lidar and are widely used to represent a range of plant and ecosystem functions (e.g. canopy architecture, ecosystem productivity, habitat quality). We estimated three functional diversity indices: richness, divergence, and evenness, combining the above mentioned three morphological traits. Functional diversity indices were calculated around each sample unit (AKA each ALS pixel) across the watershed. To calculate functional diversity around each ALS pixel, we selected a square neighborhood with a side length of 1,000 m. This spatial neighborhood was selected in order to include at least 10 GEDI footprints within each neighborhood and to be compatible with GEDI raster products . The distribution of PAI, CH, and FHD values of all the pixels within the 1,000 m square neighborhood was then analyzed to derive the functional diversity indices following the below equations from (Villéger et al., 2008;Schneider et al., 2017).
In this study, we describe functional diversity as a 3-D trait space where the axes are represented by the three functional traits. However, functional diversity can be represented as a multidimensional trait space based on the number of traits used to estimate the diversity indices (Villéger et al., 2008, Schneider et al.,2017. The values of functional diversity indices are then calculated based on the relative distribution of functional traits values extracted from each 1000 m spatial neighborhood (Supplementary Figure S2). The functional richness is the convex hull volume covered by the PAI, CH, and FHD of the pixels within the selected neighborhood mapped in the 3-D trait space. Functional divergence (Fdiv) is estimated as: (Schneider et al., 2017), where, S is the number of pixels used to map the functional divergence, dG i is the Euclidean distance between the ith pixel and the centre of gravity of 3D trait space defined by CH, FHD, and PAI as trait axes, and dG is the mean distance of all pixels to the center of gravity (Supplementary Figure S2B). The estimation of Functional evenness (FEve) requires partial weighted evenness (PEW) measurements. The Euclidian distance between nodes was estimated using the minimum spanning tree method in MATLAB (Prim, 1957). In this study, a node is a pixel defined in the 3D trait space by its PAI, FHD, and CH trait axes. The minimum spanning tree was drawn connecting all the pixels in the functional space using the minimum possible total edge weight.
where, EW l is the Euclidian distance of branch l in the minimum spanning tree, and S-1 is the number of branches. Here, S is the number of pixels used and the subscript i represents the 1,000 m square neighborhood space used to estimate the functional diversity (Supplementary Figure S2C).
The resulting values were assigned as the center pixel functional diversity to generate a seamless functional diversity map of the study area.

Statistical Analysis
We first investigated the correlation among functional diversity indices estimated from ALS and the abiotic factors. To analyze the effects of abiotic factors on functional diversity, we randomly selected 200 data points across the watershed covering the full range of functional diversity and the abiotic factors. We resampled the environmental factors (elevation, slope, aspect, topographic wetness index, distance to the nearest stream, and soil type) to the same 1,000 m spatial resolution of the functional diversity. We used the average of the continuous values (e.g. elevation). For categorical variables (e.g. soil type), we converted them into dummy variables and used the major category within each 1,000 m space. We set the minimum distance between data points greater than 1,020 m to avoid the mutual inclusion of niche spaces. A Gaussian generalized linear model with an identity link function was applied to the scaled continuous and dummy categorical variables with the diversity indices as the response variable to evaluate the effects of abiotic factors on functional diversity. The relative importance of each factor in the linear model was used (R package relaimpo, calc. relimp) to assess the capability of all factors together and each factor separately to explain the variance of each diversity index estimated from ALS and GEDI.

Functional Diversity Shift Under Fire Disturbance
To study how fire may shift the functional diversity, we applied a "space for time" approach to four different fires distributed across the study area. These are the Koke (natural fire, 2014), Whiskey mountain (prescribed fire, 2005), Break (prescribed fire, 2002) and Rabbit creek (natural fire, 1996). The burned areas varied from 0.5 to 2 km 2 . To assess the disturbance driven diversity shift, we generated circular buffer zones with 500 m distance in the unburned areas around each fire. The functional traits and diversity values of all the ALS pixels within each burned area and their surrounding unburned areas were then extracted. Trait and diversity indices of each burned pixel and the means of undisturbed values from each fire were used to predict the percent recovery. The resultant percent trait and diversity within the fires and the original traits and diversity values from the undisturbed areas were separately fitted with non-linear functions (Gompertz, Frontiers in Remote Sensing | www.frontiersin.org November 2021 | Volume 2 | Article 743320 5 polynomial, and Chapman-Richards) and the best fit models were used to interpret the post-fire recovery trajectory. Given the limited number of fires in the study area, we did not evaluate other potential factors for fire recovery rates (e.g. impact of environmental gradient on fire recovery rate). To understand the marginal spatial neighborhood size that can show the disturbance driven diversity shifts, we estimated the functional diversity with varying spatial neighborhood sizes (30-500 m).

GEDI Data and Functional Traits
To evaluate the potential of GEDI data to capture functional diversity in semi-aid ecosystems, we used GEDI 1B (v001), GEDI 2A (v001) and GEDI 2B (v001) data between April 2019 to August 2020. Data were downloaded from the Land Processes Distributed Active Archive Center (LP DAAC) using the GEDI finder tool (https://lpdaacsvc.cr.usgs.gov/services/gedifinder). In this study, we used GEDI footprints with a quality flag of 1 (Sensitivity >0.90). Horizontal geolocation accuracy of the first version GEDI data are approximately 20 m. As this may greatly impact our results, using the correlation matching algorithm (Blair and Hofton 1999) implemented in the GEDI simulator  and the ALS point cloud data, we geocorrected the GEDI data before estimating the functional diversity indices. The GEDI simulator simulates GEDI-like waveforms using the ALS point cloud and compares them with the real GEDI waveforms within a given horizontal region (∼20 m in this case). The location that has the highest correlation between the simulated and the real GEDI waveform provides the offset of the real GEDI waveform. After geocorrection, we extracted functional traits that resembles the same functional traits we derived from ALS data. The RH98 metric from the GEDI 2A product was used as the canopy height. LAI and FHD from the GEDI 2B product were used as the PAI and FHD, respectively, in this research. We used GEDI 2A and 2B products based on the default waveform processing algorithm setting. In addition, several other parameters (beam and sensitivity) were extracted to evaluate the factors affecting the accuracy of GEDI estimated functional traits. In this study, we used ALS based ground slope, ALS estimated vegetation height, ALS estimated cover, and GEDI beam (power vs coverage) and sensitivity as factors that may affect the retrieval of GEDI functional traits following (Liu et al., 2021). A relative importance matrix from a random forest regression model was used to assess the impact of the above-mentioned factors on GEDI functional traits retrievals. Given the limited number of field samples that overlap with the real GEDI footprints, we used ALS data as the reference data to assess the accuracy of GEDI predicted functional traits.

ALS Derived Functional Traits
Our field observations showed a range of functional trait values for PAI (0.19-1.88) and CH (0.31 m-2.52) with mean values of 0.76 and 1.20 m, respectively. Figures 2A,B show the correlation between the functional traits estimated from 10 × 10 m field plots and traits predicted from ALS point clouds within the field plots. Field observed functional traits show a strong correlation (r 71-88%) with ALS derived functional traits. However, the ALS underestimates the field observed vegetation heights and overestimates the PAI estimates. The ALS height underestimation is compatible with previous ALS predicted heights in similar study areas (Streutker and Glenn, 2006;Mitchell et al., 2011). This ALS height underestimation could be due to the uncertainty in bare ground detection . The PAI overestimation could be due to reflectance similarities between the canopy and bare ground.

Environmental Controls of Functional Diversity
ALS estimated functional richness values ranged between 0.0-0.11 ( Figure 3A). The evenness and divergence show relatively larger values compared to the richness. The ALS evenness and divergence range between 0.31-0.99 and 0.55-0.86 ( Figures 3B,C), correspondingly. The dense treedominant southern portion of the study area is characterized by large richness and small evenness values whereas the shrubdominant central and northern portions of the watershed are characterized by the opposite. Evenness and richness differentiate the tree-shrub ecotones and relative densities within shrub dominant areas. A variable functional divergence is observed throughout the watershed ( Figure 3C).
The trends in functional diversity indices with abiotic factors are shown in Figure 4. Functional richness slightly decreases with an increase in elevation until 1,600 m and then has a rapid increase with increased elevation ( Figure 4A). The ALS predicted functional richness increases from south to north aspects ( Figure 4B). ALS predicted functional evenness values decrease with a change from south to north aspects. In contrast, functional evenness values increase when the distance to water increases. There was no trend observed when we considered both perennial and ephemeral streams together to estimate the distance to water. Among all environmental factors studied, distance to water is the only factor that shows a trend with functional divergence predicted from the ALS data. This trend is generally positive, with functional divergence increasing with an increase in distance to water, especially at higher elevations ( Figure 4F).

Functional Diversity in Fire Disturbed Areas
We analysed ALS functional traits and diversity indices over selected burned and unburned areas to infer the fire resilience of the study area. We used a fire chronosequence spanning 1996-2014 to understand if there are shifts of post-fire traits and diversity through vegetation recovery. One way ANOVA indicated that all burned and unburned groups are significantly different (95% confidence level) (Supplementary Figure S3). Based on the tested models, a 2nd order polynomial function provided the best fit and is used hereafter to describe the post-fire recovery. Based on ALS data, CH shows an initial decrease followed by an increase towards the undisturbed state (100% recovery state) ( Figure 5A). PAI and FHD show a faster recovery  to the undisturbed state (∼18 years post-fire); however, CH does not reach the undisturbed state during the study period. Though functional richness increased initially, a drastic decrease is demonstrated ∼9 years post-fire. In contrast, functional evenness continuously increases ( Figure 5E). This post-fire increase in functional evenness indicates that fire may changes the state of the ecosystem driving it to a more functionally even ecosystem. Functional divergence shows a similar behaviour to functional richness ( Figure 5F).
We further estimated the functional diversity of the four fires at a range of spatial neighbourhood sizes with 30 m increments using ALS predicted functional indices (Supplementary Figure  S4). To calculate functional diversity, three or more measurements from each functional trait are considered. The minimum neighbourhood of 30 m was selected as it is 3 times bigger than the ALS pixel size. By selecting a square neighbourhood, we were constrained to 30 m being the minimum we could use. After the minimum neighbourhood was selected, the increment can be any size. We used 30 m intervals. At smaller neighbourhoods, functional richness from disturbed and undisturbed areas cannot be differentiated (Supplementary Figure S4). Functional richness values and their uncertainty increases with an increase in the spatial neighbourhood size. On the contrary, at smaller neighbourhoods, the differences in functional evenness and divergence between disturbed and undisturbed areas are clearly visible. These differences gradually diminish with an increase in the spatial neighbourhood (∼270 m) (Supplementary Figure S4). The uncertainty of functional evenness and divergence estimates also decrease when the spatial neighbourhood size increases. The results further indicated that functional diversity values resulted from each selected neighbourhood are significantly different (Supplementary Figures S6-S8).

GEDI Predicted Functional Traits in a Semi-Arid Ecosystem
Our random forest regression based relative importance values show that ground slope is the most important error factor for FIGURE 5 | Percent recovery of functional traits (A-C) and diversity (D-F) across the fire chronosequence (1996-2014) predicted from ALS data. The box plots represent the percent recovery data distribution within a given year. The black line and blue shaded area represent the mean percent recovery and 95% confidence interval of functional trait/diversity post-fire, respectively. The box plots show the raw data (functional traits and diversity) distribution extracted from burned areas colored based on time since fire.
Frontiers in Remote Sensing | www.frontiersin.org November 2021 | Volume 2 | Article 743320 8 all three GEDI based functional trait estimates ( Figure 6). In addition, ALS estimated canopy cover and canopy height influenced the GEDI trait retrievals. The ground slope varied between 0 and >60°in our study area. Due to this known influence of ground slope on large footprint lidar backscatter waveforms (Silva et al., 2018;Dong et al., 2019), we further investigated how GEDI predicted functional traits vary with ground slope. The distribution of ALS and GEDI estimated functional traits with slope are displayed in Supplementary Figure S9. The results indicate that GEDI predicted canopy height (r 26%) and FHD (r 51%) show a positive correlation with ground slope. This moderately strong correlation of FHD with slope could be mainly due to the shorter heights of the vegetation (<5 m) where the ground and vegetation signals likely be in mixed. We expect this slope-FHD correlation to be smaller in areas with taller vegetation. Based on our study and other related studies (Schneider et al., 2020;Liu et al., 2021), we used GEDI data only from areas with slopes less than 20°for the analysis hereafter.
To compare GEDI functional traits with those predicted from ALS, we used the same surface coverage for both ALS and GEDI (25 m). We used the GEDI simulator based geocorrection results to find overlapping ALS and GEDI footprints. We then estimated ALS functional traits at 25 m spatial resolution from those best matched ALS point clouds. GEDI functional traits show low to moderate correlation (r 11-61%) with ALS predicted functional traits (Figure 7). Our results indicate that GEDI data have difficulty estimating vegetation heights less than 2 m (red circle in Figure 7A). Further, GEDI data overestimate FHD (RMSE 2.45), canopy heights (RMSE 4.69 m) and the PAI (RMSE 1.97) in this study area. This overestimation can be attributed to its large pulse length relative to the shorter vegetation heights in this area. In addition, the differences in the algorithms used to estimate the functional traits from the two laser systems (ALS and GEDI) may further cause these overestimations.

GEDI Predicted Functional Diversity in a Semi-Arid Ecosystem
GEDI shows small functional richness values, ranging between 0.001-0.07 ( Figure 8A). The evenness and divergence show relatively larger values compared to the richness. The GEDI evenness vary between 0.55 and 0.86 ( Figures 8B,C). As in ALS diversity maps, GEDI evenness and richness differentiate the tree-shrub ecotones and relative densities within shrub  Frontiers in Remote Sensing | www.frontiersin.org November 2021 | Volume 2 | Article 743320 9 dominant areas. Compared to ALS functional divergence, discrete high functional divergence patches are observed throughout the watershed ( Figure 8C). The correlation between ALS and GEDI predicted functional richness, functional evenness, and functional divergence are 24, 20 and 5%, respectively ( Figures 8D,E). In addition, functional richness predicted from both ALS and GEDI show consistent inverse correlation (∼60%) with their respective functional evenness (Supplementary Figure S10). Though we observed moderate correlations between GEDI and ALS predicted functional traits (CH and FHD), these low correlations of functional diversity could be mainly due to the sparseness of GEDI footprints within each 1,000 m spatial neighbourhood (∼20-100 footprints compared to ∼1,600 pixels from ALS).
However, there are consistent trends between ALS and GEDI diversity with abiotic factors. As in ALS, functional richness slightly decreases with an increase in elevation until 1,600 m and then has a rapid increase with increased elevation ( Figure 9A).
The GEDI predicted functional richness shows an initial decrease on south aspects and an increase on north aspects ( Figure 9B). The functional evenness shows a consistent trend with both the ALS and GEDI data ( Figure 9C). The GEDI predicted functional evenness values decrease with a change from south to north aspects. The GEDI functional evenness and divergence resulted in with similar patterns as observed with ALS predicted evenness and divergence respectively ( Figure 9F).
Due to the limited coverage of GEDI data (April 2019-August 2020), we did not have enough (>10) GEDI footprints over the burned areas to evaluate the post-fire functional shifts at GEDI scale in this semi-arid ecosystem. However, we emphasize the importance of such a study with FIGURE 8 | Functional diversity mapped using a 1,000 m spatial neighborhood GEDI data and their correlation with ALS derived functional diversity indices. Top row: GEDI predicted (A) functional richness (GEDI FRic), (B) Functional evenness (GEDI FEve), (C) Functional divergence (GEDI FDiv). Bottom row: Correlation between (D) ALS functional richness (ALS FRic) and GEDI predicted functional richness (GEDI FRic), (E) ALS functional evenness (ALS FEve) and GEDI predicted functional evenness (GEDI FEve), (F) ALS functional divergence (ALS FDiv) and GEDI predicted functional divergence (GEDI FDiv).
Frontiers in Remote Sensing | www.frontiersin.org November 2021 | Volume 2 | Article 743320 appropriate GEDI coverage (two or potentially more years of data collection) in the future, as the GEDI mission aims to address global biodiversity and habitat change, which are directly related to the functional diversity.

ALS Based Functional Diversity in Semi-Arid Ecosystems
Our results emphasize that the relative distribution of elevation, soil, distance to water, and aspect control the trends and patterns of functional diversity in this ecosystem. Functional richness is the niche extent in the trait space. We observed large functional richness at high elevations, especially in the southern portion of the study area. This region receives the highest precipitation throughout the year, mostly as snow, which helps maintain perennial streamflow (Harrison et al., 2020). The major soil type is the Harmehl-Gabica association characterized by deep soil profiles, thick surficial A horizons, and high organic matter on north facing slopes (Seyfried et al., 2000). Due to the availability of water, soil nutrients, and solar radiation, high elevation areas provide favorable conditions for conifers, aspen and mountain sagebrush and these species have diverse morphological characteristics. The diverse topographic gradients over short distances limits the spatial extension of this functionally rich zone. Hence, the vegetation in the southern portion of the study area may have higher competition for resources leading to a lower overall functional evenness. In addition, this region shows a mean annual temperature range from 4-16°C providing favorable conditions for vegetation growth. Due to a strong dependency with elevation, we did not include precipitation or temperature as abiotic factors. Future studies including these climatic variables and their temporal variability will help us understand the temporal and spatial dynamics of functional diversity. Contrary to this study, a decrease in functional richness at high elevations due to increases in aridity and decreases in temperatures have been observed in several other studies (Schneider et al., 2017;Durán et al., 2019;Wieczynski et al., 2019). In contrast, the lower elevations in the central and north portions of the study area are relatively dry with lower and intermittent precipitation. The northern watershed consists of shallow, rocky soils with mesic soil temperatures favoring FIGURE 9 | Trends between functional diversity and abiotic controls. The red and blue curves represent the mean variation of diversity indices predicted from ALS and GEDI data, respectively, and the surrounding gray area represents the standard deviation. Functional richness with (A) elevation; (B) aspect. Functional evenness with (C) elevation; (D) aspect; (E) distance to water. Functional divergence with (F) distance to water.
Frontiers in Remote Sensing | www.frontiersin.org November 2021 | Volume 2 | Article 743320 sagebrush species. This unique and important plant community in the northern Great Basin ecoregion is referred to as the sagebrush-steppe and is co-dominated by big sagebrush and several perennial grasses and forbs. The functional evenness is the distribution of trait abundance in the occupied niche space. Our functional trait maps show that variability in traits in these regions is minimal. Rather than the abundance, the distribution of structurally similar vegetation in the niche space is reflected by the evenness. This further explains the negative correlation between functional richness and evenness. Sparsely distributed, structurally similar shrubs can effectively utilize the entire range of resources available. Hence, we observe functionally even, but functionally less-rich landscapes in the lower elevations of the study area. This low functional richness may also be a reflection of environmental tolerance (beta niches) to disturbance such as grazing (Tilman and Downing 1994). In contrast, the water availability along perennial streams favors the growth of dense, structurally variable riparian vegetation contributing to a relatively low functional evenness. Our results of this increased plant diversity with wetter environmental conditions in water limited environments are confirmed by Harrison et al. (2020). This soil and elevation dependency in functional richness and evenness in this study area are further confirmed by our generalized linear models and are consistent with both ALS and GEDI (Figure 10 A, B, D, and E and Supplementary Figure S10). Remotely sensed functional diversity indicates how the niche is grouped in trait space. Though we observe similar morphological traits at lower elevation and diverse morphological traits at higher elevations, the local grouping or the patchiness of the vegetation is mainly governed by the soils. Our relative importance values resulting from generalized linear models also confirmed elevation, soil, and distance to water as important for functional divergence (Figures 10C, F). The less relative importance of distance to water with the GEDI predicted data could be due to sparseness of the GEDI footprints across the watershed. In addition, water availability provides the necessary wetness for vegetation to grow and survive. Hence, different densities of vegetation can be observed in the trait space leading to dynamic divergence characteristics in the study area, in all ecotones. This explains why the functional divergence shows a trend only with distance to water. The variable functional diversity in the study area may resemble the maximum use of available resources. Overall, this study shows that both the fine resolution airborne lidar can predict functional richness, evenness and divergence of the semi-arid study area and the functional diversity in this study area is mainly driven by the elevation, soil, and water availability.

Functional Diversity Shifts Under Disturbance
Disturbance and time-since-disturbance show a significant influence on driving the functional diversity in the study area. Further, the functional diversity in the fire regions is highly dependent on the spatial neighborhood size as well as the size of fire. While we did not have enough GEDI footprints to characterize the relatively small fires, our GEDI simulations indicate that after 2 years of the mission, GEDI coverage is likely to be sufficient (Ilangakoon, 2020). Similar functional richness over disturbed and undisturbed areas from ALS data may be due to standing dead or burned wood in the disturbed area. Integrating other remote sensing techniques that can differentiate woody materials from leafy vegetation could resolve this (e.g. Li et al., 2018). However, the increase of functional richness of both disturbed and undisturbed areas with the expansion of neighborhood size is due to the potential merging of different traits associations. Larger neighborhood sizes include areas outside of the burn. However, this analysis provides a change in functional diversity compared to purely unburned areas. This helps to understand if an area has undergone a disturbance even if the neighborhood size is bigger than the size of disturbance. This can be especially helpful when using coarse spatial resolution data where the burned areas are only represented by a few pixels (e.g. wall to wall data from MODIS). All three functional traits and diversity indices show a clear separation between disturbed and undisturbed areas. The differences disappear at large spatial extents though, where functional evenness and divergence of burned and unburned areas converge. This convergence emphasizes the importance of using fine resolution remote sensing data such as ALS to estimate fine-scale disturbance effects. At global scales, this could be potentially accomplished fusing GEDI data with ICESat-2 and Synthetic Aperture Radar (SAR), including the upcoming NISAR (Bae et al., 2019;Dubayah, 2016, 2019). Though functional richness and divergence show relatively lower values immediately post-fire compared to the undisturbed areas, functional evenness increases immediately after fire. Disturbance results in a functionally different and disconnected landscape from its surrounding. The rate of recovery largely affects native habitat, especially for sage grouse (a threatened species), and potential changes in ecosystem processes. The ideal habitat for sage grouse requires connected mosaics of sagebrush steppe, which allow them safe migration, secure shelter, and food resources (Stiver et al., 2015). In addition, disturbance-driven fragmentation strongly affects biodiversity and resource distributions (e.g. ecological functions and processes), especially along edges of the disturbance (Collinge, 1996). Recovery from disturbance in the face of global change represents a substantial challenge to agencies that manage these lands. Our study shows that assessing functional diversity can help identify areas for restoration or other management activities to consider for treatment.
Stability and post-disturbance recovery of individual species and the community are important components of resilience to disturbance. For example, ecological resilience, an important aspect used to understand ecosystem resilience, is the estimate of post-disturbance vegetation structure recovery (Cole et al., 2014). The resilience of this ecosystem to wildfire thus can be gauged by the rate at which the vegetation recovers in structure and diversity. Our results show that this ecosystem recovers to its background state within 10-15 years, both in terms of functional traits and their diversity (except FHD). The fire return interval of sagebrush-dominated ecosystems can vary from decades to centuries (Baker, 2006;Heyerdahl et al., 2006;Mensing et al., 2006) though the invasion of non-native grass like cheatgrass can shorten fire return intervals (Balch et al., 2013). Based on historic research and our findings, the study area is likely in a stable and resilient state to wildfires. However, with increasing fire size, frequency, and severity, and nonnative grass invasion in the western US, this ecosystem may also take an alternative unstable state. Our results further emphasize the unprecedented opportunities to characterize ecosystem resilience in addition to vegetation structure to assess wildfire impacts on ecosystems across spatially continuous scales.

GEDI in Semi-Arid Ecosystems
Functional traits and their diversity in an ecosystem can indicate resource availability and response mechanisms to environmental perturbations such as fires. The potential of GEDI to predict functional traits and diversity of semi-arid ecosystems thus has important global implications. While the correlation between GEDI and ALS estimates of canopy height and FHD confirmed the potential of GEDI in these heterogeneous semi-arid ecosystems, the poor correlation with ALS predicted PAI emphasized the need to further investigate GEDI's feasibility in low-height vegetation environments with full GEDI coverage. Our study area mainly consists of short shrubs and grass (<5 m) where vegetation return is likely to be mixed in with the ground return. These mixed signals may incorrectly classify vegetation as ground leading to this PAI underestimate. We expect stronger correlation between ALS and GEDI predicted PAI in biomes where there is a clear inflection point between the ground and canopy (in areas where vegetation is taller than 5 m). However, the correlations of field estimated functional traits with those estimated from both the ALS and simulated GEDI data (Ilangakoon, 2020) provide a baseline to understand GEDI's feasibility and challenges to track trends and patterns of functional traits at regional and global scales in similar ecosystems. In this study, we used real GEDI data with the default algorithm setting for vegetation metric calculations. Our findings may improve by implementing algorithms optimized for the study area, correcting for known differences between the two instruments (e.g., pulse widths of GEDI 15.6 ns, ALS 3 ns), and with additional GEDI coverage over time. In addition, correction of ALS functional trait estimates for the reflectance differences between the bare ground and the canopy of different plant functional types may help to constrain the ALS predicted PAI overestimation compared to the field estimates, as well as GEDI estimates. GEDI data version 1 used a default global reflectance ratio (1.5) between the canopy and bare ground while we considered a value of 1 in the ALS estimates due to the observed similarity between the spectral signatures of the ground and the canopy of the dominant plant functional type (shrub). However, the study area also has grass and trees which can cause the ratio of canopy and ground reflectance to vary accordingly. Though we used the LAI profile estimated from the ALS and GEDI data to calculate FHD, the vertical resolution we used for ALS-FHD estimate is 0.20 m whereas it is 5 m for GEDI-FHD estimates . We used 0.20 m layers from the ALS point cloud in order to capture fine-scale variation of FHD and to evaluate GEDI's feasibility in this ecosystem. The differences in vertical resolution and reflectance ratios between the two datasets may help explain the low to moderate correlation and uncertainties between the ALS and GEDI trait predictions. Derivation of FHD at higher vertical resolution using original georeferenced GEDI waveforms (GEDI 1B product) may provide better FHD values for these types of low height vegetation environments. Further, we observed approximately 48% of GEDI footprints in this study area are on slopes >20°and these slopes influence our GEDI based trait estimates. Future investigations may choose to utilize a slope adaptive method to correct for the influence of the ground slope in waveform lidar predicted metrics as explained in Wang et al. (2019) or applying a physical slope correction (Park et al., 2014). While 16 months of GEDI sampling is not sufficient to provide a comprehensive description of the disturbance, GEDI's capability of estimating functional traits (CH and FHD) in this semi-arid ecosystem demonstrates its exciting potential to identify critical biophysical and ecological shifts. Footprintbased GEDI functional diversity estimates showed consistent trends with abiotic factors as with ALS data. Though the sparseness of GEDI footprints (∼20-100 footprints within a 1 km 2 area) limited our ability to map functional diversity at local and regional scales, this could be potentially accomplished fusing GEDI data with ICESat-2 and Synthetic Aperture Radar (SAR), including the upcoming NISAR (Bae et al., 2019;Qi and Dubayah, 2016;Qi et al., 2019). This information will help us understand the potential of GEDI to monitor the changes in carbon-cycle dynamics, habitats and biodiversity across the globe in semi-arid ecosystems. The errors in GEDI predicted functional traits in a semi-arid ecosystem such as our study area will result in errors of secondary products such as functional diversity, biomass, and carbon stocks. Hence, it is important to identify methods and algorithms with properly documented workflows appropriate for the ecosystem, such as our efforts here with functional traits. Ultimately, this study highlights the capabilities and opportunities for further research and refinement in order to utilize GEDI for global scale estimates of biomass, carbon, and biodiversity.

DATA AVAILABILITY STATEMENT
Publicly available datasets were analyzed in this study. This data can be found here: 10.3334/ORNLDAAC/1503.

AUTHOR CONTRIBUTIONS
NI and NG designed the study. NI, NG, HD, and LS processed and analyzed data. NI, NG, FS, HD, SH, LS, and T.G. wrote the paper.

FUNDING
This work was supported by the NASA Terrestrial Ecology (NNX14AD81G) and NASA Earth and Space Science Fellowship 2017 (NESSF 2017) (17-EARTH17F-0209) programs. The research carried out at the Jet Propulsion Laboratory, California Institute of Technology, was under a contract with the National Aeronautics and Space Administration (80NM0018D0004). This material is also based in part upon work supported by the National Science Foundation through the NEON Program. The National Ecological Observatory Network is a program sponsored by the National Science Foundation and operated under cooperative agreement by Battelle Memorial Institute.