Non-Destructive Biomass Estimation in Mediterranean Alpha Steppes: Improving Traditional Methods for Measuring Dry and Green Fractions by Combining Proximal Remote Sensing Tools

The Mediterranean region is experiencing a stronger warming effect than other regions, which has generated a cascade of negative impacts on productivity, biodiversity, and stability of the ecosystem. To monitor ecosystem status and dynamics, aboveground biomass (AGB) is a good indicator, being a surrogate of many ecosystem functions and services and one of the main terrestrial carbon pools. Thus, accurate methodologies for AGB estimation are needed. This has been traditionally done by performing direct field measurements. However, field-based methods, such as biomass harvesting, are destructive, expensive, and time consuming and only provide punctual information, not being appropriate for large scale applications. Here, we propose a new non-destructive methodology for monitoring the spatiotemporal dynamics of AGB and green biomass (GB) of M. tenacissima L. plants by combining structural information obtained from terrestrial laser scanner (TLS) point clouds and spectral information. Our results demonstrate that the three volume measurement methods derived from the TLS point clouds tested (3D convex hull, voxel, and raster surface models) improved the results obtained by traditional field-based measurements. (Adjust-R2 = 0.86–0.84 and RMSE = 927.3–960.2 g for AGB in OLS regressions and Adjust-R2 = 0.93 and RMSE = 376.6–385.1 g for AGB in gradient boosting regression). Among the approaches, the voxel model at 5 cm of spatial resolution provided the best results; however, differences with the 3D convex hull and raster surfacebased models were very small. We also found that by combining TLS AGB estimations with spectral information, green and dry biomass fraction can be accurately measured (Adjust-R2 = 0.65–0.56 and RMSE = 149.96–166.87 g in OLS regressions and Adjust-R2 = 0.96–0.97 and RMSE = 46.1–49.8 g in gradient boosting regression), which is critical in heterogeneous Mediterranean ecosystems in which AGB largely varies in response to climatic fluctuations. Thus, our results represent important progress for the measurement of M. tenacissima L. biomass and dynamics, providing a promising tool for calibration and validation of further studies aimed at developing new methodologies for AGB estimation at ecosystem regional scales.


Introduction
Aboveground biomass (AGB) constitutes one of the main terrestrial carbon pools [1][2][3] and is considered a good indicator of many ecosystem functions and services, such as water regulation, erosion control, or C fixation [4][5][6]. In addition, it is closely related to vegetation productivity [7][8][9] and functional diversity [10,11]. Thus, continuous monitoring of AGB over space and time may provide crucial information about ecosystem status and dynamics and provide information about potential shifts in ecosystem functioning due to climate or human interactions [6].
Traditional methods for AGB estimation are not very efficient as they involve the clipping, the oven-drying, and the weighting of the plant material [12]. This is an expensive, time-consuming, and destructive methodology, which hinders its applicability as a monitoring tool. For this reason, non-destructive methodologies, such as the point intercept method [13,14] or the application of allometric equations, which relate field measurements of plant structure with plant volume and biomass [15], arose to overcome these difficulties. However, they only provide punctual information that is not representative of the landscape heterogeneity and are still not easily scaled up [16].
New generations of remote sensing instruments offer the possibility of acquiring high resolution information of vegetation features over large areas in a relatively short period of time [17,18]. Optical remote sensing information, especially multispectral imagery, has been widely employed to survey vegetation over large extents [19][20][21][22], as it provides information about plant physicochemical properties and functional traits that reflect the vegetation status [23]. However, remote sensing data can be strongly affected by factors like canopy architecture, viewing geometry, and the mixing of reflectance signals from different surface components. As a result, spectral signatures often have a limited correlation with real measurements of AGB [24][25][26][27][28]. This is especially relevant in drylands where vascular plants usually appear forming patches surrounded by open areas colonized by photosynthetically active biocrusts, rocks, and bare soil that may induce errors in the analysis of the vegetation spectral response due to the mixing of surface signals [29,30].
Alternatively, LIDAR technology is used to acquire three-dimensional (3D) point clouds over large areas of the territory that allow for the characterization of the structure of vegetation canopy [31] and biomass [32,33]. Laser scanner sensors can be placed on terrestrial (terrestrial laser scanner, TLS) and aerial (airborne laser scanner, ALS) platforms. TLSs produce dense 3D point clouds that accurately describe the AGB of individual plants in a wide range of vegetation types and ecosystems [34][35][36][37][38][39][40][41][42][43][44][45] being very useful for non-destructive field measurements and monitoring and for calibration and validation of ecosystem scale measurements and models [46]. Most available LIDAR information at the ecosystem scale is provided by ALS sensors, which present low return point density and large footprints. This limits their application in heterogeneous ecosystems due to the low structure of vegetation and their patchy distribution [45,47]. To solve these limitations, several studies explored the possibility of combining structural and spectral information from different sensors to study forest structure and physiology at different spatial scales [48][49][50][51][52][53]. Nevertheless, their applicability in drylands is still very limited, mainly because of the important challenges that these heterogenous landscapes present for remote sensing studies such as the low vegetation signal-to-noise ratios, the high soil background reflectance, the presence of photosynthetic organisms in the soil, the high spatial heterogeneity at different spatial scales, the irregular growing and phenological seasons, and the low plant sizes [30]. This is the case of the steppes of Macrochloa tenacissima (L.) Kunth (=Stipa tenacissima L., alpha grass, hereafter "espartales") in which plants are often interspersed within a heterogenous matrix of biocrusts, stones, and bare soil that interferes in the overall spectral response. This feature combined with the complex morphology and phenology of M. tenacissima L., which forms dense tussocks about 0.5-1 m in height with a high and variable proportion of dead leaves and inflorescences (necromass) [54], hinders the success of the already existing methodologies for AGB monitoring based on remote sensing [27]. Indeed, to our knowledge, there are no specific methodologies combining spectral and structural information for AGB estimation and monitoring of M. tenacissima L. (or any other alpha grass species with similar morphology) at any spatial scale.
The development of simple and non-destructive strategies able to provide low-cost biomass information over M. tenacissima L. steppes is necessary due to these landscapes being among the most representative ecosystems of the Mediterranean region, covering a large fraction of the territory from Libya to the southeastern portion of the Iberian Penin-sula [55,56]. However, their extent has been dramatically reduced to less than 35,000 km 2 during the last decades, mainly due to an intensification of land use over the region (increase in grazing pressure within these ecosystems [56]) and more adverse climatic periods. Moreover, recent studies revealed that this trend may continue in the future due to ongoing global change [57]. Consequently, the ability to rapidly assess the conservation state of this important habitat, which is recognised as a priority habitat by the EU (6220-Pseudosteppe with grasses and annuals "Thero-Brachypodietea") and where M. tenacissima L. can contribute more than 90% of AGB [58], is a major challenge for land managers and conservation authorities.
This study proposes a new non-destructive methodology for monitoring spatiotemporal dynamics of AGB and green biomass of M. tenacissima L. plants by combining structural and physiological (spectral) information, measured in the field. This can be used to support ground monitoring programs and also for the rapid calibration and validation of ALS-multispectral-based AGB estimation at coarser spatial scales. To do this, we combine structural information obtained from TLS 3D point clouds with spectral information. More precisely, we aim to: (1) determine whether TLS sampling could rapidly create a high-density, three-dimensional (3D) model of M. tenacissima L. able to provide AGB more accurately than traditional field methods based on plant measurements; and (2) develop a new methodology to estimate not only biomass but also plant physiological condition (green biomass) by combining both TLS and spectral measurements. These aims are a critical issue for dryland ecosystems over the Mediterranean region, where the spatial distribution of vegetation results from interactions between climate, soil, topography, facilitation, competition, and human activities [59], configuring a wide heterogeneous landscape that makes AGB monitoring difficult.

Study Sites and Field Sampling
For this study, we used M. tenacissima L. plants from 3 representative Mediterranean steppes ecosystems from the SE of the Iberian Peninsula ( Figure 1): El Cautivo experimental area (37 • [60]. All sites are within the province of Almeria (SE Spain) and have a semiarid Mediterranean climate characterized by hot, dry summers and mild temperatures the rest of the year, with rain falling mostly in winter. The three sites show subtle differences in precipitation and temperature but contrasted soil properties and topography that induce differences in the spatial distribution and size of M. tenacissima L. plants.
El Cautivo study site is a dry badlands area (mean temperature~17.8 • C, annual precipitation~235 mm/year) in which M. tenacissima L. covers approximately 25% of the area [59,61]. Soils are characterized by low to moderate development, high salinity, and low organic carbon content (from 5 to 10 g/kg), the most frequent soil types being Endoleptic Leptosols, and Calcaric Regosols [61]. Soil properties together with steep slopes and high levels of solar radiation make the terrain very susceptible to erosion, limiting M. tenacissima L. distribution to north and northeast hillslopes [62]. At this site, M. tenacissima L. plants usually reach big sizes and form wide tussock groups in dense vegetation patches, among which lichen-dominated biocrusts often appear [63]. Balsa Blanca and Las Amoladeras are located within the Cabo de Gata Natural Park, and both show similar climate (mean temperature~18 • C, annual precipitatioñ 220 mm/year) and topography, being developed over a flat alluvial fan [60]. Dominant soils in theses ecosystems are Calcaric Leptosols (Balsa Blanca) and Haplic Calcisols (Las Amoladeras) [64] with sandy loam texture, high carbonate content, and rock fragments [65]. Soils are deeper and they have higher soil organic carbon (SOC) content at Balsa Blanca (maximum soil depth 20 cm,~15 g/kg SOC) than at Las Amoladeras (maximum soil depth 10 cm,~10 g/kg SOC) [66], mainly due to the historical land use and grazing intensity differences that lead to changes in M. tenacissima L. coverage and biomass between them. Overall Las Amoladeras presents the lowest vegetation coverage (~25%), consisting of smaller M. tenacissima L. plants, with open spaces between vegetation patches usually covered by biocrusts or bare soil. Conversely, vegetation cover in Balsa Blanca is higher (~65%), reaching larger sizes of M. tenacissima L. plants, and therefore forming dense tussocks that constitute bigger vegetation patches with smaller open areas between plants. More details about Las Amoladeras and Balsa Blanca are provided in [67] and in [60].
We selected 10 M. tenacissima L. plants on each site (total n = 30). During the plant selection and identification process, we tried to include the maximum range of variation in M. tenacissima L. volume and the largest variation in the ratio between photosynthetic (hereafter "green biomass (GB)") and no photosynthetic biomass (hereafter "dry biomass") within each site. On each plant, mean height and diameter were measured in three different positions and averaged. Then, we acquired a 3D point cloud of each plant using a Leica ScanStation 2 terrestrial laser scanner (TLS). This system operates in the green region of the spectrum (532 nm) and has a maximum range of 300 m (at 90% albedo), laser spot size of 6 mm, beam divergence of 0.15 µrad, and accuracy and precision of 4 and 2 mm [68]. To reduce shadowing effect due to the geometry of point cloud acquisition, two different scans at 1 mm of spatial resolution were done on each plant, from two different positions (azimuth 90 • , azimuth 260 • ). Finally, we measured plant canopy reflectance using an ASD FieldSpec HH field radiometer (325-1075 nm, 512 channels) with 3 nm of spectral resolution and 25 • IFOV. A total of three replicate measurements, each consisting of the average of the five individual spectra, were collected for each plant. All measurements were conducted at solar noon, and a Spectralon (ASD Inc., Boulder, CO, USA) white reference sample was measured between sampling of different plants.
Once field TLS and spectral measurements were done, individual M. tenacissima L. plants were cut at the ground level, stored in plastic bags, and transported to laboratory where green leaves (photosynthetically active biomass) were separated from dry biomass. Since M. tenacissima L. is a tussock-grass plant, plant biomass was composed of green leaves, inflorescences, and necromass. Green and dry materials were oven-dried at 70 • C for 72 h or constant weight following the procedure recommended by FAO [69].

TLS Data Processing and Volume Estimation
The different point clouds of each M. tenacissima L. plant were merged and georeferenced using five reflective targets common to both scans using Cyclone 7.1 software (Leica Geosystems). We built final 3D point clouds for each plant that were used to test three different approaches for deriving AGB: (i) a volumetric approach based on the 3D convex hull algorithm, (ii) a voxel-based model, and (iii) a raster-based model ( Figure 2). All these methods have been applied with accurate results over a wide range of vegetation types, such as shrublands [35,38,45] or forests [44], but they rely on different approaches and assumptions, the final accuracy of each one being dependent on plant density and structure. The 3D convex hull approach fits a 3D surface to the point cloud [45] by assuming that all space within the 3D surface model is occupied by plant biomass. The same assumption is presumed in the raster-based model. However, volume estimation in the raster model is based on a completely different approach, as it involves the creation of a digital surface model in which each cell represents the maximum value of all points within it. Then, plant volume is determined as the product of the cell area and the average height.
The main advantage of the voxel approach is the complete structural representation of the plant by attributing points to a voxel based on coincident 3D location. By doing this, it provides better results in vegetation with irregularly distributed canopy gaps, but it may not be appropriate in very dense vegetation where the TLS laser cannot penetrate within the canopy [44]. To solve this limitation, voxel approach was implemented according to the method developed by Vonderach [70] for filling interplant spaces where laser pulse cannot penetrate. We used the rLiDAR package [71] to fit 3D convex hull and the VoxR package [72] to build the voxel models. Raster models were built in QGIS 3.14 using the maximum height of each point in the point cloud. Spatial resolution of raster-based and voxel models may influence the final plant volume estimation [44], the best resolution being determined by point density and any occluded (non-scanned) space within the plant canopy. To determine the effect of spatial resolution on AGB estimation accuracy, 4 different resolutions (1 cm, 2 cm, 5 cm, and 10 cm spatial resolutions) were tested in both raster-based and voxel approaches.

Field Spectra Data Pre-Processing
Before the spectral analyses, all data were subjected to a pre-treatment that consisted of removing the noisy bands (350-400 nm and above 950 nm) and applying a cubic polynomial smoothing filter with a 17 nm window size (Savitzky-Golay filter) [73]. Individual spectra of each plant were then averaged to produce the plant canopy spectrum (average of 9 individual spectra). Finally, we estimated four well known vegetation indices such as the Normalized Difference Vegetation Index (NDVI; [74]), the Modified Normalized Difference Vegetation Index (mNDVI; [75]), the Total Ratio Vegetation Index (TRVI; [76]), and the Green Normalized Difference Vegetation Index (GNDVI; [75]) from ASD field spectra according to Equations (1)-(4), respectively. We chose these indices as they have been demonstrated to correctly characterize the fraction of photosynthetic active radiation absorbed by vegetation (fPAR) [77,78], being good indicators of the green vegetation fraction, green LAI fraction, or biomass [20,47,[77][78][79], and they can be easily obtained from high and medium resolution satellite images and multispectral cameras.

Plant Biomass Estimation
The total AGB of each plant was modelled by fitting linear ordinary least squares (OLS) regressions between the different volume estimation methods derived from TLS point cloud (3D convex hull, voxel, raster) (see Section 2.2) and the corresponding harvested plant biomass. In order to develop an equation with biological signification, all intercepts in regression analyses were set to zero assuming that no plant volume corresponds to no plant biomass and therefore avoiding negative values for biomass estimations in low size plants. As classical statistical regression may not effectively describe the complex nonlinear relationship between AGB and plant volume, we also tested a gradient boosting model, which is a powerful ensemble machine learning algorithm that has been proved to supply excellent performance in biomass estimation modelling [80]. More precisely, it fits boosted decision trees by minimizing an error gradient between observed and predicted values.
The model was applied by using XGBoost [81] which provides a flexible and highly scalable tree structure enhancement model that performs a second-order Taylor expansion for the objective function and uses the second derivative to accelerate the convergence speed of the model while the training is running. Optimum maximum number of boosting iterations (nrounds) and learning_rate values were set by testing the complete set of combinations between nround values from 1 to 1000 and learning rate values of 0.01, 0.02, 0.05, 0.1, 0.2, 0.3, and 0.5 in a 3-fold cross validation approach (see Supplementary Figure S1), whereas default values provided in XGBoost R package were set for all other parameters. Optimum model configuration for each TLS volumetric model (3D convex hull, voxel, and raster) was applied to the entire dataset for comparison with linear estimations.
Obtained results by both ordinary least-squares and XGBoost were compared with biomass estimation based on traditional allometric equations consisting of fitting a specific geometric form that represents plant canopy. Concretely, we use a hemisphere model based on plant diameter measured in field. Additionally, we explored the use of field measurements of plant height and diameter as direct predictors of AGB (see field sampling in Section 2.1).

Green Biomass Estimation
Green biomass fraction (onwards GBF) of each plant [%] was calculated as the ratio between GB and total harvested AGB previously separated in the laboratory. Then, we explored the relationships between GBF and the spectral indices described in Section 2.3 by fitting linear regression models. Finally, the total amount of GB of each plant was predicted by combining AGB estimation obtained by best TLS methods described in Section 2.4 and the GBF estimated based on the value of the different spectral indices, according to Equation (5): where AGB TLS methods is the predicted total biomass in grams using TLS approaches (voxel at 5 cm of resolution; 3D convex hull and raster at 2 cm of resolution) derived for point cloud, and GBF VI is the [%] of green biomass predicted using the spectral index with the best predictive capacity (see Section 2.3). In a manner similar to AGB estimation, XGBoost was used to fit a gradient boosting regression model between GB, the best VI, and plant volume obtained by the three different TLS volumetric models using the optimum values of nrounds and learning_rate obtained during the calibration of the model (Supplementary Figure S1). Performance of the different methods and equations for AGB and GB estimation was evaluated based on Adjust-R 2 , the root-mean-square error (RMSE), the relative rootmean-square error (rRMSE), and the mean bias error (MBE) values between harvested and predicted biomass.

Results
AGB of the 30 individual plants analysed in this study ranged from 184 g to 6655 g, with GBF values between~5% and~43%. As observed in Table 1, the set of plants from each site comprised a wide range of size and biomass, the mean and maximum values being higher at El Cautivo than at Balsa Blanca and Las Amoladeras, respectively. There were also differences in GBF values of the different plants within each site and between sites. M tenacissima L. plants from El Cautivo and Las Amoladeras had a mean GBF value of~14%, while Balsa Blanca had a higher GBF, with an average value of~28%. The estimated plant volume also varied between ecosystems, with the highest values obtained at El Cautivo, followed by Balsa Blanca and Las Amoladeras, respectively. Among approaches, the plant volume estimated by a hemisphere model presented lower values than the models based on TLS datasets for all plants. 3D convex hull offered the best performance while the raster and voxel models had intermediate values. There was also an effect of spatial resolution in the raster and voxel models. The bigger the pixel size, the higher the plant volume value ( Table 2).

Total Biomass Estimation
Relationships between volume estimation methods and AGB using OLS regression analyses are summarized in Table 3. Overall, volumetric methods derived from the TLS point cloud improved AGB estimation over 250 gr (~12% of total biomass) in comparison with statistical models that included field measurements and the allometric equation based on the hemisphere model. From the different TLS methods, the voxel approach at 5 cm of resolution yielded the best fit (Adjust-R 2 = 0.86, RMSE= 927.3 g, rRMSE = 46.28%, MBE = −25.4 g), but the variations in RMSE between the most and the least accurate TLS approach, considering the best resolution for voxel and raster models, were very low (~30 g; Figure 3). The regression fits obtained from voxel methods at different resolutions ranged from 0.66 to 0.86, with RMSE varying between 927.3 g (46.28%) and 1433 g (71.52%). Raster-based volume estimation showed the worst results of all TLS derived techniques, with the best estimate obtained at 2 cm resolution, but it was the least sensitive approach to spatial resolution (Table 3). Relationships between observed and predicted AGB (g) according to the best spatial resolution for voxel (5 cm) and raster (2 cm). XGBoost represents the results according to the best set of tunning parameters described in Supplementary Figure S1. Black line represents 1:1 line. Table 3. Linear regression equations derived by OLS between AGB (y, g) and field measurements, plant volume (m 3 ) obtained by the hemisphere volume, and plant volume obtained by TLS methods (3D convex hull, voxel, and raster-based models). Values in bold reflect the best fit for each approach. *** indicates a significant relationship at p-values < 0.001.

Equation
Adjust The gradient boosting improved AGB predictions in all cases with a reduction of RMSẼ 600 g in comparison to OLS linear regressions using the plant volume derived for TLS approaches as the AGB predictor (Adjust-R 2 = 0.93, RMSE~376.6 gr, rRMSE~19%). However, we found a slight overestimation of plant AGB at larger volumes with this technique, with an MBE value of~−60 g (see Figure 3).

Photosynthetic Biomass Estimation
Total photosynthetically active biomass showed weak relationships with all different spectral indices tested in this study (see Table 4 for details), but they were significantly related with GBF. From the different indices, NDVI was the best predictor of GBF (Adjust-R 2 = 0.58, RMSE = 6.33%, rRMSE = 35.2%). By combining the best estimate of AGB obtained from TLS models with the GBF prediction from NDVI, we obtained the final equations of GB of M. tenacissima L. plants (Table 5). OLS regression fits ranged from 0.56 to 0.65 for Adjust-R 2 with RMSE values between~150 and 167 g (rRMSE~50%), depending on the TLS approach. Nevertheless, there was a slight overestimation of GB at high values in all cases, being very similar between models (see Figure 4). In a manner similar to that observed for AGB, the gradient boosting analysis (XGBoost) using TLS plant volume and NDVI as predictors improved GB estimation in all cases (RMSE = 46.1-49.8 g, rRMSE~16% and Adjust-R 2 = 0.96-0.97) with a reduction of RMSE~120 g over OLS regressions. However, contrary to OLS regressions, XGBoost tended to underestimate GB as plant size increased (see Figure 4).  Table 5 for OLS regression analysis (up) and XGBoost analysis using plant volume derived for TLS and NDVI as GB predictors (down). Black line represents 1:1 line. Table 4. Linear regression equations derived by OLS between total GB (y, g) (g) and spectral indices derived from field spectra and (up) and between GBF (y, %) and spectral indices (down). Values in bold shown the best fit index for both estimations. x represents the predictor variable (spectral indices). All regressions were significant (*** p-values < 0.001).

Discussion
The allometric equations derived from field measurements and the volumetric models derived from TLS point clouds have been proven to be good methodologies for biomass estimation in a wide range of vegetation types and ecosystems [34,37,42,45,[82][83][84][85][86]. However, none of them have been widely employed in dryland ecosystems yet, and most existing examples are focused on woody shrubs [10,35,38,41] or small grasses [87]. In this study we demonstrate that, as already described for other vegetation types [50,51,88,89], the combination of structural information obtained from TLS point clouds and spectral information provides an accurate and non-destructive methodology for the measurement and monitoring of AGB and GB in M. tenacissima L. plants, which is one of the dominant species in the Mediterranean region.
Overall AGB estimation accuracies from the TLS-data-based models vary depending on the plant type. For example, in woody shrubs, the TLS methods correctly reproduce plant structure and leaves biomass, thus providing very accurate results (see for example [35,45]). AGB estimations in grasses are less accurate, even if TLS datasets are combined with spectral information [89]. This is because plants are smaller, are thinner, and present many spaces between them, making the application of other TLS data analysis approaches not based on the development of a 3D plant model necessary [90][91][92]. Macrochloa tenacissima L. represents an intermediate situation, as it is a grass that forms dense tussocks of similar size to woody shrubs and mainly composed of a mixture of green and dry (usually more than 60% of the total AGB) leaves and inflorescences [54,93,94]. This specific morphology allows M. tenacissima L. to survive under the harsh conditions that characterize the arid and semi-arid ecosystems where it inhabits [56,[94][95][96][97] but also poses a challenge for the development of an accurate plant-volume based model to estimate AGB or GB. This is the main reason explaining the differences in accuracy between our results obtained by OLS regressions and other studies focused on woody shrubs. For example, whereas our best equation has an accuracy of~0.87, [35] found that TLS methodologies can be used to estimate AGB over sagebrush shrubs with accuracies between 0.92 and 0.86 depending on the volume estimation method. Our results are closer to those obtained in studies using allometric equations to estimate AGB over heterogenic mixtures of annual grasses of low size (i.e., [87]) and tussock grass species in altitude grasslands [98]. This limitation in AGB estimation accuracy due to the structural complexity of M. tenacissima L. plants that may lead to complex relationships between plant volume and biomass can be partially solved by the application of novel and complex machine learning methods (see Figure 4).
Overall, there is an improvement in accuracy when TLS approaches are compared with traditional allometric equations based on field sampling (see Table 3 for further details). This, as well as the weak differences in accuracy between the different volumetric approaches (Figure 4), implies a breakthrough in the study of grassland ecosystems, such as the one presented in this study. In addition, our models were developed by considering a wide range of plant sizes over three different ecosystems, representing the expected range of variability for this species in the Mediterranean region. This guarantees the transferability of the method to different areas and seasons and may explain the low values of MBE in AGB (values very close to 0 in all cases; Figure 3). Model transferability is especially important in M. tenacissima L. plants since their AGB and physiological status varies in space and time due to climate and local conditions. Concretely, M. tenacissima L. tussocks show a high variability in necromass density among plants due to two main reasons: on the one hand, the dead material is unequally removed by wind, runoff, or animals over the landscape [99]; on the other hand, necromass saturation might occur, leading to the dieback of the tuft centre, involving changes in plant density with not much variation in plant volume [54,100]. These processes imply that the necromass of M. tenacissima L. plants is very sensitive to landscape position and time, leading to variations between plant volume and AGB relationships.
The photosynthetic behaviour of this species, measured as the green fraction of AGB within the canopy, also changes during the year and with landscape position in response to water availability. In fact, in the southeast of Spain, M. tenacissima L. has a growing or active period during late autumn and winter, while it remains in a latent state during spring and summer [101], leading to seasonal changes in GB that control ecosystem processes such as overall CO 2 dynamics [102]. Different studies describe good relationships between vegetation indices that represent the amount of photosynthetically active radiation such as NDVI and the GB amount [20,27,[102][103][104] and dynamics [17,18,28]. However, our results demonstrate that M. tenacissima L. GB shows weak relationships with all the vegetation indices tested at the plant scale (Table 4). This is probably due to morphological adaptations of M. tenacissima L. plants. As previously explained, the M. tenacissima L. tussock is composed of a heterogeneous mixing of green leaves and necromass [93,94] with complex structure that includes different inclination levels and the overlay of different biomass layers that promote shadows on the photosynthetic parts in order to reduce water losses (see [56] for further information about M. tenacissima L. ecology). As spectral measurements only reflect the spectral signal from the upper part of the plant canopy, they cannot correctly represent the heterogeneous organization between dry and green leaves that characterize M. tenacissima L. plants, being strongly affected by occlusion due to the leaves' geometry in the tufts.
This limitation can be solved by combining structural information describing plant volume obtained from TLS point clouds with the spectral response of the plant (Table 5). For this purpose, NDVI has been found to be the most suitable index (Table 4), as also shown in previous studies for other grasses at the plot level [89], M. tenacissima L. high steppes in Algeria [27], or even dryland vegetation in general over coarse scales [20,105]. Overall, the RMSE of GB estimates by OLS regressions ranged from~150 g to~167 g (50-55.7%) depending on the methodology used for plant volume estimation (Figure 4), with a slight overestimation of GB as plant volume increased, demonstrating the difficulty of characterizing the photosynthetic activity of water-limited ecosystems even at the plant scale. As observed for AGB models, XGBoost analysis using plant volume and NDVI as predictors increased the accuracy in GB predictions in all cases, with RMSE values about 50 g (15.38-16.60%). However, the application of XGBoost analysis also leads to a bias in GB estimations, tending to underestimate the GB of large plants (Figure 4). Despite the observed small bias in GB prediction over large plants (a slight overestimation in OLS regressions and the opposite pattern in XGBoost models), our approach offers a fast, accurate, and non-destructive methodology for the estimation of biomass dynamics in M. tenacissima L. ecosystems. This improvement in comparison with classic spectral indices represents key progress for monitoring issues, especially due to the high variability of green biomass over space and time in these Mediterranean dryland ecosystems, being a key step forward in their monitoring and the analysis of their response to ongoing climate change.
In addition, the use of volumetric approaches for biomass estimation supposes important progress for the measurement of M. tenacissima L. biomass and dynamics, providing a promising tool for calibration and validation of further studies aimed at developing new methodologies for AGB estimation at ecosystem and regional scales. The good accuracy of raster methods (both OLS and XGBoost) and their low sensitivity to spatial resolution (see Table 3) especially make them appropriate to be adapted to other remote sensing data sources obtained from aerial platforms and unmanned aerial vehicle (UAVs) with zenithal geometry of data acquisition [106,107]. For example, structural information obtained by LiDAR (ALS) or derived by photogrammetry in UAV images could be combined with optical sensors such as SENTINEL-2 or multispectral cameras to obtain GB at the ecosystem scale, as already done in other ecosystems. These techniques have been demonstrated to be powerful tools for biomass estimation in many plant types [43,44,50,51,82,[106][107][108][109][110][111] and could be useful for the monitoring of M. tenacissima L. ecosystems, where its patchy distribution leads to easy plant delimitation. In particular, UAV images offer very detailed information that can be used for generating vegetation maps in a wide type of vegetation species [112][113][114][115][116]. Furthermore, our results provide an easy way to determine M. tenacissima L. plant biomass, a fact especially important over Mediterranean tussocks steppes where this species could represent more than 90% of AGB [58]. Plant biomass monitoring would also be particularly useful for evaluating changes in ecosystem services provided by M. tenacissima ecosystems, such as carbon storage, water regulation, and the biodiversity reservoir over the Mediterranean basin, which is one of the most threatened regions by climate change [117][118][119].

Conclusions
Photosynthetic and non-photosynthetic M. tenacissima L. biomass can be accurately estimated by combining terrestrial laser scanning and spectral information. 3D point clouds are able to correctly characterize the complex morphology of M. tenacissima L., with higher accuracy than traditional methods based on field measurements of plant structure. From the different tested approaches for analysing 3D point clouds, the voxel approach at 5 cm of resolution provided the best results for AGB estimation, but differences between it and the 3D convex hull and raster models were minimal. AGB predictions based on TLS plant volume increased their accuracy when biomass models were fitted by using state-of-the-art machine learning methodologies, such as XGBoost. This is interesting, as raster models can easily be developed from data acquired from aerial platforms, thus increasing the opportunity to be easily upscaled in further studies.
Integrating TLS volumetric information with spectral measurements of plant canopy allowed for accurately estimating the green biomass of individual plants. Thus, our methodology can be used not only for AGB estimation but also for identifying both senescent and green material, NDVI being the best multispectral index for that estimation in M. tenacissima L. plants. As NDVI can easily be obtained from low resolution satellite images or by using multispectral UAV cameras, our results represent a steep forward to unlock the complex interrelationship between green fraction and biomass in Mediterranean ecosystems. Further investigation is required to refine these techniques and to develop coarse scale approaches based on information obtained from aerial platforms. The results obtained in this study can be used as an initial step in calibration and validation processes for this issue. Funding: This research was funded by the RH2O-ARID (P18-RT-5130) project founded by the Junta de Andalucia with European Union funds for regional development, the REBIOARID (RTI2018-101921-B-I00) project funded by the FEDER/Science and Innovation Ministry-National Research Agency through the Spanish National Plan for Research and the European Union, in-cluding European Funds for Regional Development; and the project "Laboratorio de espectranomica: nueva herramienta para estudiar la diversidad y funcionamiento de zonas áridas (EQC2019-006653-P)" founded by Programa Estatal de I+D+I from Ministry of Economy and Competitiveness (Grants for infra-structures and scientific-technical equipment). Field sampling was conducted within the framework of the project "Cuantificación de flujos de carbono y agua en zonas áridas a partir de información spectral" founded by Aerial Platforms for Research-ICTS (INTA) through the campaign "Investigaciones de Altura." B.R.L was funded by the FPU predoctoral fellowship from the Educational, Culture, and Sports Ministry of Spain (FPU17/01886). E.R.C. was supported by the EMERGIA program (EMERGIA20_00337) from the General Secretariat of Universities, Research, and Technology of the Council of Economic Transformation, Industry, Knowledge, and Universities.
Institutional Review Board Statement: Not applicable.

Informed Consent Statement: Not applicable.
Data Availability Statement: Data available on request from the authors.