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.

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² 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.

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 distribution of these two variables within the study area is primarily controlled by elevation and aspect (Seyfried et al., 2018).

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² 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.

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 shrub-dominated areas. Canopy height is the maximum height from all above ground returns from each 10 × 10 m pixel area. PAI was estimated as, PAI ( θ ) = L Ω ( θ ) = −   ln ( P ( θ ) ) cosθ G ( θ ) 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, P ( θ ) = ∑ σ g , θ ∑ σ all , θ 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). FHD = − ∑ p i  . ln p i 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).

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 (Dubayah et al., 2020). 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: Δ | d | =   ∑ i = 1 S 1 S   | d G i − dG ¯ |    ( Schneider   et   al .,   2017 ) , FDiv =   dG ¯ Δ | d | + dG ¯    ( 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 d G ¯ 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. PE W i = E W l ∑ i = 1 S − 1 E W l FEve =   ∑ l = 1 S − 1 min ( PE W l   , 1 S − 1 ) −   1 S − 1 1 − 1 S − 1 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.

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.

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². 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, 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).

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 (Hancock et al., 2019) and the ALS point cloud data, we geo-corrected 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 geo-correction, 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.

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 (Glenn et al., 2011). The PAI overestimation could be due to reflectance similarities between the canopy and bare ground.

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 tree-dominant southern portion of the study area is characterized by large richness and small evenness values whereas the shrub-dominant 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).

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).

Our random forest regression based relative importance values show that ground slope is the most important error factor for 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 geo-correction 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 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 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 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.

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 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.

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; Qi and 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 non-native 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.

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 (Tang and Armston, 2019). 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. Footprint-based 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² 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.