Biomass equations for rockrose ( Cistus laurifolius L.) shrublands in North-central Spain

Aim of study: To construct biomass weight equations for rockrose ( Cistus laurifolius L.) shrublands in North-central Spain comparing different methodologies and evaluating the applicability of the current Spanish open PNOA-LiDAR data. Area of study: The growing extension of Mediterranean shrublands associated with a high wildfire risk in a climate change scenario is considered a relevant source of biomass for energy use and bioproducts. Quantifying the biomass load of the shrublands provides essential information for adequate management, calling for the development of equations to estimate said biomass loads in the most extensive mo nospecific shrublands. Materials and methods: Biomass dry weight from 290 destructive sampling plots (ø4m) and 426 individual plants along with LiDAR data from PNOA were related to dasometric parameters to fit weight per surface and weight per plant equations. Main results: Three new equations improve rockrose biomass estimations in North-central Spain: a) Weight per unit area ( t DM .ha -1 ) equation (Eq. 1) based on apparent biovolume (product of crown cover in percentage by average height in meters) (R adj2 0.69, MAE 26.1%, RMSE 38.4%); b) Weight per plant ( kg DM .plant -1 ) equation (Eq. 2) from height and crown diameter (R adj2 0.87, MAE 26.5%, RMSE 45.2%) and c) Weight per unit area equation ( t DM .ha -1 ) (Eq. 3) based on LiDAR data contrasted with field data (R adj2 0.89, MAE 15.1%, RMSE 22.9%). Research highlights: Eq. 1 and Eq. 3 combined with high resolution LiDAR information offer rockrose ( Cistus laurifolius L.) biomass estimations without added field work costs that are an improvement on certain more general studies carried out in other areas of Spain.


Introduction
Wildfires are a major hazard throughout Europe, producing large environmental and economic losses and having an impact on human lives. Over 40,000 fires per year were reported between 2010 and 2016 in Greece, Spain, France, Italy and Portugal, where around 85% of the total burned area in Europe was located (EC, 2018).
According to the EU´s official soil database (LUCAS, 2018), five Mediterranean countries have over 50% of the EU28 shrublands (16 Mha), of which 8.4 Mha are located in Spain. In this country, currently ranked second in fire incidence after only Portugal (San- Miguel-Ayanz et al., 2017), between 54% and 83% of the total forest area burned from 2005 to 2015 was covered by shrublands (MAPA, 2019).
Shrub formations play an important ecological role in ecosystem restoration San Miguel et al., 2004;Maestre et al., 2009;Rey et al., 2009), in soil protection (Bochet et al., 2006;Pueyo et al., 2013), in the nutrient and carbon cycle (Yarie, 1980;Van Cleve & Alexander, 1981;Chapin, 1983), in biodiversity (Noss, 1990); Mangas et al., 2008), in carbon fixation capacity (Navarro & Blanco, 2006;Fonseca et al., 2012;Gratani et al., 2013;Pasalodos et al., 2015), in obtaining essential oils (Küpeli Akkol et al., 2012;Orhan et al., 2013;Karim et al., 2017;Mediavilla et al., 2021) and in bioenergy (Viana et al., 2012;Mediavilla et al., 2017;Esteban et al., 2019;Bados et al., 2020). However, shrublands -currently considered a resource of low economic value -usually suffer from lack of management which, together with the abandonment of traditional forestry uses, leads to a forest structure more prone to wildfires (Wessel et al., 2004, Rigueiro-Rodríguez et al., 2008. In most cases, the uncontrolled concentration of shrubs is frequently the cause of the start and spread of new forest fires in European Mediterranean countries (Baeza et al., 2002, Núñez-Regueira et al., 2004, García-Hurtado et al., 2013, Mediavilla et al., 2017. On the other hand, the EU Renewable Energy Directive 2018/2001/EC stipulates that 32% of total energy consumption by the member states should come from renewable sources by 2030. In this context, Spain has recently approved the Spanish Circular Economy Strategy (MITE-CO, 2020) to achieve a sustainable, decarbonised, efficient (when using natural resources) and competitive economy. New biomass resources are required to face the increasing demand of bioenergy in a future bioeconomy (EC, 2012;Scarlat et al., 2015;Lainez et al., 2018), especially those that do not compete with other uses of biomass, such as food, feed or products for the wood and fiber industry. Biomass derived from shrub formations is gaining importance, as a complement of the biomass derived from clearing, thinning, pruning and that obtained from plantations of fast growing species. In this context, the ENERBIOSCRUB LIFE+ project stressed the need for sustainable mobilisation of new biomass resources through the production of sustainable solid biofuels from the mechanised cleaning of shrublands to mitigate wildfire risk.
According to the Spanish Strategy for the Development of Forest Biomass (MARM, 2010), the potential of forest biomass in Spain is close to 6.6 million t DM .year -1 (tons of dry matter per year), of which 4.5 million t DM .
year -1 correspond to tree-covered forest, and 2.1 million t DM .year -1 to shrublands.
In Spain, there are studies that estimate the biomass in shrublands through allometric equations based on field inventories (Patón et al., 1999;Navarro & Blanco, 2006;Ruiz-Peinado et al., 2013) and also general methods to quantify Cistaceae bushes and its annual growth through equations based on field measurements or on parameters obtai-ned from the National Forest Map of Spain Pasalodos et al., 2015;Montero et al., 2020), but they have been built with data outside the region of Castile and León, a North-central Spain region that represents 18.6% of the country's surface area, where Cistus sp. occupies more than 326,000 ha (MFE, 2020). Besides, there are no specific biomass estimation models for rockrose (Cistus laurifolius L.) shrublands, except some estimations of shrub biomass availability along two geographical transects in the Iberian Peninsula (González-González et al., 2017a) and some equations to estimate rockrose weight per plant based on small samples from specific places (Pérez & Esteban, 2008).
These facts highlight the need to develop new equations that can improve biomass predictions in North-central Spain, where rockrose (Cistus laurifolius L.) is widely spread and forms monospecific systems of extensive coverage, occupying former pastures and marginal lands over more than 175,000 ha (MFE, 2020).
Additionally, recent studies call for deepening and broadening into forest applications from LIDAR (Light Detection and Ranging) data available for the whole Spanish territory within the framework of the National Plan for Aerial Orthography (PNOA, 2010), to generate shrubland models, as it is a scantly explored work area (Montealegre, 2017). LiDAR is an active remote sensing system based on a laser scanner that allows a better description than any other known system of soil and vegetation structures in areas with dense vegetation. LiDAR sensors can penetrate the vegetation cover and capture information from different height layers (trees, shrubs and soil), which enables this technology to characterize vegetation structures and to quantify biomass stocks (Bernal et al., 2017). LiDAR data acquisition over large areas is becoming a widespread practice in many countries thanks to the multitude of regional and national programs (countries such as Denmark, Finland, Poland, Switzerland, England, Sweden and USA acquire LiDAR data throughout their territory). In Spain, the first LiDAR coverage was completed in 2015 and the second is expected to be completed in 2021 (Fragoso-Campón et al., 2020).
The high availability of data is increasing the use of LiDAR data for forest management inventory (Domingo et al., 2018;Fernández-Landa et al., 2018). Although, the main applications of LiDAR in forest inventory have focused on forest stands (Gómez et al., 2019), its application for biomass estimation in shrublands is extremely interesting. In addition, new available LiDAR data have an increasingly higher point density, which could improve the vegetation structure estimation for low, woody vegetation (Estornell et al., 2011;Greaves et al., 2016). The development of models to predict aboveground biomass in shrublands from LiDAR data is an opportunity to quantify and understand the distribution of this important resource over large areas at a very low cost. Many studies have demonstrated the advantages of using LiDAR data for mapping surface fuel models over large areas (González-Olabarria et al., 2012, Marino et al., 2016. This also indicates the ability of this information to quantify biomass density in shrubland. However, it should be noted that the potential of LiDAR to describe shrubland structure may be seen as reduced due to poor classification of the point cloud in cases with high scrub cover and low LiDAR point density (Montealegre, 2017).
The aim of this paper is to develop new equations to predict the biomass of rockrose (Cistus laurifolius L.) shrublands in North-central Spain, based on field measurements and LiDAR data, in order to test the following hypothesis: a) New equations based on field measurements can improve rockrose biomass estimations quantified with more general methods for Cistaceae bushes. b) LiDAR information can contribute to estimate rockrose biomass stocks.

Study area
This study was carried out in three rockrose shrublands in North-central Spain, corresponding to former pasture-lands with no livestock, over an area of 110 ha in the province of Soria. The area is classified as a Mediterranean temperate climate, with milder and wetter summers than in Southern Spain, with a less extended dry season and a longer interval of cold or winter frost (Peel et al., 2007;AEMET, 2011). Fig. S1 [suppl.] shows the atmospheric characterization during 2016 in one of the study areas (Lubia), which is representative of the four locations.
The shrublands had in common: a) the presence of almost pure stands of rockrose (Cistus laurifolius L.), with an abundance of over 80% of this species; b) shrub average age between 11 and 29 years, according to a dendrochronological analysis carried out by INIA (The National Institute of Agricultural and Food Research and Technology) in 2016 (González-González et al., 2017b); c) average height between 1.1 and 1.5 m; d) low slope terrain with no stones or rocky outcrops; e) shrublands adjoining pine forests; f) no recent wildfire events.

Data sources
Field measurements and PNOA-LiDAR data were used to develop biomass estimation models. Field data collection was performed in 2016 on 3644 m 2 (290 ø4mplots) of destructive sampling in the four mentioned shrublands (Table S1 [suppl.]) and biometric data from plant community (both shrub mass and individual plants) were recorded. Airborne LiDAR information was provided by PNOA-LiDAR collected in 2010, the nearest data to the sampling date, during leaf-on conditions, at a density of 0.5 points per m 2 and a vertical accuracy of less than 0.20 m. When the study began, the only LiDAR data available was from 2010, and when data from 2017 was published, all the data processing was done. On the other hand, they were only available for the South of the province of Soria, and one of the study areas (study area 1-Torretartajo) did not have updated LiDAR data from 2017. The temporal difference between the LiDAR flight and the field data can reduce predictive capacity of LiDAR information in young scrub areas with a higher growth capacity, what causes this LiDAR to work better in mature shrublands with slower growth. This information was provided in 2 km resolution tiles covering the studied areas. It should be noted that more recent LiDAR images of all the sampled shrublands do not exist nowadays; that there were no perturbances in the scrub vegetation between the two dates (only growth) and that the average age of the vegetation is 21 years, an age from which growth slows down, so those LiDAR images are considered valid for the present study.

Weight per unit area (W, t DM .ha -1 ) equation based on field measurements
Systematic sampling was carried out to fit a Weight per unit area (t DM .ha -1 ) equation based on field measurements. The areas covered by rockrose vegetation were delimited and measured on aerial photographs (PNOA, 2010) and verified with field measurements. Sample size was based on a previous random pilot sampling carried out on 30 plots of 4 m in diameter in the same shrublands, in which a biomass average value of 15.4 t DM .ha -1 with a standard deviation of 8 t DM .ha -1 was obtained. The number of systematic sampling plots, setting a maximum error of 6%, was 288. The sampling plot centers were located at the nodes of a 55 m side square net. The corresponding UTM coordinates (Datum WGS84) were identified and located in the field with a sub-metric precision GPS. Circular plots of 4 m in diameter were marked on the terrain and the measurements of the following data were concentrated on them: shrub crown cover (CC, %), species composi-tion, number of plants per hectare (N), and average shrub height (H, m) obtained by calculating the average height of the 12 plants closest to the plot center. Subsequently, all the plants were cut at ground level, without differentiation of species composition, and were weighed with a 30 kg ± 5 g digital dynamometer. Four samples per shrubland (2.5 kg per sample), including rockrose trunk, branches and leaves, were collected and sent to the Laboratory of Biomass Characterization (LCB) at CEDER-CIEMAT in Soria (Center for the Development of Renewable Energy Sources) to determine moisture content in order to estimate dry biomass weight per plot. To analyze the shrub moisture content, the analytical sample was prepared according to the UNE 14780:2011 standard, by means of homogenization, division and drying. The analytical method, drying at 105±2 ºC, was performed in LCB following the ISO 18134-1:2015 standard.

Weight per plant (w, kg DM •plant -1 ) equation based on field measurements
Transect sampling of rockrose plants was carried out in the three study areas to fit a Weight per plant (kg.plant -1 ) equation. A total of 426 individual rockrose plants with heights between 0.2 and 2.4 meters were sampled. In this region rockrose plants rarely exceed 2 m in height, and on the other hand, small plants below 20 cm in height cannot be collected for biomass with commercial mechanization systems. To have representative plants of all sizes, at least 30 plants per 20 cm height interval were measured and evenly distributed among the three locations. The measured biometric parameters were: plant green weight (w, kg), plant height (h, m) and the average of two plant crown perpendicular diameters (d, m). Afterwards, three plants per shrubland were taken to analyze moisture content and estimate plant dry weight.

Weight per unit area (W, t DM .ha -1 ) equation based on Li-DAR data and field measurements
Airborne LiDAR point clouds were automatically classified by PNOA distinguishing ground from vegetation points. Since the study was aimed at studying shrublands, this classification was verified before data processing. Locations with shrubs of 0.4-1.0 meters high were misclassified as ground, and in order to improve the results focused on this specific forest structure type, a reclassification procedure was considered using the LASground classification method implemented in LAStools software. The point cloud classification of PNOA-Li-DAR information is automatically carried out for large wooded and treeless forest areas with very diverse structure and composition. In spite of PNOA classification works generally well, it can always be improved when it is focused on a specific forest structure type. For this reason, the LiDAR point cloud was reclassified in the formations under study. Afterwards, LiDAR tiles were processed with Fusion Software (McGaughey & Carson, 2003). A 2-meter resolution Digital Elevation Model (DEM) was generated from the reclassified ground points. This DEM was used to subtract the ellipsoidal elevation of the DEM from the Z coordinate of each Li-DAR return and normalize the LiDAR point cloud. The study area was tessellated into 20-meter pixels, computing a total of 26 LiDAR metrics for each of them. These metrics were computed using a predefined threshold height of 0.4 m, and they corresponded to a mean height (LH mean ), maximum height (LH maximum ) and minimum height (LH minimum ), mode, standard deviation, variation coefficient, variance, interquartile range, kurtosis, skewness, shrub crown cover (LCC, %) (percentage of first returns over 0.4 m), and several percentiles (ranging from the 1 st to 99 th percentile: P1, P5, P10, P20, P25, P30, P40, P50, P60, P70, P75, P80, P90, P95 and P99). Derived variables were constructed from metrics and percentiles of LiDAR heights that allowed describing the vertical distribution of biomass in shrublands. Sampling plots from the study areas were used to elaborate a Shrub weight per unit area (t DM .ha -1 ) equation based on LiDAR data and field measurements. In each plot, LiDAR mean shrub height (LH mean , m), LiDAR shrub crown cover (LCC, %) and biomass dry weight (W, t DM . ha -1 ) were considered to fit this model. Finally, the same LiDAR metrics were calculated for each field plot and used as predictive variables to estimate dry biomass.

Weight per unit area (W, t DM .ha -1 ) equation based on field measurements
Linear and non-linear (power and exponential) models based on sampling plot information were generated to fit dry biomass weight per unit area. H, CC and apparent biovolume (ABV), which is the product of H and CC (Cook, 1960), were considered as independent variables. Data analysis, statistical fits and cross-validation were implemented in R (R Development Core Team, 2008). The model was obtained with a level of significance of 0.05 for all the parameters. A cross-validation technique of 50 iterations was applied to validate the model that provided the best results in the diagnostic phase, splitting the original data into a training set of 70% of the cases and a test set of 30%. The result of the cross validation was the average of the 50 iterations. The selection of this high number of validation sets was considered highly representative for the validation process.
The following statistical metrics were used to evaluate and validate the precision and accuracy of the equations: R adj 2 , mean absolute error (MAE), root mean squared error (RMSE) and bias estimate (b).

Weight per plant (w, kg DM •plant -1 ) equation based on field measurements
Regression models were fitted to estimate dry biomass weight per plant. Data analysis, statistical adjustments and model cross-validations were implemented in R (R Development Core Team, 2008). Linear and non-linear models (power and exponential) were generated selecting the ones which provided the best results in the fitting and validation phases. The model was obtained with a level of significance of 0.05 for all the parameters. The same cross-validation procedure described in the previous section was used to evaluate and validate the biomass weight equation.

Weight per unit area (W, t DM .ha -1 ) equation based on Li-DAR data and field measurements
Linear and non-linear (power and exponential) models were generated selecting the one which provided the best results in the fitting and validation phases. A cross-validation technique was applied to evaluate the accuracy of each model comparing R adj 2 and RMSE. Data analysis, statistical adjustments and model cross-validations were implemented in R (R Development Core Team, 2008). A stepwise approach was followed to perform the variable selection through an R MASS package (Ripley et al., 2015). The model was obtained with a level of significance of 0.05 for all the parameters. The most accurate model developed, with the lowest RMSE, was used to predict the dry biomass in three of the study areas in a spatially continuous way. In all the model fitting processes the Shapiro-Wilk test was used to verify the residual normality, the Breusch-Pagan test was used to verify homocedasticity and Vif test was used to verify collinearity.

Comparison of biomass weight equations
The developed weight per unit area equations were used to predict the dry biomass weight in study areas 1, 2 and 3. Technical problems related to the storage capacity and data processing of the computer, prevented the estimation of biomass in study area 4 with the equation based on LiDAR data, so this area was not considered in the comparison. The predicted values were compared with the field biomass sampling results as another way of validation. The predicted results were also compared with two biomass equations for rockrose shrubs and Cistaceae bushes: Montero et al. (2020) and Pasalodos-Tato et al. (2015). The comparison of estimations from the different models was carried out considering accuracy statistics (mean shrub weight (W), standard deviation (SD), confidence interval (CI)) and precision statistics (MAE, bias and RMSE). Weight per unit area (t DM .ha -1 ) equation based on LiDAR data and field measurements offered a single estimation per study area, not per sampling plot, so standard deviation values when applying this model in the different sampling plots were null. Finally, as another source of comparison, the biomass values estimated with the three weight equations were compared with the field biomass estimations carried out by TRAGSA S.A., a company responsible for shrub harvesting trials in study areas 1 and 3 using a harvester-baler Biobaler WB55 (Bados et al., 2020). The procedure followed by this company was based on the sum of harvested biomass and biomass left on the terrain after mechanized harvesting. The latter included inadequate shrub clearing and losses of fine material which, after being cleared, did not enter into the baling unit and fell to the ground, or went into and out of the baling unit without being part of a bale and fell to the ground. It was estimated by a systematic sampling of random approach transects establishing square plots of 0.5 x 0.5 meters (Blasco et al., 2017).

Results
Based on field measurements and LiDAR data, three biomass equations (Eq. 1, Eq. 2 and Eq. 3) were developed to predict aboveground biomass of rockrose shrublands in North-central Spain (two equations estimating tonnes of dry matter per hectare (Eq. 1 and 3) and one equation estimating kilograms of dry matter per plant (Eq. 2) .

Weight per unit area (W, t DM .ha -1 ) equation based on field measurements
The best equation based on field measurements to predict dry biomass weight per unit area was obtained through an allometric model in which the independent selected variable was ABV (Apparent biovolume) adjusted by non-linear regression. ABV is the product of the mean shrub crown cover (CC, %) and the mean shrub height (h, m). It yielded an appropriate level of precision and accuracy explaining 69% of the variance of the estimated variable, with a MAE of 4.12 t DM .ha -1 (26.1%), bias estimate of -1.01 t DM .ha -1 (6.4%) and RMSE of 6.08 t DM .ha -1 (38.4%) ( Table 1): where W is the dry biomass weight per hectare (t DM .ha -1 ), ABV is the apparent biovolume.
Regarding multicollinearity and heteroscedasticity of this model, although residuals from the regression line of Eq. 1 corresponded to a normal distribution without autocorrelation, they showed some heteroscedasticity (studentized residual grew with the independent variable values), so a Box-cox transformation was tried. The results showed a softer growing trend of the residuals, but its behavior was worse for the lower values of the independent value (overestimating the weight). Its mean absolute error was greater than that of the original Eq. 1 and it had a 95% significant lack-of-fit statistic, so it was finally discarded.

Weight per plant (w, kg DM •plant -1 ) equation based on field measurements
The best statistical fit for dry weight per plant was obtained through an allometric model in which the independent variables were h and d fitted by non-linear regression. It yielded an appropriate level of precision and accuracy explaining 87% of the variance of the estimated variable, with a MAE of 0.46 kg DM •plant -1 (26.5%), bias estimate of -0.26 kg DM .plant -1 (-17.1%) and RMSE of 0.82 kg DM .plant -1 (45.2%) ( where w is the dry weight of the plant (kg DM •plant -1 ), h is the plant height (m) and d is the average value of the plants crown perpendicular diameters (m).
Although residuals from the regression line of Eq. 2 corresponded to a normal distribution and did not show autocorrelation, they showed a certain degree of heteroscedasticity (studentized residual grew with the estimated value or w), but as most of the w-values were low (only 12 out of the 412 values were greater than 6 kg DM •plant -1 ), the effect of such heteroscedasticity was considered irrelevant. With regard to the multicollinearity between the two explaining variables, their correlation coefficient was 0.70, what could lead to eliminating one of them. Nonetheless, the analysis of the significance of the coefficients of both variables made keeping them in the model advisable, as the removal of any of them worsened the fitting quality.

Comparison of biomass estimations with the developed equations
Eq. 1 and Eq. 3 were used to predict dry biomass weight in study areas 1, 2 and 3. These estimations, together with the field biomass sampling results and the predicted values obtained with Montero et al. (2020) model (Eq. 4) and Pasalodos-Tato et al. (2015) model (Eq. 5) are shown for comparison purposes in Table 2. Equations 1, 2, 3, 4 and 5 are described in Appendix A1 [suppl.]. The weight per plant equation (Eq. 2), in spite of its good fit, was not used to predict tonnes per hectare since it requires additional sampling of the number of plants per hectare, which increases the data collection effort and can be an additional source of error. Finally, the biomass values estimated in study areas 1 and 3 using Eq. 1, 3, 4 and 5 are shown, for comparison in Table 3. This table also includes the field biomass estimations carried out by TRAGSA after the shrub harvesting trials of both study areas.
Eq. 1, based on ABV as an explanatory variable, offers the most accurate and precise estimations in study areas 2 and 3, with a bias between -0.2% and 0.2%, MAE between 27.0% and 17.2% (4.0 and 2.0 t DM ·ha -1 ) and RMSE between 35.4% and 23.7%, respectively. The most  with an average age between 11 and 16 years) than for older ones (study area 1 with average age of 29 years) (Table S1 [suppl.]). Nevertheless, Eq.1 estimations in study areas 1 and 3 are in accordance with the biomass values estimated by TRAGSA after harvesting both shrublands (Table 3).

Discussion
The developed weight per hectare equations (Eq. 1 and Eq. 3), generated through two different methodologies, allowed appropriate precision and accuracy levels, with R adj 2 values of 0.69 and 0.89 respectively, which are of the same order or better than other more general previous studies to estimate Cistaceae shrubs (R adj 2 0.76 in Pasalodos-Tato et al. (2015) and R adj 2 0.64 in Montero et al., (2013), and RMSE values below 39%). Navarro & Blanco (2006) estimated Cistus ladanifer L. biomass as a function of shrub age ranging between 1 and 12 years, with a greater R 2 (0.97) and Cistus ladanifer L. and Erica sp. with values of the same order or lower (R 2 of 0.71). Yao et al. (2021) defined allometric models to estimate shrub biomass, including crown-related volumes as predictors, with similar R 2 (0.63-0.86).
The weight per plant equation (Eq. 2) also provided acceptable levels of precision and accuracy (R adj 2 0.87 and MAE 2.4 kg DM ˑplant -1 ), but less favorable than those published by Pérez & Esteban (2008) (R adj 2 0.96 and MAE 0.062 kg DM ˑplant -1 ), although it should be noted that the cited study was based on a sample of 30 plants located at Lubia (Soria) versus the sample of 426 plants located in study areas 1, 2 and 3. Shrub biomass estimations provided by Eq. 1 and Eq. 3 are within the average biomass accumulation confidence interval estimated by González-González et al. (2017a) for a shrub formation composed of rockrose shrubs and Cistaceae bushes in Central Spain (11.93±5.8 t DM ·ha -1 ), considering an average height of 98.2±3 cm and an average crown cover of 59.9±2 %. Montero et al. (2020) and Pasalodos-Tato et al. (2015) models overestimate shrub biomass weight between 25 and 40% with respect to the field sampling results in the three study areas, with weight per hectare values being above the systematic sampling confidence intervals. This is reasonable taking into account that Eq. 1 and 3 were fitted from plants and mass values sampled in North-central Spain, while Eq. 4 and 5 were not developed from sampling plots located in the region of Andalusia (Southern Spain). Eq. 1 and 3 estimations offer biomass values closer to TRAGSA estimations after harvesting study areas 1 and 3, being underestimated in this case between 25-43% with respect to field sampling results.
LiDAR technology provides an outstanding information source for characterizing the forest structure (Naesset, 2002) and, therefore, for predicting forest variables in a spatially continuous manner (Domingo et al., 2017;Montealegre et al., 2016). However, the great majority of the studies using PNOA-LiDAR data have focused on the quantification of the forest resources without considering the shrubs ecosystems (Fernández-Landa et al., 2018;Gómez et al., 2019). As far as the authors are concerned, this study along with the works of Estornell et al. (2011Estornell et al. ( , 2012 are unique in the employment of LiDAR data for shrubland biomass estimation in Spain. The results show that low density PNOA-LiDAR can be used as an auxiliary information to obtain biomass estimates spatially continuous. The cross-validation results, in terms of relative RMSE and R 2 , are consistent with other studies. The dry biomass models based on LiDAR data yielded a relative RMSE of 26.70% and R 2 of 0.60. Estornell et al. (2012) employed high density LiDAR data to assess the biomass of a Mediterranean forest in Valencia (Spain) achieving an R 2 of 0.67 and a relative RMSE of 28%. Greaves et al. (2016) estimated shrub biomass in Arctic Tundra using also LiDAR data and reported an R 2 of 0.61. Glenn et al. (2016) modeled the biomass of semi-arid vegetation sites in the western U.S. with LiDAR data and found a R 2 of 0.56. Zhao et al. (2021) estimated shrub aboveground biomass as a function of high-density LiDAR point cloud data acquired from an unmanned aerial vehicle with a greater R 2 (0.77).
The methodology developed in this study leveraging LiDAR could be replicated at a national scale or over other territories with well-established LiDAR programs. Further research should focus on studying how LiDAR data and other sources of information (optical or SAR satellite imagery) could be integrated to estimate biomass more accurately. It was found that the calibration of regression models using optical and LiDAR metrics improved shrub biomass estimates (Glenn et al., 2016;  Equation 5 19.6 13.6  Greaves et al., 2016;Zhao et al., 2021). Optical data could also be considered as an alternative to develop a more cost-effective biomass-inventory solution when LiDAR data is unavailable free of charge. Chen et al. (2018) evaluated the capacity of Landsat data to predict shrub biomass in a semi-arid ecosystem of China with successful results (R 2 = 0.88), although a shrub mask was initially created to limit the application of the model. In this sense, the generation of shrub cartography products along with fuel variables of the canopy, derived from Li-DAR data, could also be of great importance for the generation of fuel-model maps (Marino et al., 2016). The accurate filtering of the LiDAR points cloud into ground and shrub points hinder the use of LiDAR data for shrub biomass estimation (Montealegre, 2017). This could explain the poorest performance of Eq. 3 in the study area 3 (Table 2), as suggested by a dry biomass mean value below the confidence interval of the field sampling estimation. Study area 3 is characterized by smaller shrublands with lower mean height values which could lead to a greater confusion between ground and shrub points, and thereby, producing a lower quality DEM. In this sense, higher LiDAR data density could overcome this problem resulting in more accurate predictions (Estornell et al., 2012). The PNOA-LIDAR program is on the verge of finalizing the acquisition of the second coverage with a higher point density which could be convenient for further research to assess whether more accurate estimates are obtained. Besides point density, additional factors related to temporal differences between LiDAR data and ground data could affect the biomass estimation (McRoberts et al., 2016). Regression models calibrated with LiDAR data could be applied to obtain wall-to-wall estimates without the need for field plots. In addition, the PNOA-LiDAR program has started to develop the second coverage campaign with a higher density of points, which could be convenient to verify whether more robust biomass values could be estimated. It would be advisable to carry out new research to analyze the combination of LiDAR with other sources of information, such as optical or SAR satellite imagery, that would allow better predictions. It was found that the addition of spectral information improved shrub biomass estimations (Riaño et al., 2007;Greaves et al., 2016).

Conclusions
The study contributes to improving aboveground rockrose (Cistus laurifolius L.) shrub biomass estimations in North-central Spain. The Weight per unit area (W, t DM •ha -1 ) equation (Eq. 1) is based on general data that can be obtained from the National Forest Map of Spain, and offers more accurate and precise estimations than earlier more general studies carried out in other regions of Spain, especially in young and not highly lignified shrublands.
The Weight per unit area equation (Eq. 3) based on LiDAR and field data was developed using a methodology capable of obtaining acceptable biomass estimations at a very low cost. Based on the PNOA-LiDAR second coverage campaign, currently under development, new research is needed to analyze whether other combinations of predictive variables could improve biomass estimations.
Eq. 1 and Eq. 3, combined with high resolution Li-DAR information, could allow estimating rockrose biomass stocks in a continuous way and without additional field work costs over more than 175,000 ha of rockrose (Cistus laurifolius L.) shrublands in different regions of North-central Spain, in order to contribute to a sustainable management of shrub formations for energy use, bio-based industries, and fuel modeling for assessing fire hazard. .