Estimation of Above Ground Biomass in a Tropical Mountain Forest in Southern Ecuador Using Airborne LiDAR Data

Abstract

A reliable estimation of Above Ground Biomass (AGB) in Tropical Mountain Forest (TMF) is still complicated, due to fast-changing climate and topographic conditions, which modifies the forest structure within fine scales. The variations in vertical and horizontal forest structure are hardly detectable by small field plots, especially in natural TMF due to the high tree diversity and the inaccessibility of remote areas. Therefore, the present approach used remotely sensed data from a Light Detection and Ranging (LiDAR) sensor in combination with field measurements to estimate AGB accurately for a catchment in the Andes of south-eastern Ecuador. From the LiDAR data, information about horizontal and vertical structure of the TMF could be derived and the vegetation at tree level classified, differentiated between the prevailing forest types (ravine forest, ridge forest and Elfin Forest). Furthermore, topographical variables (Topographic Position Index, TPI; Morphometric Protection Index, MPI) were calculated by means of the high-resolution LiDAR data to analyse the AGB distribution within the catchment. The field measurements included different tree parameters of the species present in the plots, which were used to determine the local mean Wood Density (WD) as well as the specific height-diameter relationship to calculate AGB, applying regional scale modelling at tree level. The results confirmed that field plot measurements alone cannot capture completely the forest structure in TMF but in combination with high resolution LiDAR data, applying a classification at tree level, the AGB amount (Mg ha⁻¹) and its distribution in the entire catchment could be estimated adequately (model accuracy at tree level: R² > 0.91). It was found that the AGB distribution is strongly related to ridges and depressions (TPI) and to the protection of the site (MPI), because high AGB was also detected at higher elevations (up to 196.6 Mg ha⁻¹, above 2700 m), if the site is situated in depressions (ravine forest) and protected by the surrounding terrain. In general, highest AGB is stored in the protected ravine TMF parts, also at higher elevations, which could only be detected by means of the remote sensed data in high resolution, because most of these areas are inaccessible. Other vegetation units, present in the study catchment (pasture and subpáramo) do not contain large AGB stocks, which underlines the importance of intact natural forest stands.

1. Introduction

Global forests cover ~31% of the earth surface [1] and play an important role in the global carbon (C) cycle as well as for the climate system [2]. Forests are the most important C stocks and sinks because they store 70–90% of the terrestrial carbon (e.g., [3,4]). However, forested areas are globally endangered by human activities, such as deforestation, which is mainly caused by enhanced population pressure, especially in tropical countries [5]. As the United Nations Framework Convention on Climate Change [6] indicates, the principal cause of today’s deforestation is agriculture due to global food demand, as well as clear-cutting to make space for timber plantations [7].

This also holds true for the low-latitude tropical forests, covering ~10% of the Earth’s surface [8,9], where deforestation and land use changes have a more corrupting influence on ecosystems than general global warming [10]. Deforestation accounts for 9% to 11% of global anthropogenic greenhouse gas (GHG) emission [11] and therefore is considered as the second major source of anthropogenic CO₂ present in the atmosphere [3,11]. However, the tropical forests still contain ~59% of the total C stored in the global forest biomass [12] but large portions of it are released into the atmosphere due to human activities.

Therefore, the determination of the actual C stocks in the Above Ground Biomass (AGB) of tropical forests is of utmost importance because the storage is rapidly altered by the high deforestation rates observed in tropical countries (e.g., [13]). As Berenguer et al. [14] stated, disturbed forests only store 40% of the AGB compared to undisturbed forests, while the actual C stocks in primary tropical forests must be quantified to evaluate possible future emissions and to illustrate the importance of intact natural forest ecosystems as stocks and sinks for the increasing CO₂ concentration in the atmosphere [15].

For this task, the AGB must be monitored in space and time, which is particularly challenging in Tropical Mountain Forests (TMF; elevation >1000 m above sea level), because of the more complex terrain in comparison to lowland tropical forests [16]. The fast-changing forest structure in tropical mountains is hardly detectable by field measurements, where the forest AGB is estimated by means of individual tree samples taken in relatively small field plots [8,12,17]. Field observations traditionally estimate the mean forest AGB manually by means of random sampling of felled trees in the forest ecosystems [18]. The tree material is weighted in situ and the obtained values extrapolated to calculate mean AGB for the whole forest stand. More recent field methods use tree cores taken with increment borers to calculate the average Wood Density (WD) of the forest stand, besides other measurements, such as Tree Height (H) and Diameter at Breast Height (DBH) of individual trees present in the plots, to estimate the mean AGB, applying regional scale modelling (e.g., [19,20]). For AGB estimation, different empirical models (regional scale modelling) were developed, such as the models proposed by Chave et al. and Brown [21,22]. Notably, the generic allometric equation of Chave et al. [21] provides good results for an accurate AGB estimation, because the models include WD, H and DBH, in addition to specific variables for particular forest types (e.g., [19,23]). However, forest structure and thus AGB changes rapidly at fine scale in tropical mountains, due to the local topography, climate conditions, soil types and natural or anthropogenic disturbance, modifying the number of tree individuals, H and DBH (e.g., [24,25]). Therefore, plot based studies alone can estimate area-wide AGB in open forest stands or managed forests but in dense natural or unmanaged forests, these methods are limited due to the fast-changing forest structure, difficulty of access and the high cost of large field plot implementation [26]. In consequence, field plot measurements generally under- or overestimate the AGB in dense natural forest stands [3]. As Gourlet-Fleury et al. [27] clarified, for a reliable AGB estimation, the field plot distribution as well as the spatial scale must be fine enough to capture the landscape variability, which is especially complicated in tropical high mountains, because of the relation between forest structure, topography and morphometric protection of the site [16,28,29,30].

For a reliable quantification of AGB and C stocks, particularly in tropical mountain regions, high resolution data is necessary to detect the alterations in forest structure and land-use changes for smaller areas [12,31]. This can be determined by multispectral remote sensing observations [32], which facilitate information in very high spatial resolution. However, multispectral satellite data only provides two-dimensional information from the upper forest canopy, which can be used for forest type classifications [33,34] but hardly to detect the vertical structure of the forest for AGB and C stock estimations [35,36]. Furthermore, the availability of multispectral satellite data is limited, because areas covered by cloud cannot be classified, which is especially problematic in tropical regions, where cloud frequency is extraordinarily high [9,37].

An accurate estimation of AGB requires a tree level classification [23,38,39], which can be facilitated by modern instruments, like the Light Detection and Ranging (LiDAR) sensor [40,41]. From the three-dimensional LiDAR data, the horizontal and vertical forest structure can be derived as well as the dominant trees identified [42,43]. LiDAR data permits the generation of high resolution Digital Terrain Models (DTM) and high resolution Digital Surface Models (DSM), from which Canopy Height Models (CHM) can be obtained. These data can be used to derive different topographical indices and vegetation properties, such as the height of the vegetation at tree level (e.g., [36,38,44,45]). In general, the LiDAR measurements facilitate information about forest structure at tree level as well as information about the tree canopies, which can be used to detect the forest structure at fine scale and, in combination with field plot measurements, to estimate the AGB of the entire forest stand [46].

Area-wide AGB stock estimates in tropical regions are mainly obtained from forest inventory data by means of summary statistics, calculating mean parameters of WD, H and DBH and applying allometric equations to determine local and regional AGB storage (e.g., [47]). An alternative approach was presented by Asner and Mascaro [46], using Lidar data in combination with field plot measurements and applying a power-law function of mean top canopy height to estimate the AGB and C distribution for tropical landscapes. However, the mean values of WD, H and DBH vary significantly at local and regional scale [48] as well as the number of individual trees at fine scale [24,25]. Hence, specific local mean WD and height-diameter relationships must be calculated and an individual tree classification applied to improve the final results [39,49].

Nonetheless, during recent years the knowledge about AGB storage in the tropical lowland forests was extended through the synthesis of plot inventory data and remote sensed data (e.g., [50,51,52,53]) but in TMF (elevation >1000 m above sea level) information is still needed [54], because an adequate field plot distribution as well as the availability of suitable high resolution remote sensed data is limited due to the difficulty in access, the complex topography (large elevation differences, steep slopes, etc.; [55,56,57]) and the fast changing climatic conditions [24,25]. Furthermore, to generate area-wide AGB and C maps as well as to illustrate the AGB distribution at local and regional scale, the land cover variability must be considered, identifying areas of forest, pastures and disturbance adequately [54,56].

The available information about AGB stocks and distribution in TMF is principally based on field plot data installed in undisturbed forest parts, calculating the local or regional AGB stock by multiplying the mean value per plot-area with the area of the forest region [58,59]. Although attempts were carried out to correlate the AGB with topographical and climatological variables, such as elevation (e.g., [60]), slope angle (e.g., [61]), temperature (e.g., [62]) and rainfall (e.g., [50]), the results did not indicate overall trends [59]. For the Andes in South America it is generally stated that AGB decreases with increasing elevation and slope angle [24,56] but depending upon the prevailing soil type and the nutrient availability [63,64,65]. However, Spracklen and Righelato [59] synthesized field plot measurements in TMF from different parts of the world and concluded that only a low to moderate correlation between AGB and elevation as well as slope angle, temperature or rainfall exists.

The lack of correlation between AGB and topographical as well as climatological variables was also stated in other studies (e.g., [50,66]), which indicates that these parameters only explain a fraction of the AGB variability in TMF, why additional environmental parameters must be analysed [59]. Therefore, the present study estimates the AGB and C stocks of a natural TMF inside a small watershed in the south-eastern Andes of Ecuador, using high resolution LiDAR data in combination with field plot measurements, while applying regional scale modelling at tree level. From the LiDAR data the H of the individual trees as well as the forest, terrain and land cover structure can be derived, whereas from the field plot measurements the local mean WD and the specific height-diameter relationship are obtained. The AGB stock and its distribution are analysed by means of the Topographic Position Index (TPI), which identifies ridges and depressions [30,67] and the Morphometric Protection Index (MPI), which illustrates the sheltering effect of the surrounding terrain [29]. The results are compared to AGB values calculated in previous plot-based investigations within the same area. The study also intends to demonstrate the importance of intact primary forest ecosystems as C stocks considering the ongoing deforestation, especially in tropical high mountains.

2. Materials and Methods

2.1. Study Area

The research area, the San Francisco watershed, is located in the eastern escarpment of the south Ecuadorian Andes between the cities of Loja and Zamora (Figure 1), with drainage into the Amazon Basin. The catchment is part of the Cordillera Real and has an extension of ~85 km², including the “Reserva Biológica San Francisco” (RBSF) where the research station Estación Científica San Francisco (ECSF: lat. 3°58′18′′S, long. 79°4′45′′W, 1860 m above sea level; Figure 1) is situated. The altitudes range from 1600 m above sea level at the river outlet to 3200 m above sea level at the highest mountain tops [68]. The Andes between southern Ecuador and northern Peru are characterized by lower mean altitudes compared to other parts of the mountain range (Andean Depression, the Amotape-Huancabamba Depression), only reaching 3900 m above sea level at the highest mountains. Therefore, a strong exchange between the Ecuadorian costal vegetation and the Amazon vegetation exists, leading to extraordinarily high biodiversity in this region (hottest hotspot in biodiversity; for example, References [69,70]).

The vegetation in the catchment can be divided into three classes of ecosystems: Tropical Mountain Forest (TMF), subpáramo and pastures [71]. The TMF (Figure 2a) is mostly undisturbed and covers the slopes of the southern, western and the upper parts of the northern ridge up to the tree line at ~2700 m above sea level [72]. The natural forest structure is related to the topography and can be divided into lower slope (ravine) forest and upper slope (ridge) forest (e.g., [16]). The ravine forest is characterized by lower stem density but higher basal areas (diameters) and canopy heights compared to the ridge forest, where less tree species are also observed [25]. The difference in forest structure is mainly due to the climatic conditions and the prevailing soil types [73,74,75]. As identified by Homeier et al. [76], the most abundant tree genera are Miconia (Melastomatac), Ocotea (Laurac), Persea (Laurac), Ficus (Morac), Weinmannia (Cunoniac), Ilex (Aqufoliac), Hedyosmum (Chorantac), Inga (Fabac), Scheffera (Araliac) and Clusia (Clusiac).

Between ~2700 m and ~3000 m above sea level a transition zone between the TMF and a shrub-dominated subpáramo exists, where patches of Elfin Forest (Figure 2b), the so-called Ceja Andina, dominate the landscape. The Elfin forest is characterized by small trees, reaching heights up to 10 m [28]. Over 3000 m above sea level a shrub-dominated subpáramo (Figure 2c) prevails, which is the lowest and most diverse altitudinal zone of the páramo area, where aspects of the grass páramo (above) and the mountain forest (below) are mixed [77].

The pastures (Figure 2d) cover mainly the lower parts of the northern mountain ridge, where the natural mountain forest was replaced by anthropogenic slash-and-burn activities [71,78]. However, some parts of the pastures have already been abandoned and dense stands of bracken fern, which are resistant to burn activities [79], have infested these areas [68,80].

The climate in the study catchment is per-humid with marked altitudinal gradients in air temperature, relative humidity [81,82] and rainfall [78,83]. The annual mean air temperature ranges from 19.4 °C at the valley bottom to 9.4 °C at the mountain tops, whereas annual mean relative humidity varies between 77% at the pasture sites and almost 100% inside the TMF as well as at the mountain ridges. The annual rainfall amounts reach 2300 mm near the valley bottom and 6700 mm at the mountain tops. These values of annual precipitation include rain and fog water precipitation, because clouds and fog deposit water directly onto the vegetation and therefore must be considered as relevant available water input from the atmosphere [84]. The seasonal rainfall distribution of the watershed shows a clear annual cycle with the main rainy season in austral winter (between May and September) and a relative dry season in austral summer (between November and February).

2.2. LiDAR Data and Field Measurements

An airborne Light Detection and Ranging (LiDAR) survey in the study area was conducted in two campaigns in March and November 2012, due to the adverse climate conditions at the beginning and during the rainy season. The LiDAR data was obtained from a Leica Geosystem ALS50-II laser scanner installed on a Eurocopter AS350B2 Ecureil helicopter.

Table 1. Leica Geosystems ALS-50 II CM specifications.

Parameter	Value
Scanning pattern	Sinusoid
Field Of View, deg	0–75
Pulse Rate (maximum), kHz	150
Pulse Wavelength, nm d	1064
Scan Rate (maximum), Hz	90
Number of returns	4

shows the principal specifications of the equipment used. The Airborne Laser Scanning (ALS) produced point clouds with a density of at least 10 pulses per 1 m² for the entire watershed and a buffer around it, to ensure that all boundaries of the catchment were covered (total area: 108.93 km²). The point clouds of the study area were subdivided into 257 tiles to simplify the handling of the information. For more information about the LiDAR survey, please refer to Silva and Bendix [85].

The forest inventory data used includes measurements of WD (tree cores), H (ultrasonic tree height meter) and DBH (dendrometer) of the tree species within the two different mountain forest types (ravine and ridge) along an altitudinal gradient from 1900 m to 2100 m above sea level (see Figure 1) [76,86]. The data was collected between 2010 and 2014 under dense canopies without larger forest gaps (<2 m) and consists of 540 tree samples from 18 plots of 0.04 ha [24,25]. To consider the local topography as well as the general forest structure, six plots were situated at the upper slope, six plots at the mid-slope and six plots at the lower slope.

2.3. Methods

The processing chain to estimate AGB and C stocks of the TMF as well as its distribution influenced by the forest structure is shown in Figure 3. AGB and C stock were only estimated for the TMF and the Elfin Forest by means of the LiDAR data in combination with field measurements, considering only trees higher than 5 m with a DBH greater than 10 cm, as recommended by Li et al. and Gianico et al. [87,88]. For the other vegetation units (pastures and subpáramo) the AGB and respectively C stock values from previous investigations inside or next to the study area were used [89,90].

The mean AGB value for the shrub-dominated subpáramo area was taken from Eguiguren et al. [90], who evaluated the C stock in the Podocarpus National Park, next to the study area. They established 25 plots of 1 m² at different altitudes and soil types, estimating a mean C stock of 8.2 Mg C ha⁻¹ for shrub-dominated areas. The C stock can be converted to AGB, applying a factor of 2, as Tan et al. and Ota et al. [33,44] demonstrated, which results in a mean AGB of 16.4 Mg ha⁻¹ for the subpáramo area.

AGB values for recently abandoned tropical pastures (1–2 years) range between 2 Mg ha⁻¹ and 18 Mg ha⁻¹ with higher values for pasture with relict trees [8,79]. However, Knoke et al. [89] investigated C stocks for tropical pasture in southern Ecuador, based on CLM-DGVM modelling [91] and calculated values between 12.5 Mg C ha⁻¹ (low-input pastures) to 33.0 Mg C ha⁻¹ (long time abandoned pastures), which includes the whole C stock of the plant (roots and shoot). In the study catchment most of the pastures are in use and relict trees are generally isolated, in addition to a negligible pasture management (low-input pastures, own observations). Therefore, the C stock for low-input pastures was converted into AGB applying the shoot-root ratio for C4 grasses (here: mainly Setaria sphacelata; [79]), which resulted in a mean value of 10.8 Mg ha⁻¹ for the pasture site.

For an accurate AGB estimation different equations exist (e.g., [19,22]), in which one frequently used approach is the regional scale model of Chave et al. [21], who presented different equations for specific tropical forests types (dry forest, moist forest, moist mangrove forest and wet forest). Due to the climate conditions in the study area, the equation for wet tropical forest stands was selected, which was also applied in previous plot based AGB estimation studies inside the San Francisco watershed [24,25,73] and in plot-based AGB estimations in other TMF in southern Ecuador [8]. However, due to the fast-changing forest structure in the TMF [25] a calculation at tree level is necessary [23,38,39], which can be realized by means of the LiDAR data in combination with the field measurements [40,41]. The equation to estimate the AGB at tree level can be written as follows:

$${\mathrm{AGB}}_{\mathrm{tree}}=\exp (-2.557+0.940\times \ln ({\mathsf{\rho }\mathrm{D}}^2\mathrm{H}))=0.0776\times {({\mathsf{\rho }\mathrm{D}}^2\mathrm{H})}^{0.940}$$

(1)

where AGB_tree is the AGB of a specific tree (Mg), ρ is the average WD (gr cm⁻³), D is the particular DBH (cm) and H the particular height of the tree (m).

The WD was determined by the extraction of tree cores (increment borer) from the species present in the installed field plots, following the method described by [18]. This method is non-destructive and recommended for tropical forest trees. The tree cores were taken at breast height and the WD of the individual trees determined by the ratio of the oven-dry wood mass and the green wood volume [76]. However, WD varies between the tree species and during the lifetime of the tree as well as within the same tree because the wood of the branches as well as the different parts of the trunk (interior and exterior) present different values [23]. To solve this problem, mean WD for the whole forest stand was calculated, averaging all individual WD from the trees sampled [8].

Due to the selective forest inventory (plot based), H and DBH field measurements were only used to establish a height-diameter relationship [92], because obtaining these values for all individual trees inside the catchment is difficult to realize. In general, height-diameter relationships vary significantly between regions because of the species composition and local climate conditions [49]. To determine the local relationship 244 tree measurements were used and the remaining trees sampled were preserved for validation. To the selected H and DBH data a logarithmic transformation was applied, as suggested by Gianico et al. and Djomo and Chimi [88,93] and then a linear regression analysis executed (Figure 4; [94]). The obtained local height-diameter relationship for this TMF is shown in Figure 4.

The calculated coefficient of determination (R²) indicates a low positive correlation between DBH and H (R² = 0.60), which is acceptable, due to fast changing topography and soil types, modifying the number of tree individuals, H and DBH [24,25,49]. Nevertheless, the calculated p-value (<0.001) confirmed a high significance of the relationship between H and DBH for the TMF in the study catchment.

The height-diameter relationship equation (Figure 4) was transposed to calculate the DBH (Equation (2)), because H was deviated from the LiDAR data.

$$\mathrm{DBH}={\mathrm{e}}^{(\frac{\mathrm{Ln}\left(\mathrm{H}\right)-0.85}{0.56})}$$

(2)

In general, the LiDAR point clouds (raw data) facilitate information about the vertical and horizontal structure of the terrain and the vegetation [35]. From these data a DTM as well as a DSM can be generated, from which a CHM can be extracted [95]. Based on the CHM, the individual trees were detected and their specific height determined.

For this, the LiDAR point clouds were first corrected and processed by means of the software FUSION (Version 3.7), developed by the U.S. Department of Agriculture, Forest Service [42]. To assure the quality of the derived products, the point clouds were analysed with the Catalog tool, which evaluates several important metrics of the LiDAR data and the completeness of data coverage as well as the pulse return density. Then, outliers within the tiles were detected and eliminated, executing the FilterData extension of FUSION. Therefore, the mean elevation and its standard deviation within a user-defined grid (here: 10 m × 10 m) was calculated and the individual point compared to the mean value. If the point value in the established grid cell was higher or lower than the mean elevation, adding 3 times the standard deviation, the point was eliminated as an outlier [38]. This tool detects the presence of electricity towers, power lines and so forth, which cause errors in the subsequent individual tree classifications because these objects are not part of the natural surface or the vegetation.

After the outlier elimination, the GroundFilter tool of FUSION was applied to determine the lowest elevations of each pulse return (bare-earth points) as well as the highest elevation of each pulse return (canopy height). Therefore, all cloud points are screened to identify the respective extreme values [42]. By means of these data a DTM was created, using the GridSurfaceCreate extension, which averages the bare-earth points for each grid cell. However, Jochem et al. [35] indicated that for an adequate DTM generation, different resolutions should be evaluated, because the optimal resolution depends on the requirements of the final product [45]. The analysis resulted in a final cell size of 0.25 m × 0.25 m, because with lower resolution the individual trees were not classified correctly [45,96].

Then, a DSM was generated, executing the CanopyModel extension of FUSION, which determine the canopy height. Therefore, the highest elevation pulse returns were identified and averaged to obtain a single value for the established grid cell (0.25 m × 0.25 m). Subsequently, a CHM, which displays the above ground vegetation heights, was created, subtracting the generated DTM and DSM maps. For the tree detection process the CHM was smoothed, applying user defined mean and median filters (here: 5 × 5) integrated in this extension. Thereby, the local maxima are preserved while the surrounding cells (here: crowns) smoothed to adhere the highest point of a tree [42,96]. The individual trees were identified by means of the CanopyMaxima tool of FUSION, which detects local maxima by means of a variable-size evaluation window (here: 3 × 3) [97]. The tool identifies H and the location of dominant/codominant trees within the forest stand because mid-story and under-story trees are often covered by a higher canopy layer [42]. To evaluate the tree detection process for the study area, the number of detected trees as well as their H and their derived DBH were compared to the field plot measurements, which were not used for the determination of the local height-diameter relationship.

Based on the determined mean WD, the deviated DBH and the H of each detected tree, Equation (1) was applied to calculate the AGB at tree level. To calculate AGB (Mg ha⁻¹) and C stocks (Mg C ha⁻¹) of the TMF, first a land cover classification had to be performed, using the generated CHM (0.25 m × 0.25 m resolution) to determine the areas where forest is present. Therefore, a grid layer with a 1 ha × 1 ha resolution was created and overlaid with the CHM. Then, each hectare grid was analysed and classified as forest if more than 50% of the included CHM pixels were higher than 5 m [87]. For non-forest grids, the AGB values for pasture (10.8 Mg ha⁻¹) and subpáramo (16.8 Mg ha⁻¹) were assigned, depending on the altitude (threshold 2700 m above sea level). To obtain AGB for the TMF, all individual AGB at tree level within a grid cell (1 ha × 1 ha) were added up [98]. Finally, the C stock was calculated by multiplying the AGB with a factor of 0.5 [33,44].

In order to analyse the relationship of the complex topography, the forest structure and the AGB distribution in the study area, the Topographic Position Index (TPI), which identifies ridges and depressions [30,67] and the Morphometric Protection Index (MPI), which illustrates the sheltering effect of the surrounding terrain or the exposure of the cell [29], were calculated by means of the DTM (0.25 m × 0.25 m resolution), using the software SAGA GIS [99,100]. The TPI compares the elevation of each grid cell to the mean elevation within a specific distance to determine whether the position of the cell is situated on a ridge (positive values) or in a valley (negative values), which permits the separation into ravine and ridge forest [101]. The MPI evaluates, up to a given distance, whether the respective cell is protected by the surrounding relief to illustrate the topographic preference for bigger trees within the forest stands [29]. MPI-values near zero indicate unprotected locations, whereas higher positive values indicate the degree of protection of a grid cell within the catchment. In this study the distance for both indices was set to 100 m, as Detto et al. [102] recommended, due to the fast-changing topography, because canopy heights are strongly related to the proximity to depressions and the exposure of the site [29]. Furthermore, this distance also warrants the detection of smaller subsidiary valleys and ridges.

To illustrate the importance of the field plot distribution and size to capture the landscape variability for a reliable AGB estimation in natural TMF [27,103], the DBH of each detected tree was related to TPI and MPI. Therefore, a threshold of 20 cm (DBH) was defined, as Bunning et al. [104] suggested, to analyse the conditions of forest trees. The results were related to the AGB distribution inside the San Francisco catchment, because dominant or bigger trees within a forest stand represent 70% to 90% of the total AGB [41,105,106].

3. Results

To estimate the AGB at tree level (Equation (1)) WD, H and DBH of each individual tree is necessary. The mean WD of the TMF, based on the field measurements, was determined at 0.59 g cm⁻³, which was applied to every tree detected in the study catchment. The individual trees were identified by means of the CHM (Figure 5a) derived from the LiDAR data (DTM–DSM), which resulted in a total of 1,932,188 trees detected in the whole study catchment, including their position and H (Figure 5b). The H of each detected tree was used to calculate their specific DBH, applying Equation (2).

To evaluate these results, the specific H and DBH of the identified trees were compared to the field plot inventory (18 plots). As explained in Section 2.3, only the H and DBH measurements, which were not used for the determination of the local height-diameter relationship, were included for validation. During the field campaign 540 trees higher than 5 m with a DBH greater than 10 cm were sampled within the 18 plots. The tree detection process, based on the LiDAR data, found 319 dominant/codominant trees, which corresponds to 55% of the trees present in the different plots (see example Figure 5b). The H of the detected trees were compared to the H measurements, realized by means of an ultrasonic tree height meter in the field (Figure 6a) and their derived DBH (Equation (2)) to the dendrometer measurements (Figure 6b). The results indicated a very good correlation for H and DBH with a coefficient of determination (R²) of 0.92 (p-value < 0.001) and 0.93 (p-value < 0.001) respectively. However, the H of smaller trees up to 15 m were slightly underestimated as well as trees higher than 25 m, whereas trees with H between 15 m and 25 m were slightly overestimated (1:1 line; Figure 6a). The estimated DBH of the individual trees showed the same behaviour, because DBH of trees up to 25 cm as well as over 60 cm were slightly underestimated, whereas trees between 25 cm and 60 cm overestimated, caused by the applied local height-diameter relationship (Figure 6b).

To integrate the forest structure in the AGB estimation, the TPI (Figure 7a) as well as the MPI (Figure 7b) were calculated. These indices were used to identify the areas of ravine forest and ridge forest and to illustrate the protection of the site, because from these parameters the forest structure and the AGB distribution in the catchment can be deduced [29,101]. Therefore, a threshold of 20 cm (DBH) was set to depict the position of bigger trees within the forest stand [104]. All trees with respectively higher and lower DBH were marked in a different colour (Figure 7c). It is clearly visible that smaller trees with concurrently lower basal area (DBH) are mainly situated at unprotected sites, such as ridges or at slopes near the ridges (red polygons), whereas bigger trees cover the side valleys down to the valley bottom, protected by the surrounding relief (green polygons). The same distinction was made for the Elfin Forest areas above 2700 m above sea level (Figure 7c), where bigger trees are situated inside the upper side valleys and smaller, more isolated, trees at the ridges.

The ravine and ridge forest classifications were used to analyse the mean and extreme values (H, DBH and AGB) at tree level within the San Francisco watershed (

Table 2. Mean and extreme values of H, DBH and AGB for individual trees in the TMF and the Elfin Forest.

Forest Type	Tree Individuals	H [m]			DBH [cm]			AGB Tree [Mg]
TMF	1,638,767	Min	Mean	Max	Min	Mean	Max	Min	Mean	Max
Ravine	1,470,493	8.5	13.3	48.4	10.0	23.5	223.3	0.0	0.3	47.1
Ridge	168,274	8.5	11.3	26.8	10.0	17.1	77.8	0.0	0.1	3.7
Elfin Forest	293,421	Min	Mean	Max	Min	Mean	Max	Min	Mean	Max
Ravine	271,717	8.5	11.0	15.7	10.0	16.4	30.1	0.0	0.1	0.4
Ridge	21,704	8.5	8.9	11.5	10.0	12.9	17.1	0.0	0.0	0.1

). For the ravine TMF (up to 2700 m) an average H of 13.3 m was calculated, whereas the trees in ridge TMF only had an average H of 11.3 m. In the Elfin Forest (over 2700 m) average H is generally lower, reaching 11.0 m in the ravine forest and 8.9 m in the ridge forest parts. The biggest trees were found in the protected ravine TMF areas, reaching heights of up to 48.4 m, whereas maximum H in the ridge TMF reached only 26.8 m at lower and more protected sites. The maximum H in the Elfin Forest is notably lower but the biggest trees are also situated at more protected areas inside the upper side valleys, reaching heights up to 15.7 m (Table 2).

DBH was derived from H (Equation (2)), which means that bigger trees concurrently have a greater basal area. The calculated DBH of the trees in the San Francisco catchment ranges between 10.0 cm and 223.1 cm in the ravine TMF (average 23.5 cm) and between 10.0 cm and 77.8 cm in the ridge TMF (average 17.1 cm: Table 2). In the Elfin Forest, tree DBH is notably smaller, reaching maximum values up to 30.1 cm in the ravine forest (average 16.4 cm) and 17.1 cm at the upper ridges (average 12.9 cm; Table 2). The biggest trees were found in the lower ravine TMF near the river outlet as well as at more protected lower ridges.

H and DBH are directly correlated to the AGB (Equation (1), [21]) and therefore AGB per tree vary notably within the different forest types. The mean AGB per tree in the ravine TMF was estimated at 0.3 Mg, whereas the mean AGB per tree in the ridge TMF resulted in 0.1 Mg. For the ravine Elfin Forest mean AGB per tree also resulted in 0.1 Mg but for the ridge Elfin Forest a mean AGB of 0.0 Mg was determined. Maximum AGB per tree was obtained for the lower ravine TMF, indicating values up to 47.1 Mg, whereas for the biggest trees in ridge Elfin Forest only an AGB of 0.1 Mg was estimated (Table 2). To evaluate these findings, the estimated AGB per tree were compared to the AGB per tree calculated from the plot data inventory (Figure 8). It is clearly visible that the AGB per tree is generally overestimated, which is mainly due to the applied average WD value, because Werner and Homeier [25] applied the specific WD to each individual tree species and the established local height-diameter relationship, which slightly overestimates the DBH of the individual trees (see Figure 6b). However, the coefficient of determination (R² = 0.91; p-value < 0.001) indicated a very good correlation between measured and estimated AGB values per tree.

The minimum H, DBH and AGB values per tree were equal for all forest types, due to the thresholds applied in this study, because only trees higher than 5 m with a DBH greater than 10 cm were considered (see Section 2.3 and Table 2; [87,88]).

To illustrate the general forest structure within both forest types (ravine and ridge), the individual H (Figure 9a,b) and DBH (Figure 9c,d) of the trees were integrated into classes and the percentage of each class respective to the total number of trees calculated (Figure 9). The ravine TMF is generally characterized by trees with H between 10.0 m and 22.5 m, representing more than 97.2% of this forest type. The trees in the ravine Elfin Forest are notably smaller, reaching only heights up to 17.5 m, in which trees up to 15.0 m represent 97.3% (Figure 9a). H in the ridge TMF generally do not pass 20.0 m (99.2% of the forest stand), in which 91.8% of the trees only reach heights up to 15.0 m. In the ridge Elfin Forest, trees predominantly of up to 10.0 m can be found (62.5%), due to the adverse climate conditions (low temperatures [81] and strong winds [75]); biggest trees with H up to 12.5 m only exist at the more protected upper ridges (Figure 9b).

Respective to the DBH, the ravine TMF is generally characterized by trees with a DBH between 10.0 cm and 50.0 cm, representing more than 95.1% of this forest type. Most of the trees reach DBH up to 40.0 cm (89.9%), whereas the biggest trees (DBH > 50.0 cm) only occupy 4.9%. The DBH of trees in the ravine Elfin Forest do not exceed 30.0 cm, in which trees with a DBH up to 25.0 cm represent 92.7% of the forest type (Figure 9c). In the ridge TMF, trees are smaller compared to the ravine TMF, reaching mostly DBH up to 35.0 cm (97.0%). Trees with a DBH greater than 45.0 cm are generally absent in this forest type with some exceptions at the lower highly protected ridges. In the ridge Elfin Forest only trees with a DBH up to 20.0 cm exist (100.0%) but the majority reach DBH around 15.0 cm (82.6%; Figure 9d).

Figure 10 illustrates the percentage of the different tree-height classes respective to total AGB value of the two forest types. In the ravine TMF the trees with heights between 15.0 m and 30.0 m contribute 82.7% of the total AGB (Figure 10a), representing only 48.8% of the trees present there (see Figure 9a). Trees smaller than 15.0 m only add 10.6% of the total AGB, standing for 51.0% of the trees, whereas trees bigger than 30 m account for 6.7% of the total AGB, representing only 0.2% of this forest type. In the ravine Elfin Forest, the situation is similar, because trees bigger than 15.0 m contribute 53.5% of the total AGB, accounting only for 25.2% of all trees present there. Trees up to 15.0 m are most frequent in this forest type (74.8%, Figure 10a) but their contribution to the total AGB is only 46.5% (Figure 9a).

The same distribution was found in the ridge forest areas (Figure 9b). In the ridge TMF trees higher than 15.0 m account for 63.7% of the total AGB, representing 24.1% of the trees within this forest type (Figure 10b). Only 2.6% of trees are taller than 20 m but these trees contribute 20.9% to the total AGB. Trees up to 15.0 m represent 75.9% of this forest stand, which only add 36.3% to the total AGB. Also in the ridge Elfin Forest, where trees higher than 15.0 m are absent (Figure 9b), the biggest trees (taller than 12.5 m) contribute the majority of the total AGB (52.2%), representing only 37.5% of the trees; whereas trees of up to 12.5 m only account 47.8% for the total AGB in this forest type (Figure 10b).

To obtain an AGB value (Mg ha⁻¹) for areas covered by forest, a land cover classification was initially performed (Figure 11). Therefore, a grid mask of 1 ha × 1 ha was overlaid with the CHM and all grids classified as forest if more than 50% of the included CHM pixels were higher than 5 m [87]. Then, all AGB values of the individual trees within a grid cell were added up. Due to the complex terrain in the study area, it was only possible to analyse the differences between TMF and Elfin Forest, whereas a distinction between ravine and ridge forest areas was unfeasible.

The TMF covers 4608 ha of the San Francisco watershed, whereas the Elfin Forest only 1529 ha. The remaining hectares were classified as pasture (10.8 Mg ha⁻¹) or subpáramo (16.4 Mg ha⁻¹), depending on the elevation (threshold 2700 m), where the respective AGB values was set. The calculated mean AGB of the TMF was 106.2 Mg ha⁻¹, which is considerably higher compared to the Elfin Forest where only 32.8 Mg ha⁻¹ was estimated (

Table 3. AGB and C stock in the San Francisco catchment.

		AGB				C Stock
		[Mg ha−1]				[Mg ha−1]
Forest Type	N [ha]	Min	Mean	Max	SD	Min	Mean	Max	SD
TMF *	4608	10.0	106.2	664.1	94.1	5.0	53.1	332.0	47.0
Elfin Forest	1529	2.1	32.8	196.6	28.8	1.1	16.4	98.3	14.4

* TMF = Tropical Mountain Forest.

). Maximum AGB for the TMF was 664.1 Mg ha⁻¹, found at the lower ravine TMF near the river outlet, where undisturbed and dense forest stands exist, highly protected by the surrounding terrain (see Figure 6); whereas minimum AGB was 10.0 Mg ha⁻¹, displayed at the upper ridges (Figure 12). For the Elfin Forest, maximum AGB (196.6 Mg ha⁻¹) was estimated at protected areas inside the upper side valleys and minimum values (2.1 Mg ha⁻¹) at the exposed upper ridges (Table 3). The spatial distribution of AGB in the study catchment is shown in Figure 12.

As mentioned before, AGB can be converted directly into C stock, applying a factor of 0.5 [33,44]. C stock is principally stored in the TMF ecosystem (average 53.1 Mg C ha⁻¹), especially in the ravine forest inside the side valleys and at the valley bottom, where undisturbed and dense forest stands exist (max. 332.0 Mg C ha⁻¹, Table 3). The Elfin Forest (average 16.4 Mg C ha⁻¹, Table 3), the pastures (5.4 Mg C ha⁻¹; [89]) as well as the subpáramo area (8.2 Mg C ha⁻¹; [90]) in no way approximate this high value.

The total AGB for the whole study catchment was estimated at 565.3 Gg (C stock: 282.6 Gg C), in which the TMF contributes 489.3 Gg (C stock: 244.7 Gg C), 86.6% of the total amount, covering only 56.7% of the catchment area (

Table 4. Portion of the individual vegetation units respective to the total catchment area, including total AGB and C stock.

Land Cover	Area [ha]	Area [%]	AGB [Mg]	C Stock [Mg C]	C Stock [%]
TMF	4608	56.7	489,343.6	244,671.8	86.6
Elfin Forest	1529	18.8	50,117.2	25,058.6	8.9
Pasture	1200	14.8	12,960.0	6480.0	2.3
Subpáramo	785	9.7	12,874.0	6437.0	2.3
TOTAL	8122	100.0	565,294.8	282,647.4	100.0

). The Elfin Forest covers 18.8% of the catchment area and contributes 50.1 Gg (C stock: 25.0 Gg C), which is an 8.9% portion of the total AGB (C stock). The pasture sites only store 13.0 Gg (6.5 Gg C), representing 2.3% of the total AGB (C stock) amount, covering 14.8% of the catchment area. The subpáramo areas also account for 2.3% of the total AGB (12.9 Gg; C stock 6.4 Gg C) but only cover 9.7% of the watershed (Table 4). Taking only the natural forest ecosystems inside the study catchment into account, the TMF and the Elfin Forest contribute nearly 95.5% of the total AGB or C stock (Table 4).

4. Discussion

The separation into ravine forest and ridge forest within the study catchment was also described in earlier investigations (e.g., [16,24,86]), which characterized the ravine forest by larger trees with concurrently higher basal areas (DBH) compared to the ridge forest, while also in the ridge forest, fewer tree species are observed. These investigations concluded that the forest structure is mainly due to the adverse climate conditions at the ridges (stronger winds, [75]) and the different soil types [73,74,107]. However, as this study showed, bigger trees can also be found at higher elevations and steeper slopes, if the site is protected by the surrounding terrain. This could only be determined using remote sensed data in high resolution, because most of these areas are inaccessible.

The mean WD of the TMF in the study catchment was estimated at 0.59 g cm⁻³, which lies in the same range as other mean WD values reported for tropical forests in South America. For example, Chave et al. [18] presented mean WD values for Central and South America, sampling 2456 trees species and obtaining a mean value of 0.63 g cm⁻³ [15] calculated the same mean WD value (0.63 g cm⁻³) estimating Amazon forest carbon density. However, taking only the available information from Mitchard et al. [15] for Ecuador into account, the mean WD is slightly lower (0.57 g cm⁻³). This value was also applied by Gibbon et al. [108], investigating mean WD of a TMF in Peru and Spracklen and Righelato. [8], who estimated the AGB of another TMF in southern Ecuador. However, the field measurements determined WD between 0.16 g cm⁻³ and 0.92 g cm⁻³ for the different tree species [86], from which the average WD was calculated and therefore this mean value (0.59 g cm⁻³) should be more representative for the TMF in the San Francisco catchment.

The individual H in the TMF, as well as in the Elfin Forest, as obtained from the CHM and the individual tree classification, depends upon the forest type (ravine or ridge) and the degree of protection of the site. The average H values calculated in this study were similar to the values presented by Homeier et al. [86], who estimated an average tree height of 14.7 m for the ravine TMF and of 10.3 m for the ridge TMF (see Table 2). The small difference may be due to the slight underestimation of trees up to 15 m and over 25 m derived from the LiDAR data (Figure 6a). However, Werner and Homeier [25] presented notably higher average H values, estimated by means of field plot measurements at lower elevations using an ultrasonic tree height meter. They calculated an average tree height of 20.5 m for the lower ravine forest and 12.6 m for the lower ridge forest, which was confirmed by Leuschner et al. [24], who calculated 18.9 m for areas near the valley bottom and 12.0 m at higher elevations, without applying any distinction between the forest types. However, as Larjavaara and Muller-Landau [109] indicated, ultrasonic tree height meters in natural forests may produce errors of up to 20% due to the dense vegetation, which does not permit an undisturbed view up to the canopies and due to user errors. Furthermore, the mentioned variations in average H for the same forest stand are a result of the potential field plot distributions in natural TMF, which explains why the plots were installed at lower altitudes and under dense canopies without larger forest gaps (<2 m; [25]) where trees are generally taller. As Gourlet-Fleury et al. [27] clarified, the accuracy of the calculated values depends upon the plot distribution and the spatial scale, which must be fine enough to capture the landscape variability. The same is valid for the maximum H, where the published values vary between 30.0 m and 40.0 m for this TMF and between 8.0 m and 9.0 m for the Elfin Forest (e.g., [16,24,28]). However, this study identified trees of up to 48.4 m in protected lower parts of the catchment and trees up to 15.7 m in the ravine Elfin Forest. Therefore, H and forest structure in TMF are hardly detectable by small field plots because of the local topography, the difficulty in access, the climate conditions and other disturbances, which change rapidly within small scales. The forest structure is better accessed by means of an individual tree classification, which can be derived from remote sensed LiDAR data in high resolution and the derived topographic and morphometric indices (TPI and MPI; [29,67,101]).

The obtained individual H from the LiDAR data and the specific height-diameter relationship from the field data were used to calculate the DBH of each tree detected within the whole catchment (e.g., [92]). However, the accuracy of the height-diameter relationship depends on the model applied, the altitude, the climate conditions and the forest structure characteristics, as well as the location where the selected trees are sampled [49,110]. For tropical forests, only few investigations reported height-diameter relationships and the information about their performance is limited (e.g., [88,110]). In general, the height-diameter relationships vary significantly between regions because of the species composition and local climate and topographical conditions [49]. Only [93,111] presented coefficient of determination (R²) values for their height-diameter relationships, which ranges between 0.5 and 0.9. The dataset used in this study to calculate the specific local height-diameter relationship includes all tree species present in the different plots, measured within the different forest types [25]. Therefore, the calculated local height-diameter relationship (R²: 0.60; p-value < 0.001; Figure 4) is acceptable, because it is representative for the two TMF types (ravine and ridge) considering the different tree species [112]. Additionally, the logarithmic data transformation and the linear regression analysis resulted in an improved correlation between H and DBH [88,93,94]. The determined relationship slightly overestimates the DBH of the trees, especially between 25 cm and 60 cm (Figure 6b), which also leads to a small overestimation of the AGB per tree (Figure 8), in conjunction with the applied average WD.

While, comparing the estimated DBH values of this study with DBH measurements carried out by other plot-based studies in the catchment, it can be stated that the obtained results lie within the same range, because trees of up to 200 cm were identified (e.g., [16]; Table 2). The maximum DBH depends on the species and the age of the trees, as well as on the position of the trees respective to the surrounding terrain [41,113]. Average DBH in the study area, as published by Leuschner and Moser [114], are notably lower (TMF: 9.8 cm to 12.2 cm; Elfin Forest: 7.2 cm) compared to the obtained values. The difference might be explained by the plot distribution, because at inaccessible sites no data can be collected and the LiDAR data includes information of each individual tree detected. The individual tree classification applied here only considered dominant/codominant trees taller than 5 m with a DBH greater than 10 cm, in which smaller mid- and understory trees are often not detected because of higher canopies [42], which possibly increased the calculated average DBH for the forest stand. Nonetheless, the difference between TMF and Elfin Forest could be distinguished, which is due to harsher climate conditions at higher elevations, complicating the tree growth at these exposed positions (e.g., [73,74,75]). Therefore, the calculated mean values are representative for dominant/ codominant trees in the study catchment.

Highest AGB was calculated for the lower ravine TMF near the valley bottom and inside the side valleys, where protected and dense forest patches exist (664.1 Mg ha⁻¹). The AGB at the ridges are lower due to the more exposed position, which is why trees are smaller and do not accumulate much biomass and therefore do not contribute large amounts to the total AGB. In general, AGB at the ridge forests is reduced compared to the ravine forests (Table 2), because tree growth is inhibited. Additionally, the climate conditions become progressively harsher with increasing altitude [75], for which reason the Elfin Forest patches at the transition zone to the subpáramo area displayed even lower AGB (Figure 12). However, high AGB accumulations were also estimated at higher elevations inside the upper side valleys (up to 196.6 Mg ha⁻¹; Elfin Forest), where the trees are specially protected by the surrounding terrain. In general, AGB for the Elfin Forest did not exceed 50.0 Mg ha⁻¹ and mainly stayed below 20.0 Mg ha⁻¹, due to the unprotected topographic position which leads to smaller trees at these high altitudes.

The mean AGB for the TMF calculated in this study is notably lower than the values published by [16,24,25,115] for the same forest stand. These studies calculated mean AGB between 150 Mg ha⁻¹ and 200 Mg ha⁻¹, based on data from field plots of 20 m × 20 m at lower elevations under dense forest stands without larger forest gaps [25]. These field data were extrapolated to obtain a mean AGB for the whole study catchment but forest structure at higher elevations could not be included, due to the local topography, the difficulty in access, the climate conditions and natural disturbances (e.g., [3,27]). Taking into account only the trees of the lower TMF (up to 2100 m above sea level), the estimated mean AGB of this study increased to 151.6 Mg ha⁻¹, because trees are bigger at lower elevations. Above 2100 m above sea level, forest mean AGB decreased to 93.3 Mg ha⁻¹, because average H becomes lower and ridge forests dominate (Figure 7c). These values are similar to the mean AGB published by Leuschner et al. [24], who calculated 163.2 Mg ha⁻¹ at lower and 94.6 Mg ha⁻¹ at higher elevations. Therefore, Richter and Moreira-Muñoz [28] classified the TMF into Lower Mountain Forest (LMF) and Upper Mountain Forest (UMF). Nonetheless, the calculated mean AGB for the whole TMF (106.2 ha⁻¹; Table 3) is confirmed by Spracklen and Righelato [8], who estimated a similar mean AGB (104 Mg ha⁻¹) for a TMF in southern Ecuador by means of field plot measurements, which underlines the importance of an accurate field plot distribution to estimate AGB in natural forest stands.

Incidentally, the dominant/codominant trees within a forest stand represent 70% to 90% of the total AGB, as [41,105,106] indicated and therefore the inaccuracy of the present AGB estimation should be between 10% and 30%. This is verified by the comparison between the measured and calculated H, DBH and AGB values at tree level (Figure 6a,b and Figure 8), in which the coefficients of determination (R² > 0.91) indicated a high accuracy for each parameter. Therefore, the overall error of the applied AGB estimation approach should lie between 12.5% and 24.0%, as also estimated by other studies taking into consideration sampling errors [8,21].

In summary, the natural forest stands, especially the TMF, store the majority of C (over 95%; Table 4), as has also been reported [4], which highlights the importance of intact forest ecosystems as C stocks. The pastures as well as the subpáramo area are negligible C stocks taking into consideration future anthropogenic CO₂ emissions [2,11]. However, the subpáramo areas, as well as the natural forest ecosystems (TMF and Elfin forest), provide another important service, being the water supply for the local and regional population (e.g., [116]). Therefore, the primary objective should be the protection of these natural ecosystems from human disturbance [4,6,7], while considering the possibility of additional greenhouse gas emissions caused by ongoing deforestation, which is especially problematic in the tropical mountains of Ecuador [9,13]. Also, reforestation programs, such as REDD+ [117] should be expedited, however, occasionally the efforts are complicated due to the changing and more extreme climatic conditions over pasture and agricultural lands (e.g., [81,82]).

5. Conclusions

The objective of the study was a reliable quantification of AGB and its distribution in a natural TMF in southern Ecuador using LiDAR data in combination with field measurements. In general, the applied method, based on a classification at tree level, allowed a reliable AGB estimation for the whole forest stand. Using high resolution LiDAR data in combination with field measurements, avoided high under- or overestimation of the AGB [46]. Furthermore, by means of the vegetation classification at tree level and the topographic (TPI) and morphometric (MPI) indices, the AGB distribution within the whole study catchment could be detected and analysed. The results showed that high AGB (Mg ha⁻¹) also exists at higher elevations if the site is protected by the surrounding terrain (ravine forest) but these areas are generally inaccessible, which explains why remote sensed data is necessary to detect the complete forest structure.

This indicates the importance of undisturbed natural forest ecosystems as C stocks, because over 95% of the total AGB is stored in the mountain forest ecosystems, especially in the ravine TMF, where the tallest trees are located. The pasture sites, as well as the subpáramo area, do not contribute much to the total AGB, because C is basically stored in woody vegetation (trunks and branches; for example, Reference [23]) and the woody plants of these vegetation units are generally small. However, the TMF as well as the subpáramo provide another important ecosystem service: the water supply for the population, which is non-negotiable with regards to ongoing deforestation. Therefore, protection of these natural ecosystems is the main issue for the development of the local population and concurrently for facing future GHG emissions.