Aboveground Biomass Estimation in Short Rotation Forest Plantations in Northern Greece Using ESA’s Sentinel Medium-High Resolution Multispectral and Radar Imaging Missions

: Plantations of fast-growing forest species such as black locust ( Robinia Pseudoacacia ) can contribute to energy transformation, mitigate industrial pollution, and restore degraded, marginal land. In this study, the synergistic use of Sentinel-2 and Sentinel-1 time series data is explored for modeling aboveground biomass (AGB) in black locust short-rotation plantations in northeastern Greece. Optimal modeling dates and EO sensor data are also identiﬁed through the analysis. Random forest (RF) models were originally developed using monthly Sentinel-2 spectral indices, while, progressively, monthly Sentinel-1 bands were incorporated in the statistical analysis. The highest accuracy was observed for the models generated using Sentinel-2 August composites ( R 2 = 0.52). The inclusion of Sentinel-1 bands in the spectral indices’ models had a negligible effect on modeling accuracy during the leaf-on period. The correlation and comparative performance of the spectral indices in terms of pairwise correlation with AGB varied among the phenophases of the forest plantations. Overall, the ﬁeld-measured AGB in the forest plantations plots presented a higher correlation with the optical Sentinel-2 images. The synergy of Sentinel-1 and Sentinel-2 data proved to be a non-efﬁcient approach for improving forest biomass RF models throughout the year within the geographical and environmental context of our study.


Introduction
Europe's Green Deal strategy sets out the pathway for transforming the European Union (EU) towards climate neutrality by 2050. Energy transition is the cornerstone of the strategy, along with the reduction of greenhouse gas emissions, decarbonization of the energy sector, and neutralization of greenhouse gases released into the atmosphere through carbon sequestration. As part of the European Green Deal, Renewable Energy Directive II sets rules for the EU to achieve a minimum 32% share of renewables in final energy consumption by 2030 [1]. To reach such an energy target, the deployment and use of renewable energy resources, such as the bioenergy derived from organic material, have to be promoted and upgraded. Forest species have long been recognized for their potential contribution in increasing carbon sequestration and carbon storage as well as serving as an alternative bioenergy source to emission-intensive fossil fuels [2,3]. Forest plantations, despite their relatively small extent, comprising 7% of forest areas globally [4], constitute carbon sinks, add landscape diversity, support habitat biodiversity, and have a In the present study, EO passive and active data blending using a random forest (RF) modeling procedure is explored for estimating AGB over short-rotation black locust (Robinia Pseudoacacia) plantations over afforested agricultural fields in northeastern Greece. In addition, the effect of monthly variations in the canopy and leaf characteristics of black locust in AGB on overall modeling accuracy is also assessed. This objective is accomplished through (i) the analysis of the relationship of individual spectral indices and backscattering coefficients with AGB throughout the annual phenological stages and (ii) the development of monthly RF models using EO data from the Sentinel-1 and Sentinel-2 satellite duos.

Study Area
The short-rotation black locust plantations under study are located in northeastern Greece in the northern part of the Evros regional unit (Figure 1). Black locust plantations have covered approx. 44% of the total agricultural afforested area in Greece since the introduction of national Regulations 2080/92 and 1257/99 under the implementation framework of EU CAP [24]. Black locust is a fast-growing species that produces large amounts of biomass and high-density wood; it can also conserve soil and water, fixes soil nitrogen, and presents strong adaptability on low fertility soils, making it ideal for marginal land exploitation [25]. the synergetic use of passive and active data for AGB estimation in short-rotation plantations has not been previously explored in such an environmental (i.e., afforested agricultural land) research framework (i.e., optimal modeling dates).
In the present study, EO passive and active data blending using a random forest (RF) modeling procedure is explored for estimating AGB over short-rotation black locust (Robinia Pseudoacacia) plantations over afforested agricultural fields in northeastern Greece. In addition, the effect of monthly variations in the canopy and leaf characteristics of black locust in AGB on overall modeling accuracy is also assessed. This objective is accomplished through (i) the analysis of the relationship of individual spectral indices and backscattering coefficients with AGB throughout the annual phenological stages and (ii) the development of monthly RF models using EO data from the Sentinel-1 and Sentinel-2 satellite duos.

Study Area
The short-rotation black locust plantations under study are located in northeastern Greece in the northern part of the Evros regional unit (Figure 1). Black locust plantations have covered approx. 44% of the total agricultural afforested area in Greece since the introduction of national Regulations 2080/92 and 1257/99 under the implementation framework of EU CAP [24]. Black locust is a fast-growing species that produces large amounts of biomass and high-density wood; it can also conserve soil and water, fixes soil nitrogen, and presents strong adaptability on low fertility soils, making it ideal for marginal land exploitation [25]. Land use in the northern part of the Evros regional unit is dominated by agricultural activities due to the mild topography, fertile soil, and available water resources. Elevation in the area where the plantations have been established ranges from 15 to 455 m above Land use in the northern part of the Evros regional unit is dominated by agricultural activities due to the mild topography, fertile soil, and available water resources. Elevation in the area where the plantations have been established ranges from 15 to 455 m above mean sea level (m.s.l.), while the mean elevation is 99.19 m. It is noteworthy that approximately 80% of the study area can be described as lowland, according to Dikau's classification [26], with mild slopes (Figure 1). Moreover, according to Beck et al. [27], the climate of the area that contains the surveyed plots belongs mostly to the Csa (hot-summer Mediterranean climate) type, consistent with Koppen-Geiger's climate classification. Records of the meteorological station located at Alexandroupoli indicate that the months with the highest rainfall are April (53 mm), November (50.4 mm), December (48 mm), and May (47 mm), while, on the other hand, the months with the lowest rainfall are July (0 mm), September (1.2 mm), and March (6 mm).

Analysis Workflow
The methodological approach of the present study for modeling the AGB of black locust consists of five main steps ( Figure 2). In the first step, data are collected by performing field measurements on survey plots and by filtering the Sentinel-2 and Sentinel-1 image collections to create seasonal time-series. In the second step, the initial datasets undergo preprocessing by utilizing an allometric equation for AGB calculation from field measurements and by generating monthly composites from Sentinel-2 and Sentinel-1 filtered imagery. The third step involves the calculation of multispectral indices along with the value extraction from both multispectral indices and radar backscatter (VV and VH) for each survey plot. In the fourth step, Spearman rank correlation analysis is applied to the multispectral and radar variables in order to explore the correlation between each variable and the ABG for each month. Finally, in the last step, random forest models are developed using four different datasets, spectral variables, and the synergetic use of Sentinel-2 and Sentinel-1 variables. The flowchart diagram in Figure 2 presents a general description of the processing steps of our research. mean sea level (m.s.l.), while the mean elevation is 99.19 m. It is noteworthy that approximately 80% of the study area can be described as lowland, according to Dikau's classification [26], with mild slopes (Figure 1). Moreover, according to Beck et al. [27], the climate of the area that contains the surveyed plots belongs mostly to the Csa (hot-summer Mediterranean climate) type, consistent with Koppen-Geiger's climate classification. Records of the meteorological station located at Alexandroupoli indicate that the months with the highest rainfall are April (53 mm), November (50.4 mm), December (48 mm), and May (47 mm), while, on the other hand, the months with the lowest rainfall are July (0 mm), September (1.2 mm), and March (6 mm).

Analysis Workflow
The methodological approach of the present study for modeling the AGB of black locust consists of five main steps ( Figure 2). In the first step, data are collected by performing field measurements on survey plots and by filtering the Sentinel-2 and Sentinel-1 image collections to create seasonal time-series. In the second step, the initial datasets undergo preprocessing by utilizing an allometric equation for AGB calculation from field measurements and by generating monthly composites from Sentinel-2 and Sentinel-1 filtered imagery. The third step involves the calculation of multispectral indices along with the value extraction from both multispectral indices and radar backscatter (VV and VH) for each survey plot. In the fourth step, Spearman rank correlation analysis is applied to the multispectral and radar variables in order to explore the correlation between each variable and the ABG for each month. Finally, in the last step, random forest models are developed using four different datasets, spectral variables, and the synergetic use of Sentinel-2 and Sentinel-1 variables. The flowchart diagram in Figure 2 presents a general description of the processing steps of our research.

In Situ Measurements
Field data collection in the forest plantations of the study area was carried out from June to September of 2019. The parcel boundaries of the forest plantations in the whole area were manually digitized through visual interpretation of very high-resolution satellite basemaps. Since the Forest Service is responsible for the monitoring of the plantations, according to the requirements of national regulations, hard-copy records from the local forest service department were retrieved in order to populate the attributes (i.e., date of establishment, silvicultural treatments implemented, eligibility, species) of the plantation polygons through a laborious procedure. A stratified random sampling procedure was followed, considering the date of planting and silvicultural treatments for selecting the 105 square plots. The 0.09 ha plots were established within each parcel with the aid of the Global Navigation Satellite System ( Figure 1). Furthermore, in each plot, the number of trees, their diameters at breast height (d 1,3 ), and heights (H) were measured ( Table 1). The characteristics of the black locust trees are described in Table 1  AGB for individual trees was estimated by utilizing the allometric equation developed by Annighöfer et al. [28] for northern Italy. The non-linear model (Equations (1) and (2)) is based on the relationship between AGB and the diameter at breast height (DBH) and tree height and allows the calculation of the total aboveground biomass, as follows: where BM (Mg ha −1 ) represents the total dry weight of the biomass, α = −3.33 and b = 0.95 represent the fitted parameters, CF = 1.004 represents the correction factor when using the unlogged versions of the formula, D represents the diameter at breast height, and, H represents the height of the tree.

Remote Sensing Data and Preprocessing
Data acquired from the Sentinel-2 MultiSpectral Instrument (MSI) between January 2019 and December 2020 were used in this study. Clouds, cirrus, cloud shadows, and snow were masked from the Sentinel-2 surface reflectance product (L2A) within the Google Earth Engine (GEE) environment using the scene classification map resulting from the surface reflectance retrieval procedure [29]. In total, 400 Sentinel-2 tiles (granules TMF and TMG) were used. Sentinel-2 monthly composites with average values were created within the GEE. The lower spatial resolution (20 m) bands covering the red-edge and infrared part of the spectrum (i.e., 5, 6, 7, 8A, 11, and 12) were resampled to a 10-m spatial resolution. In order to compress the spectral information of Sentinel-2, to limit data redundancy and spectral information, 10 spectral indices ( Table 2) were calculated for each monthly composite. The selection of spectral indices was based on their merits for AGB estimation, identified in earlier studies [20,[30][31][32].
A Sentinel-1 dataset was also retrieved and processed by the GEE, consisting of 29 ground range detected (GRD) C-band SAR images of interferometric wide (IW) swath mode and VV and VH dual-polarization acquired between January 2020 and December 2020 ( Figure 3). The images were of descending pass, and all of them belonged to the same relative orbit (i.e., 109). Regarding preprocessing, each 10-m scene was initially calibrated into σ 0 backscatter coefficient values and despeckled using a 5 × 5 Lee filter, while terrain correction involved the use of SRTM 30 [30]. The final step involved the conversion of pixel values into decibels (dB), according to the following Equation (3): Forests 2021, 12, x FOR PEER REVIEW 6 of 17 relative orbit (i.e., 109). Regarding preprocessing, each 10-m scene was initially calibrated into σ 0 backscatter coefficient values and despeckled using a 5 × 5 Lee filter, while terrain correction involved the use of SRTM 30 [30]. The final step involved the conversion of pixel values into decibels (dB), according to the following Equation (3): Sentinel-1 time series backscattering coefficients for VH and VV polarizations at 10 m resolution were also generated using a monthly average composition method. VV and VH variables were chosen because they are the most useful and widely used variables in the literature; they involve the estimation of AGB with the aid of C-band SAR imagery [19,20,23,30]. Specifically, VH polarization is generally preferred due to the fact that it is less influenced by soil moisture [20], and VV returns are linked to crown biomass [31]. Moreover, VV polarization, although influenced by soil moisture and surface roughness, when combined with VH polarization, can possibly reduce the influence of surface scattering and foliar water content [33], thus making the S1 signal more sensitive to AGB [19]. Table 2. Sentinel-2-and Sentinel-1-based features evaluated as black locust AGB predictors.

Index
Formula Reference  Sentinel-1 time series backscattering coefficients for VH and VV polarizations at 10 m resolution were also generated using a monthly average composition method. VV and VH variables were chosen because they are the most useful and widely used variables in the literature; they involve the estimation of AGB with the aid of C-band SAR imagery [19,20,23,30]. Specifically, VH polarization is generally preferred due to the fact that it is less influenced by soil moisture [20], and VV returns are linked to crown biomass [31]. Moreover, VV polarization, although influenced by soil moisture and surface roughness, when combined with VH polarization, can possibly reduce the influence of surface scattering and foliar water content [33], thus making the S1 signal more sensitive to AGB [19]. Table 2. Sentinel-2-and Sentinel-1-based features evaluated as black locust AGB predictors.

Index
Formula Reference

Random Forest Modeling and Assessment
Initially, only pixels located completely within the plot boundary were used to calculate the mean plot values. Plot values were extracted for each survey plot from backscatter coefficients and spectral indices, and then they were used in the statistical analysis. The relationship among spectral indices/backscattering coefficients and AGB throughout a year was explored through a pairwise correlation (Spearman rank correlation) analysis. This correlation approach does not require assumptions of the distribution of the data being normal or linear [40].
RF analysis was performed using the randomForest package [45]. Concerning RF tuning parameters, the number of trees in the forest (ntree) was optimized based on the out-of-bag (OOB) estimated error; as for the minimum number of observations in a node (nodesize) and the number of random variables used in each tree (mtry), they were set to default values [45]. Based on this optimization procedure, ntree was set to 500, nodesize was set to 5, and mtry was set to 1/3 of the predictor variables.
RF models were fitted using all the data, while model performance was accessed using the OOB sample. RF as a bagging technique considers a random bootstrap sample of two-thirds of available data to train each tree; the remaining one-third of the original data, called the out-of-bag-sample, is used to calculate the mean prediction error (out-ofbag error) get unbiased out-of-bag estimates [41]. In other words, self-testing is possible even if all data are used for training since the RF model can be fitted into one sequence, with cross-validation being performed along the way [46]. A 10-fold cross-validation was performed along the way [46], using two-thirds of the available training data to grow each tree, and the remaining one-third was used to calculate the OOB error [47]. Finally, the calculations of the coefficient of determination (R 2 ) and the root mean square error (RMSE) based on OOB were used as indicators for the model's accuracy and the verification of the best dataset for AGB estimation.

Results
From June to October, during the phenophases of flowering to defoliation [48], all spectral indices presented positive, statistically significant correlations with black locust biomass. Figure 4 illustrates the average temporal NDVI and NDWI profiles of black locust plantations (2019 and 2020). The symmetrical seasonal change demonstrated from the NDVI values records the opening of the buds and mid-spring leaf development, with peak foliage development occurring in late May-early June. Fall leaf coloring begins in approximately mid-September, and defoliation begins in October, followed by leaf-off in November [48]. The leaf yellowing in summer (starting in July) is related to drought stress and high temperatures and is also noted in the NDVI decrease over this period. These findings, related to black locust phenology, are also confirmed by NDWI seasonal variation ( Figure 4) since positive values have been linked to healthy vegetation [33].
The strongest relationship between spectral indices and AGB is noted in July except for WET and MCARI indices, which present a higher correlation in August. In all months except May and January, RSR values present statistically significant correlations with measured black locust biomass throughout the year. findings, related to black locust phenology, are also confirmed by NDWI seas variation (Figure 4) since positive values have been linked to healthy vegetation [33] The strongest relationship between spectral indices and AGB is noted in July ex for WET and MCARI indices, which present a higher correlation in August. In all mo except May and January, RSR values present statistically significant correlations measured black locust biomass throughout the year. In January, none of the indices are correlated with AGB (p > 0.05), while in Febr and March, only the RSP index presents statistically significant correlations (p < 0 Furthermore, in April, four indices (NDVI, NDI45 PSSRa, and RSR) present a correla with p < 0.05. Additionally, in May, further to RSR, NDI45 is not correlated with A while in the same month, three indices (MCARI, GEMI, and WET) present their wea correlation. In November and December, a few indices present statistical, albeit w correlations with AGB.
Spearman's correlation analysis ( Table 3) also indicated that Sentin backscattering coefficients of the VV band in March, May, July, and October statistically correlated (p < 0.01) with AGB. The March composite presents the maxim correlation (R = 0.35). On the contrary, the measured VV polarization backscatte coefficients do not present a significant correlation in January, April, June, August, November. Table 3. Spearman's rank correlation coefficients (*: 0.01 < p < 0.05; **: p < 0.01) for Sentinel-2 spectral indices with black locust AGB measured in the 105 plots.  In January, none of the indices are correlated with AGB (p > 0.05), while in February and March, only the RSP index presents statistically significant correlations (p < 0.01). Furthermore, in April, four indices (NDVI, NDI45 PSSRa, and RSR) present a correlation, with p < 0.05. Additionally, in May, further to RSR, NDI45 is not correlated with AGB, while in the same month, three indices (MCARI, GEMI, and WET) present their weakest correlation. In November and December, a few indices present statistical, albeit weak, correlations with AGB.

Month
Spearman's correlation analysis ( Table 3) also indicated that Sentinel-1 backscattering coefficients of the VV band in March, May, July, and October are statistically correlated (p < 0.01) with AGB. The March composite presents the maximum correlation (R = 0.35). On the contrary, the measured VV polarization backscattering coefficients do not present a significant correlation in January, April, June, August, and November. Correlation analysis between field-measured AGB and monthly VH values indicates a significant correlation for all months except June and October ( Table 4). The highest correlation is noted in January (R = 0.47). During March, May, and July, both VV and VH backscattering coefficients provide high significant correlations (p < 0.01) with AGB.
Subsequently, four classes of RF models were trained and evaluated. RF models were developed using solely monthly spectral index values (S2) coupled with VV (S2/S1_VV), VH (S2/S1_VH), and both VV and VH backscattering coefficients (S2/S1) ( Table 5). For the models consisting only of Sentinel-2 data, the highest R 2 (0.52) and lowest RMSE (11.26 Mg ha −1 ) are noted when August spectral indices ( Figure 5) are used (R 2 = 0.379, RMSE = 12.825 Mg ha −1 ). In July and September, models achieved nearly equal accuracy (R 2 = 0.37, RMSE = 12.87 Mg ha −1 and R 2 = 0.37, RMSE = 12.99 Mg ha −1 , respectively). On the other hand, during the leaf-off season in winter (December, January, February) and early spring (March), the Sentinel-2 based models performed poorly (R 2 below 0). The synergistic use of passive and active Sentinel data slightly decreased the accuracy compared to the (best performing) Sentinel-2 August model (R 2 = 0.50, R 2 = 0.51, and R 2 = 0.50 for S2/S1_VV, S2/S1_VH, and S2/S1, respectively). The S2/S1_VH models presented equal or slightly better performance than S2/S1_VV models. The S2/S1 models with VV or  The synergistic use of passive and active Sentinel data slightly decreased the accuracy compared to the (best performing) Sentinel-2 August model (R 2 = 0.50, R 2 = 0.51, and R 2 = 0.50 for S2/S1_VV, S2/S1_VH, and S2/S1, respectively). The S2/S1_VH models presented equal or slightly better performance than S2/S1_VV models. The S2/S1 models with VV or VH or both VV and VH values were slightly improved over the spectral models of low performance (January to April and October to December). However, S2/S1 models using May to September data provided similar levels of accuracy to the Sentinel-2 indices models, implying that backscattering information had no significant effect on spectral-only models ( Figure 6).

AGB Models Using Sentinel-1 and Sentinel-2 Data
The regular assessment of tree biomass within plantations is increasingly important, and remote sensing offers a good alternative to costly, time-consuming, and difficult-toimplement field measurements [49]. In this study, remote-sensing-based AGB models were developed over short-rotation black locust forest plantations in northeastern Greece using open-access medium-high resolution multispectral and radar data. Black locust is the most common forest species used for afforestation purposes in agricultural and marginal lands in the region as well as on a national scale in Greece [24].
The results of this analysis indicate that the complementary use of data from the Figure 6. Seasonal variation of R 2 (a) and RMSE (b) of the AGB models using only monthly spectral index values (S2), coupled with VV (S2/S1_VV), VH (S2/S1_VH), and both VV and VH backscattering coefficients (S2/S1).

AGB Models Using Sentinel-1 and Sentinel-2 Data
The regular assessment of tree biomass within plantations is increasingly important, and remote sensing offers a good alternative to costly, time-consuming, and difficult-toimplement field measurements [49]. In this study, remote-sensing-based AGB models were developed over short-rotation black locust forest plantations in northeastern Greece using open-access medium-high resolution multispectral and radar data. Black locust is the most common forest species used for afforestation purposes in agricultural and marginal lands in the region as well as on a national scale in Greece [24].
The results of this analysis indicate that the complementary use of data from the Sentinel-2 and Sentinel-1 sensors provide marginal or no improvement in AGB modeling accuracy. Backscatter intensity slightly improved the poor modeling accuracy during the leaf-off period, implying that VV and VH values derived from the leaf-off season might carry better information for AGB estimation than spectral values. On the other hand, during the leaf-on period, the integration of Sentinel-1 bands with e Sentinel-2 spectral indices in the modeling procedure slightly decreased the models' accuracy. The inclusion of noisy or irrelevant predictors may reduce the predictive accuracy of the RF models, as noted in previous research findings [50]. It seems that since SAR C-band wavelength penetration is limited in the foliage of black locust due to the presence of leaves, the backscatter intensity provided no information on top of the spectral response recorded in the Sentinel-2 indices. Earlier research exploring the synergy between the active and passive sensors employed in this study presented diverse results related mainly to the wavelength of the radar data employed, which affected its sensitivity to forest biomass. In a similar research context in West Africa, Forkuor et al. [51] also found that Sentinel-2 was a better predictor of AGB than Sentinel-1; however, their combined use increased the accuracy by 0.07 (R 2 = 0.90%). Identical improvement in modeling accuracy was noted by Li et al. [52] in China. On the other hand, Nuthammachot et al. [53] identified that R 2 increased by only 0.02 when Sentinel-1 data was fused to an already highly accurate (R 2 = 0.84%) Sentinel-2 model over a forest area in Indonesia. On the other hand, Navarro et al. [20] observed a significant decrease of accuracy for a fused Sentinel-2/Sentinel-1 model, especially over low AGB density.
The RF spectral-only monthly models during winter (December, January, February) and early spring (March) presented poor performance (negative R 2 ), indicating that the models' predictions are worse when considering the mean AGB for prediction. In line with this pattern, the results of the pairwise correlation between the individual spectral indices and black locust indicated very low, non-statistically significant correlations for the majority of indices over the same period. The reflectance recorded by optical, passive sensors during winter months over deciduous forest stands [54]-as in the case of the dense black locust plantations-is a mixture of the reflectance from the woody structural elements and the heterogeneous (spatially and temporally) vegetation cover of the plantation floor. The recorded signal is further influenced by tree-to-tree self-shading and the shadows cast onto the ground from dense branches and brown seed pods of black locust during the leaf-off period.
With regard to the low performance of the May and June models, it must be noted that during this period, a peak in vegetation greenness and canopy foliage density is also observed (Figure 4). Zhu and Liu [55], in their study, also highlighted that the use of imagery acquired during vegetation peak may result in low accuracy in AGB models due to saturation. The highest accuracy was observed for the Sentinel-2 models generated from the July, August, and September spectral indices. The earlier research of our team over natural deciduous forest stands distributed within the same climatic type region confirmed that the dry season, prior to the beginning of the senescence phase, was the optimal period for forest parameter estimation [56].
During July, eight out of ten used spectral indices presented the highest correlation with black locust AGB compared to their (statistically significant) counterparts of August. However, the August RF model attained much higher accuracy. Previous studies also demonstrated that the strength of the RF modeling framework lies in taking into account the not-so-evident interactions among features [57], as in the case of RS-based variables, which carry information from different parts of the electromagnetic spectrum.

Information Content of Individual Variables
Spectral indices have been used for describing phenological changes in short-rotation plantations [58]. While NDVI is considered effective for modeling AGB [55], in our case, the pairwise correlation did not identify this index among the ones presenting the highest correlations with black locust biomass. Our results also indicated that NDVI during peak plantation growth in May did not present the highest correlation with AGB. NDVI shows stronger correlations to AGB in July and at the beginning of the fall season in September. This weaker relationship in the peak season can be expected due to the saturation effect of dense canopies [55] and the non-linear relationship between NDVI and green biomass [59]. Similar to NDVI, MCARI and GNDVI (including information from the green part of the spectrum) were not among the top-ranked indices in terms of correlation with AGB. Barati et al. [60] also noted that the use of the green band may further decrease the ability of the index to provide information on vegetation biomass.
The reduced simple ratio (RSR) presents the highest significant correlation during the period of leaf emergence (March to April). RSR, considering spectral information from the red, infrared, and short wave infrared parts of the spectrum, has a linear correlation with leaf area and sensitivity to moisture content [61], leading to a statistically significant (albeit low) correlation with AGB.
NDI45 presents the highest, statistically significant correlation with black locust AGB during the maximum leaf senescence and leaf-fall periods (September to November). Prior research findings have also indicated that NDI45 derived from Sentinel-2 imagery is sensitive to AGB [53,62]. NDI45, a modified version of the NDVI, is estimated upon spectral response in the red-edge part of the spectrum, replacing the near-infrared reflectance of the original NDVI formula. In the study of Mariën et al. [63], information from the red-edge part of the electromagnetic spectrum was also found to be the most relevant to the chlorophyll degradation noted during the leaf senescence process in deciduous forest species.
NDWI shows a highly significant correlation to AGB in the leaf-on season. Specifically, in July, NDWI has a better relationship to AGB than all the other indices. Other studies have used NDWI for the Leaf Area Index (LAI), biomass estimation, and gross primary production [64] in a range of ecosystem types and observed the sensitivity of NDWI to canopy leaf water content.
Previous studies using single-time spectral indices for AGB estimation have quite often reported unsatisfactory accuracy results due to saturation, especially when using foliage peak imagery, demonstrating that seasonal time series of a spectral index can offer a viable alternative for increasing modeling accuracy [55]. However, despite the current availability of high spatial-temporal satellite observations, the use of a (single) spectral index time-series to identify the optimal AGB modeling point is not a trivial task. Rapid transitions in canopy phenology add further obstacles to the challenges associated with the variability in canopy-related information during leaf development and senescence by different spectral indices [65]. Seasonal spectral variations over deciduous plots are driven by complex structural and biochemical changes occurring at tree level. Quite often, significant interannual variation in these responses might exist due to climate forcing [65]. Factors such as forest floor plant cover [66] and seasonal differences in shadow influences [54,67] add further constraints to the identification of the optimal modeling date based on a single index. Robust modeling procedures (i.e., RF), exploring the use of complex interactions and non-linear relationships between several indices on a monthly basis, seem to provide a viable solution.
With regard to the dataset of Sentinel-1 images, the best results in the current study were provided by those of VH polarization, which is in agreement with previously published research findings [19,20]. Furthermore, concerning deciduous trees, Sentinel-1 images are strongly affected by seasonal phenological behavior, making the use of timeseries highly recommended in such AGB estimations [11,20]. In the case of black locust, the results indicate that the optimal time period for the use of SAR imagery in AGB estimations is from late winter to the beginning of spring, aligned with the findings of Mauya et al. [11]. This observation can be attributed to the fact that during the leaf-off period, C-band wavelength saturation is at a minimum, and part of the signal backscatters onto the trunk of the tree [19,20]. On the other hand, soil wetness during the leaf-off period introduces variability in the backscatter that influences the robustness of the modeling procedure [22]. ESA's Earth Explorer Biomass mission, planned for launch in 2021, will be carrying the first-ever P-band SAR sensor, which is expected to contribute to AGB estimation processes.

Conclusions
Biomass and carbon stock in forest plantations can significantly contribute to mitigating climate change effects and a carbon-neutral Europe, providing sustainable bioenergy for substituting fossil fuel usage. Forest plantations of fast-growing species can also restore degraded, marginal land to productive use and act as a phytoremediator to absorb heavy metals in contaminated land. Regular assessment of the quantitative and qualitative status of plantation biomass through remote sensing is important for silvicultural management, identification of infections and health monitoring, mitigation of any negative environmental impacts, and prioritization of biomass harvesting and energy planning.
In this study, Sentinel-2 and Sentinel-1 time series data were employed for estimating AGB in black locust forest plantations in northeastern Greece. Spectral indices and backscattering coefficients were employed to examine monthly changes in AGB modeling accuracy. The modeling framework was based on RF regression, while individual pairwise correlations between remote sensing variables and AGB were also explored.
The highest accuracy was observed for the Sentinel-2 models generated from July, August, and September spectral indices. In particular, the August model achieved a rather satisfactory prediction accuracy (R 2 = 0.52). However, the indices used for the development of the July and September models presented higher, statistically significant correlations with the AGB. This finding confirms the well-known complex mechanisms driving forest species reflectance, from leaf to landscape, as well as the efficiency of RF modeling to handle these complex, hidden interactions among indices containing information from different parts of the spectrum. It also underlines the caution that should be applied when employing variable selection procedures in RF models. Furthermore, the negative impact on modeling accuracy from the inclusion of the Sentinel-1 bands into the spectral indices models should also be noted. Both polarizations demonstrated statistically significant correlations with AGB in July and September. Again, these findings indicate that individual variable correlations might not be a robust approach in RF and that noisy predictors can have a negative effect on the predictive accuracy of RF.
Single-date indices acquired during peak foliage growth presented both low pairwise correlations with AGB as well as poor prediction ability when integrated into the RF model. The relatively flat terrain allowed us to study the temporal changes induced in the correlations between spectral indices and the AGB of a forest species, with minimal cloud obstruction and distortions resulting from shadowing or terrain normalization procedures. One index (i.e., RSR) presented a stable, positive correlation with AGB for many months throughout the year, while MCARI presented the lowest correlation throughout the year. The comparative performance of the spectral indices in terms of correlation with AGB varied among the phenostages and different ranges of biomass.
Overall, the high temporal and spectral resolution of the Sentinel satellites facilitated the development of monthly black locust AGB models and contributed to the identification of the optimal phenostage in terms of modeling accuracy. The results of the study indicate that such a monthly modeling approach, based only on Sentinel-2 data, can be transferred across the country for biomass mapping over black locust plantations.