Calibration of GEDI footprint aboveground biomass models in Mediterranean forests with NFI plots: A comparison of approaches

Research article Calibration of GEDI footprint aboveground biomass models in Mediterranean forests with NFI plots: A comparison of approaches Adri´ an Pascual a,* , Paul B. May b , Aar´ on C´ ardenas-Martínez c , Juan Guerra-Hern´ andez d , Neha Hunka a , Jamis M. Bruening a , Sean P. Healey e , John D. Armston a , Ralph O. Dubayah a a Department of Geographical Sciences, University of Maryland, College Park, MD, USA b Department of Mathematics, South Dakota School of Mines and Technology, Rapid City, SD, USA c Departamento de Geografía Física y An´ alisis Geogr´ afico Regional, Universidad de Sevilla, 41004, Seville, Spain d Forest Research Centre, School of Agriculture, University of Lisbon, Instituto Superior de Agronomia (ISA), Tapada da Ajuda, 1349-017, Lisboa, Portugal e US Forest Service Rocky Mountain Research Station, Ogden, UT, 84401, USA ARTICLE INFO Keywords: Spaceborne lidar GEDI calibration Biomass modelling Inventory plots ABSTRACT Observations from the NASA Global Ecosystem Dynamics Investigation (GEDI) provide global information on forest structure and biomass. Footprint-level predictions of aboveground biomass density (AGBD) in the GEDI mission are based on training data sourced from sparsely distributed field plots coincident with airborne laser scanning surveys. National Forest Inventories (NFI) are rarely used to calibrate GEDI footprint biomass models because their sampling and positional accuracy prevent accurate colocation with GEDI or ALS. This omission can limit the harmonization of jurisdictional biomass estimates from NFI’s and GEDI; however, there are methods available to improve the colocation of NFI plots with GEDI footprints. Focusing on Mediterranean forests in Spain, we compared different approaches to the collocation of NFI and GEDI data: (i) simulated waveforms from ALS; (ii) nearest-neighbor on-orbit GEDI waveforms; and (iii) imputed GEDI waveforms imputed to NFI plot locations using a novel geostatistical method. These methods are potential solutions to improve the local performance of biomass models and address potential local systematic deviations between GEDI and NFI estimates. We assess the advantages and limitations of these methods to locally calibrate GEDI biomass models and quantify the impact of geolocation errors in reference NFI plot data. The new biomass models from each method were used to predict footprint level AGBD, which were then gridded for a province in the North-West of Spain. It was found that the imputation approach is not sensitive to common errors in NFI plot geolocation, but it can outperform ALS-based simulation in some cases, highlighting the benefit of information from multiple GEDI footprints proximate to NFI plots for improving biomass predictions. This research provides users with benchmark of available techniques to locally-calibrate GEDI footprint biomass models. 1. Introduction One of the critical aspects of terrestrial ecosystems monitoring is the estimation of above-ground biomass (AGBD) using three-dimensional (3D) vertical vegetation structure metrics (Bastos et al., 2022; Li et al., 2023). Climate policies, monitoring programs and environmental management agencies use biomass stocks and fluxes to inform on the effects of climate change and human disturbances (Dubayah et al., 2022a; Friedlingstein et al., 2022; Ma et al., 2023). The Global Ecosystem Dynamics Investigation (GEDI) mission has made available billions of biomass predictions across temperate and tropical ecosystems during the last four years (Kellner et al., 2022), dramatically expanding the ability to estimate carbon stocks and fluxes (Patterson et al., 2019). The GEDI estimation of biomass is globally consistent to reference estimates (Armston et al., 2023) despite gaps across biomes and continents in the GEDI Forest Structure and Biomass Database (FSBD) used to build GEDI biomass models (Duncanson et al., 2022), and considering the substantial heterogeneity in biomass stocking (Bruening et al., 2023; Crockett et al., 2023). GEDI spatial sampling density is unmatched by any existing * Corresponding author. E-mail addresses: [email protected] (A. Pascual), [email protected] (J. Guerra-Hern´ andez), [email protected] (N. Hunka), [email protected] (J.M. Bruening), [email protected] (S.P. Healey), [email protected] (J.D. Armston), [email protected] (R.O. Dubayah). Contents lists available at ScienceDirect Journal of Environmental Management journal homepage: www.elsevier.com/locate/jenvman https://doi.org/10.1016/j.jenvman.2025.124313 Received 25 October 2024; Received in revised form 20 January 2025; Accepted 21 January 2025 Journal of Environmental Management 375 (2025) 124313 Available online 31 January 2025 0301-4797/© 2025 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY-NC-ND license ( http://creativecommons.org/licenses/bync-nd/4.0/ ). spaceborne sensor, but still sparse enough that the chances of spatial colocation between GEDI footprints and field measurements is minimal. GEDI biomass estimates are not harmonized to the National Forest Inventories (NFI) plot-level data that countries use for monitoring and reporting. A harmonization between NFIs and GEDI data is important to understand forest biomass stocking (Knapp et al., 2020; Araza et al., 2022; Hunka et al., 2023 Pascual et al., 2023) and to tailor global GEDI estimates to local conditions, paving the way to provide solutions to mission, reduction and verification (MRV) projects over territories lacking activity data and where GEDI is the only source to sustain projects and forest management through time (Neeff, 2021; Marino and Bautista, 2022). Airborne laser scanning (ALS) data, when available and scheduled combined with field inventories, is an incumbent method to calibrate biomass. Indeed, GEDI biomass models were built following the approach (Duncanson et al., 2022; Kellner et al., 2022), but the approach relies on ALS and when available, it rarely matches the acquisition of biomass data in most NFI programs. Avoiding the ALS dependency can make many thousands of inventory plots actionable for biomass calibration. The scale of 12.5-m radius footprints is substantially smaller than NFI plots and this exacerbate disparities between GEDI estimates and NFI biomass values. Despite scaling mismatches between plots and footprints, several strategies have been pursued to calibrate biomass using on-orbit and simulated GEDI data. When ALS is available, users can recreate the GEDI biomass model building using NFI biomass as response variable (Pascual et al., 2023) and a single ALS-simulated footprint as predictor information. However, ALS-noncontingent solutions are more applicable. To date, two ALS-independent methods have been applied to calibrate biomass using GEDI. One is the matching of plots to GEDI data using assignment methods such as nearest neighbor (e.g., Zhang et al., 2022; Bullock et al., 2023). Alternatively, there is a novel probabilistic imputation approach to infer GEDI canopy height at any given location presented in May et al. (2024). These two methods rely on the distribution of multiple GEDI measurements to create a pair between the best match to an NFI plot (neighbors) and to impute canopy height estimates using footprints around NFI locations (imputation). We explored the advantages and limitations of these biomass calibration alternatives based on using 25-m GEDI footprint estimates of canopy height as predictors. Our study area is in Spain, where high-quality ALS surveys and GEDI estimates are temporally-coincident to an NFI program that is improving geolocation accuracy (Pascual et al., 2021; Pascual and Guerra-Hern´ andez, 2023). Geolocation errors in GEDI footprint data exist (e.g., Roy et al., 2021; Li et al., 2023; Tang et al., 2023; Holcomb et al., 2024). However less attention is paid to field plots positioning accuracy (Johnson and Barton, 2004; Mauro et al., 2011). Geolocation errors in NFI surveys can be problematic when inference relies on remotely-sensed metrics (Babcock et al., 2018; Andersen et al., 2022; Pascual et al., 2020; Silveira et al., 2023), and these have not been addressed in the few studies that have tried to calibrate biomass using GEDI metrics. For instance, the study from Bullock et al. (2023) paired GEDI footprints to <300 plots in Paraguay that were <200 m apart. How far GEDI data to NFI plots can range to ensure improvements in biomass calibrations requires more analyses (Zhang et al., 2022). Following with Bullock et al. (2023), geolocation errors in the inventory plots were not acknowledged despite practical uncertainties related to use of non-survey grade positioning instruments under dense tree canopies. We know enhancing plot geolocations improves the relationships between laser metrics and area-based forest attributes (Fassnacht et al., 2014; Lu et al., 2016; Pascual et al., 2018; Liu et al., 2023). Large offsets between nominal and real plot positions - and edge effects from forest fragmentation - can reduce the availability of suitable NFI data and compromise efforts in biomass calibration especially over highly-fragmented and diverse forests such as in the Mediterranean (Pascual et al., 2020). Hence, using as many scarce field plots as possible in calibrations is fundamental. The waveform imputation approach presented in May et al. (2024) ensures all NFI plots are used to calibrate biomass, but at the cost of some level of errors in the imputed GEDI canopy height metrics for user-specified locations. The approach uses Bayesian inference, with traceable error values that can be subsequently incorporated into new biomass calibration models. Trade-offs between data characteristics (simulated versus on-orbit data), scale mismatches between NFI plots and sources of GEDI canopy height, or the quality of new estimates must be thoughtfully evaluated before generalizing an approach to produce locally-tailored NFI biomass estimates from spaceborne GEDI metrics. Using a province in the North-West of Spain with contrasting forest management regimes as testing area, we implemented the three following calibration strategies to calibrate GEDI biomass estimates using NFI data. i) Use of simulated GEDI canopy height from ALS data coregistering the NFI locations - following GEDI L4A model building; ii) Create pairs of nearest neighbors between plots and on-orbit GEDI data, and iii) Implement the waveform imputation method from May et al. (2024). Standard and enhanced-geolocated positions for >1000 plots were used to disentangle the performance of these methods towards frequently-ignored but common geolocation errors present in NFI data. Our benchmarking study is a comparison of alternatives to calibrate global GEDI biomass estimates to local conditions - more alike to existing field data and more aligned to forest management and operational monitoring needs. This a relevant research topic with multiple applications such as helping NFI programs to fill data gaps, improving the estimation of forest carbon to be optimize through management, empowering managers with methods to correct biases and uncertainties in GEDI estimates, but also towards the assimilation of GEDI data into MRV projects that requires locally-calibrated estimates to meet verification standards. 2. Material and methods 2.1. Study area The study area is Leon province in the North-West of Spain. Located in the Southern border of the Cantabric mountains, Leon is the seventh largest province in Spain and home of emblematic forest reserves such as The Picos de Europa National Park (Fig. 1). The 586,000 ha of forest area are dominated by hardwood species (71%) followed by coniferous and mixed forests (26% and 5, respectively). Some of the dominant species are Quercus pyrenaica (dominant in 36.2% of the forest area), Pinus sylvestris (10.9%) or sparse Quercus ilex (7.3%). The province is foreststructurally rich and species-diverse: the Spanish NFI (SNFI) quantifies 39 forest types, from low-land sparse oak forests and restoration forests, to dense pine stands and shade-tolerant species such as old-growth beech and chestnut forests in mountainous areas. The province is a suitable to disentangle the performance of GEDI biomass estimates considering climate zones, transition areas, presence of mountain oldgrowth forests and large extensions of conifer forests from past reforestation efforts. As of today, only one model in the GEDI L4A product is scoped for all these forest types (Kellner et al., 2022). Recent ALS surveys and coincident field measurements of known and enhanced geolocation (Pascual et al., 2021; Pascual and Guerra-Hern´ andez, 2023; Hunka et al., 2023) create an optimal scenario to improve forest management information through biomass calibrations using GEDI canopy height data and high-quality NFI data. 2.2. National forest inventory data The SNFI is designed as concentric plots ranging up to 25-m radius following a 1-km grid for the systematic sampling design A. Pascual et al. Journal of Environmental Management 375 (2025) 124313 2 Fig. 1. The province of Leon in NW Spain is the study region for the research. The NFI sampling grid where plot geolocation has been improved is show in the map over the Spanish Forest Map used to retrieve the main tree species. Scenes of well-represented forest strata in the subset of NFI plots are shown. Fig. 2. Distribution of geolocation offset correction in SNFI-4 plots used for the study. We showed three examples in different forest types using a lidar-based canopy height map in the background. A. Pascual et al. Journal of Environmental Management 375 (2025) 124313 3 (´ Alvarez-Gonz´ alez et al., 2014, (Guerra-Hern´ andez et al., 2021). Species or tree structural measurements are part of the treeand stand-level array of variables monitored documented including biomass, which is estimated using SNFI species-specific allometries (Pascual and Guerra-Hern´ andez, 2023). We used the latest survey available for Leon province collected between August–November 2019. There are 1360 plots available but for this research we used 1161 plots for which upgraded geolocation information was provided by the SNFI headquarters. Correcting the geolocation offset took an average of 2 h 21’ for each plot. This allowed us to have a pair of coordinates, enhanced-geolocated and NFI standard NFI positioning, which is more inaccurate (Pascual et al., 2020). The average geolocation offset in the NFI samples was 22.7 m and it was below 50 m for 90% of the plots, some of which exceeded 100 m (Fig. 2). 2.3. Airborne lidar to estimate biomass and simulate GEDI Airborne lidar data is publicly available in Spain thanks to the Land Survey Institute (Plan Nacional de Orthophotography A´ erea, PNOA). Point cloud density increased from 1 to 2 points m −2 in early PNOA surveys to ~5 points m −2 in recent projects, improving quality in 3D structural mapping and the suitability of these surveys to support operational forest management. The schedules of the SNFI and PNOA airborne lidar surveys usually are not synchronized, but Leon province is an exception. The ALS mapping using the RIEGL LMS-Q1560 sensor was conducted in 2021 from June 21st to October 21st in leaf-on conditions. Percentile height distributions of airborne lidar echoes were computed using a 2.5-m height-break to filter out non-forest vegetation as shown in the 10-m resolution canopy height map produced (Fig. S1). Metrics from ALS point clouds co-registering well-geolocated NFI plots were used to calibrate and predict biomass. The ALS-based biomass map was used to assess the effect of new calibrations compared to existing onorbit estimates. Details on the available 25-m biomass map can be found in Pascual and Guerra-Hernandez (2023) and in the Supplementary (Fig. S2). The ALS point cloud data was used to simulate GEDI canopy height metrics at NFI locations using the GEDI simulator (Hancock et al., 2019). The ALS data pass the requirements on seasonality (leaf-on) and minimum point cloud density (>4 points m −2 ) to ensure quality in the simulated GEDI waveforms (Hancock et al., 2019). The simulator operates by generating waveforms for given locations from discrete ALS data and it was used to build GEDI biomass models (Duncanson et al., 2022; Kellner et al., 2022). We simulated GEDI metrics twice on each plot: for enhanced-geolocated locations and for uncorrected positions available from previous inventories. We added waveform noise in the GEDI simulations to mimic on-orbit conditions. 2.4. Forest structural metrics and biomass estimates from GEDI In this section, we describe how we match NFI and GEDI footprint data: we used 3 calendar years of GEDI data from 2019 to 2021 on canopy height (L2A product, Dubayah et al., 2020) and AGBD estimates (L4A product, Dubayah et al., 2022b; Kellner et al., 2022). The GEDI data quality filtering was as follows (more details in Table S1): waveform fidelity in the selected footprints exceeded 0.95 and it was greater than canopy cover and the algorithm selection in GEDI L2A was optimized for each footprint. Quality flags in L2A and L4A were considered as well as maximum allowed difference of 50 m between GEDI elevation lowest mode and TanDEM-X terrain elevation (Krieger et al., 2007) available in the GEDI L2A. Biomass in GEDI is only predicted for the leaf-on season (L4A quality flag =1) and this has relevant implications for broadleaf forest types. A total of ~236, ~203 and ~239 thousand footprints passed the filters (Fig. 3). For each pair of NFI locations, we searched in 25-m intervals for GEDI footprints ranging 100 m or less from center locations. Waveform acquisition time provided in the GEDI L2A product was used to control the pairing of GEDI footprints to NFI Fig. 3. Impact of GEDI data sparsity in the paring of NFI plots to the most adjacent GEDI footprint. Data sparsity is a function on the density of GEDI tracks, the proportion of cross-tracks and the effect of quality filters at reducing the number of high-quality GEDI footprints along each ground track. A. Pascual et al. Journal of Environmental Management 375 (2025) 124313 4 plots. Unmatched NFI plots in year 2019 were subsequently explored in 2020 and 2021 until a match was found, or otherwise the candidate plot was removed from the training data. In the 0–25 m distance range between plots and footprints, we computed 134 and 148 plot-waveform pairs for enhanced and standard geolocation, respectively. Below the 100-m threshold, we listed 538 and 535 data pairs, which shrinks the training data actionable for re-calibrations by ~50% when the option is to use nearest-neighbor as in Bullock et al. (2023). The novel approach on waveform imputation presented in May et al. (2024), however, maximizes all reference data available for GEDI biomass estimation. 2.5. Waveform imputation The waveform imputation approach is based on interpolation of GEDI spatially-incomplete observations around NFI locations and it uses forest masks to account for edges in the process (May et al., 2024). The method uses a multivariate spatial model (Taylor-Rodriguez et al., 2019; May et al., 2022) to acknowledge the spatial correlation of GEDI data before generating probabilistic predictions of GEDI-like canopy height metrics at the specified NFI geolocations. The predicted attribute of interest (canopy height – GEDI L2A relative height 98, RH 98 ) is given by the expected value of the predictive distribution and the prediction uncertainty is given by the dispersion of the distribution around the expected value. The density of GEDI footprints around plot locations and the distance search to pair RH metrics to plot attributes are important factors in the imputation method (May et al., 2024). The imputed outcomes are predicted potential values of RH metrics for a given NFI plot at the specified location and plot size. Prediction errors are incorporated in imputed RH 98 and subsequently into biomass calibrations that can be customized to jurisdictions or forest types by grouping NFI data. We imputed RH 98 for each pair of NFI geolocations accuracies – enhanced and standard. 2.6. Local calibrations of GEDI biomass data The impact of data sparsity is more evident when cross-intersection of GEDI tracks is minimal or none-existing, leading to large data gaps around sampling plots especially for broadleaf forest types affected by biomass quality filters (Table S3). By using only leaf-on GEDI data, we narrowed the subset of GEDI measurements to the interval between March 16th to November 3rd corresponding to the 2000–2021 median leaf-on season period in Leon province. Mountain forests in northern latitudes such as in Picos de Europa NP have a shorten growing season: 90 days less approx. (Caparros-Santiago and Rodriguez-Galiano, 2024). Deciduous broadleaf trees (DBT), evergreen broadleaf (EBT), evergreen needleleaf (ENT) and deciduous needleaf (DNT) are the four PFT that the GEDI mission uses for biomass mapping. The strata were grouped for Europe into one due to the few GEDI FSBD sites used to build biomass models. One model is used for all biomass predictions across European forest while grassland, shrubs & woodlands (GSW) have their global model (Table 1). Acquisition date and PFTs are used to action a leaf-on quality flag that excludes leaf-off conditions in EBT or DNT forests and woodlands (71% of Leon province forest area is dominated by broadleaved tree species). We used the following three sources of GEDI canopy height as predictor of biomass for two NFI plot geolocation series. 1) ALS-simulated GEDI-alike estimates for canopy height (RH 98 ), 2) waveform imputation with associated uncertainties on imputed RH 98 values following the method presented in May et al. (2024), 3) pairs of NFI plots and adjacent GEDI measurements for canopy height. We used GEDI data acquired during the first three years of mission (2019/21). These pairs met the following criteria: labeled as high quality according to the GEDI L2A and L4A quality flag (a combination of multiple waveform quality metrics), less than 100 m from the NFI plot center position - half the search window used in i. e., Bullock et al. (2023) - and undisturbed between the time of NFI measurements and the GEDI acquisition. Plots harvested or partially-harvested as interpreted from changes in basal area and dominant height from comparing consecutive measurements (SNFI-3 versus SNFI-4) were removed from the training data. The associated uncertainty (i.e., variance of the covariates) for predictor RH 98 was assumed zero when using ALS-simulated (1) and onorbit GEDI estimates (3). For imputed RH 98 , we incorporated the predicted variance in the calibration. For these three cases, we trained a simple regression model (Eq. (1)) between squared-root (sqrt) transformed GEDI RH 98 and reference NFI biomass - following GEDI L4A model building approach the approach but using only one predictor (Table 2; Kellner et al., 2022). The sqrt-sqrt transformation helps induce a linear relationship between the covariate and the response.  AGBD √=βo+β1 RH98 +100 √+ ε [Eq. 1] where AGBD is the reference estimate in Mg ha −1 for a given NFI plot and RH 98 is the GEDI RH 98 imputed, simulated or paired on-orbit value to the specific plot. We trained a re-calibration model for each method using SNFI-4 data for enhanced and for standard geolocations. Forest type-specific estimates for four region-dominant tree species well NFIrepresented were assessed in detailed. The species included were Quercus pyrenaica (QP, 323 plots), actively managed Pinus sylvestris (PS, 174), Quercus ilex (QI, 68), and Fagus sylvatica that dominates northern valley (FG, 87). To assess model performance, we used the absolute and relative root mean squared error (RMSE) between fitted values and the NFI biomass data in the training set. 2.7. Gridded biomass products The new NFI-calibrated GEDI estimates were applied to on-orbit Table 1 The GEDI L4A biomass models for Europe. Strata Model equation DBT EU AGBD =0.963 × (−96.531 +7.175 × RH50 +100 √+2.921  RH98 +100 √)2 ENT EU EBT EU DNT EU GSW AGBD =1.118 ×(−124.832 +12.426 × RH98 +100 √)2 Note: EBT, Evergreen Broadleaf Trees; GSW, Grasslands & Shrubs and Woodlands; EU, Europe: AGBD, aboveground biomass density (Mg ha −1 ); RH 50 and RH 98 , relative height 50 and 98, respectively, available in the GEDI L2A product. Table 2 Performance of aboveground biomass calibration models (AGBD, Mg ha −1 ) using on-orbit GEDI, waveform imputation and simulation of GEDI metrics using airborne lidar. The standard deviation (SD) of the model error, the root mean square error (RMSE) and model bias are presented for two classes of geolocation accuracies for the NFI plots used for model building. Geolocation accuracy Canopy height (RH 98 , m) Model β 0 /β 1 RMSE (Mg ha ¡1 ) RMSE (%) Bias (Mg ha ¡1 ) Enhanced Waveform Imp. −101.65/ 1.04 73.05 79.91 −2.46 ALS simulation −111.60/ 1.13 64.38 70.42 0.00 GEDI L2A −53.26/ 0.58 73.72 82.19 0.00 Standard Waveform Imp. −100.11/ 1.02 73.00 80.18 −2.68 ALS simulation −15.32/ 0.18 64.46 70.80 0.00 GEDI L2A −45.67/ 0.51 75.47 84.52 0.00 A. Pascual et al. Journal of Environmental Management 375 (2025) 124313 5 footprint passing quality filters for both L2 and L4 products, and collected during GEDI mission’s first epoch, from April 2019 to March 2023. We used hexagonal grids to summarize the three NFI-calibrated GEDI biomass distributions and to create a baseline by gridding onorbit GEDI data. Specifically, we used the H3-scale of 10 (1.5 ha in grid resolution, 75.9 m hexagon edge) to size 8 (73.7 ha and 531.3 m edge, respectively). Biomass was predicted within two forest extent masks as forests represent only 37% out of 15,581 km 2 . We used the strata used by GEDI mission (PFT) for calibration and prediction of biomass and class ‘forest’ in the 10-m ESA land cover product (Zanaga et al., 2022). We did not consider GEDI footprints ranging within the GSW class or within land cover classes in the ESA product other forest areas. To assess the benefits of calibrations towards forest management planning in the region we used inter-comparisons of quantile distributions and maps between the three tested calibration alternatives and existing baselines based on GEDI and ALS-based inference. 3. Results 3.1. New calibrations of aboveground biomass The ability of waveform imputation to detect changes in the geolocation of NFI plots had a minimal influence in the distribution of canopy height estimates when accounting for the correction in the known geolocation offset (Fig. 4). Differences in RH 98 canopy height distributions were almost identical: as low as 7 cm while for ALSsimulated GEDI canopy height the value reached 1.40 m for the entire dataset (1152 plots). The bias in the predicted biomass is higher than for ALS-simulated metrics and for on-orbit GEDI canopy height but it is a low value (<5 Mg ha −1 ). Biomass calibrations using ALS-simulated GEDI Fig. 4. Distribution of relative squared mean errors (RMSE) and model bias in the estimation of aboveground biomass density (AGBD, Mg ha −1 ) considering geolocation accuracy in NFI plot locations. Comparisons using SNFI-4 data as reference are presented for on-orbit GEDI, waveform imputation and simulation of GEDI metrics using airborne lidar (above). Geolocation offset correction in NFI plots over the gradient of canopy height data. The y-axis shows the difference between relative height 98 (RH 98 , m) computed from lidar-simulated waveforms and from the waveform imputation approach (below). A. Pascual et al. Journal of Environmental Management 375 (2025) 124313 6 canopy height showed the lowest model error and the most similar match to the distribution of reference biomass (Table 2). The largest RMSE was observed for the pairing of GEDI footprints to NFI plots; 10% RMSE compared to using ALS-simulated canopy height as model predictor. 3.2. Forest type effects and window search threshold Biomass was modelled for four relevant forest types with respect of management regimes (>650 plots, Table 3) and of different stocking patterns: mean values ranged from ~30 Mg ha −1 (sparse Quercus oaks) to ~ 200 Mg ha −1 (dense beech stands in mountain forests). Errors in biomass models for sparse oaks (~90% RMSE) tripled the values for Mediterranean pine forests regardless of the calibration method applied or the geolocation accuracy in the training data. These results for individual forest types showed that simulated GEDI canopy height was not always the most efficient solution in calibrations. Imputation improved the scores in ~5% from ALS-simulation for pine forests (P. sylvestris) while the improvement reached 15% compared to data pairs between adjacent GEDI data and NFI plots. Estimates based on adjacent on-orbit GEDI canopy height data systematically showed the worst model performance and higher values were observed when increasing the search window to include more GEDI footprints as candidate for the pairing to NFI plots (Fig. S5). In this approach, the proportion of on-orbit GEDI footprints passing L4A filters was lower for broadleaved forests compared to evergreen forests, which impacts the results for imputation and the plot-to-footprint nearest neighbor approach. For instance, the RMSE was 30% higher compared to imputation and ALS-simulation in Quercus pyrenaica forests, the most dominant forest type. 3.3. Distributions of GEDI footprint and gridded estimates Calibrated biomass estimates remained consistent regardless of the forest mask used to apply the new models to GEDI on-orbit data (Fig. S7). Differences mostly occurred in the 10–40 Mg ha −1 interval and all distributions similarly captured nearly zero-biomass conditions (Fig. S6). The bi-modal shape predictions showed lower maxima for ALS-simulation compared to on-orbit estimates both operational values (GEDI L4A) and the new calibration based on pairing on-orbit GEDI to NFI samples. For waveform imputation, predictions were similar to ALSsimulation when using PFTs as prediction strata and to on-orbit data when using the ESA forest map. The overall improvement from the NFIbased calibration was low in the nearest neighbor approach (Fig. 5). Our results correspond to size 10 of the H3-hexagonal grid, which is ~0.75 the scale of the GEDI L4B product (Fig. S3) and comparable to the grids we generated. The new calibrations based on imputation and ALSsimulated canopy height estimates were more similar to biomass estimates from ALS-based inference. There is a resolution mismatch between the later (25 m) and new GEDI biomass grids, but its use for validation helps to visualize the correction around the 90th percentile of the distribution (90–120 Mg ha −1 ). The selected masks impact the layout of biomass maps, especially for non-forest and low-biomass conditions in the central sections of Leon province (Fig. S9). Differences between gridded L4A and new calibrations were higher over mountainous forest in the Northern range of the province. The average change in mean biomass for the 20,610 hexagons (Fig. S10) showed high mean and median values for GEDI L4A gridding compared to the three NFI-calibrated biomass distributions i.e., 6.85 and 2.44 Mg ha −1 , respectively, for waveform imputation, or 6.72 and 4.12 Mg ha −1 , respectively, when using on-orbit RH 98 as predictor. Waveform imputation showed lower estimates (Fig. 6, red scale) than GEDI L4A predictions, while the use of ALS simulation and on-orbit GEDI produced higher values than the gridded estimates from GEDI L4A (Fig. 6, blue scale). Opposite tendencies are especially notorious in clusters over mountainous forests and steep terrain such as in the Leonside of ‘Picos de Europa’ National Park where biomass estimates from ALS-based inference successfully excluded non-forest patches in prediction strata. Biomass distributions in these areas looked more aligned between on-orbit GEDI and imputation methods. The discrepancies towards the ALS-based simulation approach showcased the need to carefully evaluate predictions in challenging environments and also the different properties between simulated and on-orbit data. Biomass recalibrations in these complex environments with respect to topography and phenology are especially relevant for the conservation and sustainable management of these old-growth high-biodiversity forests. 4. Discussion 4.1. Biomass estimates at plot-level Our study showcased an integration of global spaceborne data from GEDI into regional biomass calibration to produce locally-tailored estimates at jurisdictional level. Among the tested methods, ALS-simulation produced the overall best model fitting results, but the approach is not operational globally as it requires discrete lidar and assumes that GEDI waveform simulation errors are negligible compared to on-orbit GEDI properties. Hence, ALS simulation should be regarded as a theoretical upper bound to benchmark operational and scalable options: waveform imputation and data pairs between GEDI footprints and inventory plots based on geolocation. Simulated GEDI canopy height from airborne lidar yielded the lowest error when using all NFI plots available for the 39 forest types included in the training database, but fitting by homogenous groups of NFI plots (forest types) resulted in better scores for waveform imputation on pine strata – actively managed for timber production - and comparable results for other hardwood specie. Our calibration methodologies have fundamental principles: both ALS-simulation and nearest neighbor approaches use a single 25-m footprint waveform as predictor of canopy height, while imputation uses neighboring footprints to infer canopy height metrics. Waveform imputation and the NFI-to-GEDI neighbor approach rely on the distribution and density of GEDI metrics around specific targets, and this can be more informative of biomass distribution compared to relying on a single well-geolocated small-size footprint to summarize the structural conditions of a considerably larger plot. The mismatch in scale between response and explanatory variables exist: we calibrated 50-m footprint biomass estimates using 25-m Table 3 Performance of aboveground biomass calibration models (AGBD, Mg ha −1 ) by means of the percentual root mean squared error (RMSE) using on-orbit GEDI, waveform imputation and simulation of GEDI metrics using airborne lidar in five dominant forest types. The standard deviation of the model error, the root mean square error (RMSE) and model bias are presented for two classes of geolocation accuracies for the NFI plots used for model building. Forest Type RMSE (%) Mean AGBD (Mg ha ¡1 ) NFI plots Ratio GEDI shots/Forest area Imputed RH 98 Simulated RH 98 On-orbit RH 98 Enhanced Standard Enhanced Standard Q. pyrenaica 53.63 56.78 67.04 89.22 89.00 64.41 323 1.08 P. sylvestris 28.84 33.72 38.17 50.91 51.59 75.32 174 1.22 Q. ilex 89.76 100.16 99.45 107.55 110.32 29.93 68 1.32 F. sylvatica 44.03 39.76 41.65 47.94 46.20 190.1 87 0.67 A. Pascual et al. Journal of Environmental Management 375 (2025) 124313 7 footprint GEDI data. Despite scaling inconsistencies, the calibration of biomass by forest type was effective to lower estimation errors as shown in previous studies (Bruening et al., 2023; Crockett et al., 2023). We found that, regardless of the approach, biomass calibration was far from optimal for low-stature sparse forests (e.g., Quercus ilex), whose horizontal structural conditions present a challenge to the GEDI signal processing especially when canopy cover is low and much of the ground within the 25-m footprint is exposed (Dubayah et al., 2020; Dorado-Roda et al., 2021; Li et al., 2023). Considering sparse forests are not well suited for the area-based estimation of stocking (Pascual et al., 2023; Guerra-Hern´ andez et al., 2024), imputation equaled the performance of ALS-simulation according the low bias observed. The sparse oak forests we evaluated are evergreen and had higher coverage of GEDI data compared to other broadleaved forests assessed that were more affected by phenology conditions and filtering criteria. The phenology flag used in GEDI products lowered the proportion of footprints actionable for biomass calibrations in broadleaved forests when using on-orbit data. Results for waveform imputation might have been even better if we added more observations: data flagged using canopy height (L2A product) instead of the more-restrictive biomass flag (L4A product). Higher GEDI coverage for evergreen forests compared to broadleaved trees restricted the full capacities of waveform imputation as the density and the distribution of GEDI footprints play a role (May et al., 2024). The GEDI mission uses both leaf-off and leaf-on data for canopy height products but for consistency with GEDI L4A models, we applied the phenology flag used for GEDI L4A product. This can be revisited for i.e., small jurisdictions where data restriction can limit improvements in biomass estimates, or for specific strata that management identifies as priorities. We dropped >50% of on-orbit canopy height estimates for the broadleaved forest types assessed, but at more local scales in northern latitudes and higher elevations, e.g., Picos de Europa NP, the value is higher as the growing seasons shortens compared to lower elevations where most of the active forest management (silviculture) in the region takes place, mainly over pine plantations in mature stages and evolved to irregular stand structures. We calibrated biomass at jurisdictional level for each method, but we could have stacked predictions for individual forest types prioritized by management needs and economy in the region (i.e., incorporating pine-dominated areas or mixed forests) or by ecoregions as shown by Bullock et al. (2023) in Paraguay or by Bruening et al. (2023) in the United States. Indeed, we could have aggregated forest types using fewer classes at higher NFI hierarchies (i.e., genus) to ensure enough training data exist in Leon province i.e., by genus or by age class as in Jha et al. (2024). However, most of the represented forest types had <30 samples of enhanced geolocation, and that is the reason we present biomass map at jurisdictional-level for the three alternatives. Further research will expand this biomass calibration effort nationwide and for the entire Iberian Peninsula over the next years making use of incoming ALS acquisition in Spain and in Portugal – the first country-wide survey in the country. 4.2. Biomass predictions and gridding Our footprint-level biomass calibrations were applied to forest-only areas using a land cover product that is accurate for forests and bareearth topography, but less accurate for shrublands and grasslands (Zanaga et al., 2022). Mixed conditions between forests, shrublands and grasslands over mountainous areas and in the proximity of ridges can produce unrealistic high estimates for GEDI canopy height (Liu et al., 2021), but there are more factors to assess the presence of outliers in biomass maps (Knott et al., 2023). To lower the impact of those GEDI outliers over sensitive land cover types, the GEDI L4B product assigns zero mean and uncertainty to these footprints ranging over bare-earth terrain using the same product we used (Dubayah et al., 2022a; Armston et al., 2023). Managers interested in turning actionable biomass estimates and recalibrated values from GEDI should carefully evaluate the distribution of the data as these environments present challenging conditions for canopy height mapping (Planet Labs, 2024). At the prediction stage, we found biomass calibrations based on ALSsimulation produced higher values, on average, compared to GEDI L4A gridded estimates (Fig. S11), which was the opposite as for imputation and on-orbit data. Properties between real and simulated GEDI data are different and these can trigger unexpected differences in the distribution and values of biomass estimates in final maps as observed using quantile comparisons (Fig. S8): differences not captured by the fitted biomassheight calibrations at plot/footprint level. The quantile distribution of biomass for waveform imputation and for ALS-simulation remained similar below the 90th quantile. However, at more local scales and in presence of mixed land cover types that lower the proportion of forests within a given jurisdiction such as within the boundaries of Picos de Europa NP, biomass estimates can disagree more as the new calibration models were forest-specific. For instance, a substantial proportion of change estimates were above 100 Mg ha −1 in Picos de Europa NP but of different sign: models from ALS-simulation above-predicted while the imputation approach produced lower average estimates compared to gridded GEDI L4A baseline. 4.3. Support of NFI programs and conservation finance The sparsity of GEDI tracks around inventory plots affects both the NFI-to-GEDI pairing approach and the waveform imputation. The latter ensures all inventory data are used in biomass calibrations at the cost of Fig. 5. Distribution quantiles of the three biomass calibrations strategies compared to a 25-m biomass map built with airborne lidar (ALS) and NFI data. A. Pascual et al. Journal of Environmental Management 375 (2025) 124313 8 some prediction error in the canopy height predictor (May et al., 2024). Using nearest neighbor with a 100-m maximum search for data pairing heavily reduced the training data available (>50%). Waveform imputation method lacks sensitivity towards geolocation errors in the reference plots and this has more pros than cons from the optics of biomass calibrations and MRV activities. Geolocation errors are considerable in scale and distribution in many NFI programs and environmental monitoring networks (Pascual et al., 2018; Andersen et al., 2022). The scale of geolocation offsets in many forest sampling programs is likely larger than the 23-m offset average in Leon province, especially where NFI programs are in early stages and/or where there is limited access to high-end field positioning equipment (Silveira et al., 2023). Fig. 6. Change map between three methods to generate newly-calibrate estimates of aboveground biomass density (AGBD, Mg ha −1 ) gridded at ~ 500-m grid compared and to subtracted from GEDI L4A gridding at the same scale. The prediction stratum to apply the new biomass calibration equations was class ‘Forests’ in the ESA V200 landcover product. A. Pascual et al. Journal of Environmental Management 375 (2025) 124313 9