Growing Stock Volume Retrieval from Single and Multi-Frequency Radar Backscatter

: While products generated at global levels provide easy access to information on forest growing stock volume (GSV), their use at regional to national levels is limited by temporal frequency, spatial resolution, or unknown local errors that may be overcome through locally calibrated products. This study assessed the need, and utility, of developing locally calibrated GSV products for the Romanian forests. To this end, we used national forest inventory (NFI) permanent sampling plots with largely concurrent SAR datasets acquired at C- and L-bands to train and validate a machine learning algorithm. Different conﬁgurations of independent variables were evaluated to assess potential synergies between C- and L-band. The results show that GSV estimation errors at C- and L-band were rather similar, relative root mean squared errors (RelRMSE) around 55% for forests averaging over 450 m 3 ha − 1 , while synergies between the two wavelengths were limited. Locally calibrated models improved GSV estimation by 14% when compared to values obtained from global datasets. However, even the locally calibrated models showed particularly large errors over low GSV intervals. Aggregating the results over larger areas considerably reduced (down to 25%) the relative estimation errors.


Introduction
Forests are among the most biodiverse terrestrial ecosystems and a key element for carbon sequestration as they store large amounts of organic matter (i.e., biomass). Therefore, forest above ground biomass (AGB) estimation is a sensitive research topic, as information on AGB levels and dynamics is needed to estimate greenhouse gases flux and thus to shape policies development, implementation, and monitoring [1]. This importance is highlighted by the countless forest inventory programs aimed at evaluating, monitoring, and reporting, among others, AGB or Growing Stock Volume (GSV) levels. Such programs use an array of data sources from in situ measurements to information acquired by earth observation (EO) platforms. In situ surveys, based on systematic sampling grids, are the backbone of traditional national forest inventory (NFI) programs sometimes stretching back centuries [2]. However, NFI programs are expensive and time consuming while not providing for a synoptic view across the landscape [3].

SAR Data Processing and Extraction
The images acquired on the same orbit were first co-registered in the radar geometry. For Sentinel-1 products, the orbital state vectors were updated using the Precise Orbit Ephemerides, with the co-registration process being carried out by relative orbit after slice assembly (i.e., concatenation) as the Sentinel-1 image frames (i.e., slices) acquired from the same orbital path are provided in slices of variable ground footprints. For the ALOS PAL-SAR-2 data, the SAR processing was carried out on a path/frame basis. Image co-registration was based on a cross-correlation algorithm [33] with the first image of the temporal The in situ samples were collected within the second cycle (2013-2018) of the Romanian NFI. The NFI was established in 2008 using a 4 × 4 km grid, although a denser 2 × 2 km grid is used in plains where forest cover is low. The NFI is a two stage (aerial photography followed by in situ surveys) continuous forest inventory with a five-year cycle that covers the entire national territory. The field surveys comprise, at the end of the five-year cycle, about 24,000 permanent sample plots. To increase field work efficiency, the measurements are realized on four permanent sampling plots (PSP) at each grid node (Figure 1b). Each PSP contains three concentric circles (7.98 m, 12.62 m, and 25 m in radius) where different forest characteristics are assessed, including forest type, tree species, diameter at breast height (DBH), height (H), lying deadwood, and ground vegetation. DBH and height are assessed in the first two circles, 7.98 m (5.6 cm ≤ DBH ≤ 28.5 cm) and 12.62 m (DBH ≥ 28.5 cm). The GSV estimation is based on species-specific (43 main species) allometry based on DBH and height. For conifers, the volume of branches (>5 cm diameter) was established through specific equations as percentage of stem volume and was added to the stem volume [32]. The volume measurements are scaled (by the circle area), aggregated, and reported per hectare. NFI estimates GSV with a sampling error of 1.79%, with target differences between measurements and control surveys less than 1%. The NFI objective is to derive statistics of forest properties over the national territory. As such, individual plot measurements are not optimized to calibrate or evaluate remote sensing products due to the relatively small sampling area.
For each PSP, the data provided included the main forest species, forest inventory date, percentage forest cover, mean diameter at breast height, and height, as well as per hectare GSV. Of the 1815 PSPs available for this study, 1153 were inventoried between 2015 and 2016, coinciding with the EO data acquisition period. The remaining plots were inventoried in 2014. As forest stands in Romania consist of homogeneous tracts (species composition, age, and DBH classes), the NFI plots are representative of larger areas were not falling on stand border, as stands range between 6-10 ha in hilly and 10-15 ha in mountainous regions. The DBH, H, and GSV over the in-situ data, reached 96 cm, 50 m, and 1577 m 3 ha −1 , respectively, depending on the sample plot ( Table 1). As small pockets of large trees influence average values, when relatively small areas are assessed, the PSP-wise GSV was high (>1000 m 3 ha −1 ) for some plots ( Figure 2).

Earth Observation Data
The NFI plots used in this study were covered within 13 ALOS PALSAR-2 frames and six ascending and descending Sentinel-1 relative orbits (Figure 1). 104 L-band ALOS PALSAR-2 dual-polarized (HH and HV polarizations) datasets were provided by JAXA

Earth Observation Data
The NFI plots used in this study were covered within 13 ALOS PALSAR-2 frames and six ascending and descending Sentinel-1 relative orbits ( Figure 1). 104 L-band ALOS PALSAR-2 dual-polarized (HH and HV polarizations) datasets were provided by JAXA (2014-2017 period) as single look complex (SLC) images, while 1034 Sentinel-1 A and B granules were downloaded from open repositories for the years 2015-2016, the period when most (83) ALOS PALSAR-2 datasets were acquired. Sentinel-1 ground range detected (GRD) dual polarized (VV and VH polarizations) C-band datasets (A and B satellites) were downloaded for the six relative orbits (29,102, and 131 ascending and 7, 80, and 109 descending passes) covering the selected study areas.

SAR Data Processing and Extraction
The images acquired on the same orbit were first co-registered in the radar geometry. For Sentinel-1 products, the orbital state vectors were updated using the Precise Orbit Ephemerides, with the co-registration process being carried out by relative orbit after slice assembly (i.e., concatenation) as the Sentinel-1 image frames (i.e., slices) acquired from the same orbital path are provided in slices of variable ground footprints. For the ALOS PALSAR-2 data, the SAR processing was carried out on a path/frame basis. Image co-registration was based on a cross-correlation algorithm [33] with the first image of the temporal data series being used as reference. Each co-registered image was multilooked to obtain a ground pixel spacing of approximately 30 m (4 × 8 pixels in range and azimuth for ALOS PALSAR-2 and 20 × 3 pixels in range and azimuth for Sentinel-1). The backscatter was converted to gamma0 using an ellipsoid-based area as reference, to account for the influence of topography [34]. The backscatter coefficient (γ•) was subsequently topographically normalized using the real scattering area, derived pixel-wise [35], from the one arc-second Shuttle Radar Topography Mission (SRTM) digital elevation model (DEM). Each image was then orthorectified to the Universal Transverse Mercator (UTM) coordinate system (zone 35 North, datum WGS84) using a look-up table that related the coordinates of each pixel in the radar geometry with the coordinates of the same pixel in the map geometry. The look-up table was generated using image orbital information and the digital elevation model [36].
At each PSP, the backscatter coefficient of all images acquired before the forest inventory date was extracted for the pixel containing the center coordinates. The backscatter values were averaged for all available dates by polarization and by relative orbit (Sentinel-1). The standard deviation of the backscatter temporal series (i.e., temporal stability) was also computed at each PSP location for every sensor, polarization, and relative orbit (Sentinel-1). Out of the 1815 available NFI samples, 704 (665 forest and 39 non-forest) were used for modelling and validation of the retrieved GSV. The reduced number of useful samples was determined by the intersection between ALOS PALSAR-2 and Sentinel-1 orbits, the correspondence between NFI inventory date and ALOS PALSAR-2 and Sentinel-1 image acquisition dates, and topography, as pixels affected by geometric distortions were masked out during SAR data processing. As both ascending and descending Sentinel-1 orbits were tested, the number of PSPs affected by layover and shadow was high.

Growing Stock Volume Retrieval
Several modeling approaches can be used to retrieve forest GSV (or AGB) from radar backscatter coefficients, including parametric and nonparametric models [3,[8][9][10]13,16,24,[37][38][39][40][41]. Within a previous study, some of these models were evaluated using the Romanian NFI and ALOS PALSAR-2 datasets [42]. The results suggested that non-parametric models provide the lowest errors and bias, regardless of polarization or forest species over the selected study area. Therefore, we used a non-parametric modelling approach, based on Random Forests [43], to assess the synergies between the C-and L-band dataset for GSV retrieval. As non-parametric models offer the opportunity to include non-linearly related variables, and have no assumptions regarding the statistical properties of the data, such models are often preferred when a sufficiently large dataset of samples is available for model parameterization. The models use ensemble learning methods to improve the overall predictive power with respect to any of the constituent models by aggregating their predictions. In random forest regression, each tree is built using a deterministic algorithm by selecting a random set of variables and a random sample from the training dataset [43]. Although RF models provide high fit statistics there are known drawbacks, including difficulties in interpreting the results, potential overfitting, and high computational demands. A total of 85 RF models were trained and evaluated by gradually increasing the independent variables, starting from single polarized to multi-polarization multi-sensor data (see Section 4 and Table A1 in Appendix A for detailed information of the predictor variables used). Each sensor was individually tested to generate a reference baseline and allow for cross-sensor comparisons. The use of additional variables (e.g., forest type, temporal backscatter variability, and local incidence angle) was also evaluated to ascertain the opportunity for GSV retrieval improvements. For each sensor, we increased the independent variables, starting with the average backscatter, temporal stability (i.e., the standard deviation of the time series, sd), local incidence angle (LIA), and the forest type (Ft). For the C-band we analyzed the ascending and descending passes separately, as well as their combination. The models were calibrated using the 704 PSP samples. In this study, TreeBagger from MATLAB ® software package (v. 2020b) was used to construct the RF classifier. The number of decision trees was set to 200 and a curvature test was used to select the best split predictor and grow unbiased trees [44]. The remaining parameters were kept as default, with surrogate splits being allowed to account for missing values in the data.

GSV Retrieval Accuracy
Model performance was evaluated at both pixel and grid levels. At pixel (i.e., plot) level, the retrieval was evaluated using repeat random sub-sampling with cross-validation for each individual sensor and polarization, as well as for dual polarized and dual frequency configurations. To reduce variability due to random sampling effects, 100 rounds were performed by randomly splitting the 704 PSPs into training (75%) and validation (25%) samples. During each round, the models were calibrated using the training samples and subsequently used to estimate the GSV for the validation samples. The observed and predicted GSV for the validation samples was accumulated over the 100 rounds and used to compute four accuracy error metrics: the root mean squared error (RMSE, Equation (1)), the relative RMSE (RelRMSE, Equation (2)), the bias (Equation (3)), and the Pearson's correlation coefficient (r).
where: P j = predicted values, O j = in situ observed values, n = number of samples, and O and P = mean values for the observed and predicted values, respectively. As GSV estimation errors are inherently large when using small NFI plots for SAR models calibration [45,46], accuracy metrics (as per Equations (1)-(3)) were also derived over predefined grids of 10 × 10, 20 × 20, and 30 × 30 km. The grid size was a compromise between having enough samples in each tile to obtain a more reliable estimate of the NFI average stem volume and having enough tiles to compute the accuracy metrics by volume intervals. As the Sentinel-1 strips covered the 1815 available PSPs (Figure 1), the analysis was focused on the GSV maps derived from C-band VH polarized data (see Section 4.2) to maximize the number of tiles containing at least 4 in situ plots (i.e., one NFI grid node). The errors were estimated for both ascending and descending satellite passes. Models using the VH backscatter and its temporal standard deviation (see Table 2, in bold) were calibrated using the reduced set (704) of NFI for cross-comparison purposes. At each tile, the average NFI GSV was compared to the average GSV for the corresponding predicted pixels.

Local vs. Global GSV Retrieval
The GSV overall accuracy estimates were compared against values derived at continental to global levels. We used the GlobBiomass dataset [20], the only global product currently available at comparable pixel spacing (100 m) and derived from comparable remote sensing datasets (C-and L-band SAR data). The GlobBiomass product provides GSV estimates for the reference year 2010 based on the combined use of Envisat ASAR (C-band) and ALOS PALSAR (L-band) sensors and ancillary information from Landsat imagery. GSV values from the first NFI cycle (2008-2012) were used as reference to match the dates of the GlobBiomass product. As NFI data collection protocol and location did not change between the first and second NFI cycles and we have used the same SAR wavelengths to derive the local product, input datasets influence on the results were minimized. Please note that Envisat ASAR and ALOS PALSAR sensors were decommissioned in 2012 and 2011, respectively while data from Sentinel-1 and ALOS PALSAR-2 were not available until 2014. The GlobBiomass GSV values were extracted at the location of the available NFI plots to compute the same error metrics.

Results and Discussions
The higher backscatter variability at L-band when compared to C-band, due to the limited number of multi-temporal datasets available is evident for all the main species in the Carpathians ( Figure 3). Additionally, notice the different GSV ranges, significantly larger for beech forests when compared to the coniferous and oak forests. The backscatter coefficient shows the specific raise with GSV up to a SAR wavelength specific saturation point. By main species, the average backscatter values were slightly higher at L-band over the coniferous forests ( Figure 3b). For display reasons, a logarithmic model was fitted to the data. The fit was very similar for both SAR wavelengths over the high biomass beech forests (Figure 3c). Different saturation points between the two wavelengths are apparent over the coniferous and oak forests (Figure 3b,d), suggesting changed sensitivity to GSV.

Pixel-Wise GSV Estimation
Over all sensors, polarizations, and predictor variables combinations, the RMSE, Rel RMSE, and r ranged between 238-278 m 3 ha −1 , 55-65%, −5-6 m 3 ha −1 , and 0.14-0.48, respectively (Tables 2 and A1). As showing results for all possible combinations of the predictor variables (19) is problematic, Table 2 shows the error metrics for a selection of models to understand their variations when expanding the number of predictor variables. Notice that adding more variables does not translate into more accurate predictions. Additional configurations were provided in Table A1 for completeness. The data is presented by SAR wavelength, multi-frequency configurations and descending C-band passes combinations.

Pixel-Wise GSV Estimation
Over all sensors, polarizations, and predictor variables combinations, the RMSE, Rel RMSE, and r ranged between 238-278 m 3 ha −1 , 55-65%, −5-6 m 3 ha −1 , and 0.14-0.48, respectively (Table 2 and Table A1). As showing results for all possible combinations of the predictor variables (19) is problematic, Table 2 shows the error metrics for a selection of models to understand their variations when expanding the number of predictor variables. Notice that adding more variables does not translate into more accurate predictions. Additional configurations were provided in Table A1 for completeness. The data is presented by SAR wavelength, multi-frequency configurations and descending C-band passes combinations.
Within the single sensor single polarization configuration, the L-band data provided similar relative RMSE errors when compared to the C-band data. The small difference was explained by the much denser C-band time series, and the use of average backscatter over 3 years prior to the forest inventory date. This allowed for a considerable reduction in the environmental induced noise (i.e., backscatter variation due to changes in soil and vegetation water content) at C-band, where over 150 images were averaged over each relative orbit. In contrast, fewer images (5-10) were available for the ALOS PALSAR-2 sensor over each acquisition frame. Improved GSV retrieval accuracy at C-band was observed when adding as an independent variable either the local incidence angle (LIA) or the forest type (FT). However, such improvements were marginal for all accuracy metrics. The RMSE, RelRMSE, and r improved, on average, by 15 m 3 ha −1 , 3%, and 0.1, respectively. Dual polarized models were slightly more accurate (2-5%) when compared to single polarized   Table 2. GSV accuracy as a function of the independent variables used for the single polarized models. Bold numbers show overall results for the models analyzed by GSV intervals in Figures 4 and 5. Abbreviations as follows: RMSE-root mean squared error, RelRMSE-relative RMSE, r-Pearson's correlation coefficient, L-stands for L-band, C-stands for C-band, H-horizontal, V-vertical, a-ascending pass, d-descending pass, sd-standard deviation, Ft-forest type, and LIA-local incidence angle.  Within the single sensor single polarization configuration, the L-band data provided similar relative RMSE errors when compared to the C-band data. The small difference was explained by the much denser C-band time series, and the use of average backscatter over 3 years prior to the forest inventory date. This allowed for a considerable reduction in the environmental induced noise (i.e., backscatter variation due to changes in soil and vegetation water content) at C-band, where over 150 images were averaged over each relative orbit. In contrast, fewer images (5-10) were available for the ALOS PALSAR-2 sensor over each acquisition frame. Improved GSV retrieval accuracy at C-band was observed when adding as an independent variable either the local incidence angle (LIA) or the forest type (FT). However, such improvements were marginal for all accuracy metrics. The RMSE, RelRMSE, and r improved, on average, by 15 m 3 ha −1 , 3%, and 0.1, respectively. Dual polarized models were slightly more accurate (2-5%) when compared to single polarized models at both C-and L-bands, while the addition of co-to cross-polarized backscatter ratios further decreased the estimation error, albeit marginally. Adding the time series standard deviation to backscatter and co-to cross-polarized ratios improved the accuracy metrics slightly, suggesting that some information can be added to that already contained in the backscatter coefficient. Similarly, adding LIA and/or FT to the models only resulted in marginal improvement, most likely being related to the increased number of variables available in the model. No synergies were observed in terms of retrieval accuracy, when combining information from both ascending and descending C-band passes with model performance not improving over the use of either pass. However, GSV estimation from both ascending and descending passes has merit over the rough Carpathian topography, as it limits the extent of 'not-observed' areas, i.e., areas with topographic distortion due to shadows and layover which are masked out during SAR processing. Further, the similar error level from both passes provides for a homogeneous GSV estimation over the entire landscape. Little synergy was also observed when simultaneously using C-and L-band for GSV retrieval. Improvements in the RMSE, RelRMSE, and r were noticed with respect to using single pass C-band data, but were marginal, i.e., 10 m 3 /ha, 3%, and 0.05, respectively. In addition, such improvements were observed when adding the LIA to the various combinations of SAR frequencies, polarizations, and passes. The information provided by LIA and Ftype was largely interchangeable in all configurations (single, dual, multi-frequency, and multi-pass), as simultaneously using both variables did not improve the results. In fact, the opposite was observed for some of the tested combinations. It seems that LIA partitions data in a similar way to Ftype as oak forests are mostly found on lesser slopes in the hilly region, with conifers being found towards the mountain tops on steep terrain. This explanation is further supported by the fact that adding the LIA from different sensors or satellite passes provided no tangible improvement of the accuracy metrics.

Independent
Error analysis by intervals (Figure 4 for models in bold in Table 2) shows much higher relative RMSE errors over forests supporting low (<100 m 3 ha −1 ) GSV values at both C-and L-bands, which is consistent with previous findings [14,47]. The high relative RMSE are related to the high signal variability due to the influence of local surface conditions (i.e., soil surface roughness and moisture) as direct scattering from the ground dominates the signal in forests supporting low GSV levels [48]. In addition, the relative RMSE metric is unstable at lower GSV values, as the denominator approaches zero the relative error approaches infinity. Over the remaining GSV intervals the relative RMSE decreased significantly at both wavelengths with the minimum values (around 20%) being observed for the 400-600 m 3 ha −1 GSV range. The estimation error was largely dominated by bias, as the coefficient of variation (CV) of the error was below one (data not shown), for most GSV intervals except the 400-600 m 3 ha −1 where the random error dominated, and the 0-100 m 3 ha −1 interval where bias and random error contribution was balanced (0.8 < CV < 0.9). Note that the bias was positive over low GSV ranges and negative over the high GSV ranges (Figure 5b), with a crossover at the GSV range with the lowest errors, 400-600 m 3 ha −1 . The standard deviation (SD) of the predicted values did not vary across GSV intervals (75-85 m 3 ha −1 ) except for the first bin, where the SD was twice as much (Figure 5a). Similar values were observed over temperate forests in Poland in previous studies [14]. Overall, the accuracy assessment shows that GSV is over-and under-estimated at low and high GSV values, respectively as the models fitting aims at the point defined by the average of the observed and predicted values [14]. ha −1 . The standard deviation (SD) of the predicted values did not vary across GSV intervals (75-85 m 3 ha −1 ) except for the first bin, where the SD was twice as much (Figure 5a). Similar values were observed over temperate forests in Poland in previous studies [14].
Overall, the accuracy assessment shows that GSV is over-and under-estimated at low and high GSV values, respectively as the models fitting aims at the point defined by the average of the observed and predicted values [14].  The relatively large errors observed for all combinations of predictor variables may be related to potential mismatches between the Earth Observation data and the in situ samples, as the coordinates at each PSP were reconstructed from the node grid coordi- Rel. RMSE (%) GSV reference (m 3 ha -1 ) d) Figure 5. GSV retrieval metrics for 10 × 10 km grid cells for single polarization C-band models. Panels (a-d) show the same information as in Figure 4. The relatively large errors observed for all combinations of predictor variables may be related to potential mismatches between the Earth Observation data and the in situ samples, as the coordinates at each PSP were reconstructed from the node grid coordinates using the specified distance and azimuth. This may introduce displacements with respect to the real location of the NFI plot where field crews were not able to precisely measure the required azimuth and distance, or the exact location of the grid node (e.g., GPS location error under dense forest canopy may reach several meters). Such positioning errors coupled with the relatively small area inventoried may induce errors particularly for border stands. For the remaining NFI samples, the positioning errors should have limited effects as Romanian forests stands, averaging between 3.5 and 15 ha from plains to mountains, have homogeneous structure and species composition.

Grid Based GSV Accuracy
Over the 10 × 10 km grid (244 cells), the RMSE and RelRMSE decreased to 158 m 3 ha −1 and 35%, respectively (average values), with marginal differences being observed between values from ascending (35.2%) and descending (33.8%) satellite passes. Averaging GSV values from ascending and descending passes did not result in improved accuracy metrics. Further aggregation through the larger 20 × 20 (125 cells) and 30 × 30 km (78 cells) grids improved RMSE (119 and 105 m 3 ha −1 , respectively) and RelRMSE (29% and 25%, respectively) values, but also increased the bias from 13 to 30 and, 35 m 3 ha −1 , respectively, as only two or three bins were available. The RelRMSE for the 30 × 30 km aggregation level was similar to that observed over the Swedish forests (21.4%) at the same scale. However, our retrieval statistics indicate some residual bias, which was not the case in Sweden, where a longer wavelength (L-band) was used for retrieval and more samples were aggregated [49].
When compared to the pixel-based assessment, similar patterns were observed over all biomass intervals except for the 0-100 m 3 ha −1 range for which no grid cell was available ( Figure 5 for models in bold in Table 2). The improvement in retrieval statistics seemd related to the absence of low GSV plots (<100 m 3 ha −1 ) where large estimation errors were observed at pixel level (Figures 4 and 5). Notice that the lower number of available cells at larger grids precluded a similar analysis for the 20 × 20 and 30 × 30 km grids.

Comparison to Global Products
Our results support previous findings which demonstrated increased retrieval accuracies for locally calibrated models [14,[23][24][25]. When compared to the GlobBiomass values ( Figure 6), locally derived estimates showed potential to improve the RelRMSE by about 15% (55.4 vs. 70.5%), decrease the RMSE by about 54 m 3 ha −1 (236 vs. 290 m 3 /ha), and increase the correlation between predicted and observe values from 0.2 to 0.5. Such improvements, although significant and at a higher (30 vs. 100 m) spatial resolution, are still far from the required accuracies needed for operational forest management which target 5-10% accuracies for GSV at stand level. Nevertheless, one should notice that aggregating the pixel wise GSV to lower spatial resolution (e.g., stand level) has the potential to significantly decrease the RMSE and RelRMSE values [13,50] as also demonstrated by the improvements observed at grid level. The incentive in using SAR based GSV estimation from locally calibrate models resides in the opportunity to better fine tune the model to the local conditions when compared to using globally available estimate. In addition, SAR-based maps provide wall to wall cover, yearly updates, and the opportunity to derive regional statistics.
nificantly decrease the RMSE and RelRMSE values [13,50] as also demonstrated by the improvements observed at grid level. The incentive in using SAR based GSV estimation from locally calibrate models resides in the opportunity to better fine tune the model to the local conditions when compared to using globally available estimate. In addition, SAR-based maps provide wall to wall cover, yearly updates, and the opportunity to derive regional statistics.

Conclusions
This study assessed the utility of C-and L-band data for GSV estimation over high growing stock volume forests in the Carpathians. The study included representative areas covering the three main forest species, oak, beech, and coniferous, making up the bulk of the Romanian forests. The in situ NFI data and the SAR imagery were acquired between 2015 and 2016. GSV retrieval was carried out at 30 m pixel size using a machine learning (Random Forests) algorithm and various configurations of the independent variables from single sensor single polarization to multi-sensor multi-polarization. The observed RelRMSE varied between 55-66% depending on the input predictor variables. For similar predictor variables, the mapping accuracy was slightly higher at L-when compared to Cband.
The availability of significantly more datasets at C-band reduced speckle as the much denser Sentinel-1 time series seemed to have compensated for the inherent limitations of a shorter wavelength. Such a finding suggests that using denser L-band time series, as those available from the soon to be launched NISAR mission, have the potential to further

Conclusions
This study assessed the utility of C-and L-band data for GSV estimation over high growing stock volume forests in the Carpathians. The study included representative areas covering the three main forest species, oak, beech, and coniferous, making up the bulk of the Romanian forests. The in situ NFI data and the SAR imagery were acquired between 2015 and 2016. GSV retrieval was carried out at 30 m pixel size using a machine learning (Random Forests) algorithm and various configurations of the independent variables from single sensor single polarization to multi-sensor multi-polarization. The observed RelRMSE varied between 55-66% depending on the input predictor variables. For similar predictor variables, the mapping accuracy was slightly higher at L-when compared to C-band.
The availability of significantly more datasets at C-band reduced speckle as the much denser Sentinel-1 time series seemed to have compensated for the inherent limitations of a shorter wavelength. Such a finding suggests that using denser L-band time series, as those available from the soon to be launched NISAR mission, have the potential to further improve GSV estimation accuracy in high volume forests. The GSV retrieval errors decreased by 5-10% when adding temporal information (i.e., backscatter coefficient temporal standard deviation) as an independent variable eliminating the need for a second polarization. This suggests a high correlation between the information of co-and crosspolarized channels, which may be explained by the overall high GSV levels and thus the relative low correlation between GSV and backscatter. Little synergies were observed when using both C-and L-bands, as well as when jointly using C-band ascending and descending passes. While the utility of the pixel level GSV values is somewhat reduced, their aggregation at grid level decreased the estimation errors (25-35%) and provided more reliable estimates which may be useful for regional and national assessments. Future work should focus on reducing GSV errors by fusing the Lidar data acquired by NASA's Global Ecosystem Dynamics Investigation (GEDI) mission to obtain a GSV gridded product that provides more accurate estimates for stock levels at regional to national levels. In addition, dense L-band temporal series from the future NISAR mission should be tested to ascertain the potential to improve the locally derived GSV estimates, as datasets from the P-band BIOMASS mission may not be available over the Carpathians.
As for the limitations of this study, the most important is related to the design of NFIs, which are not optimized for calibrating and validating remote sensing products [49], the large differences between the amount of data collected by the two SAR sensors, which limited a like-for-like comparison, and the available DEM (SRTM DEM) used for terrain normalization as a more precise DEM (e.g., Lidar based or Tandem-X DEM) allows for improved scattering area estimation reducing the effect of topography on the backscatter and thus improving the retrieval of the target biophysical characteristic [34,51]. Table A1. GSV accuracy as a function of the selected independent variables for a range of configurations. Abbreviations as follows: RMSE-root mean squared error, RelRMSE-relative RMSE, r-Pearson's correlation coefficient, L-stands for L-band, C-stands for C-band, H-horizontal, V-vertical, a-ascending pass, d-descending pass, sd-standard deviation, Ft-forest type, and LIA-local incidence angle.