Next Article in Journal
Physical Activity in Forest and Psychological Health Benefits: A Field Experiment with Young Polish Adults
Next Article in Special Issue
Linking Forest Management Practices to the Functional Composition of Plant Communities
Previous Article in Journal
Analysis of Spatial-Temporal Changes and Driving Factors of Vegetation Coverage in Jiamusi City
Previous Article in Special Issue
Is the Invasiveness of Pittosporum undulatum in Eucalypt Forests Explained by the Wide Ranging Effects of Its Secondary Metabolites?
 
 
Search for Articles:
Title / Keyword
Author / Affiliation / Email
Journal
Article Type
 
 
Advanced Search
 
Section
Special Issue
Volume
Issue
Number
Page
 
Journals
Forests
Volume 14
Issue 9
10.3390/f14091903
forests-logo
Submit to this Journal Review for this Journal Propose a Special Issue
► Article Menu

Article Menu

share Share announcement Help format_quote Cite question_answer Discuss in SciProfiles thumb_up
...
Endorse textsms
...
Comment
first_page
settings
Order Article Reprints
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Mapping Above-Ground Carbon Stocks at the Landscape Scale to Support a Carbon Compensation Mechanism: The Chocó Andino Case Study

by
Francisco Cuesta
1,*,
Marco Calderón-Loor
2,*,
Paulina Rosero
1,
Noam Miron
2,
Andrei Sharf
2,
Carolina Proaño-Castro
3 and
Felipe Andrade
3
1
Grupo de Investigación en Biodiversidad Medio Ambiente y Salud-BIOMAS, Universidad de las Américas, Vía a Nayón, Quito 170124, Ecuador
2
Albo Climate, Ahad Haam, 9, Tel Aviv 6525101, Israel
3
Fundación Futuro, Av. De Las Palmeras N45-74 y de las Orquídeas, Quito 170502, Ecuador
*
Authors to whom correspondence should be addressed.
Forests 2023, 14(9), 1903; https://doi.org/10.3390/f14091903
Submission received: 3 June 2023 / Revised: 28 August 2023 / Accepted: 5 September 2023 / Published: 19 September 2023
(This article belongs to the Special Issue Biodiversity and Ecosystem Functioning in Forests)
Browse Figures
Versions Notes

Abstract

:
(1) Background: Tropical Mountain forests (TMF) constitute a threatened major carbon sink due to deforestation. Carbon compensation projects could significantly aid in preserving these ecosystems. Consequently, we need a better understanding of the above-ground carbon (AGC) spatial distribution in TMFs to provide project developers with accurate estimations of their mitigation potential; (2) Methods: integrating field measurements and remote sensing data into a random forest (RF) modelling framework, we present the first high-resolution estimates of AGC density (Mg C ha−1) over the western Ecuadorian Andes to inform an ongoing carbon compensation mechanism; (3) Results: In 2021, the total landscape carbon storage was 13.65 Tg in 194,795 ha. We found a broad regional partitioning of AGC density mediated primarily by elevation. We report RF-estimated AGC density errors of 15% (RMSE = 23.8 Mg C ha−1) on any 10 m pixel along 3000 m of elevation gradient covering a wide range of ecological conditions; (4) Conclusions: Our approach showed that AGC high-resolution maps displaying carbon stocks on a per-pixel level with high accuracy (85%) could be obtained with a minimum of 14 ground-truth plots enriched with AGC density data from published regional studies. Likewise, our maps increased precision and reduced uncertainty concerning current methodologies used by international standards in the Voluntary Carbon Market.

1. Introduction

Human-induced land use change is one of the main contributors to the current atmospheric concentration of greenhouse gases [GHG; [1] ]. Reducing GHG emissions is of transcendent importance, as curving them will directly impact the global climatic system [2] and can help achieve the goals of the Paris Agreement [3]. To meet the ambitious 1.5 °C goal set by the Paris Agreement, substantial reductions in deforestation and associated CO2 emissions are required, with targets of approximately 70% by 2030 and 95% by 2050 [4,5]. Moreover, forest restoration and regrowth, leading to CO2 sequestration, must continue as a primary mitigation strategy [6].
However, the extent and distribution of climate mitigation opportunities available from preserving, restoring, or enhancing tropical mountain forest cover have been insufficiently explored [7,8]. One major challenge in designing carbon compensation projects at subnational and national levels is the uncertainty associated with the spatial distribution of carbon stored in above-ground biomass (AGB) [9,10,11]. This uncertainty is particularly pronounced in tropical mountain forests (TMFs) due to sharp environmental gradients, limited data availability, and knowledge gaps [12,13]. Likewise, in many cases, tropical mountain forests’ remnants are embedded in human-dominated landscapes at different successional stages of recovery or degradation [11,14], hence increasing the spatial variation of above-ground carbon (AGC) [13].
Accurate estimation of AGC density (Mg C ha−1) at the landscape scale depends mainly on the availability and quality of ground-truth and remote sensing (RS) data together with the modeling techniques [15,16,17]. While local forest plot measurements provide accurate information at the plot level, their scarcity in tropical mountainous regions poses a challenge when extrapolating their values to the landscape scale [18]. Integrating local plot data with existing high-quality ground-truth AGC datasets offers a promising approach to improve AGC stock predictions at a pixel scale and to estimate model uncertainty [19,20,21].
The integration of RS data is crucial for extrapolating plot-based AGC density estimates and deriving pixel-level estimates of AGC stocks and changes [9,22]. Recent advances in RS technology have provided comprehensive coverage, high spatial resolution, and short return cycles. The Sentinel-1 C-band Synthetic Aperture Radar (SAR) and the Sentinel-2 multispectral instrument by the European Space Agency (ESA), are examples of these advances [23]. The combination of different sensors has the potential to overcome limitations associated with single remote sensing techniques for AGC density estimates [24].
The choice of a modelling framework is also critical to improving the accuracy of AGC predictions beyond field data [25,26]. While parametric models such as linear regressions have been widely used [27,28,29], non-parametric algorithms, such as machine learning models, offer a better ability to capture complex relationships and deal with skewed data [30]. Machine learning algorithms have demonstrated high accuracy in producing spatially explicit AGC estimates at the pixel level [15,31]. However, previous AGC mapping exercises in tropical regions have often been limited by coarse spatial resolution (~500 m to 1000 m) and high uncertainties, rendering them less suitable for local-scale carbon compensation projects [32,33,34,35].
Natural landscapes in Ecuador have changed rapidly in the last four decades [36,37,38,39]. The lowland and mountain forests of north-western Ecuador, the Chocó-Andes region, is one of the three regions of the country where most of the deforestation has occurred [40,41]. One of the largest remaining forest tracts is embedded in the UNESCO Chocó Andino Biosphere Reserve—RBCA hereafter [42]. From 1991 to 2017, forest loss within the RBCA was estimated at an annual gross loss rate of 0.66%. However, forest regrowth was observed for the same period, equal to 8.5% of the forest coverage in 1991 [43]. Moreover, the RBCA has witnessed a high density of forest conservation and restoration efforts in the last decade, offering a unique opportunity to establish a compensation mechanism to reduce deforestation and enhance carbon sequestration.
This research presents an innovative approach to map AGC forest stocks and their temporal variation using RS coupled with advanced statistical methods at high spatial resolution (10 m-pixel resolution) where local spatiotemporal ground-truth data is limited. Our approach leverages machine learning techniques to establish complex non-linear relationships between observed AGC values, satellite imagery channels, and ancillary environmental data. Specifically, we developed a set of high-resolution maps of AGC, with geospatially explicit uncertainty estimates for the RBCA to inform an ongoing local carbon compensation mechanism, the Zero Carbon program (https://nftree.com.ec/, accessed on 2 June 2023). This study focused exclusively on the AGC of standing trees ≥ 5 cm in diameter, not on the belowground or necromass carbon pools.

2. Materials and Methods

2.1. Study Area

This study was conducted in the western versant of the Ecuadorian Andes (0°11′17.8656″ N, 0°7′36.21″ S; −78°34′8.8566″ W, −78°54′46.242″ W), which spans over 286,800 ha from the outer foothills of the western cordillera of the Andes to its inner foothills, along 4000 m of elevation gradient (Figure S1). Forests in the study area are humid, evergreen ecosystems with limited or no seasonality characteristic of tropical mountain systems. Forest composition gradually changes along the elevation gradient driven by temperature conditions [44,45] defining four ecosystems with different structure, composition, and above-ground biomass [46]. Mean annual air temperature diminishes from 21.6 °C at 653 m a.s.l. to 7.2 °C at 3507 m a.s.l., decreasing ca. 0.56 °C per 100 m a.s.l., reaching values of 0.7 °C at the upper end (3507 m a.s.l.). By 2020, native forest covered nearly 68% (i.e., 195,000 ha) of the study area, 3% (i.e., 10,000 ha) was mountain grasslands located above the forest line (above 3900 m a.s.l.), and 9% was mountain shrublands (i.e., 25,000 ha). The remaining 20% (i.e., 57,000 ha) constituted agricultural lands. Agriculture and cattle raising constitute an essential livelihood, with over 80% of the Biosphere Reserve’s productive land use being dedicated to extensive cattle raising [43].

2.2. Ground-Truth Data

The ground-truth data inputs for this study come from an in situ AGB dataset from the Pichincha Forest plot network, established in 2015 [45,47], enriched with a published regional AGB dataset: https://doi.org/10.5061/dryad.59zw3r26f, accessed on 8 November 2022 [48] coming from a regional Andean forest plot network—Red Bosques Andinos—RBA (https://redbosques.condesan.org/, accessed on 8 November 2022) and a subset of an extensive AGB pantropical dataset [19]. Both the local and regional databases provided AGB estimations based on permanent plots to study forest dynamics. From the larger pantropical dataset [19], we generated a subset based on the following criteria: (1) AGB records with tree cover < 70%, based on [49], were removed from the analysis. Albeit, removing these plots might limit the applicability of our AGC model to less mature forests (i.e., forests between 30% and 70% of tree cover), more than 90% of the forest area inside the study area has a tree cover > 70% (Figure S2); (2) the geographic extent included only AGB data from the north-western Amazonia, the Andes (Perú and Bolivia), and the Chocó forests (Colombia, Ecuador, and Panama); (3) latitudinal and longitudinal coordinates were not degraded or incomplete; (4) AGB field measurements were sampled from 2000 onwards; (5) extreme AGB values falling below the 5th percentile or exceeding the 95th percentile were removed.
The resulting dataset with 494 AGB plots was randomly split into a training and validation dataset, 80% and 20%, respectively (Table 1), covering the tropical lowland and montane forests of Colombia, Ecuador, Perú, Bolivia, and Panamá (Figure 1). The dataset covers an elevation gradient of 3240 m a.s.l. with a mean value of 431 m a.s.l. (±531 m a.s.l.), with a skewed distribution towards the lower section of the elevation range (i.e., 75% of the plots correspond to an altitudinal range bounded between 100 and 560 m a.s.l.; Figure 2a). All AGB plot locations had a mean forest cover of 97% ± Sd = 4.4 (Figure 2b). The mean plot size was 0.63 ha (±0.32 ha). Lastly, to obtain AGC values, we multiplied the individual plot AGB value by a conversion constant of 0.465 [50].

2.3. Remote Sensing Data Assembling

Multispectral Sentinel-2 (S-2, hereafter) and Sentinel-1 Synthetic Aperture Radar (S-1, hereafter) imagery from the European Space Agency were accessed from the Google Earth Engine (GEE) platform [51]. We used atmospherically corrected S-2 scenes (bands B2-B8, B11, and B12) to predict AGC in the RBCA. We created a stack of images for the target year, i.e., 2018–2021 ± one year (e.g., images from 2017 to 2019 were used for producing the year 2018 mosaic) to minimize missing data and cloud interference [52]. We then applied the s2cloudless algorithm implemented in GEE for masking clouds and cloud-shadows to every stack. Each multiyear stack contained around 5 to 20 valid observations at pixel-level (year and location dependent). In addition, we created yearly S-1 mosaics. S-1 images were retrieved from the ‘Descending orbit’ (to eliminate potential discrepancies in the backscatter) and pre-processed by applying an additional border noise correction, a multi-temporal boxcar speckle filtering and a radiometric terrain normalization ([53]; Table S1).
Different spectral indices were calculated and appended to the cloud-free stacks (Table S1). Likewise, elevation and slope data (30-m) from the Shuttle Radar Topography Mission (SRTM) dataset [54] were included as predictor variables. Lastly, X (longitude) and Y (latitude) coordinates were generated for each pixel of the study area and integrated with the RS dataset. The geographical coordinates were included to capture the environmental gradients that influence our study region, e.g., harsher climatic conditions towards the east and a marked carbon stocks variation from north to south. All the aligned and registered layers were resampled to a 10 m resolution if necessary. The final composites were produced by computing the median value of every variable/band at the pixel level in the 2017–2019 period.

2.4. Above-Ground Carbon (AGC) Regression Models

We used a random forest (RF) algorithm [55] to develop a regression-based machine learning model. RF is a computationally efficient tree-based machine learning model robust to random and systematic noise [56,57]. RF builds several trees, previously specified by the user, iteratively from random training data samples (n = 394). The ensemble of the predictions of all trees constitutes the final model. We parametrized the model to generate 1000 trees with 15 variables per split; a grid search algorithm was used to find the optimum number of trees. We used the scikit-learn library for fine tuning the RF model in Python [58]. A total of 23 variables were used as predictors (Table S2).
To increase the robustness of our regression and to provide a measure of the model’s estimation uncertainty at the pixel level we used the mean of ten distinct RF models as the AGC value for every pixel and every time step. Our dataset, n = 494, was divided into ten folds [59]; for each fold, an 80% random sample (i.e., n = 394) was used to train an RF model. For every time step, ten AGC maps were generated and averaged to produce AGC density estimates at the pixel level.

2.5. Error and Uncertainty Estimation

We assessed the error estimation and precision for the predicted AGC models by comparing the predicted AGC (i.e., mean of the 10 maps) against the observed ground-truth AGC. Error and precision were calculated for the entire training dataset (n = 393), validation dataset (n = 100) and the Pichincha Forest plot network exclusively (n = 23). To evaluate the performance of the AGC models, we used the adjusted R2 metric, which considers the model’s goodness of fit while adjusting for model complexity to avoid overfitting. Additionally, we employed three widely used performance indicators: Root Mean Square Error (RMSE), Mean Absolute Error (MAE), and the Mean Absolute Percentage Error (MAPE) [23,60]. These metrics allowed us to quantify the magnitude of differences between the observed and predicted AGC values and provided a comprehensive evaluation of model accuracy.
Lastly, we assessed the generalization capabilities of our AGC model by comparing the outputs of independent random forest (RF) models against localized ground-truth data. Four distinct areas of interest (AOI) were defined, each containing a different number of local plots (Table 2). In an iterative process, 600 different RF models were trained for each AOI. For constructing the training and validation datasets for each AOI, the following procedure was followed: (1) plots belonging to the specific AOI under consideration were removed from the entire dataset; (2) the remaining dataset was randomly split into training and validation subsets, with proportions of 80% and 20%, respectively, in each iteration; (3) local plots were then randomly selected and incrementally added to both the training and validation datasets; (4) the maximum number of local plots that were added to the training subset was limited to 80% of the total local plots available within the AOI. This iterative process allowed us to progressively evaluate the impact of adding more localized data on the model’s performance. The validation dataset, consisting of 20% of the data, including local validation points, was used to estimate the model’s performance for each AOI. This approach ensured a comprehensive assessment of the model’s ability to generalize across different areas with varying degrees of local data availability.

2.6. Carbon Maps’ Post-Processing

A post-processing was applied to all the generated AGC maps. First, non-forest areas were masked-out from the results using the Global Forest Change V1.10 [49]. Then, a temporal coefficient of variation (CV) at the pixel-level was calculated and used for improving the temporal consistency of the AGC maps. Pixels with a temporal CV > 1.5, i.e., high temporal uncertainty, and high slope (>25°) were smoothed by using a focal window. We used a circular kernel of 80 m radius to calculate the average AGC inside the moving window. In each time-step the AGC in the selected pixels was replaced with the smoothed AGC value. Then, a temporal correction was applied to all the previously smoothed pixels. AGC values in these pixels were replaced by their minimum temporal AGC value to reduce temporal inconsistencies [49].

2.7. Carbon Stocks Spatial Patterns, Gains, and Losses

As elevation is the single most important factor controlling AGC stocks in tropical mountain ecosystems [61,62,63], AGC maps were decomposed into six elevation ranks to assess the spatial distribution of the above-ground carbon across the study area. Elevation ranks were set based on natural breaks in the histogram of the elevation grid. We also generated a profile along the elevation gradient and summarised the AGC value per pixel for each elevation rank and generated AGC averages per each rank. Lastly, we compared our AGC model against pre-existing AGC global datasets [19,64,65] across the elevation gradient in the RBCA.
We estimated carbon stock change between 2018 and 2021 at the pixel scale to assess carbon gains and losses. We defined the AGC net change (Mg ha−1 yr−1) as the difference between AGC stock in 2021 and AGC stock in 2018 divided by the elapsed time (t; in years) [AGC net change = (AGC2021 − AGC2018)/4]. Positive values were defined as forest carbon stock gain, whereas negative values as carbon loss.

3. Results

In 2018, the forests of the RBCA stored a mean of 70.26 (±8.1) Mg C ha−1, ranging from 55.6 to 102.8 Mg C ha−1 (Table 3, Figure 3). In 2021 a slight increase in AGC density was observed at 70.64 (±8.1) Mg C ha−1 at the pixel scale. In 2021, total carbon storage at the landscape scale was 13.66 Tg (million megagrams) in its above-ground compartment, evidencing an annual net gain of 71,886 Mg C, equalling a productivity of 0.37 Mg C ha−1 yr−1 (Table 4). This net increase resulted from a greater increase in pixels that gained biomass (mean = 0.60 ± 0.66 Mg C−1 ha−1 yr−1) than pixels that lost biomass (mean = −0.49 ± 0.61 Mg C−1 ha−1 yr−1) between 2018 and 2021.
We found a broad regional partitioning of standing carbon stocks mediated primarily by elevation, longitude, and latitude (Figure 4). The mean AGC of the piedmont forest (340–1200 m a.s.l.) was 80 (±6.2) Mg C ha−1 and decreased monotonically along the gradient to as low as 59.8 (±3.8) Mg C ha−1 in the upper montane forests (>3500 m a.s.l.; Figure 5).
Based on the analysis of the entire validation dataset (n = 100), we successfully predicted above-ground carbon (AGC) values across the study region with high precision (R2 = 0.6; p < 0.0001) and accuracy at a fine resolution of 0.01 ha (RMSE = 44.3 Mg C ha−1 in 2021; Table 5). Subsequently, we focused solely on the Pichincha plots (n = 23), which resulted in improved precision (adj. R2 = 0.7); however, it also led to an increase in the error, accounting for nearly 28% of the mean AGC density (RMSE = 23.8 Mg C ha−1; Table 5). Likewise, the RF models with different training subsets for the four AOIs depict that our AGC model can achieve an overall mean accuracy of ~85%, when at least 14 local sampling plots are combined with regional plots to reach a relatively high mapping accuracy (Table 6). Table 5 and Table 6 report differences in RMSE, MAE, and MAPE at the RBCA level. Main reasons for these differences are the number of RF models that were employed in each analysis. While we used the average of 10 maps for calculating the accuracy and error metrics in Table 5, the tabular output of 600 RF models was employed for calculating and assessing the generalization capabilities of our approach in Table 6.
When comparing our AGC model with the global maps available for the tropics [19,22,65], we found that all global maps tended to overestimate AGC density across the study area, and the variation reported by all global maps was at least two times higher than our results (Table 7). Further, the three global datasets failed to represent the AGC/elevation trend reported for our study area. All global AGC maps predicted the montane forests between 1800 and 2500 m a.s.l. to contain the higher AGC density followed by the range between 1200 and 1800 m a.s.l. (Figure S3, Table 7).

4. Discussion

Our predictive random forest (RF) models for above-ground carbon (AGC) revealed significant carbon storage in the montane forest of the RBCA, with a total of 13,662,109 Mg C stored in the above-ground compartment in 2021. The piedmont (600–1200 m a.s.l.) and lower montane (1200–1800 m a.s.l.) forests on the northwestern versant of the study area were identified as the primary sinks, accounting for nearly 60% of the total AGC stored in the RBCA. Above this elevation range, the AGC stocks showed a steady decrease (refer to Table 7 and Figure 5). This observation aligns with previous studies that reported a decline in AGC with increasing elevation along elevation gradients in the RBCA and the tropical Andes [48,61,62,66], reinforcing the capability of our AGC models to accurately capture the spatial variation of AGC influenced by elevation.
The observed increase in AGC over time was found to be lower than that reported for mountain forests in the Andes using permanent plots (0.67 ± 0.08 Mg C ha−1 yr−1; [48]). Additionally, between 2015 and 2019, 16 out of the 21 local plots in our study area exhibited an AGC increase at an annual rate of 1.20 ± 0.89 Mg ha−1 yr−1 [46]. However, at larger spatial scales environmental heterogeneity increases, leading to varied AGC trends. As a result, many pixels across the RBCA indicated a decrease in AGC between 2018 and 2021. These pixels showed a mean AGC loss of −0.49 (±0.61) AGC ha−1 yr−1 (Table S3). On the other hand, pixels reporting a consistent increase in carbon stocks showed a mean AGC gain of 0.60 (±0.66) AGC ha−1 yr−1. To better understand the trajectory of these pixels over time, further analyses are needed. For example, conducting a time series analysis for a subset of pixels (e.g., 15% of all pixels exhibiting AGC loss and gain) surrounding field plot locations using spectral–temporal segmentation algorithms (e.g., LandTrendr https://github.com/eMapR/LT-GEE, accessed on 2 June 2023), can provide valuable insights into forest cover changes based on the pixel’s spectral history, utilizing moderate or high-resolution imagery [67]. Such an approach would complement our understanding of AGC pixel trends over shorter timeframes.

4.1. Uncertainty and Its Impact on the Calculation of Carbon Credits

Our approach has proven highly effective in estimating above-ground carbon (AGC) densities at the landscape scale, providing predictive values (Mg ha−1) at the pixel level that accurately reflect AGC spatial variation within complex landscape mosaics, such as the western versant of the Ecuadorian Andes. The increased AGC resolution achieved through our methodology has significant implications for the design of local Payment for Ecosystem Services (PES) schemes, particularly in determining the size of carbon benefits that can be generated for any project area. Notably, our high-resolution maps offer enhanced precision and reduced uncertainty compared to current methodologies used by international standards in the Voluntary Carbon Market (e.g., Plan Vivo, Verra), where average AGC density per hectare is estimated for each forest stratum [68,69].
Significant differences emerged when comparing our approach with available global AGC datasets (Table 7). Global datasets consistently estimated higher AGC contents at the pixel scale, with considerable dispersion around the means, ranging from 31% to 53% variation. In contrast, our model demonstrated a considerably lower dispersion, with only 12% variation. On the landscape scale, on average, the three datasets estimated 7.4 million Mg more carbon stored in the RBCA than our estimates, potentially leading to an overestimation of the mitigation potential (refer to Table 6). Several factors contribute to these discrepancies: (1) The inability of the three global datasets to capture the AGC/elevation trend reported for tropical mountain ecosystems (e.g., [12,48,61]). All global AGC maps predicted that montane forests between 1200 and 2500 m a.s.l. contained the highest carbon density in our study area (Table 7, Figure S3). As this elevation band constitutes 60% of the forest cover in the RBCA, it significantly increases the estimated total carbon at the landscape scale. (2) There is limited plot coverage in tropical mountain regions of the Andes, which forms the basis for global products. Plot distribution tends to be skewed toward lowland areas [18,70] leading to a lack of representation in higher elevation regions. This limitation hinders the applicability of global datasets and increases the uncertainty surrounding the quantified AGC estimates. (3) Globally available forest masks at a 30 m resolution (i.e., [49]), may overestimate forest cover remnants in agricultural mosaics. A 30 m pixel resolution product might inaccurately depict forest cover remnants in complex topographic regions, where forest fragments exist in areas of high slope and are embedded within diverse land use mosaics (e.g., palmetto plantations and cacao agroforestry systems). This can lead to overestimation of carbon stocks in such regions.

4.2. Robustness of the AGC Model

Our approach demonstrated a remarkable reduction in the number of local plots required to achieve a statistically valid sample for obtaining highly accurate AGC estimates. Current methodologies for REDD+ projects typically recommend a minimum of 45 local plots for training (70%) and validating remote-sensing-based AGC models [68]. However, our approach showcased that AGC high-resolution maps, displaying carbon density at a per-pixel level with high accuracy (≥85%), can be obtained with as few as 15 local plots enriched with AGC data from published regional studies and existing forest and carbon datasets. This significant reduction in the number of required local plots can lead to substantial cost savings and improved efficiency in REDD+ project designs. Furthermore, our proposed approach allows for near real-time monitoring of AGC dynamics across the entire study region. Although we utilized S-1 and S-2 imagery, our approach is sensor-agnostic and can be applied to different satellite imagery sources, whether public or commercial, at various spatiotemporal resolutions. These findings can have a meaningful impact on REDD+ project implementation, particularly for projects covering extensive areas with diverse land-use/land-cover types and limited accessibility, where reaching a statistically valid sampling size using traditional plot-based biomass measurements can be challenging.
The inclusion of environmental covariates in the RF model significantly improved the accuracy of our AGC estimates. Elevation, longitude, and latitude were identified as the most influential factors in the RF model (refer to Figure 4). This aligns with previous research that highlights the strong influence of environmental variables on AGC spatial distribution, especially in complex topographic regions [61,63,71]. Interestingly, an average pixel-level error of 20.6 Mg C ha−1 compares favourably with past studies that reported uncertainties ranging from 21 Mg C ha−1 [15] to 40 Mg C ha−1 in tropical forests [60,72,73,74], see Figure S4. The accuracy achieved by our AGC model is a crucial aspect for the success of carbon sequestration schemes and can support more informed decision-making in climate mitigation efforts.

4.3. Implications for Forest Conservation and Restoration in the RBCA

The AGC maps for the RBCA were generated to be integrated into the Zero Carbon Compensation Program. The Zero Carbon Program emerges as a mechanism to promote, on the one hand, the decarbonization of Grupo Futuro (GF) companies (https://www.futuro.com.ec/, accessed on 2 June 2023) and, on the other, as a financing mechanism for the conservation of biodiversity. The Zero Carbon Program is structured into three components. The first focuses on quantifying greenhouse gas (GHG) emissions from GF companies, also known as carbon footprint (based on international quality standards, i.e., ISO 14064 [75], 14067 [76] standards; IPCC Protocol for the AFOLU sector). The second component centres on implementing strategies to reduce GF companies’ carbon footprint. The third component focuses on designing and implementing compensation strategies through a trust fund, the Carbon Neutral Fund, formed in 2019. Funds from this trust are invested in the conservation and restoration of forests via conservation agreements between GF and landowners in the RBCA.
The third component, related to the compensation mechanism, is implemented on a technological platform called NFTrees (https://nftree.com.ec/, accessed on 2 June 2023). Our carbon maps were uploaded in blockchain technology to allow interested buyers to purchase digital tokens to compensate for their emissions. The purchase is made through a credit card or wire transfer directly to the Carbon Neutral Trust Fund. The NFTrees platform allows transactional traceability and avoids double accounting (i.e., different parties cannot claim the same carbon removal or reduction credit). Likewise, blockchain technology generates a unique and unalterable record (a Token). The token guarantees that any hectare of forest in which the carbon has been compensated cannot be sold again. Currently, the GF’s Zero Carbon Program model generates USD 250,000 per year, equivalent to the carbon offset of around 13,000 Mg of CO2-eq of 12 companies of the GF, which translates into the conservation of 3000 ha of forests within the RBCA. With the AGC maps integrated with blockchain technology and digital contracts to gain transparency and traceability, the Zero Program envisions generating USD 1.9 million in the next five years to finance the conservation of 22,600 hectares of forest and benefit at least 700 families within the RBCA.

4.4. Limitations and Further Research

RS mapping of above-ground carbon (AGC) stocks relies on various factors that can influence the quality of the outputs. While we introduced a method to reduce the need for extensive local ground-truth data, the quantity and quality of regional/global AGC data remain critical for accurately representing AGC in areas with high environmental heterogeneity, such as complex topography and diverse climates. The adequate representation of such heterogeneity requires a sufficient number of well-distributed regional records, obtained from forest plots in regional/global datasets. In the case of the RBCA, AGC stocks were partly underestimated due to the scarcity of up-to-date AGC data from mountainous areas, particularly above 1000 m a.s.l. (Figure 2a). The lack of updated regional ground-truth data from regions with similar environmental conditions can increase the Root Mean Square Error (RMSE) of models and reduce their prediction capabilities. To address this, initiatives releasing global AGB/AGC datasets that can be used for training and validating RS-based models must be properly supported with fair compensation mechanisms, ensuring the maintenance of long-term plot monitoring programs [77]. Additionally, the incorporation of ground-truth data from the Global Ecosystem Dynamics Investigation (GEDI) Light Detection and Ranging (LiDAR) sensor [78] as a source of ground-truth data can increase the prediction capabilities of RS-based AGC models (e.g., [79]).
Accurate carbon credits and potential mitigation estimates need the combination of removals and emissions (i.e., AGC stocks) and activity data (i.e., land cover change). Locally generated land change information needs to be developed alongside AGC data for improving wall-to-wall AGC estimates while ensuring the generated data’s transparency. Here we have used a previously generated global forest cover dataset at 30 m [49] for generating a forest vs. non-forest mask for each of the four AGC maps (2018–2021). Although the Global Forest Change [49] product has been validated and extensively used worldwide, it presents limitations when used at local scales, particularly in landscape mosaics of forest remnants embedded in agricultural lands.
Similarly, further research is needed to produce temporally consistent AGC stock maps, especially in tropical mountainous regions. The atmospheric and topographic characteristics of the regions add an extra layer of complexity to the consolidation of quality multitemporal satellite images. Here, we implemented and applied different pre-processing algorithms to the multitemporal stacks. However, some remnant noise can influence the AGC stocks at the pixel level. Time series analyses to decompose the raw data can be beneficial for removing seasonality and noise signals from the stacks [80]. Furthermore, adding the AGC time series in the training dataset can also improve the temporal consistency of AGC estimates. Although different authors have used statistical approaches for generating historical AGC maps (e.g., validating these products requires updated and multitemporal local and regional AGB/AGC data (e.g., [81]).

5. Conclusions

In this study, we have addressed a critical knowledge gap by providing comprehensive insights into the spatial distribution of above-ground carbon (AGC) density in the TMFs of the western Ecuadorian Andes. Our findings hold significant implications for carbon compensation projects aimed at preserving these threatened ecosystems. By integrating field measurements and remote sensing data within an RF modeling framework, we have achieved a pioneering advancement in accurately estimating AGC density at high spatial resolutions. By employing this approach, we have successfully produced the first high-resolution estimates of AGC density (Mg C ha−1) across the western Ecuadorian Andes. Notably, our approach yields a remarkable accuracy of 85%, with RF-estimated AGC density errors of only 15% (RMSE = 23.8 Mg C ha−1) on any 10 m pixel along a 3000 m elevation gradient, encompassing diverse ecological conditions. The resulting AGC density maps provide a detailed overview of carbon stocks at a per-pixel level, offering insights into the factors controlling the spatial distribution of AGC density at the landscape scale. The outcomes of our study revealed a pronounced regional partitioning of AGC density predominantly driven by elevation variations. This nuanced understanding of AGC distribution enriches the foundation for carbon compensation mechanisms, enabling project developers to make more precise estimations of the mitigation potential inherent in TMFs. Furthermore, our investigation established the feasibility of obtaining highly accurate AGC density maps by utilizing a minimum of 14 ground-truth plots augmented with AGC density data from published regional studies. This approach not only enhances the precision of carbon stock assessments but also contributes to reducing the uncertainty associated with existing methodologies employed within the Voluntary Carbon Market and international standards.

Supplementary Materials

The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/f14091903/s1, Figure S1: Location of the 23 permanent plots across an elevation gradient on the Choco Andino Biosphere Reserve (RBCA for its Spanish acronym); Figure S2: Forest cover density (>70%) in 2019 on the Chocó Andino Biosphere Reserve. Forest cover density derived from [49]; Figure S3: AGC stock estimates by four different models. Maps A-C correspond to global carbon models for the tropics at different spatial resolutions; Figure S4: Variation in AGC uncertainty at pixel scale based on the coefficient of variation of 10 random forest models. White areas within the RBCA represent no-forest cover (i.e., pixels with less than 70% of forest cover, based on [49] for the year 2021). Polygons inside de RBCA delimit Quito municipality conservation and sustainable use areas (ACUs for its Spanish acronym); Table S1: Description of satellite bands used to build the Chocó Andino Biosphere Reserve AGC model; Table S2: Variable relative importance of the random forest model performed to predict AGC in the RBCA for the reference year 2019. Number of Trees: 1000; Table S3: Pixels that increase and decrease in above-ground carbon for the period of 2018–2021 in the Chocó Andino Biosphere Reserve.

Author Contributions

Conceptualization, F.C., M.C.-L., C.P.-C. and F.A.; methodology, M.C.-L., F.C., P.R. and A.S.; formal analysis, M.C.-L., F.C. and P.R.; writing—original draft preparation, F.C., M.C.-L. and P.R.; writing—review and editing, F.C., M.C.-L., P.R., N.M., A.S., C.P.-C. and F.A.; project administration, F.C. and C.P.-C.; funding acquisition, F.C. and C.P.-C. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by Universidad de las Américas, Ecuador (Grant No. DIV.FCC.22.02), and Fundación Futuro (Grant No. UDLA-002).

Data Availability Statement

Not Applicable.

Acknowledgments

We want to thank Hedi Zisling for his help with retrieving and preparing the RS data. We thank Reut Keller, Huanhuan Huang, and Ira Bar for their inputs on an earlier version of the manuscript. We express our gratitude to Esteban Pinto, Antonella Bernardi, Andrea TeránValdez, Laura Cifuentes, and Emilia Ayala for their invaluable work in leading and collecting the AGB in the Pichincha long-term monitoring transect. Likewise, we would like to express our gratitude to the owners and administrators of the private reserves where the forest dynamics permanent plots are established: Alejandro Solano, María Emilia Arcos, and Oliver Torres from Mashpi Shungo and Pambiliño; Mateo Roldán from Mashpi Lodge/Fundación Futuro, Sonia Saltos from Río Bravo, Juan Diego Ortíz, and Sol Acosta from Sacha Urcu, Inti Arcos and Nina Duarte from Intyllacta, Richard Parson and Favio Fernández from Bellavista Cloud Forest, Germán Toasa from El Cedral Ecolodge, María Elisa Vizuete and Enrique Maldonado from Verdecocha and the Jocotoco foundation from the Yanacocha Reserve. We thank Universidad de las Americas for financially supporting this research (Grant No. DIV.FCC.22.02) and Fundación Futuro (Grant No. UDLA-002). We acknowledge ALBO Climate for their continuous support throughout this project.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Pörtner, H.-O.; Roberts, D.C.; Adams, H.; Adler, C.; Aldunce, P.; Ali, E.; Begum, R.A.; Betts, R.; Kerr, R.B.; Biesbroek, R. Climate Change 2022: Impacts, Adaptation and Vulnerability; IPCC Sixth Assessment Report; IPCC: Geneva, Switzerland, 2022. [Google Scholar]
  2. Bastin, J.-F.; Finegold, Y.; Garcia, C.; Mollicone, D.; Rezende, M.; Routh, D.; Zohner, C.M.; Crowther, T.W. The Global Tree Restoration Potential. Science 2019, 365, 76–79. [Google Scholar] [CrossRef]
  3. IPCC Summary for Policymakers. In Climate Change 2022: Mitigation of Climate Change. Contribution of Working Group III to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change; Shukla, P.R.; Skea, J.; Slade, R.; Al Khourdajie, A.; van Diemen, R.; McCollum, D.; Pathak, M.; Some, S.; Vyas, P.; Fradera, R.; et al. (Eds.) Cambridge University Press: Cambridge, UK, 2022; pp. 1–53. [Google Scholar]
  4. Roe, S.; Streck, C.; Obersteiner, M.; Frank, S.; Griscom, B.; Drouet, L.; Fricko, O.; Gusti, M.; Harris, N.; Hasegawa, T.; et al. Contribution of the Land Sector to a 1.5 °C World. Nat. Clim. Chang. 2019, 9, 817–828. [Google Scholar] [CrossRef]
  5. Streck, C. REDD+ and Leakage: Debunking Myths and Promoting Integrated Solutions. Clim. Policy 2021, 21, 843–852. [Google Scholar] [CrossRef]
  6. Goetz, S.J.; Hansen, M.; Houghton, R.A.; Walker, W.; Laporte, N.; Busch, J. Measurement and Monitoring Needs, Capabilities and Potential for Addressing Reduced Emissions from Deforestation and Forest Degradation under REDD+. Environ. Res. Lett. 2015, 10, 123001. [Google Scholar] [CrossRef]
  7. Cook-Patton, S.C.; Leavitt, S.M.; Gibbs, D.; Harris, N.L.; Lister, K.; Anderson-Teixeira, K.J.; Briggs, R.D.; Chazdon, R.L.; Crowther, T.W.; Ellis, P.W.; et al. Mapping Carbon Accumulation Potential from Global Natural Forest Regrowth. Nature 2020, 585, 545–550. [Google Scholar] [CrossRef] [PubMed]
  8. Aragón, S.; Salinas, N.; Nina-Quispe, A.; Qquellon, V.H.; Paucar, G.R.; Huaman, W.; Porroa, P.C.; Olarte, J.C.; Cruz, R.; Muñiz, J.G.; et al. Aboveground Biomass in Secondary Montane Forests in Peru: Slow Carbon Recovery in Agroforestry Legacies. Glob. Ecol. Conserv. 2021, 28, e01696. [Google Scholar] [CrossRef]
  9. Herold, M.; Carter, S.; Avitabile, V.; Espejo, A.B.; Jonckheere, I.; Lucas, R.; McRoberts, R.E.; Næsset, E.; Nightingale, J.; Petersen, R.; et al. The Role and Need for Space-Based Forest Biomass-Related Measurements in Environmental Management and Policy. Surv. Geophys. 2019, 40, 757–778. [Google Scholar] [CrossRef]
  10. González-Jaramillo, V.; Fries, A.; Zeilinger, J.; Homeier, J.; Paladines-Benitez, J.; Bendix, J. Estimation of above Ground Biomass in a Tropical Mountain Forest in Southern Ecuador Using Airborne LiDAR Data. Remote Sens. 2018, 10, 660. [Google Scholar] [CrossRef]
  11. Christmann, T.; Menor, I.O. A Synthesis and Future Research Directions for Tropical Mountain Ecosystem Restoration. Sci. Rep. 2021, 11, 23948. [Google Scholar] [CrossRef] [PubMed]
  12. Cuni-Sanchez, A.; Sullivan, M.J.P.; Platts, P.J.; Lewis, S.L.; Marchant, R.; Imani, G.; Hubau, W.; Abiem, I.; Adhikari, H.; Albrecht, T.; et al. High Aboveground Carbon Stock of African Tropical Montane Forests. Nature 2021, 596, 536–542. [Google Scholar] [CrossRef]
  13. Christmann, T.; Palomeque, X.; Armenteras, D.; Wilson, S.J.; Malhi, Y.; Oliveras Menor, I. Disrupted Montane Forest Recovery Hinders Biodiversity Conservation in the Tropical Andes. Glob. Ecol. Biogeogr. 2023, 32, 793–808. [Google Scholar] [CrossRef]
  14. Mathez-Stiefel, S.L.; Peralvo, M.; Báez, S.; Rist, S.; Buytaert, W.; Cuesta, F.; Fadrique, B.; Feeley, K.J.; Groth, A.A.P.; Homeier, J.; et al. Research Priorities for the Conservation and Sustainable Governance of Andean Forest Landscapes. Mt. Res. Dev. 2017, 37, 323–339. [Google Scholar] [CrossRef]
  15. Asner, G.P.; Brodrick, P.G.; Philipson, C.; Vaughn, N.R.; Martin, R.E.; Knapp, D.E.; Heckler, J.; Evans, L.J.; Jucker, T.; Goossens, B. Mapped Aboveground Carbon Stocks to Advance Forest Conservation and Recovery in Malaysian Borneo. Biol. Conserv. 2018, 217, 289–310. [Google Scholar] [CrossRef]
  16. Rodríguez-Veiga, P.; Quegan, S.; Carreiras, J.; Persson, H.J.; Fransson, J.E.S.; Hoscilo, A.; Ziółkowski, D.; Stereńczak, K.; Lohberger, S.; Stängel, M.; et al. Forest Biomass Retrieval Approaches from Earth Observation in Different Biomes. Int. J. Appl. Earth Obs. Geoinf. 2019, 77, 53–68. [Google Scholar] [CrossRef]
  17. Bustamante, M.M.C.; Roitman, I.; Aide, T.M.; Alencar, A.; Anderson, L.O.; Aragão, L.; Asner, G.P.; Barlow, J.; Berenguer, E.; Chambers, J. Toward an Integrated Monitoring Framework to Assess the Effects of Tropical Forest Degradation and Recovery on Carbon Stocks and Biodiversity. Glob. Chang. Biol. 2016, 22, 92–109. [Google Scholar] [CrossRef] [PubMed]
  18. Blundo, C.; Carilla, J.; Grau, R.; Malizia, A.; Malizia, L.; Osinaga-Acosta, O.; Bird, M.; Bradford, M.; Catchpole, D.; Ford, A. Taking the Pulse of Earth’s Tropical Forests Using Networks of Highly Distributed Plots. Biol. Conserv. 2021, 260, 108849. [Google Scholar]
  19. Avitabile, V.; Herold, M.; Heuvelink, G.B.M.; Lewis, S.L.; Phillips, O.L.; Asner, G.P.; Armston, J.; Ashton, P.S.; Banin, L.; Bayol, N.; et al. An Integrated Pan-Tropical Biomass Map Using Multiple Reference Datasets. Glob. Chang. Biol. 2016, 22, 1406–1420. [Google Scholar] [CrossRef] [PubMed]
  20. Araza, A.; De Bruin, S.; Herold, M.; Quegan, S.; Labriere, N.; Rodriguez-Veiga, P.; Avitabile, V.; Santoro, M.; Mitchard, E.T.A.; Ryan, C.M. A Comprehensive Framework for Assessing the Accuracy and Uncertainty of Global Above-Ground Biomass Maps. Remote Sens. Environ. 2022, 272, 112917. [Google Scholar] [CrossRef]
  21. Réjou-Méchain, M.; Tymen, B.; Blanc, L.; Fauset, S.; Feldpausch, T.R.; Monteagudo, A.; Phillips, O.L.; Richard, H.; Chave, J. Using Repeated Small-Footprint LiDAR Acquisitions to Infer Spatial and Temporal Variations of a High-Biomass Neotropical Forest. Remote Sens. Environ. 2015, 169, 93–101. [Google Scholar] [CrossRef]
  22. Baccini, A.; Walker, W.; Carvalho, L.; Farina, M.; Sulla-Menashe, D.; Houghton, R.A. Tropical Forests Are a Net Carbon Source Based on Aboveground Measurements of Gain and Loss. Science 2017, 358, 230–234. [Google Scholar] [CrossRef]
  23. Chen, L.; Ren, C.; Zhang, B.; Wang, Z.; Xi, Y. Estimation of Forest Above-Ground Biomass by Geographically Weighted Regression and Machine Learning with Sentinel Imagery. Forests 2018, 9, 582. [Google Scholar] [CrossRef]
  24. Duncanson, L.; Armston, J.; Disney, M.; Avitabile, V.; Barbier, N.; Calders, K.; Carter, S.; Chave, J.; Herold, M.; MacBean, N. Aboveground Woody Biomass Product Validation: Good Practices Protocol; Committee on Earth Observation Satellites: Reston, VA, USA, 2021.
  25. Musthafa, M.; Singh, G. Improving Forest Above-Ground Biomass Retrieval Using Multi-Sensor L-and C-Band SAR Data and Multi-Temporal Spaceborne LiDAR Data. Front. For. Glob. Chang. 2022, 5, 822704. [Google Scholar] [CrossRef]
  26. Li, Y.; Li, M.; Li, C.; Liu, Z. Forest Aboveground Biomass Estimation Using Landsat 8 and Sentinel-1A Data with Machine Learning Algorithms. Sci. Rep. 2020, 10, 9952. [Google Scholar] [CrossRef]
  27. Nelson, R.; Jimenez, J.; Schnell, C.E.; Hartshorn, G.S.; Gregoire, T.G.; Oderwald, R. Technical Note: Canopy Height Models and Airborne Lasers to Estimate Forest Biomass: Two Problems. Int. J. Remote Sens. 2000, 21, 2153–2162. [Google Scholar] [CrossRef]
  28. Lu, D. Aboveground Biomass Estimation Using Landsat TM Data in the Brazilian Amazon. Int. J. Remote Sens. 2005, 26, 2509–2525. [Google Scholar] [CrossRef]
  29. Rodrigues de Moura Fernandes, M.; Fernandes da Silva, G.; Quintão de Almeida, A.; Marques Fernandes, M.; Ribeiro de Mendonça, A.; Guimarães Gonçalves, F.; Henrique Breda Binoti, D.; Pereira Martins Silva, J. Aboveground Biomass Estimation in Dry Forest in Northeastern Brazil Using Metrics Extracted from Sentinel-2 Data: Comparing Parametric and Non-Parametric Estimation Methods. Adv. Sp. Res. 2023, 72, 361–377. [Google Scholar] [CrossRef]
  30. Fassnacht, F.E.; Hartig, F.; Latifi, H.; Berger, C.; Hernández, J.; Corvalán, P.; Koch, B. Importance of Sample Size, Data Type and Prediction Method for Remote Sensing-Based Estimations of Aboveground Forest Biomass. Remote Sens. Environ. 2014, 154, 102–114. [Google Scholar] [CrossRef]
  31. Mascaro, J.; Asner, G.P.; Knapp, D.E.; Kennedy-Bowdoin, T.; Martin, R.E.; Anderson, C.; Higgins, M.; Chadwick, K.D. A Tale of Two “Forests”: Random Forest Machine Learning Aids Tropical Forest Carbon Mapping. PLoS ONE 2014, 9, e85993. [Google Scholar] [CrossRef]
  32. Cartus, O.; Kellndorfer, J.; Walker, W.; Franco, C.; Bishop, J.; Santos, L.; Michel Fuentes, J.M. A National, Detailed Map of Forest Aboveground Carbon Stocks in Mexico. Remote Sens. 2014, 6, 5559–5588. [Google Scholar] [CrossRef]
  33. Rodríguez-Veiga, P.; Saatchi, S.; Tansey, K.; Balzter, H. Magnitude, Spatial Distribution and Uncertainty of Forest Biomass Stocks in Mexico. Remote Sens. Environ. 2016, 183, 265–281. [Google Scholar] [CrossRef]
  34. Saatchi, S.; Malhi, Y.; Zutta, B.; Buermann, W.; Anderson, L.O.; Araujo, A.M.; Phillips, O.L.; Peacock, J.; ter Steege, H.; Lopez Gonzalez, G.; et al. Mapping Landscape Scale Variations of Forest Structure, Biomass, and Productivity in Amazonia. Biogeosci. Discuss. 2009, 6, 5461–5505. [Google Scholar] [CrossRef]
  35. Zhang, Y.; Liang, S. Fusion of Multiple Gridded Biomass Datasets for Generating a Global Forest Aboveground Biomass Map. Remote Sens. 2020, 12, 2559. [Google Scholar] [CrossRef]
  36. de Koning, G.H.J.; Veldkamp, A.; Fresco, L.O. Land Use in Ecuador: A Statistical Analysis at Different Aggregation Levels. Agric. Ecosyst. Environ. 1998, 70, 231–247. [Google Scholar] [CrossRef]
  37. Sierra, R. Dynamics and Patterns of Deforestation in the Western Amazon: The Napo Deforestation Front, 1986–1996. Appl. Geogr. 2000, 20, 1–16. [Google Scholar] [CrossRef]
  38. Messina, J.P.; Walsh, S.J.; Mena, C.F.; Delamater, P.L. Land Tenure and Deforestation Patterns in the Ecuadorian Amazon: Conflicts in Land Conservation in Frontier Settings. Appl. Geogr. 2006, 26, 113–128. [Google Scholar] [CrossRef]
  39. Cuenca, P.; Arriagada, R.; Echeverría, C. How Much Deforestation Do Protected Areas Avoid in Tropical Andean Landscapes? Environ. Sci. Policy 2016, 56, 56–66. [Google Scholar] [CrossRef]
  40. Sierra, R.; Stallings, J. The Dynamics and Social Organization of Tropical Deforestation in Northwest Ecuador, 1983–1995. Hum. Ecol. 1998, 26, 135–161. [Google Scholar] [CrossRef]
  41. Lepers, E.; Lambin, E.F.; Janetos, A.C.; DeFries, R.; Achard, F.; Ramankutty, N.; Scholes, R.J. A Synthesis of Information on Rapid Land-Cover Change for the Period 1981–2000. Bioscience 2005, 55, 115–124. [Google Scholar] [CrossRef]
  42. Kleemann, J.; Zamora, C.; Villacis-Chiluisa, A.B.; Cuenca, P.; Koo, H.; Noh, J.K.; Fürst, C.; Thiel, M. Deforestation in Continental Ecuador with a Focus on Protected Areas. Land 2022, 11, 268. [Google Scholar] [CrossRef]
  43. Wiegant, D.; Peralvo, M.; van Oel, P.; Dewulf, A. Five Scale Challenges in Ecuadorian Forest and Landscape Restoration Governance. Land Use Policy 2020, 96, 104686. [Google Scholar] [CrossRef]
  44. Van Manen, M.L.T.; Jansen, B.; Cuesta, F.; León-Yánez, S.; Gosling, W.D. From Leaf to Soil: N-Alkane Signal Preservation, despite Degradation along an Environmental Gradient in the Tropical Andes. Biogeosciences 2020, 17, 5465–5487. [Google Scholar] [CrossRef]
  45. Llerena-Zambrano, M.; Ordoñez, J.C.; Llambí, L.D.; van der Sande, M.; Pinto, E.; Salazar, L.; Cuesta, F. Minimum Temperature Drives Community Leaf Trait Variation in Secondary Montane Forests along a 3000-m Elevation Gradient in the Tropical Andes. Plant Ecol. Divers. 2021, 14, 47–63. [Google Scholar] [CrossRef]
  46. Pinto, E.; Cuesta, F.; Bernardi, A.; Llerena-Zambrano, M.; Pérez, A.J.; van der Sande, M.T.; Gosling, W.D.; Burgess, K.S. Determinants of Above-Ground Carbon Stocks and Productivity in Secondary Forests along a 3000 m Elevation Gradient in the Ecuadorian Andes. Plant Ecol. Divers 2023. submitted. [Google Scholar]
  47. Pinto, E.; Cuesta, F. Monitoreo de Biodiversidad, Contenidos de Carbono, Productividad y Rasgos Funcionales En Bosques Montanos; CONDESAN: Quito, Ecuador, 2019. [Google Scholar]
  48. Duque, A.; Peña, M.A.; Cuesta, F.; González-Caro, S.; Kennedy, P.; Phillips, O.L.; Calderón-Loor, M.; Blundo, C.; Carilla, J.; Cayola, L. Mature Andean Forests as Globally Important Carbon Sinks and Future Carbon Refuges. Nat. Commun. 2021, 12, 2138. [Google Scholar] [CrossRef]
  49. Hansen, M.C.; Potapov, P.V.; Moore, R.; Hancher, M.; Turubanova, S.A.; Tyukavina, A.; Thau, D.; Stehman, S.V.; Goetz, S.J.; Loveland, T.R.; et al. High-Resolution Global Maps of 21st-Century Forest Cover Change. Science 2013, 342, 850–853. [Google Scholar] [CrossRef]
  50. Martin, A.R.; Doraisami, M.; Thomas, S.C. Global Patterns in Wood Carbon Concentration across the World’s Trees and Forests. Nat. Geosci. 2018, 11, 915–920. [Google Scholar] [CrossRef]
  51. Gorelick, N.; Hancher, M.; Dixon, M.; Ilyushchenko, S.; Thau, D.; Moore, R. Google Earth Engine: Planetary-Scale Geospatial Analysis for Everyone. Remote Sens. Environ. 2017, 202, 18–27. [Google Scholar] [CrossRef]
  52. Scarth, P.; Armston, J.; Flood, N.; Denham, R.; Collett, L.; Watson, F.; Trevithick, B.; Muir, J.; Goodwin, N.; Tindall, D. Operational Application of the Landsat Timeseries to Address Large Area Landcover Understanding. Int. Arch. Photogramm. Remote Sens. Spat. Inf. Sci. 2015, 40, 571. [Google Scholar] [CrossRef]
  53. Mullissa, A.; Vollrath, A.; Odongo-Braun, C.; Slagter, B.; Balling, J.; Gou, Y.; Gorelick, N.; Reiche, J. Sentinel-1 Sar Backscatter Analysis Ready Data Preparation in Google Earth Engine. Remote Sens. 2021, 13, 1954. [Google Scholar] [CrossRef]
  54. Farr, T.G.; Rosen, P.A.; Caro, E.; Crippen, R.; Duren, R.; Hensley, S.; Kobrick, M.; Paller, M.; Rodriguez, E.; Roth, L. The Shuttle Radar Topography Mission. Rev. Geophys. 2007, 45. [Google Scholar] [CrossRef]
  55. Breiman, L. Random Forests. Mach. Learn. 2001, 45, 5–32. [Google Scholar] [CrossRef]
  56. Denisko, D.; Hoffman, M.M. Classification and Interaction in Random Forests. Proc. Natl. Acad. Sci. USA 2018, 115, 1690–1692. [Google Scholar] [CrossRef]
  57. Pelletier, C.; Valero, S.; Inglada, J.; Champion, N.; Marais Sicre, C.; Dedieu, G. Effect of Training Class Label Noise on Classification Performances for Land Cover Mapping with Satellite Image Time Series. Remote Sens. 2017, 9, 173. [Google Scholar] [CrossRef]
  58. Pedregosa, F.; Varoquaux, G.; Gramfort, A.; Michel, V.; Thirion, B.; Grisel, O.; Blondel, M.; Prettenhofer, P.; Weiss, R.; Dubourg, V.; et al. Scikit-Learn: Machine Learning in Python. J. Mach. Learn. Res. 2011, 12, 2825–2830. [Google Scholar]
  59. Bengio, Y.; Grandvalet, Y. No Unbiased Estimator of the Variance of K-Fold Cross-Validation. Adv. Neural Inf. Process. Syst. 2004, 5, 1089–1105. [Google Scholar]
  60. Asner, G.P.; Mascaro, J.; Anderson, C.; Knapp, D.E.; Martin, R.E.; Kennedy-Bowdoin, T.; van Breugel, M.; Davies, S.; Hall, J.S.; Muller-Landau, H.C. High-Fidelity National Carbon Mapping for Resource Management and REDD+. Carbon Balance Manag. 2013, 8, 7. [Google Scholar] [CrossRef]
  61. Girardin, C.A.J.; Farfan-Rios, W.; Garcia, K.; Feeley, K.J.; Jørgensen, P.M.; Murakami, A.A.; Cayola Pérez, L.; Seidel, R.; Paniagua, N.; Fuentes Claros, A.F.; et al. Spatial Patterns of Above-Ground Structure, Biomass and Composition in a Network of Six Andean Elevation Transects. Plant Ecol. Divers. 2014, 7, 161–171. [Google Scholar] [CrossRef]
  62. Phillips, J.; Ramirez, S.; Wayson, C.; Duque, A. Differences in Carbon Stocks along an Elevational Gradient in Tropical Mountain Forests of Colombia. Biotropica 2019, 51, 490–499. [Google Scholar] [CrossRef]
  63. Slik, J.W.F.; Aiba, S.; Brearley, F.Q.; Cannon, C.H.; Forshed, O.; Kitayama, K.; Nagamasu, H.; Nilus, R.; Payne, J.; Paoli, G. Environmental Correlates of Tree Biomass, Basal Area, Wood Specific Gravity and Stem Density Gradients in Borneo’s Tropical Forests. Glob. Ecol. Biogeogr. 2010, 19, 50–60. [Google Scholar] [CrossRef]
  64. Baccini, A.; Goetz, S.J.; Walker, W.S.; Laporte, N.T.; Sun, M.; Sulla-Menashe, D.; Hackler, J.; Beck, P.S.A.; Dubayah, R.; Friedl, M.A. Estimated Carbon Dioxide Emissions from Tropical Deforestation Improved by Carbon-Density Maps. Nat. Clim. Chang. 2012, 2, 182–185. [Google Scholar] [CrossRef]
  65. Spawn, S.A.; Sullivan, C.C.; Lark, T.J.; Gibbs, H.K. Harmonized Global Maps of above and Belowground Biomass Carbon Density in the Year 2010. Sci. Data 2020, 7, 112. [Google Scholar] [CrossRef] [PubMed]
  66. Girardin, C.A.J.; Malhi, Y.; AragÃO, L.E.O.C.; Mamani, M.; Huaraca Huasco, W.; Durand, L.; Feeley, K.J.; Rapp, J.; Silva-Espejo, J.E.; Silman, M.; et al. Net Primary Productivity Allocation and Cycling of Carbon along a Tropical Forest Elevational Transect in the Peruvian Andes. Glob. Chang. Biol. 2010, 16, 3176–3192. [Google Scholar] [CrossRef]
  67. Kennedy, R.E.; Yang, Z.; Gorelick, N.; Braaten, J.; Cavalcante, L.; Cohen, W.B.; Healey, S. Implementation of the LandTrendr Algorithm on Google Earth Engine. Remote Sens. 2018, 10, 691. [Google Scholar] [CrossRef]
  68. Terra Global Capital. VCS Methodology for Carbon Accounting for Mosaic and Landscape-Scale REDD Projects, V2.2—VM0006; 2017. Available online: https://verra.org/methodologies/vm0006-methodology-for-carbon-accounting-for-mosaic-and-landscape-scale-redd-projects-v2-2/ (accessed on 2 June 2023).
  69. Berry, N.; Baumann, D.; Foster, B.; Kaonga, M.; Langan-Farmer, J.; Wheeler, C. Agriculture and Forestry Carbon Benefit Assessment Methodology—PM001. Version 0.1. The Landscapes and Livelihoods Group and Plan Vivo Foundation Technical Advisory Committee. 2022, p30. Available online: https://www.planvivo.org/Handlers/Download.ashx?IDMF=9362bb39-5dc5-45c1-a240-600148494ae9 (accessed on 2 June 2023).
  70. Malizia, A.; Blundo, C.; Carilla, J.; Acosta, O.O.; Cuesta, F.; Duque, A.; Aguirre, N.; Aguirre, Z.; Ataroff, M.; Baez, S.; et al. Elevation and Latitude Drives Structure and Tree Species Composition in Andean Forests: Results from a Large-Scale Plot Network. PLoS ONE 2020, 15, e0231553. [Google Scholar] [CrossRef] [PubMed]
  71. Santoro, M.; Cartus, O.; Carvalhais, N.; Rozendaal, D.; Avitabile, V.; Araza, A.; De Bruin, S.; Herold, M.; Quegan, S.; Rodríguez-Veiga, P. The Global Forest Above-Ground Biomass Pool for 2010 Estimated from High-Resolution Satellite Observations. Earth Syst. Sci. Data 2021, 13, 3927–3950. [Google Scholar] [CrossRef]
  72. Asner, G.P.; Powell, G.V.N.; Mascaro, J.; Knapp, D.E.; Clark, J.K.; Jacobson, J.; Kennedy-Bowdoin, T.; Balaji, A.; Paez-Acosta, G.; Victoria, E. High-Resolution Forest Carbon Stocks and Emissions in the Amazon. Proc. Natl. Acad. Sci. USA 2010, 107, 16738–16742. [Google Scholar] [CrossRef]
  73. Mascaro, J.; Asner, G.P.; Muller-Landau, H.C.; van Breugel, M.; Hall, J.; Dahlin, K. Controls over Aboveground Forest Carbon Density on Barro Colorado Island, Panama. Biogeosciences 2011, 8, 1615–1629. [Google Scholar] [CrossRef]
  74. Asner, G.P.; Knapp, D.E.; Martin, R.E.; Tupayachi, R.; Anderson, C.B.; Mascaro, J.; Sinca, F.; Chadwick, K.D.; Higgins, M.; Farfan, W. Targeted Carbon Conservation at National Scales with High-Resolution Monitoring. Proc. Natl. Acad. Sci. USA 2014, 111, E5016–E5022. [Google Scholar] [CrossRef]
  75. ISO 14064-1:2018; Greenhouse Gases—Part 1: Specification with Guidance at the Organization Level for Guantification and Reporting of Greenhouse Gas Emissions and Removals. ISO: Genève, Switzerland, 2018. Available online: https://www.iso.org/standard/66453.html (accessed on 5 August 2022).
  76. ISO 14067:2018; Greenhouse Gases—Carbon Footprint of Products—Requirements and Guidelines for Quantification. ISO: Genève, Switzerland, 2018. Available online: https://www.iso.org/standard/71206.html (accessed on 5 August 2022).
  77. de Lima, R.A.F.; Phillips, O.L.; Duque, A.; Tello, J.S.; Davies, S.J.; de Oliveira, A.A.; Muller, S.; Honorio Coronado, E.N.; Vilanova, E.; Cuni-Sanchez, A. Making Forest Data Fair and Open. Nat. Ecol. Evol. 2022, 6, 656–658. [Google Scholar] [CrossRef]
  78. Liang, M.; González-Roglich, M.; Roehrdanz, P.; Tabor, K.; Zvoleff, A.; Leitold, V.; Silva, J.; Fatoyinbo, T.; Hansen, M.; Duncanson, L. Assessing Protected Area’s Carbon Stocks and Ecological Structure at Regional-Scale Using GEDI Lidar. Glob. Environ. Chang. 2023, 78, 102621. [Google Scholar] [CrossRef]
  79. Rishmawi, K.; Huang, C.; Zhan, X. Monitoring Key Forest Structure Attributes across the Conterminous United States by Integrating Gedi Lidar Measurements and Viirs Data. Remote Sens. 2021, 13, 442. [Google Scholar] [CrossRef]
  80. Zhang, X.; Jia, W.; Sun, Y.; Wang, F.; Miu, Y. Simulation of Spatial and Temporal Distribution of Forest Carbon Stocks in Long Time Series—Based on Remote Sensing and Deep Learning. Forests 2023, 14, 483. [Google Scholar] [CrossRef]
  81. Powell, S.L.; Cohen, W.B.; Healey, S.P.; Kennedy, R.E.; Moisen, G.G.; Pierce, K.B.; Ohmann, J.L. Quantification of Live Aboveground Forest Biomass Dynamics with Landsat Time-Series and Field Inventory Data: A Comparison of Empirical Modeling Approaches. Remote Sens. Environ. 2010, 114, 1053–1068. [Google Scholar] [CrossRef]
Forests 14 01903 g001
Figure 1. Above-ground biomass (AGB) reference dataset for the tropical lowland and montane forests of Colombia, Ecuador, Perú, Bolivia, and Panamá.
Figure 1. Above-ground biomass (AGB) reference dataset for the tropical lowland and montane forests of Colombia, Ecuador, Perú, Bolivia, and Panamá.
Forests 14 01903 g001
Forests 14 01903 g002
Figure 2. Above-ground biomass data distribution along (a) elevation gradient and (b) tree cover values for the year 2000 based on [49].
Figure 2. Above-ground biomass data distribution along (a) elevation gradient and (b) tree cover values for the year 2000 based on [49].
Forests 14 01903 g002
Forests 14 01903 g003
Figure 3. Variation in above-ground carbon density (AGC, Mg C ha−1) at 10 m pixel resolution throughout the Chocó Andino Biosphere Reserve (RBCA for its Spanish acronym) in the western Ecuadorian Andes, derived from an integrated model of remote sensing, local, regional, and global field plot datasets integrated in a random forest algorithm.
Figure 3. Variation in above-ground carbon density (AGC, Mg C ha−1) at 10 m pixel resolution throughout the Chocó Andino Biosphere Reserve (RBCA for its Spanish acronym) in the western Ecuadorian Andes, derived from an integrated model of remote sensing, local, regional, and global field plot datasets integrated in a random forest algorithm.
Forests 14 01903 g003
Forests 14 01903 g004
Figure 4. Variable importance in a random forest model based on the Gini impurity index for the Chocó Andino Biosphere Reserve (RBCA for its Spanish acronym).
Figure 4. Variable importance in a random forest model based on the Gini impurity index for the Chocó Andino Biosphere Reserve (RBCA for its Spanish acronym).
Forests 14 01903 g004
Forests 14 01903 g005
Figure 5. The variation with elevation of the predicted above-ground carbon density (AGC, Mg C ha−1) at a pixel size of 10 m in the remnant forests of the Chocó Andino Biosphere Reserve (RBCA for its Spanish acronym).
Figure 5. The variation with elevation of the predicted above-ground carbon density (AGC, Mg C ha−1) at a pixel size of 10 m in the remnant forests of the Chocó Andino Biosphere Reserve (RBCA for its Spanish acronym).
Forests 14 01903 g005
Table 1. Training and validation dataset’s composition.
Table 1. Training and validation dataset’s composition.
DatasetTrainingValidation
Pichincha Forest plot network—RBCA158
Regional37992
Total394100
Table 2. Metadata of four areas of interest (AOI) where model’s generalization capabilities were assessed.
Table 2. Metadata of four areas of interest (AOI) where model’s generalization capabilities were assessed.
AOIAOI Area
(km2)		Coordinates Extent (°)Median
Elevation (Interquartile Range)
(m a.s.l.)		Local Plots Size (n)
LatitudeLongitude
RBCA286879.381,
78.36		−0.643,
−0.395		1267 (1380)23
Bolivia1,165,573−69.5,
−68		−15.5,
−14		1494 (1462)23
Peru1,804,27372.146,
−70.986		−13.368,
−11.781		1519 (2188)28
Bolivia and Peru7,761,764−72.206,
−68		−15.5,
−11.67		1200 (1973)67
Table 3. Estimated Above-Ground Carbon density (AGC, Mg C ha−1) and AGC stocks (landscape) in the Chocó Andino Biosphere Reserve (RBCA) for the period of 2018–2021.
Table 3. Estimated Above-Ground Carbon density (AGC, Mg C ha−1) and AGC stocks (landscape) in the Chocó Andino Biosphere Reserve (RBCA) for the period of 2018–2021.
YearAGC Pixel-Level Mean (±Sd)AGC Pixel-Level Range (Min–Max)Landscape Total (Tg)
201870.26 (±8.1)55.6–102.813.45
201970.65 (±8.2)55.6–102.813.71
202070.64 (±8.2)55.7–102.813.76
202170.64 (±8.1)55.6–101.913.66
Table 4. Forest cover, AGC stocks, and AGC net change between 2018 and 2021 at the pixel and landscape scale in the Chocó Andino Biosphere Reserve (RBCA).
Table 4. Forest cover, AGC stocks, and AGC net change between 2018 and 2021 at the pixel and landscape scale in the Chocó Andino Biosphere Reserve (RBCA).
Variable20182021Productivity
AGC (Mg)13,446,45113,662,109
Forest cover (ha)198,727194,795
Landscape (Mg C yr−1) 71,886
Pixel (Mg C ha−1 yr−1) 0.37
Table 5. Models’ error assessment based on different datasets for the year 2021 based on the linear fit between the observed above-ground carbon values against the random forest modelled estimates at 10 m pixel resolution, where: AGC: above-ground carbon; RMSE: Root Mean Square Error; MAE: Mean Absolute Error; and MAPE: Mean Absolute Percentage Error.
Table 5. Models’ error assessment based on different datasets for the year 2021 based on the linear fit between the observed above-ground carbon values against the random forest modelled estimates at 10 m pixel resolution, where: AGC: above-ground carbon; RMSE: Root Mean Square Error; MAE: Mean Absolute Error; and MAPE: Mean Absolute Percentage Error.
Total Plots (n)TypeObserved AGCAdj. R2RMSEMAEMAPE (%)
(Mg ha−1)(Mg C ha−1)(Mg C ha−1)
494Entire dataset85.70.8712.18.290.12
393Training dataset85.80.958.36.070.09
100Validation dataset85.10.620.615.710.24
23RBCA entire dataset84.50.723.817.390.23
8RBCA validation dataset810.434.628.620.39
Table 6. Models’ error assessment in four areas of interests (AOIs) with different number of local plots, where: AGC: above-ground carbon; RMSE: Root Mean Square Error; NRMSE; Normalized RMSE; and MAPE: Mean Absolute Percentage Error. Minimum plots refer to the minimal number of local plots added to the dataset that ensured a mean accuracy of ≥85%.
Table 6. Models’ error assessment in four areas of interests (AOIs) with different number of local plots, where: AGC: above-ground carbon; RMSE: Root Mean Square Error; NRMSE; Normalized RMSE; and MAPE: Mean Absolute Percentage Error. Minimum plots refer to the minimal number of local plots added to the dataset that ensured a mean accuracy of ≥85%.
AOITotal Plots (n)Minimum Plots (n)Observed
AGC (Mg ha−1)		RMSE
(Mg C ha−1)		NRMSEMAPE (%)
RBCA232384.5 (40.7)18.4 (3.7)0.22 (0.044)16.4 (0.01)
Bolivia231982.9(44)19.9 (4.7)0.24 (0.06)14.8 (0.02)
Peru281487.6(34.3)16.3 (2.7)0.19 (0.03)15.4 (0.03)
Bolivia and Peru673590 (36.3)21.5 (2.6)0.24 (0.03)15.2 (0.01)
Table 7. Estimated above-ground carbon (AGC) stocks in the Chocó Andino Biosphere Reserve along the elevation gradient by the three global datasets and this study, using 2019 as a reference midpoint of the period studied (2018–2021). The six elevation bins used resemble the ecosystems with distinct carbon contents (based on [46]).
Table 7. Estimated above-ground carbon (AGC) stocks in the Chocó Andino Biosphere Reserve along the elevation gradient by the three global datasets and this study, using 2019 as a reference midpoint of the period studied (2018–2021). The six elevation bins used resemble the ecosystems with distinct carbon contents (based on [46]).
Elevation (m a.s.l.)[22][19][65]This Study
AGC (±Sd)
Mg ha−1		Landcape
(Mg C) 		AGC (±Sd)
Mg ha−1		Landcape
(Mg C) 		Mean (±Sd)SumAGC (±Sd)
Mg ha−1		Landcape
(Mg C)
340–120079.05 (±12.9)2,692,02388.57 (±23.2)4,008,94291.7 (±36.9)3,895,62380.01 (±6.19)3,359,412
1200–180094.02 (±25.3)5,909,56697.53 (±32.4)6,996,586100.2 (±42.4)6,391,24372.98 (±6.27)4,514,733
1800–2500122.23 (±30.7)6,507,414125.74 (±41.9)7,270,708113 (±42.2)6,040,24665.87 (±4.18)3,273,434
2500–3000106.3 (±35.9)2,920,852106.93 (±52.9)3,549,32399.1 (±49.3)2,815,01663.82 (±4.31)1,641,145
3000–350094.6 (±34.1)1,352,74481.42 (±45.3)1,463,63882.5 (±50.2)1,245,40162.45 (±4.16)856,181
3500–477771.5 (±27)91,70045.87 (±29.1)89,00043.51 (±39.2)43,37359.78 (±3.75)59,682
All100.7 (±31.7)19,474,30089.13 (±47.2)23,378,19799.87 (±44.1)20,430,90370.65 (±8.2)13,706,235
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Cuesta, F.; Calderón-Loor, M.; Rosero, P.; Miron, N.; Sharf, A.; Proaño-Castro, C.; Andrade, F. Mapping Above-Ground Carbon Stocks at the Landscape Scale to Support a Carbon Compensation Mechanism: The Chocó Andino Case Study. Forests 2023, 14, 1903. https://doi.org/10.3390/f14091903

AMA Style

Cuesta F, Calderón-Loor M, Rosero P, Miron N, Sharf A, Proaño-Castro C, Andrade F. Mapping Above-Ground Carbon Stocks at the Landscape Scale to Support a Carbon Compensation Mechanism: The Chocó Andino Case Study. Forests. 2023; 14(9):1903. https://doi.org/10.3390/f14091903

Chicago/Turabian Style

Cuesta, Francisco, Marco Calderón-Loor, Paulina Rosero, Noam Miron, Andrei Sharf, Carolina Proaño-Castro, and Felipe Andrade. 2023. "Mapping Above-Ground Carbon Stocks at the Landscape Scale to Support a Carbon Compensation Mechanism: The Chocó Andino Case Study" Forests 14, no. 9: 1903. https://doi.org/10.3390/f14091903

APA Style

Cuesta, F., Calderón-Loor, M., Rosero, P., Miron, N., Sharf, A., Proaño-Castro, C., & Andrade, F. (2023). Mapping Above-Ground Carbon Stocks at the Landscape Scale to Support a Carbon Compensation Mechanism: The Chocó Andino Case Study. Forests, 14(9), 1903. https://doi.org/10.3390/f14091903

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop