Estimating Aboveground Biomass and Carbon Stocks in Periurban Andean Secondary Forests Using Very High Resolution Imagery

Periurban forests are key to offsetting anthropogenic carbon emissions, but they are under constant threat from urbanization. In particular, secondary Neotropical forest types in Andean periurban areas have a high potential to store carbon, but are currently poorly characterized. To address this lack of information, we developed a method to estimate periurban aboveground biomass (AGB)—a proxy for multiple ecosystem services—of secondary Andean forests near Bogotá, Colombia, based on very high resolution (VHR) GeoEye-1, Pleiades-1A imagery and field-measured plot data. Specifically, we tested a series of different pre-processing workflows to derive six vegetation indices that were regressed against in situ estimates of AGB. Overall, the coupling of linear models and the Ratio Vegetation Index produced the most satisfactory results. Atmospheric and topographic correction proved to be key in improving model fit, especially in high aerosol and rugged terrain such as the Andes. Methods and findings provide baseline AGB and carbon stock information for little studied periurban Andean secondary forests. The methodological approach can also be used for integrating limited forest monitoring plot AGB data with very high resolution imagery for cost-effective modelling of ecosystem service provision from forests, monitoring reforestation and forest cover change, and for carbon offset assessments.


Introduction
Urban areas and forests play an important role in the carbon (C) cycle and the global climate because they are main drivers in the increase and regulation, respectively, of atmospheric carbon dioxide (CO 2 ), a major greenhouse gas [1,2].Cities in the developing world are often documented as sources of high levels of carbon dioxide emissions [3].It is also reported that without forests and their regulating ecosystem process and services, sensu [4], carbon dioxide concentration in the atmosphere would be about 43% higher than current concentrations [5].Indeed, the International Panel for Climate Change (IPCC) has highlighted the importance of developing baseline estimates of carbon stocks in terrestrial ecosystems, which require mapping vegetation biomass distribution from the local to the global scale [6,7].Precise carbon stock estimates in and around urban areas, therefore, represent key and necessary information to define carbon emission mitigation strategies and programs at the local and regional level (e.g., REDD+, [3,8]).Latin America and the Caribbean regions have about 80% of their population residing in urban areas.Most periurban areas in Neotropical Latin America have historically undergone land use change [9].Recently this urbanization trend has increased and led to the concern of deleterious effects to the hydrology and to increased carbon emissions [10,11].In particular, periurban forests in Neotropical America are providers of key ecosystem services to the highly urbanized population in the region, but they are increasingly being lost [6,12,13].Among terrestrial ecosystems, tropical forests are possibly the most important in regulating global carbon cycle.Lewis et al. [14] estimated that tropical forests are sinks for about 40% of terrestrial carbon on Earth, and responsible for almost half of the terrestrial net primary productivity [15].Pan et al. [7] calculated that tropical forests store about 475 ˆ10 9 Mg of carbon, highlighting the paramount role these ecosystems have in regulating climate and mitigating climate change.Thus, information is needed on how to rapidly map, monitor, and estimate carbon stocks and other ecosystem services in periurban, Neotropical forests in a rapid, cost-effective manner.
In South America, estimations of aboveground biomass (AGB), a key proxy for carbon stocks [16], have been derived especially for lowland tropical rainforests [17][18][19] and for semi-arid savannas and dry forests [20].One of the less studied Neotropical types of forest cover, but with a high potential contribution to the terrestrial C cycle, are early secondary forests whose development follows highly intervened forest land covers.These high Andean periurban secondary forest formations, typical of mountain areas, are present at elevation between ca.2600 and 3200 m [21].This vegetation type is characterized in general terms by a closed, dense canopy and low-to-medium woody plants that develop via ecological succession after land abandonment-typically cattle, pastures, or crops.In Colombia, it is estimated that in 1995 they covered about 7% of the country's terrestrial area [22].Due to the dynamics of rural land abandonment and migration to urban areas and subsequent urbanization characterizing Colombia in recent decades [9,23], it is very likely that Colombian secondary forests are currently following an expansion trend.More importantly, Neotropical secondary forests are still not well-characterized in terms of composition, successional dynamics, carbon cycle, and even biodiversity [24].In this context, the estimation of AGB and carbon in these ecosystems would represent important information to better characterize such a little studied successional forest type from the Neotropics.
Ground-based methods for forest AGB estimation are generally based on permanent plots, forest inventory, and monitoring methods that use species/family-specific or general allometric equations obtained from destructive methods and in situ dendrometric measurements (e.g., diameter at breast height (DBH), wood density, height, volume, and crown area (CA)), [19,[25][26][27].Although these approaches are common and can reach acceptable levels of accuracy, they are not adequate to map AGB distribution at regional scales in the Neotropics, because of the costs, access, and safety issues, and time necessary to perform extended field-based forest inventories.
Remote sensing, both active and passive, provide some of the most time-efficient and cost-effective approaches to derive AGB estimation at regional and national scale.Radar, optical, and LiDAR data have been extensively used to estimate AGB with a variety of methods (see [28,29] for comprehensive reviews).The use of optical data (visible, near-infrared) has been largely exploited for vegetation biomass estimation, due to the large amount and type of optical imagery available.Especially, sensors from the Landsat series (i.e., TM, ETM+, OLI) have been historically used to map biomass and carbon in a variety of ecosystems, for the relevance of their spectral bands, the continuity of the program, and the suitability of the 30 m spatial resolution for regional mapping (e.g., [30][31][32][33]).Several research discussed AGB modelling using Landsat imagery in tropical regions [34][35][36][37].Most of these studies reported difficulties in modelling AGB in complex forest stands such as successional forest.Lu [38] demonstrated that the complexity of forest vegetation results in highly variable standing stocks of AGB and an even more variable rate of AGB accumulation following a deforestation event such as the case in successional forests.Most available optical-based methods are established based on the relationships between field-derived AGB and different vegetation indices (VI).In doing so the aim is the identification of an analytical relationship between on-the-ground AGB estimates and spectral indices sensitive to photosynthetically active radiation and scaling with amount of biomass [39][40][41][42].A vegetation index is generally derived from the spectral reflectance of two or more bands, and proportional to the value of biophysical parameters like leaf area index (LAI), net primary productivity (NPP), and absorbed photosynthetically active radiation (APAR) [43].The choice of adequately performing indices (VIs) and satellite data type depends on the scale of analysis, type of ecosystem and environmental conditions, vegetation density, spectral information available, and the nature of the field information available.
In this work-due to the difficulty in AGB measurements in often inaccessible, dense Andean forests, as well as a limited spatial extent and availability of secondary forest plot data-we hypothesize that the use of very high resolution (VHR) data from commercial satellites, such as the GeoEye-1 and Pleiades-1A satellites, coupled with limited field data can be used to model AGB in complex secondary forests of the Colombian Andes.Previous work with VHR imagery includes the use of IKONOS by Thenkabail et al. [44] to estimate AGB and carbon stock of oil palm plantations in African savannas.Zhou and colleagues [45] exploited VIs and texture images derived from Quickbird and topographic information to derive AGB of black locust plantations in the Loess Plateau of China.Multispectral high-resolution Worldview-2 images were also used by Zhu et al. [46] to derive AGB for mangrove ecosystems using neural networks algorithms.Pereira et al. [47] used GeoEye-1 imagery to estimate biomass and carbon stock of coffee crops in Brazil exploiting a series of VIs.
The aim of this study is to develop an approach to model AGB and carbon stocks in complex secondary periurban Neotropical forests using very high resolution satellite imagery and limited field data from periurban early and late secondary forest monitoring plots.Specific objectives of this research are to: (i) identify an exploitable relationship between VIs derived from VHR satellite imagery and AGB estimated in the field in Andean periurban secondary forests; (ii) assess the effect of data processing on the performance of the biomass estimation models.
Such an approach should yield a low cost and rapid method for forest inventory and monitoring activities, and assessing the regulating ecosystem services and carbon dioxide offset potential of secondary periurban Neotropical forests in the Colombian Andes [8,19,24].

Study Area
The study area is located in the Cundiboyacense high plain, in the Eastern Colombian Andes, at an average altitude of 2800 m (Figure 1).The annual temperature is 14 ˝C, and precipitation varies from 600 mm¨year ´1 in the center of the region and 1200 mm¨year ´1 in the western part.The high plain is characterized by two rainy periods, occurring from April to June and from September to November [48].The high elevation plain environment is highly influenced by the adjacent city of Bogotá and contiguous urbanizations, with about 10 million inhabitants and a very high population density [49].The region is one of the more important agro-industrial centers of the country, with agriculture, pastures, and mining as dominant activities in the rural non-urban areas [10,50].
Secondary periurban Andean forests are present throughout the eastern and northern peripheries of Bogotá [12].The five most common vegetation families found in 20 permanents plots (20 m ˆ20 m)-established in early and late secondary forests in the study are-are Ericaceae, Melastomataceae, Cunoniaceae, Primulaceae, and Asteraceae, representing 56% of individuals with a basal diameter above 5 cm (Norden, Posada, et al., unpublished data).The five dominant genera are Miconia, Weinmannia, Cavendisha, Myrsine, and Myrcianthes, representing 51% of individuals.Eighty species of shrubs and trees have been identified in the area and the five dominant are Other studies have reported that below 2700 m common forest species are Ilex kunthiana Triana 1872, and Vallea stipularis Linnaeus filius 1782, while in drier areas shrubs of the genera Myrcianthes and Morella are more frequent.Other studies report that in the mountains between 2700 and 3000 m the genera that dominate forests are Clusia, Miconia, Weinmania, and Cedrela [21].Other studies report that in the mountains between 2700 and 3000 m the genera that dominate forests are Clusia, Miconia, Weinmania, and Cedrela [21].
In this study, we used eight available 20 m × 20 m permanent plots of young and old secondary vegetation, established in Bogotá city´s northern limits in the suburban municipalities of Guatavita, Guasca, Torca and Tabio, Colombia (Figure 1), at altitudes varying between 2600 and 3100 m, in moderate slope hills.Most plots are dominated by species of low to medium height, as a result of an In this study, we used eight available 20 m ˆ20 m permanent plots of young and old secondary vegetation, established in Bogotá city´s northern limits in the suburban municipalities of Guatavita, Guasca, Torca and Tabio, Colombia (Figure 1), at altitudes varying between 2600 and 3100 m, in moderate slope hills.Most plots are dominated by species of low to medium height, as a result of an early stage of succession, with the exception of two plots of late secondary forests.We identified 57 species in these plots and the most common species were Cavendishia bracteata, Miconia squamulosa, Myrcianthes leucoxyla, Myrsine guianensis, and Weinmannia tomentosa, representing 51% of all the individuals (Norden, Posada et al., unpublished data).

In situ Biomass Estimation
Aboveground biomass was estimated in these eight plots using allometric equations for lower and upper montane Andean forests [51,52].For estimation of AGB in plots located in lower montane Andean forests, we used the following allometric equation [51]: AGB " 0.107314 ¨DBH 2.422   where AGB refers to the aboveground biomass and DBH refers to stem diameter at breast height (1.3 m).For AGB estimation in upper montane Andean forests, we used, for a pool of 10 species, the following allometric equation [52]: The reported biomass values correspond to an average of biomass for both equations (Table 1).We measured all stems with a basal diameter (BD) higher than 5 cm, at 5 cm above the soil surface, for a total of 1354 individuals (Table 1).In five plots we also had measures of DBH while in the remaining three plots DBH was estimated with species-specific or genus-specific equations between BD and DBH.

Number of Species
Above-Ground Biomass (Mg¨ha ´1)

Remote Sensing-Derived Biomass Estimation
We selected VHR GeoEye-1 and Pleiades-1A imagery in order to cover all the eight permanent plots without cloud presence, and acquired in a period of time close to the field census of the vegetation (Table 2).Both VHR sensors cover the blue, green, red, and near infra-red spectral ranges, while the spatial resolution of multispectral imagery is 1.65 m for GeoEye-1 and 2 m for Pleiades-1A.Every scene covers a 5 km ˆ5 km area.

Data Pre-Processing
Data pre-processing for the VHR imagery involved different steps, and led to the production of five groups of data.First, we performed standard radiometric conversion of digital numbers to radiance using sensor specific calibration coefficients provided in the images' metadata (e.g., [53]).In a Forests 2016, 7, 138 6 of 17 second step we applied geometric correction to minimize topographic and sensor geometry image distortions [54].We used a series of field points taken with a sub-meter resolution GPS (Spectra Precision MobileMapper 120) around the 20 m ˆ20 m plots, and post-processed using the nearest local geodetic station BOGA.The geometric correction processing was carried out using first order polynomic equations and nearest neighbor resampling [55] in ENVI 4.3.The root mean square error (RMSE) varied from 0.60 to 1.12 m for the four scenes.The first group of VHR data (hereafter Group A) did not follow any further pre-processing.
Four additional data groups were prepared using additional pre-processing steps to account for atmospheric and topographic effects.These steps are very crucial since different VHR scenes captured at different dates and locations were used in the analysis.Therefore, the images were subjected to either: (i) top-of-atmosphere (TOA), atmospheric correction for data Group B; (ii) radiative transfer model atmospheric correction for Group C; ((iii) relative normalization for Group D; and (iv) topographic correction, for Group E.
Top-of-atmosphere (TOA) reflectance (ρ λ ) was computed using the ENVI v5.2 software: where, L λ is the Radiance in Watts/(m 2 ˆSr); d is the Earth-Sun distance (astronomical distance units); SE λ is the solar irradiance in Watt/m 2 ; and θ is the sun elevation angle in degrees.Atmospheric correction was carried out using the Fast Line-of-sight Atmospheric Analysis of Spectral Hypercubes (FLAASH) method [56].The FLAASH model corrects wavelengths in the visible-mid infrared spectra and it includes an updated MODTRAN4 radiation transfer module [57].The FLAASH atmospheric correction used the average height at each of the study sites within each image, and day visibility values derived from the nearest meteorological station (airport of Bogotá El Dorado, about 30 km distance from the plots) as inputs.Historical meteorological data were derived from the Colombian Institute for Hydrology, Meteorology and Environmental Studies-IDEAM (Table 3).According to the climatic data characteristics, we used the sub-Arctic summer (SAS) atmospheric model and the rural aerosol model within FLAASH.For data group D relative normalization was performed using the dark-object subtraction method followed by the flat bed calibration approach, which does not require pre-knowledge about the atmosphere and/or atmospheric radiative transfer modelling.Consistent dark objects (mostly water bodies) and bright objects (mostly man-made structures and roofs) were used to implement the algorithm [58]:

Relative reflectance values "
Radiance values ´XtDark1u where X tDark1u is the average value of a set of selected dark pixels and X tBright1u is the average value of a set of selected bright pixels.First, the radiation scattering was corrected as the dark-object subtraction method suggests [59], using the Dark Subtract module of ENVI 4.3 to choose a set of dark pixels, compute their average value, and subtract it from the radiance scale image.Then, using the flat-field calibration method [60], the corrected radiance values were turned into relative reflectance values by dividing the radiance values by the average value of a set of chosen bright pixels.Topographic effects on illumination can have a direct effect on the recorded images, especially in the case of large regions with steep slopes [61,62].We examined the use of two topographic correction methods.The first method is the C-correction algorithm [61] applied on the FLAASH atmospherically corrected data set, and the second method is the empirical algorithm suggested by Meyer and colleagues [63] and implemented on high resolution Quickbird imagery by [64].We used the most precise digital surface model (DSM) available for the study area, i.e., Intermap Technologies ® NEXTMap ® World 30™ DSM, with a 30 m ground sampling distance and 10 m LE95 vertical accuracy.The empirical method did not provided meaningful results in two of the study sites; hence only the results of the C-correction method were included in the subsequent analyses.
Finally, at the end of data pre-processing, we had five different data groups as follows: 1. Data group A: data radiometrically calibrated to radiance units.
Data group C: data atmospherically corrected to reflectance units.

4.
Data group D: data normalized to relative reflectance values.

5.
Data group E: data corrected from topographic effect.

Aboveground Biomass Estimation and Carbon Mapping
Five vegetation indices were selected based on studies of remote sensing-derived biomass estimates, and also according to their applicability with the type of imagery used in this study [30,36,[65][66][67][68].Table 4 details the analytic expression of the selected indices and related reference source.

Vegetation Index Equation Reference
Normalized Difference Vegetation Index (NDVI) NDVI " pNIR´REDq pNIR`REDq [69] Vegetation index number (VIN) VIN " NIR RED [70] Ratio Vegetation Index (RVI) RVI " RED NIR [70] Normalized Difference Greenness Index (NDGI) NDGI " pGREEN´REDq pGREEN`REDq [71] Transformed Vegetation Index (TVI) TVI " pNDVI`0.5q|NDVI`0.5| a |NDVI `0.5| [72] All the permanent plots were georeferenced with sub-meter precision using the Spectra Precision MobileMapper 120 ® GPS to derive precise plot boundary polygons.We selected all pixels from inside the plot boundaries (Figure 2), while excluding a few pixels where negative values were observed, a phenomenon commonly observed when radiative-transfer-model-based atmospheric correction is applied.We calculated a single VI value per plot for the five indices by averaging the values of all pixels selected in every plot.
For each plot, the five vegetation indices were then related to the correspondent AGB values estimated from the field observations (Table 1), for all five data groups.To this end, we tested three different models (AGB as a function of VI): linear, semi-log, and log-log.We used both in situ AGB values proportional to the total area of the VHR pixels selected per each plot (Figure 2), and the total AGB of the plot.These two options produced very similar results, hence we only reported the results related to the total plot biomass.The coefficient of determination (R 2 ) and RMSE were calculated as indications of the models goodness of fit and accuracy, respectively [29]. source.
All the permanent plots were georeferenced with sub-meter precision using the Spectra Precision MobileMapper 120 ® GPS to derive precise plot boundary polygons.We selected all pixels from inside the plot boundaries (Figure 2), while excluding a few pixels where negative values were observed, a phenomenon commonly observed when radiative-transfer-model-based atmospheric correction is applied.We calculated a single VI value per plot for the five indices by averaging the values of all pixels selected in every plot.For each plot, the five vegetation indices were then related to the correspondent AGB values estimated from the field observations (Table 1), for all five data groups.To this end, we tested three different models (AGB as a function of VI): linear, semi-log, and log-log.We used both in situ AGB values proportional to the total area of the VHR pixels selected per each plot (Figure 2), and the total To better illustrate the application of the best fit VI model, we derived a VHR carbon map for secondary forests in one specific test area around an existing monitoring plot.First, we performed a parallelepiped supervised classification [55] of the VHR images of the Tabio area's secondary forests vegetation, using field data and training points delineated on-screen using the GoogleEarth interface.Pre-processing was previously applied to the raw data (radiometric, geometric and FLAASH atmospheric correction).The result was a secondary forests vegetation mask image with 92% user accuracy and 90% producer accuracy using 50 randomly selected and visually validated points analyzed in GoogleEarth.In a second step we applied, to the masked pre-processed image, our best model equation to derive a secondary forest aboveground biomass layer for the region.
Carbon stocks can vary substantially according to vegetation type (e.g., broadleaf, conifers) and biome.Accordingly, we assumed that aboveground carbon = 0.471 ˆAGB, based on average C content value calculated using 134 tropical angiosperms species from Thomas and Martin's review [73].Using this relationship, we finally developed a carbon stock map for the Tabio test area.

Results
The goodness of fit and accuracy were calculated using R 2 and RMSE, for the best linear, semi-log, and log-log models of the five VIs (Table 5).All p-values were reported in Table 5 so as to indicate the process used for final model selection.For four out of the five vegetation indices, model goodness increased with the implementation of additional pre-processing steps (e.g., applying atmospheric correction (Group B, C, D) on radiance (Group A) images and further apply topographic correction (Group E)).For the normalized difference greenness index (NDGI)-the only index exploiting the Green band-no evident improvement in model goodness was noticeable among the data groups.The importance of the red wavelength as a most critical band in providing vegetation biophysical information is well-supported in the literature (e.g., [44]).Nevertheless, NDGI shows the best R 2 in the radiance and TOA data groups (A, B).Considering both R 2 and RMSE parameters, the best performing models involved the ratio vegetation index (RVI), transformed vegetation index (TVI), and NDVI indices in Group E, while the Vegetation Index Number (VIN) was the worst performing.In absolute terms, the AGB model which performed better was derived from the combination of the topographic corrected data (calculated on FLAASH atmospherically corrected images) and the ratio vegetation index (Figure 3) in a semi-log form (R 2 = 0.582, p = 0.028): log AGB " ´3.208 ˆRVI `2.185 Forests 2016, 7, 138 9 of 17 Table 5. Goodness of fit and accuracy parameters (R 2 , Root Mean Square Error, p-value) for all vegetation indices and data processing groups in the study.Group A: imagery radiometrically calibrated to radiance units; Group B: top-of-atmosphere (TOA) reflectance data; Group C: data atmospherically corrected to reflectance units; Group D: data normalized to relative reflectance values; Group E: data corrected from topographic effect.Root mean square error (RMSE) expressed in Mg ha ´1.* Significance level p =0.05.For the purpose of better synthesizing our results, we report AGB model equations only for the semi-log form for all indices and data processing groups (Table 6).In general terms, we found how the synergy betweeen the atmospheric and topographic corrections notably improved model goodness of fit in the AGB estimation models.The best performing model for AGB estimation (1) was used to derive maps of aboveground carbon stocks (Mg•ha −1 ) for the Tabio test area.Aboveground biomass values of secondary forests were calculated for the corresponding VHR scenes (Pleiades-1A and GeoEye-1).Using the relationship of C stocks = 0.471 × AGB from Thomas and Martin [73], accordingly we mapped the spatial distribution of aboveground C stocks (Mg•ha −1 ) (Figure 4).For the purpose of better synthesizing our results, we report AGB model equations only for the semi-log form for all indices and data processing groups (Table 6).In general terms, we found how the synergy betweeen the atmospheric and topographic corrections notably improved model goodness of fit in the AGB estimation models.The best performing model for AGB estimation (1) was used to derive maps of aboveground carbon stocks (Mg¨ha ´1) for the Tabio test area.Aboveground biomass values of secondary forests were calculated for the corresponding VHR scenes (Pleiades-1A and GeoEye-1).Using the relationship of  The best performing model for AGB estimation (1) was used to derive maps of aboveground carbon stocks (Mg•ha −1 ) for the Tabio test area.Aboveground biomass values of secondary forests were calculated for the corresponding VHR scenes (Pleiades-1A and GeoEye-1).Using the relationship of C stocks = 0.471 × AGB from Thomas and Martin [73], accordingly we mapped the spatial distribution of aboveground C stocks (Mg•ha −1 ) (Figure 4).

Discussion
The present study investigated the use of remote sensing to cost-effectively estimate AGB biomass and carbon stocks for Andean successional forests in periurban areas of Bogotá, Colombia.The AGB and carbon stock mapping estimates for a mid-successional stage secondary forest in the test area (ca.6-51 Mg¨ha ´1) were within the range reported by other studies of secondary forests in the Neotropics (e.g., [74]).We identified significant analytical relationships between vegetation indices extracted from multiple Geoeye-1 and Pleiades-1A VHR imagery and AGB values estimated using field measurements and allometric equations.A basic underlying assumption of the modelling approach is that larger tree densities (i.e., canopy density) are associated with larger aboveground biomass estimates [65].Among the different preprocessing workflows used, the study confirmed the importance of the synergy between atmospheric and topographic correction when estimating biomass using satellite-based vegetation indices.At an average height of 2800 m our Andean study area is characterized by the frequent presence of dense aerosols and rugged terrain.As such, these are key factors when deriving the spectral characteristics of Andean Neotropical vegetation, and thus represent significant elements to account for when calculating vegetation indices in these environments.Also, the correction model we applied (FLAASH) has the advantage of being easily adaptable to specific types of environmental conditions, allowing for the proper mitigation of the effects of site-specific atmospheric conditions.
We tested the use of empirical atmospheric correction technique through image-to-image normalization and compared the results with the FLAASH model-based atmospheric correction and at sensor radiance results.We found that image-to-image normalization produced results comparable to the FLAASH module results.The manual selection of less-optimal dark and bright objects in each image could explain the slightly inferior results compared to the FLAASH model output results.
However, we believe image-to-image normalization has the potential of providing better results in the case of the existence of variations in illumination conditions due to sparse clouds [75], which is often the case in the Andeans plains..More research is needed to test other statistically-based empirical atmospheric correction as well as different topographic correction models [76,77].
In this study we also tested the use of topographic illumination variation correction algorithm [64].Our results show slight improvement in the AGB model goodness of fit and RMSE when the topographic correction is applied to the images.Although spectrally normalized vegetation indices such as the NDVI are less affected by the topographic-induced illumination variations [78,79], we believe that robust modeling of bidirectional distribution reflection function (BDRF) should provide for a more robust approach to handle the sun-object-sensor geometry variations common in forested land covers [80,81].However, modeling the BDRF requires multi-view and multi-date images, which can be hard to achieve in our study area due to the persistent cloud cover found throughout the year.We also believe that the presence of a high resolution digital elevation model could have improved the topographic correction process and possibly also the AGB model results.
Overall, the VHR imagery provided satisfactory results for AGB and carbon modelling, considering the limitations posed by the number and size of the plots (20 m ˆ20 m) available for this analysis.The ratio vegetation index-the simple ratio between the NIR and the RED bands [70]-provided the best results among our AGB estimation models.It is not possible, however, to generalize about the goodness of the model relating RVI to in situ AGB, that is, if this particular index also has an optimal performance in different successional forests and environmental conditions.However, it does represent a significant relationship for high Andean secondary forest vegetation.
The study was limited by the small number of available permanent plots.Further research will involve the establishment of a larger network of regional and urban to rural gradient plots, to better take into account the heterogeneity of diverse species assemblages and their ecological successional stages, together with stochastic small-scale disturbances.Additional sources of error in our approach are also introduced by the selection of specific allometric equations for estimating AGB [82].Similarly, we did not account for the presence and contribution to total biomass and carbon stocks from nonarboreal vegetation, such as mosses, lianas, ferns, and epiphytes, which can contribute a considerable amount of biomass in tropical secondary forests [83].

Conclusions
Secondary forests like those in this study represent a dominant land cover type, to the extent that in several tropical countries their surface area already exceeds that of mature forest cover [84].In addition, periurban forests that develop after land clearing or reforestation activities or, similarly, after designation as a conservation or protected area, represent key reservoirs for terrestrial biodiversity [85,86], as well as providing a large array of other provision, regulation, and cultural ecosystem services [12,24].For these reasons, we believe the results of the present study added a relevant piece of information for the characterization of these little known successional forests.
Overall, this study's approach and findings provide key baseline information for landscape level analyses in the Colombian Andes relating to local to regional AGB and carbon stocks estimates for little studied secondary forest types.Also, it provides a comprehensive, rapid, and cost-effective framework for AGB estimation using small field plots coupled with available VHR imagery.This approach can specifically be applied to facilitate periurban and regional mapping assessments of forest ecosystem service provision.It can also be used to better understand the effects of direct and indirect drivers of forest cover change, urbanization rates, and ecosystem service provision, as well as for forest inventory and monitoring of REDD+ (Reducing Emissions from Deforestation and forest Degradation) projects.Methods and mapping framework findings could also be exploited to develop carbon offset protocols to mitigate anthropogenic urban carbon emission from adjacent Bogotá and other cities in the Andes.

Figure 2 .
Figure 2. Overlap between plot boundaries (green line) and selected very high resolution pixels, in yellow (GeoEye-1 imagery) for the Guasca, Colombia, study site.

Figure 2 .
Figure 2. Overlap between plot boundaries (green line) and selected very high resolution pixels, in yellow (GeoEye-1 imagery) for the Guasca, Colombia, study site.

Figure 3 .
Figure 3. Best performing aboveground biomass (AGB) estimation model based on the ratio vegetation index (RVI) from Group E.

Figure 3 .
Figure 3. Best performing aboveground biomass (AGB) estimation model based on the ratio vegetation index (RVI) from Group E.

Figure 4 .
Figure 4. Aboveground carbon stock distribution map for secondary Andean in the Tabio test area near Bogotá, Colombia in Mg•ha −1 (based on Pleiades-1A imagery).

Figure 4 .
Figure 4. Aboveground carbon stock map for secondary Andean in the Tabio test area near Bogotá, Colombia in Mg¨ha ´1 (based on Pleiades-1A imagery).

Table 1 .
Plot ID, forest successional status, number of individuals, number of species, and average above-ground biomass (Mg¨ha ´1).

Table 2 .
Remote sensing images and acquisition configurations used in the present study.

Table 3 .
Parametrization of the Fast Line-of-sight Atmospheric Analysis of Spectral Hypercubes (FLAASH) atmospheric correction model.SAS is the sub-Arctic summer model.

Table 4 .
Selected vegetation indices and reference source.GREEN, RED and near infrared (NIR) relate to the pixel values for the corresponding bands.

Table 6 .
Aboveground biomass (AGB) semi-log models and model goodness of fit (R 2 ) and root mean square error (RMSE) for all vegetation indices and five data processing group.* Significance level p = 0.05.(Best performing model in bold.)