Long-term annual estimation of forest above ground biomass, canopy cover, and height from airborne and spaceborne sensors synergies in the Iberian Peninsula

Environmental Research 259 (2024) 119432 Available online 27 June 2024 0013-9351/© 2024 The Authors. Published by Elsevier Inc. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/bync-nd/4.0/). Long-term annual estimation of forest above ground biomass, canopy cover, and height from airborne and spaceborne sensors synergies in the Iberian Peninsula M.A. Tanase a , * , M.C. Mihai a , S. Miguel a , A. Cantero b , J. Tijerin c , P. Ruiz-Benito c , D. Domingo d , g , A. Garcia-Martin e , g , C. Aponte f , M.T. Lamelas e , g a Universidad de Alcal´ a, Environmental Remote Sensing Research Group, Departamento de Geología, Geografía y Medio Ambiente, Colegios 2, 28801, Alcal´ a de Henares, Spain b HAZI Fundazioa, Vitoria-Gasteiz, Spain c Universidad de Alcal´ a, Grupo de Ecología y Restauraci´ on Forestal, Departamento de Ciencias de la Vida, Facultad de Ciencias, 28805, Alcal´ a de Henares, Spain d iuFOR, EiFAB, Universidad de Valladolid, 42004 Soria, Spain e Centro Universitario de la Defensa de Zaragoza, Academia General Militar, Ctra. de Huesca s/n, Zaragoza 50090, Spain f Instituto de Ciencias Forestales ICIFOR-INIA, CSIC, Madrid, Spain g GEOFOREST-IUCA, Departamento de Geografía, Universidad de Zaragoza, Pedro Cerbuna 12, 50009 Zaragoza, Spain ARTICLE INFO Keywords: Forest attributes Temporal trends Spanish forests Chronosequence ABSTRACT The Mediterranean Basin has experienced substantial land use changes as traditional agriculture decreased and population migrated from rural to urban areas, which have resulted in a large forest cover increase. The combination of Landsat time series, providing spectral information, with lidar, offering three-dimensional insights, has emerged as a viable option for the large-scale cartography of forest structural attributes across large time spans. Here we develop and test a comprehensive framework to map forest above ground biomass, canopy cover and forest height in two regions spanning the most representative biomes in the peninsular Spain, Mediterranean (Madrid region) and temperate (Basque Country). As reference, we used lidar-based direct estimates of stand height and forest canopy cover. The reference biomass and volume were predicted from lidar metrics. Landsat time series predictors included annual temporal profiles of band reflectance and vegetation indices for the 1985–2023 period. Additional predictor variables including synthetic aperture radar, disturbance history, topography and forest type were also evaluated to optimize forest structural attributes retrieval. The estimates were independently validated at two temporal scales, i) the year of model calibration and ii) the year of the second lidar survey. The final models used as predictor variables only Landsat based metrics and topographic information, as the available SAR time-series were relatively short (1991–2011) and disturbance information did not decrease the estimation error. Model accuracies were higher in the Mediterranean forests when compared to the temperate forests (R 2 =0.6–0.8 vs. 0.4–0.5). Between the first (1985–1989) and the last (2020–2023) decades of the monitoring period the average forest cover increased from 21 ±2% to 32 ±1%, mean height increased from 6.6 ±0.43 m to 7.9 ±0.18 m and the mean biomass from 31.9 ±3.6 t ha −1 to 50.4 ±1 t ha −1 for the Mediterranean forests. In temperate forests, the average canopy cover increased from 55 ±4% to 59 ±3%, mean height increased from 15.8 ±0.77 m to 17.3 ±0.21m, while the growing stock volume increased from 137.8 ±8.2 to 151.5 ±3.8 m 3 ha −1 . Our results suggest that multispectral data can be successfully linked with lidar to provide continuous information on forest height, cover, and biomass trends. 1. Introduction Forests, the most biodiverse terrestrial ecosystems, are key for observing and monitoring Earth condition as they can either tie up atmospheric carbon in long-lasting wood and soil stores (carbon sink) or contribute as carbon source (e.g., biomass loss due to fire, logging, land use change). Forests are also a major component of rural development, provide species habitat and protective functions, and contribute goods * Corresponding author. E-mail address: [email protected] (M.A. Tanase). Contents lists available at ScienceDirect Environmental Research journal homepage: www.elsevier.com/locate/envres https://doi.org/10.1016/j.envres.2024.119432 Received 2 May 2024; Received in revised form 4 June 2024; Accepted 15 June 2024 Environmental Research 259 (2024) 119432 2 and services (Ojea et al., 2010). During the past century, environmental and anthropogenic factors have caused serious threats to forest ecosystems integrity leading to habitat degradation and the related loss of services (Michel and Seidling, 2014; Senf et al., 2017). Many forest ecosystems have experienced an increase in disturbance rate, magnitude, and frequency, with recent trends reaching unprecedented levels (Millar and Stephenson, 2015). Globally, approximately 67 million ha of forests are lost through fires each year while further 10 million hectares are affected by biotic agents including insects and diseases (Lierop et al., 2015). In Europe alone, 17% of forests were affected by some type of disturbance during the past 20 years (Senf and Seidl, 2020). In 2022, the last year with full records, the European Union experienced the second-worst year for wildfires, with about 0.84M ha being scorched out of which 38% (0.315M ha) were recorded in Spain alone. Indeed, the 2022 Spanish fire season was characterized by a threefold increase in burned area with respect to decadal averages (San-Miguel-Ayanz et al., 2023). Further projected changes in global climate due to increasing greenhouse gas concentrations will have severe effects on terrestrial biomes, particularly forests (Anderegg et al., 2020; Seidl et al., 2017). The rise of temperature and alteration of precipitations patterns are accelerating affecting both local and broad-scale ecosystem processes, including disturbance regimes (Senf and Seidl, 2020). Indeed, temporal trends observed in forest demography suggest a worldwide increase of climate-induced tree mortality (Brienen et al., 2015; Senf et al., 2018) with similar trends being observed for the Spanish forests (Astigarraga et al., 2020). According to future climate projections disturbances will further intensify (Seidl et al., 2014) with potentially major impacts on global carbon sequestration (Kurz et al., 2008) as forest conversion contributes to a large share of greenhouse gas emissions. Therefore, understanding forest dynamics by monitoring forest spatial extent and structural attributes (i.e., forest height, forest canopy cover, above ground biomass) is fundamental to understanding their current state, predict future trajectories, and develop adequate management plans and policies. Reliable and consistent monitoring of forest height, canopy cover, and biomass are essential for estimating forest-related carbon emissions, analyzing forest degradation, and evaluating the effectiveness of forest management strategies in response to current environmental challenges (Lang et al., 2023). Forest canopy height correlates with many biophysical parameters including basal area, above ground biomass, net primary productivity allowing understanding the structure and condition of forest ecosystems (Bolton et al., 2017; Goetz et al., 2007; Lang et al., 2019; Potapov et al., 2021). In addition, canopy height directly characterizes habitat heterogeneity, being ranked as a high-priority biodiversity variable (Skidmore et al., 2021). Similarly, information on forest canopy cover enables the estimation of forest area and thus carbon stocks monitoring and associated greenhouse gas emissions and removals (Estoque et al., 2021). Further, monitoring changes in canopy cover is essential for various applications in ecology, hydrology, and forest management, including establishing restoration priorities, understanding wildlife species habitat, or estimating fire impacts (Gast´ on et al., 2019; Yin et al., 2020). Land managers invest large amounts of effort and resources to estimate forest structural attributes for multiple management objectives (Hudak et al., 2020). Remote sensing technologies are the pillar of many successful programs, such as the National Aeronautics and Space Administration (NASA) Carbon Monitoring System and the European Space Agency (ESA) Climate Change Initiative, focused on land use management, resource allocation, and environmental monitoring. In this context, the spatially and temporally sparse in-situ information, usually acquired within national forest inventories, is often used to estimates wall-to-wall forest attributes through modelling or imputation methods. Imputation methods, such as the k-Nearest Neighbour (kNN), fill in the missing data points using spatially explicit predictor variables derived from remote-sensing sensors based on the similarity between the values of the predictors associated with a target pixel and those associated with the k nearest training sample units (Nguyen et al., 2018). Their main advantage consists in assigning a set of measured attributes that actually occur (from the in-situ data) thus ensuring prediction of realistic values (Hudak et al., 2008). However, the use of imputation methods across large spatial and temporal spans may be problematic due to the oversimplification of forest structure distribution, challenges in validation, and limitations in capturing spatial and temporal variations (Johnson et al., 2020; Nguyen et al., 2018; Paik, 1997; Riley et al., 2016). Methods based on modelling and remote sensing data are increasingly used to estimate forest structural attributes (Hudak et al., 2020; Lang et al., 2019; Tanase et al. 2014, 2021) as they provide more accurate estimates by incorporating a wide range of parameters to develop complex relationships between the target variable and its predictors. Nevertheless, such methods may distort marginal distributions and covariation between the predictor variables which may result in unrealistic predictions (Matasci et al., 2018b). Forest structure and its changes may be identified and monitored by means of remote sensing using optical (Bolton et al., 2020; Hansen et al., 2013; Senf et al., 2018), radio detection and ranging or radar (Becker et al., 2023; Belenguer-Plomer et al., 2019; Tanase et al., 2015) light detection and ranging or lidar (Domingo et al., 2019; Goetz et al., 2010) sensors. Optical sensors are suitable for mapping land-cover since large changes in forest canopy cover can be indirectly observed. Lidar sensors provide the most accurate information to estimate vertical forest structure (Nelson et al., 2007; Sexton et al., 2009) while synthetic aperture radar (SAR) sensors are more sensitive to changes in forest structure when compared to optical data and less limited by spatial and temporal availability when compared to lidar. However, working with SAR data has its own challenges as the signal is affected by unrelated factors (e.g., precipitation, topography) which often hinder the retrieval of forest structural information (Tanase et al., 2019). Therefore, integrated frameworks that allows for harnessing the strengths of different sensor technologies to improving the information recovered from each sensor alone are increasingly developed (Belenguer-Plomer et al., 2021; Hudak et al., 2020; Joshi et al., 2016; Tanase et al., 2022; Zhang et al., 2019). As remote sensing technologies advance, new ways to collect, distribute, and analyze data become available. The digital revolution was accompanied by a ‘sensing’ revolution that provided unprecedented amounts of data on most Earth processes (Wilson, 2013). These new data streams offered advanced possibilities to further understanding an integrated earth system. In particular, advances in lidar technology coupled with the opening of the Landsat archive in 2008, allowed for the development of wide-ranging applications including monitoring disturbance, vegetation phenology, and ecosystem dynamics or the estimation of forest structural attributes (Bolton et al., 2020; Chen et al., 2021; Hudak et al., 2020; Pasquarella et al., 2017). An increasing number of studies leverage the accurate lidar estimated forest structural attributes to extend the information over broader areas and larger time spans (Coops et al., 2021). Such studies highlight the importance of Landsat time series in understanding forest ecosystem dynamics and changes across time and space. The integration of Landsat archives with newer lidar datasets has allowed obtaining backdated estimates of forest structural attributes (Caughlin et al., 2021; Hudak et al., 2020; Matasci et al., 2018b) facilitating synoptic views at landscape scales across the past decades (Matasci et al., 2018a). Using similar approaches, we aimed to evaluate historical trends (1985–2023) of the Spanish forest’s structural attributes (i.e., forest height, forest canopy cover, and above ground biomass) at regional level. The selected regions correspond to distinct Mediterranean (Madrid region) and temperate (Basque Country) biomes (Olson et al., 2001). The specific objectives of the study were to: i) leverage the systematic Spanish national lidar surveys to estimate reference forest structural attributes including above ground biomass, forest canopy cover, and forest height; M.A. Tanase et al. Environmental Research 259 (2024) 119432 3 ii) evaluate active/passive predictor variables synergies and transfer the models across time to produce annual, pixel-wise estimates for each forest attribute and region; and iii) analyze the temporal trend of each forest attribute at regional level. The evaluated predictor variables included optical reflectance, radar backscatter, and pixel-based disturbance history. To our knowledge this is the first study appraising synergies between spaceborne optic and radar data to predict chronosequence of forest attributes spanning several decades. 2. Material and methods 2.1. Study area The selected areas are in the peninsular Spain and comprise two autonomous communities (from herein referred to as regions), Madrid and the Basque Country covering 8028 km 2 and 7234 km 2 , respectively (Fig. 1). The natural vegetation (forests and shrublands) covers about 55% of the surface in the Madrid region and 56% in the Basque Country. The selection of administrative divisions is related to the self-governing status and the subsequent broad-scope competences which include agriculture and forestry. As forest related data are estimated by regional Agencies differences may arise with respect to the methods and temporal scales at which such data are estimated and made available. In the context of this study, such data include national and regional forest inventories and lidar surveys. The Madrid region, in the centre of the Iberian Peninsula, extends across elevations between 400 m a.s.l. in the south and 2400 m in the north. The relief changes from lowlands to mountains which results in variable ecological conditions. According to data from the Spanish meteorological agency (AEMET), the climate is Mediterranean with hot summers and an average annual rainfall of 450 mm, which occurs mainly in spring and autumn, with an important gradient from southern lowlands (<500 mm year −1 ) towards northern mountains (1500 mm year −1 ). The average monthly temperature ranges between 1 ◦C in winter and 32 ◦C in summer with a gradient from south (hotter) to north (colder). According to the forest species map (see section 2.8), the main tree species are oaks (i.e., Quercus ilex, and Q. pyrenaica) and pines (i.e., Pinus halepensis, P. sylvestris, P. pinea, and P. pinaster) which represent 26% and 11.5%, respectively, of the natural vegetation areas (Fig. 1a). An additional 12% (by area) of the vegetation consists of mixed forest species. Grasslands and riparian vegetation represent 16.5% and 2%, respectively, with the remaining areas (33%) being covered by shrublands (i.e., Cistus sp., Q. coccifera, Retama sphaerocarpa L., Thymus sp., Erica sp., and Macrochloa tenacissima L.). Forest height and biomass distribution is heavily skewed as small values dominate (Fig. 1b and 2). The Basque Country, on the Bay of Biscay, extends across largely mountainous terrain with elevations ranging from the sea level to over 2600 m a.s.l. It has a mild oceanic climate in the north that gradually changes toward a Mediterranean climate with semi-continental characteristics in the south. Average precipitations, typically exceeding 1000 mm year −1 in the north, drop to about 450 mm year −1 in the south. The average monthly temperature ranges between 5 ◦C in winter and 20 ◦C in summer (AEMET). The mountainous terrain modulates the temperature with hotter summers being common in the south. According to the forest species map, the main tree species are conifers (43%) with Pinus radiata representing about 26% of the forested area which stands at about 400 000 ha. Beech forests (Fagus sylvatica) occupy about 14% of the area while species from Quercus genus (Q. robur, Q. petrea, Q. pyrenaica, Q. fgaginea, and Q. ilex) occupy about 22% of the forested area (Fig. 1c). Eucalyptus plantations occupy about 6% while the mixed Atlantic forests occupy about 10% of the area. When compared to the Madrid region these forests are characterized by higher height and canopy cover (Fig. 1d and 2). 2.2. Modelling framework The study builds upon existing frameworks (Bolton et al., 2020; Fig. 1. Study area with the location of forest inventory plots (4th survey), the main tree species (a, c), and the normalized height estimated from the first lidar survey in Madrid (a, year 2010) and the Basque Country (d, year 2012). Main tree species according to the regional forest map (see section 2.8). M.A. Tanase et al. Environmental Research 259 (2024) 119432 4 Hudak et al., 2020) with the main changes related to the use of multi-temporal lidar surveys, standardized in-situ data from national forest inventories, a focus on largely homogeneous landscapes, pixel-wise disturbance history estimated with the Continuous Change Detection and Classification – Spectral Mixture Analysis (CCDC-SMA) model (Chen et al., 2021), and the use of annual statistics (i.e., percentiles) estimated from satellite time series as predictor variables (Fig. 3). Briefly, lidar data from the first Spanish national lidar survey (2008–2015) were used to estimate reference forest attributes directly (i.e., forest height and forest canopy cover) or indirectly, through non-parametric models trained with in-situ data (above ground biomass and growing stock volume). Additional non-parametric models (i.e. random forests) were trained to predict the lidar-based reference forest attributes using satellite imagery. A sensitivity analysis was carried out to select the optimum combination of predictor variables (optical, radar, topographic, disturbance history) in each region and for each forest attribute. The trained models were subsequently used to generate chronosequence (at annual step) for each forest attribute of interest. Temporal model transferability was independently tested against reference attributes estimated from the second Spanish national lidar survey (2015–2021). For the Basque Country additional independent datasets were used for validation including partial forest inventories and forest height estimated from stereo-orthophotography. A detailed description of each step is provided in the following sections. 2.3. Field data We used field data collected during the fourth Spanish national forest inventory (NFI) in 2011 in the Basque Country (n =1617) and 2012–2013 in Madrid (n =1055). The NFI sampling grid follows a regular pattern with one sample plot every square kilometre within forested areas (Sanz and Soto, 1990). All trees with a diameter at breast height (DBH) greater than 7.5 cm and height greater than 130 cm are recorded for their species and tree size (DBH and height). Tree height is measured indirectly with hypsometer or laser range finders while forest canopy cover is estimated visually. The plot radius depends on the individual tree size with a 5-m radius when the recorded DBH was below 12.5 cm, a 10-m radius when the recorded DBH was between 12.5 and Fig. 2. Distribution of forest canopy cover (a) and height (b) based on the first lidar survey (same years as Fig. 1). Fig. 3. Flowchart of the methodology to estimate annual maps of forest structural attributes (forest height, forest canopy cover and above ground biomass) from airborne and spaceborne data archives. RF stands for random forests, NFI stands for National Forest Inventory, CCDC-SMA stands for Continuous Change Detection and Classification-Spectral Mixture Analysis, AGB stands for Above Ground Biomass, while x and y stand for the last digit of the lidar survey year. M.A. Tanase et al. Environmental Research 259 (2024) 119432 5 25.5 cm, a 15-m radius when the recorded DBH was between 25.5 and 42.5 cm, and a 25-m radius when the recorded DBH was above 42.5 cm. For the Madrid region allometric models (Montero et al., 2005) of each species were used to estimate tree level above ground biomass (AGB) based on tree size. Tree level data were subsequently summed to obtain plot estimates (t ha −1 ). For the Basque Country the growing stock volume (GSV, m 3 ha −1 ) was estimated for each tree using species specific allometric models (Sanz and Soto, 1990) which were subsequently aggregated to obtain plot level estimates. GSV was used for the Basque Country to take advantage of the additional forest inventories carried out in 2005 and 2022 for a subsample of the national forest inventory grid. These datasets were available as per plot aggregates (mean height and total GSV) as opposed to the per tree measurements provided in the NFI database. Based on the 16 051 trees measured in the Madrid region the plot average above ground biomass is 101.1 t ha −1 , ranging between 0 and 618.9 t ha −1 . The plot average GSV in the Basque Country, based on the 38 475 measured trees, reaches 157.4 m3 ha −1 ranging from 0.6 to 877 m3 ha −1 . In both regions 98% of the measured trees had a height estimate. 2.4. Airborne lidar data Classified point cloud data in compressed format (LAZ) were downloaded from the Spanish National Center for Geographic Information for the available lidar flights: August–November 2010 and August–September 2016 (Madrid region) and July–August 2012 and May–October 2017 (the Basque Country) for the first and second Spanish national lidar survey, respectively. The data were examined for extent, point density, consistency, overlapping areas or gaps, and the accuracy of the existing classification. Using the ground class, digital elevation models (DEM) were interpolated for each flight and used to estimate the height above ground (i.e. normalized height) for the remaining classes. The normalized point cloud data were used to produce generic metrics (e.g., cover and density, height percentiles) representing proxies of forest structural characteristics. The lidar metrics were estimated using the Fusion Area Processor, a free and open source lidar processing software developed and maintained at the USDA Forest Service (McGaughey, 2021), as it allows for parallel processing of large datasets. All lidar metrics were estimated at 30 m pixel spacing due to the relatively low lidar point densities (0.5-to 2 points m −2 ), the spatial resolution of the Landsat data, and the size of the NFI plots. The cut-off height was 2 m to limit returns from understory vegetation. 2.5. Remotely sensed data Landsat Collection 2 Level 2 (1985–2023) for path 201 rows 032 (Madrid region) and 030 (the Basque Country) were used. Image processing was carried out in Google Earth Engine and include filtering out images with cloud cover above 90%, standardising band names of Landsat TM, ETM+, OLI, and OLI-2 sensors, filtering noisy pixels (i.e., cloud contamination, cloud shadows, snow, and radiometric-saturation) using the available quality bands, and applying topographic corrections (Soenen et al., 2005) to account for variations in reflectance values due to illumination effects from terrain (Young et al., 2017), a particularly relevant issue for the rugged terrain of both study areas (Shepherd and Dymond, 2003). In addition, we used spectral unmixing (Souza et al., 2005) to remove residual cloud and shade as the quality bands provided with the Landsat imagery are affected by frequent errors (Chen et al., 2021; Souza et al., 2013). Pixels with a cloud fraction above 0.3 in the Madrid region and above 0.2 in the Basque Country or a shade fraction equal to one were masked. The corrected images were used to estimate annual statistics for individual bands and selected vegetation indices (as described in section 2.7) as well as pixel-wise forest disturbance history, all used as predictor variables when generating forest attributes chronosequence. The forest disturbance metrics were estimated using the CCDC-SMA algorithm (Chen et al., 2021) adapted to the specific conditions of the Iberian Peninsula. Adaptation included the generation of new spectral libraries for the main biomes, Mediterranean and Atlantic, and tuning model parameters. Parameter tuning was carried out for each biome by cross-checking accuracy result for each parameter combination. The minimum number of observations used to mark a change in the harmonic models was six for both biomes. The CCDC-SMA database was used to calculate statistics of forest disturbance considering only disturbances detected within 10 years of the target year as follows: number of disturbances; number of years since last disturbance; the magnitude (last, maximum, and accumulated - sum) of the soil, photosynthetic vegetation (PV), and non-photosynthetic vegetation (NPV) fractions and the Normalized Difference Fraction Index (NDFI). The NDFI index, derived from the PV, NPV, soil, and shade fractions, is a synthetic index sensitive to subtle canopy variations (Bullock et al., 2020). The second European Remote Sensing (ERS) satellite (ERS-2) and the Environmental Satellite (Envisat) Advanced Synthetic Aperture Radar (ASAR) were used to estimate the backscatter intensity for the Madrid and respectively the Basque Country. All available SAR images were downloaded as multi-look (speckle-reduced), ground range precision images (IMP, 12.5 m pixel spacing) from the open-source European Space Agency (ESA) repository. Data from two orbital strips (20 images from 00051 to 00323 relative orbits) were downloaded, processed, and mosaicked for the Madrid region (seven dates from April to October 2010). Five images (May to September 2011 from relative orbit 00166) were downloaded for the Basque Country. Different sensors were used as the ERS-2 sensor data acquisition ended in June 2011. The IMP images were imported into the commercial GAMMA software and the relative state vectors of each image were updated using the Precise Orbit Data (POD) downloaded from the ESA repository. The IMP images acquired from the same relative orbit were transformed to slant range and concatenated to form the original strip. The strips were subsequently coregistered using as reference (master) the first date in each strip temporal series, automatically collected offsets, and a lookup table generated from the master orbital parameters and a DEM (Frey et al., 2013; Werner et al., 2005). The DEM was obtained from the first lidar flight available in each region. The offsets were modelled using least-squares regression and a third-degree polynomial function that was subsequently used to correct the lookup table needed to co-register each image to the reference image. After co-registration, each strip was multi-looked in range (1) and azimuth (4) to obtain a ground pixel spacing of approximately 30 m. The SAR intensity was transformed to the radar backscatter coefficient (γ◦) after applying the absolute calibration factors available in the metadata. The images were orthorectified, by strip, to the Universal Transverse Mercator (UTM) coordinate system using a lookup table estimated from the DEM and the orbital information of the master strip. To correct for possible inaccuracies in the input data, a refinement of the lookup table was applied in form of offsets between the master strip and a DEM-based simulated SAR image transformed to the radar geometry (Werner et al., 2002). The backscatter intensity was topographically normalized using a DEM-based scattering area (Frey et al., 2013). 2.6. Reference forest structural attributes Forest height (the 99th percentile of the normalized elevation, m) and cover (first returns above 2 m divided by the total number of first returns within the resolution cell, %) were estimated directly from the lidar point clouds. Indirect estimation, using NFI plots for model calibration and validation, was used to predict above ground biomass (Madrid region) and volume (the Basque Country). Average values of lidar-based metrics together with topographic (e.g., elevation, aspect) information were extracted at the location of each NFI plot using a 25 m radius. Regional forest maps (see section 2.8) were used to extract the main forest species at each location. random forests (RF) non-parametric modelling applied to regression (Breiman, 2001) was used to estimate M.A. Tanase et al. Environmental Research 259 (2024) 119432 6 above ground biomass and growing stock volume. RF is a flexible approach that accounts for non-linear variable interactions. RF generates a virtual forest of regression trees while employing an out-of-bag sampling with replacement strategy to limit overfitting. To further limit overfitting, few lidar metrics were used after eliminating the highly correlated (r >0.7) ones. Among correlated metrics those with higher correlation with the dependent variable (e.g., AGB) were kept. The data was randomly split between training (90%) and validation (10%) using 10 independent runs. The mean error (ME, eq. (1)), mean absolute error (MAE, eq. (2)), root mean squared error (RMSE, eq. (3)), R 2 (observed vs. predicted) were calculated for each run using the validation data. The ME quantifies systematic bias, where a negative ME indicates predictions systematically smaller than the reference values (Lang et al., 2023). The overall prediction performance were calculated by averaging all runs. The entire NFI dataset was subsequently used to calibrate a full model and generate the reference cartographic products. ME =1 n∑ n i=1 (Predictedi−Observedi)(eq. 1) MAE =1 n∑ n i=1 |Predictedi−Observedi|(eq. 2) RMSE =∑ n i=1 1 n(Predictedi−Observedi)2(eq. 3) 2.7. Chronosequence For each forest attribute, the lidar-based reference values were used to train RF models using many combinations of predictor variables starting from individual sensors (i.e., optic and radar) and progressively adding auxiliary predictors including forest type (by grouping forest species from the same family), topographic variables (elevation, slope, orientation, aspect, and roughness), and disturbance statistics. Random sampling was used to extract values for reference and predictor variables in each region (i.e., Madrid or the Basque Country). Only sample units falling in areas classified as forest (n ~ 40 000) were kept and divided in training and validation sets. To limit possible variations due to sampling, the average error and R 2 were estimated across 10 independent runs after splitting the data into training (90%) and validation (10%). The entire dataset was subsequently used to calibrate a full model and generate the annual cartographic products forming the chronosequence of each forest structural attribute (i.e. tree height, forest canopy cover, above ground biomass). An independent validation was carried out using as reference the forest attributes estimated from the second lidar survey and the corresponding annual maps. To this end, random sampling was used to match the predicted and the observed values at pixel level and estimate the RMSE, ME, MAE, and R 2 . For the Basque Country further comparisons using in situ estimates from 2005 (n =3103) and 2022 (n =420) and forest gridded metrics estimated from a stereo point cloud (n =450 000, 1 ha grid, 2021 survey) were also carried out. For each region and forest attribute, a sensitivity analysis was carried out to determine the number of predictor variables at which the decrease of MAE flattens out. The analysis was carried out by using all predictor variables to determine variable importance, determined as the mean decrease in error when including the variable in the model, followed by runs that used an increasing number of variables (5–60) in steps of five after sorting the importance from large to small, i.e., the most important variables were added first. RMSE, ME, MAE, and R 2 were estimated for each model to determine the cut-off threshold. The selection of the final model, for each region and attribute, was carried out considering the trade-off between estimation error, model complexity, and the temporal nature of the output cartographic product. The aim was to use a reduced number of predictor variables draw from few types (e.g., optic, SAR, disturbance, topography) within a common framework across regions and forest attributes. Predictor variable selection for each forest structural attribute was carried out by iteratively assessing each sensor type and their combination with an increasing number of auxiliary predictor types. As predictor variables, annual statistics (10th, 25th, 50th, 75th and 90th percentiles, p) and the standard deviation (stdDev) were estimated for the green, red, near-infrared (NIR), and short-wave infrared bands (SWIR1, SWIR2) using cloud masked optical images. Images acquired from all operational Landsat sensors available each year were used. For Landsat 7 images acquired after May 31, 2003 when the Scan Line Corrector (SLC) failed (i.e. SLC-off data) the provided SLC gap mask was applied. The same statistics were estimated for two vegetation indices (VIs), the Normalized Burn Ratio (NBR) and the Enhanced Vegetation Index (EVI), and two tasselled cap components (wetness and greenness). Statistics (by band and VIs) were also estimated for three vegetation seasons, spring, summer, and autumn, using the following date of year (doy) intervals: 60–151, 152–243 and 244–334, respectively. As the number of valid pixels used to estimate the predictor variables (percentiles and standard deviation) varies due to cloud cover and the number of active Landsat sensors, four additional variables were always included: valid pixels count, and minimum, maximum, and median date of year (doy). Similarly, annual statistics were estimated from the C-band radar backscatter images acquired by the ERS-2 (VV polarization, Madrid region) and the Envisat’s ASAR (HH polarization, the Basque Country). Eight Grey-Level Co-Occurrence Matrix (GLCM) textural features (i.e., mean, variance, homogeneity, contrast, dissimilarity, entropy, second moment, correlation) (Haralick et al., 1973) were estimated for all directions (0◦, 45◦, 90◦, and 135◦) using the 50th percentile annual statistic and a 5 ×5 kernel size as it provides more balanced omission and commission errors (Balling et al., 2023). Such textural features provide neighbourhood information from adjacent pixels potentially improving image analysis (Balling et al., 2023; Hethcoat et al., 2021). Due to the limited number of images available the seasonal statistics could not be estimated for the radar sensors. 2.8. Forest/non-forest mask The forest structural attributes were estimated only in areas classified as forest in each regional forest map. For the Madrid region the forest/shrubland species map was generated in 2008 by the General Directorate of Biodiversity and Natural Resources, Ministry of the Environment, Territorial Planning and Sustainability with the latest update being applied in April 2022. Notice that yearly maps are not available as the map is continuously updated on the official repository (https://datos.comunidad.madrid/catalogo/). For the Basque Country, forest species maps are available for different years (2005, 2010, 2016, 2018, 2020, 2021, 2022) on the official repository (https://www.geo. euskadi.eus/). Information from the forest species map closest in time with the lidar data acquisition was used when modelling the reference forest attributes (i.e., AGB and GSV) for each lidar survey. For the chronosequence estimation the latest forest map was used to mask non-forest areas in both regions. According to the Spanish national forest inventory data, the increase in forest area between the first (1970‵s) and the fourth (2010‵s) surveys was about 2% for the Basque Country and 11% for the Madrid region. The use of the latest forest species map to mask nonforest areas seemed appropriate as it considers the gains during the last decades. It also allows consistent comparison across years avoiding oscillations due to potentially variable forested areas. 3. Results 3.1. Reference data Except for forest canopy cover in the Basque Country (R 2 =0.30), M.A. Tanase et al. Environmental Research 259 (2024) 119432 7 strong relationships were observed between the in-situ data and the lidar estimated tree height and forest canopy cover (0.61<R 2 <0.74, Fig. 4). The root mean absolute error (MAE) for forest height was 1.3 m (13% relative MAE) for the Madrid region and 2.8 m (20%) for the Basque Country while for the forest canopy cover (FCC) was 12% and 15%, respectively. The modelling error for AGB in the Madrid region was 20.3 t ha −1 (30%) while for the Basque Country the GSV was modelled with a MAE of 49.2 m 3 ha −1 (33%). The most important predictor variables for modelling above ground biomass and stock volume were the canopy average height, the proportion of returns in the 5–10 m strata, and the main forest species, all within the top five variables in both regions (Fig. 5). 3.2. Chronosequence of forest structural attributes Overall, differences between models including optical predictor variables were marginal, except when information related to topography and species was added (Table 1). Out of bag errors for the final models (Table 1 last column) were similar with those observed during the evaluation of different combinations of attributes. Notice that Table. 1 shows error metrics only for selected combinations of predictor variables. Results for models based on SAR data (except for the reference one) as well as those based on seasonal statistics or combination of seasonal and yearly statistics of optical data are not shown as their utility was limited as explained below. The importance of the predictor variables varied across models, regions and forest attributes (Fig. S1). When considering the first 30 most important variables across the four model types based on Landsat data (e.g. optic, optic-disturbance, optic-disturbance-topography, and opticdisturbance-topography-species), three forest attributes (AGB, FCC, H) and two regions and (i.e., 24 models), two predictor variables (GREENNESS_stdDev and NIR_p10) were common to all models while further 10 (NIR_p25, NIR_p50, NIR_stdDev, NBR_stdDev, GREEN_p75, SWIR1_stdDev, SWIR2_stdDev, GREENNESS_p10, GREENNESS_p50, WETNESS_stdDev) were common to 75% of the models. The auxiliary information related to Landsat images availability, count and minimum date of year, was important for 80% of the models. For the 12 models that included topographic information, elevation was always present within the 30 most important variables while species was always present for the six models that included this variable. For the Madrid region, 18 of the selected predictor variables were common between the three forest attributes modelled using the combination of optic and topographic features while for the Basque Country 11 were common (Fig. S2) when only considering the model selected to estimate the chronosequence (i.e. optic-disturbance-topography). Across regions, seven variables were common to all six models (elevation, GREENNESS_p10, GREENNESS_p50, GREENNESS_stdDev, NBR_stdDev, NIR_p10, and NIR_p25) while further 14 were common to at least four models (EVI_p10, GREEN_p75, GREEN_stdDev, GREENNESS_p75, GREENNESS_p90, NIR_p50, NIR_p75, NIR_p90 NIR_stdDev, SWIR1_stdDev, WETNESS_stdDev count, doy_minimum, and roughness). The effect of modelling parameters, i.e., number sampling units (10 000, 15 000, 20 000, 40 000) and number of trees in the regression model (250, 350, 500) on model accuracy was also marginal with median difference across models below 5% for both regions. Modelling errors flatten out when reaching between 20 and 30 variables depending on the forest attribute and region (Fig. 6). Although the number of predictors at which estimation accuracy flattened varied, a fixed number (30, see Fig. S2) was used across both regions when training the final models to simplify the process, since adding more predictors did not change the results. The use of seasonal statistics allowed for marginally smaller MAE (Table S1, supplementary material). However, the few image acquisitions available at the beginning (1984–1989) and towards the end (2011–2013) of the Landsat 5 mission coupled with the Landsat 7 scan line error (since June 2003) and the setbacks of Landsat 8 mission development resulted in no statistics for many pixels for some seasonsyear partitions particularly in the Basque Country where cloud cover is frequent. Therefore, seasonal metrics were discarded to avoid incomplete mapping products (Fig. 7). Nevertheless, as improved temporal frequencies are becoming available, most notably through image collection from virtual constellations such as the Harmonized Landsat Sentinel product (Claverie et al., 2018), the use of seasonal metrics may become a viable option albeit over shorter monitoring periods. 3.3. Cross comparisons and trends analysis Cross comparisons of the estimated chronosequence maps with independent reference data obtained from the second lidar flight showed large agreement for the Madrid region (0.7 <R 2 <0.8, Table 2). The relative MAE was 30% for heights, 42% for AGB, and 45% for FCC while the ME was 6.8 t ha −1 , 5 %, and 0.15 m for above ground biomass, forest canopy cover, and forest height, respectively. For the Basque Country lower agreement was observed (0.6<R 2 <0.7, Table 2) with relative MAE reaching 35% for stock volume, 23% for forest canopy cover and 26% for forest height. All forest metrics were underestimated with ME reaching −30.3 m 3 /ha −1 , -8% and −0.9 m for height. For height and biomass (or volume), the highest MAE by species was observed for pine forests (Fig. S3) in both regions while errors for FCC were similar across species. For forested areas, relative MAE (data not shown) reached the highest values for small conifers (Madrid region) and eucalypt forests (the Basque Country). For shrublands, relative errors of all forest Fig. 4. Comparison of above ground biomass (AGB, a) forest canopy cover (b) and forest height (c), and above ground biomass (AGB, c) from in-situ field measurements (4th Spanish National Forest Inventory plots) against the corresponding lidar-based reference layers. Height and forest canopy cover were estimated directly form the lidar data. AGB and growing stock volume (GSV) were modelled using random forests and lidar-based metrics based on the first lidar survey. For visualization, a conversion factor (0.55) was applied for the Basque Country to estimate above ground biomass from volume. The correction was estimated as the mean wood density for the main species (IPCC, 2006). Dashed and dotted lines show linear regression for Madrid (black circles) and Basque Country (grey circles) regions, respectively. Black and grey text show errors for Madrid region and the Basque Country, respectively. M.A. Tanase et al. Environmental Research 259 (2024) 119432 8 attributes were twice as higher when compared to forests. Except for FCC in the Basque Country, the estimation error increased with increasing forest height (Fig. S4). Agreement with additional in-situ estimates (2005 and 2022) and forest gridded metrics was lower than with the 2nd lidar survey, with R 2 values decreasing (down to 0.24–0.38) and mean absolute errors increasing (Table 2). On average, the stock volume was underestimated by 10 m 3 ha −1 and the forest height was overestimated by 2.4 m when compared to 2005 and 2022 in-situ datasets. When compared to the gridded metrics (2021 dataset) the errors reached 80.7 m 3 ha −1 and 4.9 m for GSV and height, respectively (Table 2). Notice that similar errors (81.9 m 3 ha −1 for GSV and 3.7 m for height) were also observed between Fig. 5. Selected predictor variables for models used to generate the reference above ground biomass in the Madrid region (a) and growing stock volume in the Basque Country (b). Variables sorted in descending order by the importance (i.e., decrease in mean squared error from inclusion of the variable across all random forest trees). Table 1 Modelling errors for above ground biomass (AGB), forest canopy cover (FCC), tree height (H), and growing stock volume (GSV) as a function of sensors and auxiliary variables combinations. R 2 represents the coefficient of determination between observed and predicted values while MAE stands for mean absolute error. Region Forest attribute Optic (n a =58) SAR (n = 14) Optic Disturbance (n = 75) Optic Disturbance Topography (n = 79) Optic Disturbance Topography Species (n = 80) Optic Disturbance Topography Species SAR (n = 93) Out of bag errors for the selected predictors (n =30, Fig. S1) R 2 /MAE Madrid AGB (t ha −1 ) FCC (%) H (m) 0.77/13.8 0.77/0.09 0.61/2.0 0.1/30.7 0.16/0.2 0.06/3.6 0.77/13.8 0.76/0.09 0.61/2.0 0.80/13.3 0.79/0.09 0.64/1.9 0.83/11.9 0.83/0.08 0.70/1.7 0.84/12.0 0.84/0.08 0.71/1.8 0.77/13.6 0.77/0.09 0.62/1.94 Basque country GSV (m 3 ha −1 ) FCC (%) H (m) 0.39/61.6 0.48/0.2 0.44/4.4 0.03/82.2 0.11/0.2 0.05/6 0.38/60.8 0.48/0.2 0.44/4.3 0.41/60.2 0.5/0.1 0.47/4.3 0.42/59.2 0.52/0.1 0.48/4.2 0.43/58.9 0.52/0.1 0.5/4.1 0.39/61.2 0.47/0.16 0.44/4.38 a n - number of predictor variables used in the model. Fig. 6. Decrease of MAE (mean absolute error) and R 2 (predicted vs. observed values) as a function of the number of predictor variables for above ground biomass (AGB), forest canopy cover (FCC) and tree height (H) in the Madrid region (a, b, c) and the Basque Country (d, e, f). For visualization, a conversion factor (0.55) was applied for the Basque Country to estimate above ground biomass from volume (see caption in Fig. 4). M.A. Tanase et al. Environmental Research 259 (2024) 119432 9 the 2022 in-situ data and the 2021 gridded metrics. For all forest variables, ascending trends (positive slope) were observed since 1985 (Fig. 8). Year to year estimates were rather stable for the Madrid region while for the Basque Country higher variation was observed particularly at the beginning of the study period. For the final years (2019–2023) the trend shows decreasing values in the Basque Country in agreement with the large area (40 000 ha) of Pinus radiata forests affected by fungi. (Mycosphaerella sp.) in 2017 and 2018. Within the analysed period, the average forest height increased from 6.6 ±0.43 m to 7.9 ±0.18 m in Madrid region and from 15.8 ±0.77 m to 17.3 ±21 m in the Basque Country while the average forest canopy cover increased by 11% (from 21 ±2% to 32 ±1%) and 4% (from 54 ±4% to 59 ±3%), respectively. Total above ground biomass increased by about 62% in the Madrid region while total volume increased by 14% in the Basque Country with mean values increasing from 31.9 ±3.6 to 50.4 ±1 t ha −1 and respectively from 137.8 ±8.2 to 151.5 ±3.8 m 3 ha −1 . These values were estimate as the difference between the modelled data for the years falling in the first decade (1985–1989 average) and the years falling in the last decade (2020–2023 average). As a static forest mask was used, the average values seem related to changes in forest structure across the studied period. 4. Discussions We retrieved relevant forest structural variables (i.e., height, forest canopy cover and, biomass) from long-term archives of remotely sensed data to understand the evolution of forested areas for the last four decades at regional levels for temperate and Mediterranean forests in the Iberian Peninsula. The study dealt with a range of issues including sensors synergies and limitations, modelling approaches, and predictor variables selection, to ensure a temporally and spatially consistent dataset. We used a bottom-up approach starting with the in-situ measurements carried out on gridded variable-radius field-based plots. To avoid random variations in the extracted lidar metrics a fixed radius was used (25 m following the National Forest Inventory) under the assumption that expansion of field measurements holds for small areas usually occupied by homogeneous forest patches. This also allowed binning many lidar hits as the first Spanish national lidar flight point density was low (0.5 points m −2 ). The reference forest variable Fig. 7. Number of mapped forest pixels as a function of the type of statistics estimated from optical data. Table 2 Cross-comparison of chronosequence maps against the 2nd lidar survey (n =45 000) and auxiliary datasets in 2005, 2021, and 2022 (only in Basque Country) for the above ground biomass (AGB), growing stock volume (GSV), forest canopy cover (FCC) and height (H). Region Forest attribute Cross 2nd lidar flight comparison In situ 2005 In situ 2022 Gridded metrics from stereo point cloud (2021) a R 2 /Mean Absolute Error Madrid AGB (t ha −1 ) FCC (%) H (m) 0.76/18. 0.81/0.12 0.71/2.3 -/- -/- -/- Basque Country GSV (m 3 ha −1 ) FCC (%) H (m) 0.60/66.9 0.67/0.16 0.62/5.0 0.22/94.3 -/−0.23/5.8 0.31/103.8 -/−0.28/6.0 0.35/80.7 0.40/0.19 0.38/4.93 a Available at 100 m pixel spacing on the Basque Country spatial data repository (geo.euskadi.eus). Fig. 8. Forest structure trend between 1985 and 2023 for the Madrid region and the Basque Country for above ground biomass (AGB) and growing stock volume (GSV), mean forest canopy cover (FCC), and mean forest height. Notice that the trend for total biomass (and volume) value is affected by year-to-year variation in the number of mapped pixels, related to cloud masking, albeit such variation is marginal (standard deviation <5000 ha) M.A. Tanase et al.