Carbon Stock Estimations in a Mediterranean Riparian Forest: A Case Study Combining Field Data and UAV Imagery

: This study aims to estimate the total biomass aboveground and soil carbon stocks in a Mediterranean riparian forest and identify the contribution of the di ﬀ erent species and ecosystem compartments to the overall riparian carbon reservoir. We used a combined ﬁeld and object-based image analysis (OBIA) approach, based on unmanned aerial vehicle (UAV) multispectral imagery, to assess C stock of three dominant riparian species. A linear discriminator was designed, based on a set of spectral variables previously selected in an optimal way, permitting the classiﬁcation of the species corresponding to every object in the study area. This made it possible to estimate the area occupied by each species and its contribution to the tree aboveground biomass (AGB). Three uncertainty levels were considered, related to the trade-o ﬀ between the number of unclassiﬁed and misclassiﬁed objects, leading to an error control associated with the estimated tree AGB. We found that riparian woodlands dominated by Acacia dealbata Link showed the highest average carbon stock per unit area (251 ± 90 tC ha − 1 ) followed by Alnus glutinosa (L.) Gaertner (162 ± 12 tC ha − 1 ) and by Salix salviifolia Brot. (73 ± 17 tC ha − 1 ), which are mainly related to the stem density, vegetation development and successional stage of the di ﬀ erent stands. The woody tree compartment showed the highest inputs (79%), followed by the understory vegetation (12%) and lastly by the soil mineral layer (9%). Spectral vegetation indices developed to suppress saturation e ﬀ ects were consistently selected as important variables for species classiﬁcation. The total tree AGB in the study area varies from 734 to 1053 tC according to the distinct levels of uncertainty. This study provided the foundations for the assessment of the riparian carbon sequestration and the economic value of the carbon stocks provided by similar Mediterranean riparian forests, a highly relevant ecosystem service for the regulation of climate change e ﬀ ects.


Introduction
Riparian forests and floodplain areas have been referred to as important carbon sinks, and thus are crucial systems to mitigate the effects of climate change and to provide a regulation ecosystem service [1][2][3][4]. Riparian forests as subjected to frequent floods and resulting geomorphologic processes, exhibiting high spatial and temporal variability with consequent high potential for long-term C storage. Patterns of sediment and nutrient changes combined with distinct levels of water availability promote Forests 2020, 11 rapid changes in species composition, distribution, and density. These factors can stimulate higher rates of biomass accumulation in riparian zones when compared with terrestrial forest systems [1,[3][4][5]. Carbon storage and sequestration vary along environmental gradients and for different riparian species [4,6]. Carbon reservoirs in riparian forests seem to be positively correlated with annual precipitation and negatively correlated with maximum temperature [4,7]. However, no clear patterns concerning the effects of climate and geological setting in C sequestration have yet emerged from the available literature, as most of the studies are on temperate regions [3,4,6]. In addition to climate, several local factors may influence the carbon storage capacity of riparian ecosystems like the floodplain width, the gradient of inundation, the valley geometry, the channel complexity and the flow regime [3,4,8].
Wide floodplains with complex channel forms and strong lateral connectivity promoting complex vegetation structures are expected to have high carbon stock per unit area [4]. Additionally, C stocks in the different compartments seem to be ruled by different processes, although manifold drivers remain uncertain [3,9]. While the carbon stock in aboveground woody biomass and roots appears to be driven by site-specific factors, such as distance to the groundwater table, stand age, vegetation structure, successional stage and species composition, the carbon stock in the understory litterfall and soil are additionally sensitive to short-term changes in environmental variables such as mean annual temperature and precipitation and hydrological disturbance [3]. A detailed review concerning the physical factors that influence the retention of organic carbon in floodplains and riparian systems can be found in Sutfin et al. [4].
The most important gap concerning C stocks of riparian forests is the lack of knowledge across diverse climates and related vegetation types. Research on riparian carbon stock is very scarce and most studies came from South and North America, with limited studies in Europe [1][2][3][10][11][12], and particularly for the Mediterranean areas. As far as we know, only Cabezas and Comín [11] addressed the carbon stock capacity of a Mediterranean riparian forest, located in the Ebro River (NE Spain) and only at the topsoil level and in the river sediments. So, quantifying the carbon stocks and sequestration in these ecosystems may help to define the baseline conditions which can be further used to support restoration designs or for addressing the impact of climate change on C stocks in riparian ecosystems [9]. Field data are required to quantify the carbon stocks of riparian ecosystems, and the relative contribution of the different compartments across distinct ecoregions.
Field-based estimations can be combined with GIS techniques and remote sensing approaches to obtain large-scale maps of vegetation, soil and total carbon stock distributions. Although in a very limited number, these spatial modelling approaches have been successfully applied in the USA [13,14] and Mexico [15], but also in Europe [12,16,17]. The most common methodologies rely on calibrating satellite measurements (spectral, textural and geometric variables) to in situ estimates of aboveground biomass (AGB) at field study plots. AGB is usually determined using a combination of allometric equations that establish the relationships among field plot inventory data (stem diameter, tree height, tree density,) and AGB. These allometric equations are, in turn, developed from trees that have been dissected, oven-dried and weighed [18] and mostly directed for commercial woody species [19,20]. Thus, field surveys are an urgent need to develop appropriate allometric biomass models across a wide range of riparian ecosystems, which may be important to assess the global significance of riparian systems on the overall carbon cycle.
There is a vast literature concerning the use of remote sensing approaches to estimate AGB for terrestrial forest ecosystems [18]. However, research on AGB and carbon stock in riparian and floodplain forests is extremely limited [13,14]. The information is even scarcer for riparian systems located in the Mediterranean region, where estimates of AGB and carbon sequestration are extremely needed. Multispectral and hyperspectral data with high spatial resolution have been successfully used in floodplain areas to perform remote-sensing AGB-retrieval models, using in situ AGB as training and validation data [13][14][15][16]21]. However, in riparian habitats, and especially in Mediterranean regions, very high spatial resolution is mandatory to map vegetation due to the limited width and high structural complexity of riparian vegetation patches [22][23][24]. Unmanned aerial vehicle (UAV) imagery Forests 2020, 11, 376 3 of 21 may overcome this issue since the sub-decimeter spatial resolution and high spatial accuracy may be achieved with these sensors [25]. Although rare, some studies have demonstrated the potential of UAV imagery to finely describe and monitor riparian systems thought the classification and mapping of riparian species, but also by assessing the vegetation health condition, at a local scale [26,27] and using multi-temporal approaches [28]. Recent studies have also concluded that AGB can be accurately estimated using tree height and crown diameter variables derived from airborne high spatial resolution images [19] or LIDAR data. However, this "tree-centric" remote sensing approach cannot be applied for riparian forests, since it is not possible to individualize trees or measure the vertical structure of the riparian forest accurately, in the amalgam of vegetation that characterizes Mediterranean riparian systems. Thus, in this study we applied an object-based image analysis (OBIA) approach to a UAV multispectral imagery, to identify, not individual trees, but clusters of similar neighbouring pixels from the same riparian species, hereafter referred to as riparian objects (RO). In very high spatial resolution imagery, riparian objects are composed of several adjacent pixels. Object-based image analysis (OBIA) initially aggregates these adjacent pixels into spectrally homogeneous regions-the objects-so that the variance within objects is less than the variance between objects [24,29,30].
This study primarily aims to estimate the carbon stock of Mediterranean riparian forest and assess the potentiality of remote sensing methods for C stocks assessments in riparian ecosystems. For this, we used in situ data and the spatial distribution of the tree AGB is obtained using a remote sensing OBIA approach applied to the UAV-acquired imagery. Four underlying questions were addressed: (i) what are the differences in the aboveground C stock on tree biomass of the dominant species (i.e., Alnus glutinosa (L.) Gaertner, Salix salviifolia Brot.and Acacia dealbata Link )?; 9ii) What is the contribution of the different ecosystem compartments (tree + understory + soil) to the overall riparian carbon reservoir? (iii) What are the key spectral variables most suitable for the remote sensing classification of these riparian trees' species? (iv) What is the confidence level of the assessment of the riparian C reservoir when using remote sensing data?
The results are expected to contribute to understanding and evaluating the importance of riparian forests as providers of a climate regulation ecosystem service and to support the management and restoration of these multifunctional ecotones.

Study Area and Riparian Typologies
The study was conducted in a 3km-long stretch of River Alcolobre, a tributary of Tagus River, located in centre-west Portugal, between 39 • 24 39" N, 8 • 15 41" W and 39 • 23 44"N, 8 • 15 07" W ( Figure 1, Video 1-DOI:10.5281/zenodo.3667672). The studied stretch is spread over calcareous Mesozoic formations in a typical Mediterranean climate characterized by mild winters and hot dry summers and frequently has irregular interannual fluctuations of precipitation [31]. Flood peaks usually occur in early winter, followed by a slow decline in flow and a consequent drying during late spring and summer. The riparian area has been managed and protected since 1983 by The Navigator Company, a pulp and paper company, under responsible and sustainable forestry practices in compliance with recognized international forest management certification systems.
The study area was selected due to the presence of a black alder (Alnus glutinosa) woodland with high conservation status (priority habitat 91E0-Alluvial forests with Alnus glutinosa according to the Annex I of Directive 92/43/EEC) and the simultaneous occurrence of other types of riparian woodland. The latter includes an almost monospecific forest composed by an Iberian endemic shrub -Salix salviifolia (borrazeira-branca)-and an area highly invaded by silver wattle (Acacia dealbata), an invasive species native to southeastern Australia.
Specifically, the study area, at its northern part, is composed of a mature and well-preserved alder woodland with elderberries (Sambucus nigra L.) on the edge, grey-willows (Salix atrocinerea Brot.), black poplars (Populus nigra L.), alder buckthorn (Frangula alnus Miller) and common hawthorn (Crataegus monogyna Jacq.) in the understory. The central area is composed of a dense forest of old mimosa trees with alders mixed with narrowed-leaved ash (Fraxinus angustifolia Vahl), and elm-leaf blackberry (Rubus ulmifolius Schott) on the understory. In the southern area, there is a young but dense willow woodland of borrazeira-branca with small width as it has been constrained by umbrella pine (Pinus pinea L.) plantations. The study area was selected due to the presence of a black alder (Alnus glutinosa) woodland with high conservation status (priority habitat 91E0-Alluvial forests with Alnus glutinosa according to the Annex I of Directive 92/43/EEC) and the simultaneous occurrence of other types of riparian woodland. The latter includes an almost monospecific forest composed by an Iberian endemic shrub -Salix salviifolia (borrazeira-branca)-and an area highly invaded by silver wattle (Acacia dealbata), an invasive species native to southeastern Australia.
Specifically, the study area, at its northern part, is composed of a mature and well-preserved alder woodland with elderberries (Sambucus nigra L.) on the edge, grey-willows (Salix atrocinerea Brot.), black poplars (Populus nigra L.), alder buckthorn (Frangula alnus Miller) and common hawthorn (Crataegus monogyna Jacq.) in the understory. The central area is composed of a dense forest of old mimosa trees with alders mixed with narrowed-leaved ash (Fraxinus angustifolia Vahl), and elm-leaf blackberry (Rubus ulmifolius Schott) on the understory. In the southern area, there is a young but dense willow woodland of borrazeira-branca with small width as it has been constrained by umbrella pine (Pinus pinea L.) plantations.

Field Data Collection
We collected field data on 14 sampling plots (SP) across the riparian forest to estimate in situ the total aboveground biomass (AGB), in the understory component and the organic and mineral soil

Field Data Collection
We collected field data on 14 sampling plots (SP) across the riparian forest to estimate in situ the total aboveground biomass (AGB), in the understory component and the organic and mineral soil layer, for the three dominant riparian woodlands: alder, willow, and acacia woodlands.

Tree Inventory
The tree inventory was conducted in sampling plots (SP) of 314 m 2 (10m of radius) using a stratified random design to represent each riparian woodland type ( Figure 1). Four sampling plots were established with a radius of 5m (SP10 and SP11) or using a semi-circular design (SP08, SP13) due to the small width of the riparian corridor. Sampling plots were separated by at least 50m to avoid spatial autocorrelation in the data. The central point of each sampling plot was geo-referenced by applying a post-processing differential correction to the data recorded through a global positioning system (GPS) with an estimated accuracy of less than 1m (Figure 2a). Within each SP, and to estimate the total biomass in the aboveground component, we measured diameter at breast height (DBH), total height (H), and identified all trees at the species level, with DBH ≥ 7cm following the methodological approach of Ravindranath and Ostwald [32] (Figure 2b). The position of each tree inside each plot was identified by recording the azimuth and the distance to the central point of the sampling plot. We considered pure and dominated plots for each riparian species with more than 75% of the basal area tree cover. The number of stems, as well as the basal area (calculated as the sum of the DBH cross-sectional area), of all the trees in the plot were reported to the hectare unit area.

Allometric Models for Biomass Quantification
To obtain the in situ AGB for each riparian woodland type, we used two methods: 1) speciesspecific biomass allometric models applied to the trees inventoried in each SP and 2) the biomass harvesting method for the construction of a non-linear biomass model for the species where no allometric models exist: the endemic Salix salviifolia. Regarding the first method, we carried out a literature review of all the species-specific allometric models and organized the database by dendrometric ranges of DBH and H of the trees used to build the models. The models selected were the ones that fitted the DBH and H ranges of the trees under study and in analogous climatic conditions. Regarding the second methodology, and only for the endemic willow (Salix salviifolia), where no allometric equation was found in the literature with the referred combination criteria, we develop a new allometric model using a destructive approach. For that, 15 stems were harvested from 5 to 30 cm DBH and the total aboveground fresh weight measured in the field with a digital scale (Figure 2c and 2d). A sample of approximately 300g was brought to the laboratory for dry weight calculations and used to convert fresh to dry weight biomass. A simple non-linear regression model was developed by relating the total aboveground biomass (in kg) with DBH (cm) and/or H (m). The best model used only DBH: (1)

Understory Biomass
To estimate the biomass in the understory component, we collected all green material (shrubs and herbaceous vegetation) within 50 × 50 cm square plots placed at five locations in each sampling plot: the central point and in the extreme N, S, W and E coordinates ( Figure 2e). All the material was kept in identified bags, taken to the laboratory and placed in the oven for 48 h at 60 • C until constant weight. The understory biomass was converted to the respective plot unit area and to carbon stocks by applying a conversion factor of 0.47 [5,33].

Organic and Mineral Soil Layer
For the quantification of the biomass in the organic soil layer, we collected all the partially decomposed material at the soil surface within the 50 × 50cm square plots in the center of each sampling plot. All the material was collected and placed in identified bags. The samples were taken to the laboratory and placed in the oven for 48 h at 70 • C for dry weight calculation and the carbon stock calculated by multiplying the biomass dry weight by a conversion factor of 0.47. For the quantification of the carbon in the mineral soil layer, we used a soil probe and randomly collected five soil cores up to 20 cm within the plot area ( Figure 2f). We also collected the respective soil bulk density to convert the organic material in carbon stock per unit area using the methodology described in FAO [34]. Field data were collected between 4 May and 16 July 2018.

Allometric Models for Biomass Quantification
To obtain the in situ AGB for each riparian woodland type, we used two methods: (1) species-specific biomass allometric models applied to the trees inventoried in each SP and (2) the biomass harvesting method for the construction of a non-linear biomass model for the species where no allometric models exist: the endemic Salix salviifolia. Regarding the first method, we carried out a literature review of all the species-specific allometric models and organized the database by dendrometric ranges of DBH and H of the trees used to build the models. The models selected were the ones that fitted the DBH and H ranges of the trees under study and in analogous climatic conditions. Regarding the second methodology, and only for the endemic willow (Salix salviifolia), where no allometric equation was found in the literature with the referred combination criteria, we develop a new allometric model using a destructive approach. For that, 15 stems were harvested from 5 to 30 cm DBH and the total aboveground fresh weight measured in the field with a digital scale (Figure 2c,d). A sample of approximately 300g was brought to the laboratory for dry weight calculations and used to convert fresh to dry weight biomass. A simple non-linear regression model was developed by relating the total aboveground biomass (in kg) with DBH (cm) and/or H (m). The best model used only DBH: The biomass equations collected and developed were then applied to each of the stems measured in each plot and the total species-specific biomass calculated on an area basis.
We used the conversion factor of 0.5 for the conversion of biomass into carbon and the calculations were reported on an area basis (tC ha −1 ) [33,35].

UAV Imagery Acquisition and Processing
The study area was surveyed on 11th August 2018 using a fixed-wing UAV SenseFly eBee platform, operated by Terradrone Co. The UAV weighs 700g with payload and battery and has a maximum operational flight time of approximately 45 min, depending on wind speed and direction. The flight control system on-board is a SenseFly autopilot. Installed in a laptop, the software eMotion 3 plans and controls the flight in real-time, using a telemetry radio link up to 3 km. The flight mission was prepared to cover the entire study area using parallel flight lines with an overlap of 80% (both lateral and longitudinal). The altitude of each flight line was around 120 m from the terrain heights, i.e., referenced to the terrain topography to obtain a ground distance as constant as possible.
The RGB sensor was used to produce a true-colour composite orthomosaic, while the MSP sensor was used to produce reflectance maps with spectral information for each band. Since the UAV platform can only operate with one sensor at a time, and due to battery and weather limitations, four flights were conducted to cover the entire study area. A total of 520 images were acquired with the RGB sensor, at an altitude of 120 m, while 2854 images were acquired with the MSP sensor at an altitude of 100 m. With these flight settings, the acquired multispectral images had a spatial resolution of 10 cm while the RGB orthophoto was obtained with 3 cm of spatial resolution. An 80%-80% of forwarding and side overlap between neighbouring images was set, taking into consideration the terrain orography and the vegetation characteristics (vegetation density and land use). This is a critical aspect of mission planning. A reduced overlap between neighbouring images may inhibit the process of creating an accurate image mosaic by producing zones without radiometric information [36,37].
The flights were conducted under sunny and cloudless conditions, with a wind speed of 5 m/s, between 11:00 AM and 02:00 PM. The MSP camera has an upward-facing irradiance sensor that captures information for lighting conditions and sun angle at the exact time of capture of each image. The irradiance sensor can measure the arbitrary incoming irradiance for each image frame so that the irradiance measurement can be used for radiometric normalisation. Image quality evaluation was performed by visual examination and mainly involves the detection of vignetting effects and the presence of blurry images and oblique scenes. All images were included in the dataset since none of these problems was detected. The Pix4DMapper software (Pix4D S.A., Lucerne, Switzerland) was used for the radiometric correction in order to convert sensor radiance into normalised surface reflectance values for each band [38,39]. The algorithm uses the illumination parameters such as sun angle and sun irradiation collected before each flight-using a radiometric target-and the irradiance measurements to generate a 'surface reflectance map', which will allow the comparison of spectral measurements and derive calibrated vegetation indexes for the entire study area.
As for the geometric correction, 12 GPS RTK points were deployed as ground control points (GCPs) across the entire study area, using the real-time network correction based on the nearest base station (Entroncamento), located within 20 km of the study site. The coordinates and heights of the points were collected and subsequently used in the photogrammetric process to enhance horizontal and vertical precision. The geo-referencing root-mean-square error (RMSE) of the dataset was 0.16 m. The overall projection error was 1.27 pixels in all datasets.

Image Segmentation
An object-based approach was used to identify clusters of similar neighbouring pixels from the same riparian species, hereafter referred to as riparian objects (RO). The object-based approach is particularly useful for species classification using very high spatial resolution imagery, when compared with the pixel-based approach, since it promotes the statistical separability between classes, by initially aggregating adjacent pixels into spectrally, texturally and geometrically homogeneous regions-the objects [29,40]. The object-oriented classification includes two steps: image segmentation and object classification. The segmentation process should produce semantically and ecologically meaningful objects, which means having a translation in the real-world [29]. In this study, segmentation should ideally produce polygons that correspond to individual trees, in order to estimate AGB using variables that can be measured remotely, such as the crown diameter [19]. However, in our riparian forest, it was not possible to individualize trees due to its high spatial variability, vegetation density and strata complexity. Thus, to remotely estimate the total tree AGB and the respective spatial distribution along the study area, we first segmented the UAV image and classified the objects in the three dominant riparian classes. Then, we applied the in situ carbon stocks estimations, per unit canopy tree cover, for each class, to the entire riparian forest.
We started by digitalized manually, in the RGB-geometrically corrected imagery, the riparian patches for the three dominant classes: Alnus glutinosa (RO_AlnGlu), Acacia dealbata (RO_AcaDea) and Salix salviifolia (RO_SalSal). The geographic location of the species patch is known from the field surveys. We termed this dataset as RO_Manual and used it as ground truth data (Figure 3). In order to obtain the total canopy cover area for each species, we summed the areas of the RO_Manual for corresponding species. The remote tree average AGB estimate per unit area was obtained for each dominant species by summing the AGB values obtained with the field data (described in 2.

and 2.3)
Forests 2020, 11, 376 8 of 21 in all the RO_Manual and dividing the value by the total canopy cover area. It should be noted that the remote estimation of the tree average AGB was reported to the unit area of canopy cover, which is an image-based measure. However, the soil unit area was used for the characterization of the in situ carbon stocks of the riparian forest (described in 2.2). Then, the "Mean-shift" algorithm [41,42] was used in the segmentation phase using the Segmentation toolbox (ESRI ArcGIS10.5). The optimizable parameters in the "Mean-shift" algorithm include the spatial detail, the spectral detail and the minimum segment size. The parameters were established following the criteria of previous studies [24,43,44] by using a systematic trial and error Then, the "Mean-shift" algorithm [41,42] was used in the segmentation phase using the Segmentation toolbox (ESRI ArcGIS10.5). The optimizable parameters in the "Mean-shift" algorithm include the spatial detail, the spectral detail and the minimum segment size. The parameters were established following the criteria of previous studies [24,43,44] by using a systematic trial and error approach, validated by visual image inspection, which means for this case, how well the objects match with the manually digitalized patches, i.e., the RO_Manual (Figure 3). In our study, we selected a spatial detail of 16, a spectral detail of 5 and the minimum segment size of 20. The segmentation algorithm was then applied to two image-composite bands, namely Green-Red-NIR and Green-Red-RedEdge. These two datasets were termed as GRNIR and GRRedEdge, respectively. Image-composite bands with NIR and RedEdge in simultaneous were not considered since they showed high redundancy (correlation coefficient > 0.95).

Classification and Accuracy Assessment
Linear discriminant analysis (LDA) was used to separate the riparian objects according to the corresponding species. The linear separator permits us to classify any riparian object based on its mean spectral variables. Before the separation process, we identified the best set of spectral variables able to perform species classification.
The selected variables describe the spectral reflectance characteristics of each riparian object class (Table 1). We calculated the mean and standard deviation of each UAV band and tested several well-referenced Vegetation Indices (VI), which assesses the greenness level of the objects and its associated primary productivity and vegetation growth [21,[45][46][47][48]. We did not select any spectral indices developed to minimize the effect of soil or the atmospheric effects in the reflectance of the vegetation, since all the study area showed a large vegetation cover and the images were collected from a low altitude platform. Training objects were selected to represent about 0.1% of the total objects, for each class, i.e., 31 RO_AlnGlu, 34 RO_AcaDea and 32 RO_SalSal. The full dataset of training objects was composed of three matrices (RO_GRNIR, RO_GRRedEdge, and RO_RGB_Manual) with eleven spectral variables in each one.
Correlation coefficients (r) between all pairs of variables were initially used to avoid redundancy in the data. RedEdgeCI, RedEdgeNDVI and GNDVI were eliminated since they show a high correlation (r>0.99) with RedEdgeMSR and NDVI, respectively, in both RO_GRNIR, RO_GRRedEdge and RO_RGB_Manual datasets.
An optimization algorithm, that implements an adaptation of the Leaps and Bounds Algorithm, was then applied to each matrix, returning the best variable subset of a given cardinality that is optimal for linear discrimination. The computations were performed by the function eleaps from R package subselect (v.014; Cerdeira et al. [49,50]), using 'ccr12' as criterion. The resulting subset of spectral variables was then used to devise a linear discriminator for each UAV image-composite band, namely GRNIR and GRRedEdge. Since three main riparian species were considered in the training set, the discriminant is a 2-dimensional space. For each UAV image-composite band, the accuracy of the obtained discriminator was achieved through 100 bootstrap replicates, where, in each replicate, 99 random points were selected with replication as training data and the remaining points (± 32 to 41 points) were used to test. The accuracy was measured by the percentage of misclassified objects, summarized in a confusion matrix. All computations related to LDA were performed by the functions lda and predict from R package MASS (v.7.3-51.5; Ripley et. al., [51]). The use of bootstrap to access accuracy is implemented in caret R package (v. 6.0-84; Kuhn M, [51]).
A bivariate normal density was then fitted to the training data points, in the 2-dimension discriminant space, corresponding to each species. We used the ellipsis associated to the regions of probabilities p = 95%, 99%, and 99.9% to classify each object as follows: if the object falls in an ellipse of probability p, it is classified in one of the three riparian classes, AlnGlu, AcaDea or SalSal; if it falls outside every ellipses, it is considered as "unclassified". The "unclassified" objects include the "canopy gap" and other riparian species and vegetation-related classes found in the gallery, but to a lesser extent.
Combining the area correspondent to the classified objects, with the AGB per unit area of canopy cover, for each of the three species, we obtained the remote estimations of the total amount of tree AGB in the riparian forest, for the three main riparian species. Repeating the classification process for each probability level, it was possible to get an accuracy interval for the remote tree AGB estimation. We additionally characterize the riparian objects classified in the three species using landscape metrics, namely the mean patch size (MPS), the number of objects (Nump), and the percentage of class area (CA) relative to the total study area, using the three probability levels (95%, 99%, and 99.9%) [52].

Carbon Stocks in Situ
The field results reveal that the compartment with the highest contribution to the total carbon stock in the riparian forest is the woody tree compartment (79%), followed by the understory vegetation (herbaceous, shrubs and partially decomposed litterfall) with 12%, and lastly the soil mineral layer with 9% ( Figure 4).
The plots dominated by Acacia dealbata showed the highest total carbon stock (tree + understory + soil), averaging 251 ± 90 tC ha −1 (mean ± standard error) (Figure 4). The diameter at breast height in these plots averages 11 ± 3 cm and the average number of trees is 6102 ± 3682 trees/ha. The carbon fraction in the aboveground tree compartment is 192 ± 84tC ha −1 , followed by the understory, with 45 ± 10 tC ha −1 , and the mineral soil, with 15 ± 1t tC ha −1 (Figure 4).
Pure and dominated plots of Alnus glutinosa showed an average DBH of 26 ± 2 cm and an average tree density of 743 ± 165 trees/ha. The total carbon stock (tree + understory + soil) averages 162 ± 12 tC ha −1 , having the highest fraction in the aboveground tree biomass component (140 ± 13 tC ha −1 ). The understory vegetation and litterfall accounts for 9 ± 3 tC ha −1 and the mineral soil layer for 13 ± 1 tC ha −1 .
Pure and dominated plots of Alnus glutinosa showed an average DBH of 26 ± 2 cm and an average tree density of 743 ± 165 trees/ha. The total carbon stock (tree + understory + soil) averages 162 ± 12 tC ha −1 , having the highest fraction in the aboveground tree biomass component (140 ± 13 tC ha −1 ). The understory vegetation and litterfall accounts for 9 ± 3 tC ha −1 and the mineral soil layer for 13 ± 1 tC ha −1 . The plots dominated by Salix salviifolia are quite heterogeneous regarding density, displaying an average of 6436 ± 2793 stem/ha and a DBH of 7 ± 0 cm. The total carbon stock (tree + understory + soil) averages 73 ± 17tC ha −1 . The aboveground tree compartment has the lowest average among the three species, with 48 ± 18 tC ha −1 , followed by the soil, with 13 ± 1 tC ha −1 , and the understory, with 12 ± 1 tC ha −1 .
Two plots with a mixture of species, other than the abovementioned, were also measured. In addition to Acacia dealbata, Alnus glutinosa and Salix salviifolia, species like Fraxinus angustifolia, Populus nigra and Salix atrocinerea were also inventoried. These species DBH range was 13 ± 7, 26 ± 6 and 17 ± 6 cm and occupying an area of 10%, 12% and 72%, respectively. The biomass was not estimated for these species as they represented only a small fraction of the riparian forest.

OBIA Species Classification and Tree AGB Maps
The discriminant confusion matrix showed high and similar average accuracy with the three band combinations: GRNIR image-composite band (99.5%) followed by the RGB_Manual (99.0%) and the GRRedEdge (98.3%) ( Table 2).  The plots dominated by Salix salviifolia are quite heterogeneous regarding density, displaying an average of 6436 ± 2793 stem/ha and a DBH of 7 ± 0 cm. The total carbon stock (tree + understory + soil) averages 73 ± 17tC ha −1 . The aboveground tree compartment has the lowest average among the three species, with 48 ± 18 tC ha −1 , followed by the soil, with 13 ± 1 tC ha −1 , and the understory, with 12 ± 1 tC ha −1 .
Two plots with a mixture of species, other than the abovementioned, were also measured. In addition to Acacia dealbata, Alnus glutinosa and Salix salviifolia, species like Fraxinus angustifolia, Populus nigra and Salix atrocinerea were also inventoried. These species DBH range was 13 ± 7, 26 ± 6 and 17 ± 6 cm and occupying an area of 10%, 12% and 72%, respectively. The biomass was not estimated for these species as they represented only a small fraction of the riparian forest.

OBIA Species Classification and Tree AGB Maps
The discriminant confusion matrix showed high and similar average accuracy with the three band combinations: GRNIR image-composite band (99.5%) followed by the RGB_Manual (99.0%) and the GRRedEdge (98.3%) ( Table 2). Misclassifications mainly occurred between Alnus glutinosa and Acacia dealbata objects (Table 2, Figure 5) and reached a maximum of 1.4% in the GRRedEdge image-composite band. The efficient complete search algorithm identified some differences in the selection of the best set of variables, according to the band combination (File S2). MSR, RedEdgeMSR and GCI and the average values of the red and NIR bands were typically considered as significant variables to distinguish the distinct RO classes. On the contrary, the NDVI band was consistently avoided in the best sub-set models. Nevertheless, and since the GRNIR image achieved the highest classification accuracy (99.5) (Table 2), NIR, GCI, RedEdgeMSR, Red and MSR variables were selected for the subsequent species classification.
Misclassifications mainly occurred between Alnus glutinosa and Acacia dealbata objects (Table 2, Figure 5) and reached a maximum of 1.4% in the GRRedEdge image-composite band. The efficient complete search algorithm identified some differences in the selection of the best set of variables, according to the band combination (File S2). MSR, RedEdgeMSR and GCI and the average values of the red and NIR bands were typically considered as significant variables to distinguish the distinct RO classes. On the contrary, the NDVI band was consistently avoided in the best sub-set models. Nevertheless, and since the GRNIR image achieved the highest classification accuracy (99.5) (Table  2), NIR, GCI, RedEdgeMSR, Red and MSR variables were selected for the subsequent species classification. The regions corresponding to each species in the space of linear discriminant are shown in File S3, while Figure 6 shows the respective OBIA classification maps for an exemplification area dominated by each of the studied woodlands. As expected, the objects located in the northern part of the riparian forest were mainly classified as Alnus glutinosa, although some Acacia dealbata and Salix The regions corresponding to each species in the space of linear discriminant are shown in File S3, while Figure 6 shows the respective OBIA classification maps for an exemplification area dominated by each of the studied woodlands. As expected, the objects located in the northern part of the riparian forest were mainly classified as Alnus glutinosa, although some Acacia dealbata and Salix salviifolia stands were recognized (Figure 6a). In the central area of the riparian forest, the majority of the objects were identified as Acacia dealbata, while some Salix salviifolia stands and large Alnus glutinosa patches were are also predicted (Figure 6b). The southern part was mainly mapped as Salix salviifolia, although some intermingled small Acacia dealbata and Alnus glutinosa patches were forecasted (Figure 6c). The unclassified objects were mainly related to other species-vegetated patches, bare soil areas mostly located in the external limits of the study area, and to low albedo features, like shadows, typically located inside the gallery.
Concerning the tree cover area, Acacia dealbata showed the highest percentage of predicted occupation in the riparian forest (CA), ranging from 22.3% to 32.1%, followed by Alnus glutinosa (CA from 20.9% to 29.8%) and by Salix salviifolia (CA from 9.9% to 14.2%) ( Table 3). Salix salviifolia objects showed the young character of the woodland area with the smallest mean patch size (MPS from 3.0 m 2 ± 0.2 to 2.3m 2 ± 0.1) (mean ± standard error), while Alnus glutinosa and Acacia dealbata exhibit larger riparian patches (MPS from 4.3 ± 0.4 to 3.8 ± 0.3 m 2 and MPS from 4.2 ± 0.3 to 3.7 ± 0.2 m 2 , respectively) ( Table 3). Although the majority of the unclassified objects represent small patches (MPS from 1.9 ± 0.1 to 1.5 ± 0.0 m 2 ) the predicted cover area occupies 46.9% to 23.9%, according to the probability levels of 95% and 99.9%, respectively. Table 3. Number of objects (Nump), mean patch size (MPS) ± standard error and percentage of class area (CA) of the classified riparian objects in the three dominant species (Acacia dealbata, Alnus glutinosa and Salix salviifolia) and the unclassified riparian objects (Unclass) relative to the total study area, using the three regions of classification probability (95%, 99% and 99.9%).  salviifolia stands were recognized (Figure 6a). In the central area of the riparian forest, the majority of the objects were identified as Acacia dealbata, while some Salix salviifolia stands and large Alnus glutinosa patches were are also predicted (Figure 6b). The southern part was mainly mapped as Salix salviifolia, although some intermingled small Acacia dealbata and Alnus glutinosa patches were forecasted (Figure 6c). The unclassified objects were mainly related to other species-vegetated patches, bare soil areas mostly located in the external limits of the study area, and to low albedo features, like shadows, typically located inside the gallery.

Riparian Carbon Stocks
The in situ average estimation of the total aboveground, understory and soil carbon stocks in our Mediterranean riparian forest revealed values ranging from 122.4 to 201.2 tC ha −1 . Acacia dealbata show the highest average for the total carbon stocks (250.7 ± 89.5 tC ha −1 ), followed by Alnus glutinosa (161.6 ± 11.8 tC ha −1 ) and Salix salviifolia (73.1 ± 16.9 tC ha −1 ). Carbon reservoirs increase with forest age, successional stages and vegetation development [15]. The contrasting carbon stocks estimates found for each riparian species mainly reflect the capacity for tree growth and the biomass accumulation in the aboveground woody component, which is substantially different between Alnus glutinosa and Salix salviifolia. Acacia dealbata showed a high capacity to grow under high inter-and intraspecific competition levels, showing a high number of stems with considerable dimensions per unit area. Higher values for total C stocks (soil and aboveground carbon stocks) were reported in a study conducted on the Danubian floodplains [1]. In the referred study, the highest total C stocks (474 tC ha −1 ) were observed in mature riparian hardwoods dominated by Fraxinus excelsior L., Acer campestre L., Alnus incana (L.) Moench, Carpinus betulus L. and Quercus robur L., followed by cottonwoods dominated by Populus sp. (403 tC ha −1 ). Softwoods dominated by Salix alba L. and Acer negundo L. attained 356 tC ha −1 , while meadows and reeds achieved 212 tC ha −1 . Nevertheless, Cartisano et al. [2] estimated lower values (88 tC ha −1 ) of aboveground woody biomass in poplar-dominated riparian forests located in Paglia River (central Italy). Nevertheless, the reported values for total riparian carbon stocks, including trees, shrubs, herbaceous, litterfall and soil are highly variable in the available literature. In a recent synthesis developed mostly in temperate riparian forests located in North and South America, Dybala et al. [7] estimated that riparian forests stored an average of 68-158 tC ha −1 of biomass at maturity. In our study, we found a high stock variability, especially in the tree component (average range from 48.4 to 191.5 tC ha −1 ) and in the understory compartment (average range from 8.6 to 44.6 tC ha −1 ) and less variability in the soil component (average range from 12.9 to 14.6 tC ha −1 ). Nevertheless, those values were mentioned as competing with the highest estimates from other terrestrial forest ecosystems in Portugal, ranging values up to 150 tC ha −1 in stone pine (Pinus pinea) stands with varying ages and tree densities [53], 33 tC ha −1 in a cork oak (Quercus suber L.) forests (with 177 tree ha −1 and 50 years old) [54] and approximately 200 tC ha −1 in a eucalyptus (Eucalyptus globulus Labill.) stands [55].
Human disturbance, like land-use and land-cover changes, streamflow regulation and channel alterations may also influence carbon storage dynamics in riparian ecosystems. Worldwide patterns of deforestation and agricultural land-use intensification are likely to cause significant carbon emissions and a reduction in the amount of carbon stored by riparian and floodplain areas [9,56,57]. In a study conducted on riparian meadows of Sierra Nevada, California, USA, the authors observed that degraded riparian systems, due to logging and livestock grazing, stored a lesser amount of carbon in soils when compared with undisturbed riparian meadows [58]. In our study area, although we observed a well-preserved riparian forest, with an interesting dynamic of natural regeneration, especially in the northern part, all the study area is surrounded by an intensive forest production system, mainly composed by Eucalyptus and Pinus plantations (Video 1-DOI:10.5281/zenodo.3667672). Land-use and land-cover changes may promote changes in the transport of sediments and nutrients causing alterations on the organic carbon fluxes along the floodplain. Thus, higher values of carbon stocks could be obtained in more near-natural riparian systems, although these pristine ecosystems are relatively rare in Mediterranean systems. Additionally, the presence of old-stage successional Acacia dealbata stands in our study area reflects riparian degradation patterns, resulting from the direct long-standing human disturbance within the river system [59]. Although these stands achieved the highest carbon stocks, Acacia dealbata is an aggressive alien invasive species in Mediterranean streams, reducing biodiversity, causing water depletion and habitat loss and thus its control should be promoted [59].
Concerning the relative contribution of the different ecosystem compartments to the overall carbon reservoir we found that, like in other studies conducted in riparian and floodplain areas, the woody tree compartment shows the highest inputs (79%), followed by the understory vegetation (12%) and lastly by the soil mineral layer (9%). Nevertheless, carbon stocks and the relative importance of the distinct carbon pools are driven by large-scale climate and geological variables and also by local-scale factors like valley geometry, channel complexity, soil type, hydrological connectivity and microbial activity [1,3,4]. In our study, the higher fraction of biomass accumulated under Acacia dealbata reflects the higher productivity of the species, which influences leaf fall turnover rates and its accumulation in the understory. In Mediterranean regions, carbon reservoirs are highly influenced by the strong seasonality of precipitation and flooding, which affects productivity, local transport, and redistribution of organic matter. Soils in Mediterranean riparian areas are seasonally flooded and subjected to strong water stream events. Soil structure and all the surface soil layer itself is lost, and thus the soil carbon stock is likely highly spatially and temporally variable, and usually lower than in temperate regions [3,4]. Also, the soil texture, characterized by coarse gravel with little or no aggregation elements, severely reduces the capacity to retain any carbon.

AGB Remote-Sensing Approach
This study represents a "combine and assign remote-sensing" approach sensu Goetz et al. [18] to estimate and mapping the carbon storage of a Mediterranean riparian forest. In the proposed methodology, in situ estimates of AGB are associated with distinct vegetation classes, derived from a remote-sensing approach, to provide a finer-grained spatial distribution of tree AGB in the study area. Different AGB weights are assigned according to the riparian species and a detailed tree AGB map was produced with the geographic location of the carbon stock values across the riparian forest. These maps can be used for management proposes for the identification of conservation areas for long-term carbon reservoirs or to prioritize zones to improve or restore carbon stock accumulations. In addition, it can be used to assess the economic value of the carbon ecosystem service provided by the riparian forest. Carbon storage can be seen as a co-benefit of riparian restoration, in addition to biodiversity improvement and conservation [5,20].
Despite the advantages, the proposed approach is not able, for a given species, to discriminate AGB gradients according to the successional stage, DBH or stem density variability. A unique average AGB per unit area value is assigned to a vegetation class and assumed to be representative for a given spatial unit, i.e., the riparian object. Moreover, spectral variables were used to priorly discriminate and classify the riparian species, but were not able to directly discriminate AGB variability between the riparian species. In highly dense canopy cover conditions, which characterize our study area, similar VI spectral behaviour and saturation effects may explain the low performance of vegetation indices to discriminate AGB gradients [21,60]. Although highly variable AGB could be found in mature and advanced successional forests, as observed in our riparian system, the complicated vegetation stand structures often result in similar reflectance behaviour, generating poor correlation between biomass and spectral signatures [60].
The majority of the "combine and assign remote-sensing" studies in floodplain areas made use of allometric equations relating field measurements with remote-sensing derived variables to estimate AGB but varies in the way that those field data are combined. Usually, AGB is modelled for spectrally dissimilar land-use classes or habitat types which represent distinct biomass levels (e.g., high water-related herbaceous zone, vs low water-related woody zone) by relating in situ values, at the field plot level, to remote-sensing derived variables at the pixel level [14]. In other cases, object-based approaches (OBIA) are used and spectral information is combined with ancillary geo-information (topography maps, groundwater models, object size) to derive the classification schemes rules to discriminate the AGB for the distinct vegetation classes [12,16]. However, in our study area, all the riparian vegetation classes represent co-existence species with identical canopy closure subject to similar biophysical parameters. The distinct riparian objects (RO) had similar patch sizes (MPS) and were evenly distributed along the riparian forest. Thus, no lateral water-gradient was observed on its spatial distribution, nor another geometric attribute that could be used to improve AGB accuracy in our OBIA approach.
Our results also support the direct ability of multispectral UAV imagery data to classify riparian species and to indirectly estimate tree AGB. Other studies have also demonstrated the importance of UAV for species classification in riverine ecosystems [26,27,38]. Additionally, high classification accuracy was obtained with the distinct image bands' combinations showing that for practical proposes riparian objects may be obtained, during the segmentation phase of the OBIA approach, with the Green-Red-NIR, but also with the Green-Red-RedEdge colour-composite bands without loss of classification accuracy. However, in similar ecosystems, with high canopy cover conditions, vegetation indices (VI) developed to suppress saturation effects, like MSR, RedEdgeMSR, and GCI, should be selected to improve the classification accuracy.
The proposed methodology represents a conservative approach in the estimation of the carbon stocks in riparian forests. To minimize the misclassification of the dominant species, only objects with high confidence levels (95% to 99.9%) were classified in the three dominant species. Thus, high values of carbon stocks could be expected in a riparian forest if all vegetation patches would be classified. In addition, total carbon stocks were naturally underestimated using the exclusively remote sensing approach, since understory and soil components were not included due to the superimposition of canopies and the inherent associated uncertainty. Accurate field-remote temporal monitoring of the carbon stocks is particularly needed to improve our process-based knowledge of how riparian and riverine ecosystems dynamics contributes to carbon storage throughout the various stages of succession. However, no single approach, whether exclusive field or remote sensing, can be expected to provide consistent and especially large-scale estimates of biomass, but the use of different methods in a synergistic way may overcome the restrictions of each one.

Conclusions
Carbon stocks from riparian forests may have an uneven contribution to the total carbon pools of Mediterranean systems given the small area occupied in the landscape. In our study, we found an in situ average of above and belowground of carbon stocks ranging from 122.4 to 201.2 tC ha −1 . Although variable according to the riparian species, stem density and vegetation development stage, the woody tree compartment represents the large majority of the carbon reservoirs in our riparian forest (79%). The combined field and remote sensing approach adopted in this study allowed us to evaluate and map, at a very high spatial resolution, the spatial distribution of the tree AGB, which varies from 734.1 to 1052.8 tC, according to the 95%, 99%, and 99.9% probability levels. Nevertheless, these values represent a conservative approach for the carbon stocks estimates, since up to 46.9% to 23.9% of the total study area (according to the distinct probability levels) was not assigned to any vegetation class under study. Our conclusions also support the suitability of multispectral UAV imagery data to indirectly estimate tree AGB via a priori riparian species classification. Spectral vegetation indices developed to suppress saturation effects were consistently selected as important variables for the discrimination of co-existence woody riparian species.
Riparian forests and similar floodplain and wetlands areas represent highly dynamic and challenging ecosystems for ground survey measurements due to their restricted accessibility and thus required advanced but comprehensive methods to estimate biomass and carbon stocks such as remote sensing approaches. Future prospective research should be able to use remote-sensing approaches to discriminate AGB gradients according to the vegetation species, successional stage, DBH or stem density variability. In some cases, LIDAR-derived information may help to achieve high ABG accuracy estimates since it has the unique capability of measuring the three-dimensional vegetation structure [28]. Nevertheless, in highly dynamic riverine ecosystems the extraordinary strata complexity of the uneven-aged forests may impose serious challenges also for LIDAR sensors with limited spatial coverage and relatively large ground footprint.