Determinants of Aboveground Biomass across an Afromontane Landscape Mosaic in Kenya

Afromontane tropical forests maintain high biodiversity and provide valuable ecosystem services, such as carbon sequestration. The spatial distribution of aboveground biomass (AGB) in forest-agriculture landscape mosaics is highly variable and controlled both by physical and human factors. In this study, the objectives were (1) to generate a map of AGB for the Taita Hills, in Kenya, based on field measurements and airborne laser scanning (ALS), and (2) to examine determinants of AGB using geospatial data and statistical modelling. The study area is located in the northernmost part of the Eastern Arc Mountains, with an elevation range of approximately 600–2200 m. The field measurements were carried out in 215 plots in 2013–2015 and ALS flights conducted in 2014–2015. Multiple linear regression was used for predicting AGB at a 30 m × 30 m resolution based on canopy cover and the 25th percentile height derived from ALS returns (R2 = 0.88, RMSE = 52.9 Mg ha−1). Boosted regression trees (BRT) were used for examining the relationship between AGB and explanatory variables at a 250 m× 250 m resolution. According to the results, AGB patterns were controlled mainly by mean annual precipitation (MAP), the distribution of croplands and slope, which explained together 69.8% of the AGB variation. The highest AGB densities have been retained in the semi-natural vegetation in the higher elevations receiving more rainfall and in the steep slope, which is less suitable for agriculture. AGB was also relatively high in the eastern slopes as indicated by the strong interaction between slope and aspect. Furthermore, plantation forests, topographic position and the density of buildings had a minor influence on AGB. The findings demonstrate the utility of ALS-based AGB maps and BRT for describing AGB distributions across Afromontane landscapes, which is important for making sustainable land management decisions in the region.


Introduction
Afromontane tropical forests are globally important ecosystems, which maintain high biodiversity and provide valuable ecosystem services, such as carbon sequestration [1].Tropical forests store large amounts of carbon, but due to forest degradation and deforestation, they can emit carbon to the atmosphere and boost global warming [2].It tropical countries, such as Kenya, forests are under constant pressure because of increasing demand for food and shelter [3].In East Africa, between 2002 and 2008, wooded vegetation cover decreased in area by 5.1%, 15.8% and 19.4% from forests, woodland and shrubland, respectively [4].Hence, mapping above ground biomass (AGB) is an essential task for monitoring carbon stocks and dynamics across tropical African landscapes.Measuring AGB is also required for implementing payments for ecosystem services schemes and the United Nations Reducing Emissions from Deforestation and Forest Degradation (REDD+) mechanism [5].Furthermore, AGB maps can be used for studying the determinants that control AGB distribution [6], which is important for the management of carbon stocks and understanding how carbon stocks might change in the future.
Field measurements provide accurate information about AGB, but sampling can be spatially limited due to difficult terrain or lacking resources.Remote sensing provides a cost-efficient alternative to mapping the wall-to-wall distribution of AGB across landscapes.Nonetheless, the use of optical satellite imagery for AGB mapping is often hindered by frequent cloud coverage.Furthermore, over dense crown covers, the data obtained by optical sensors may have limited capabilities in describing vegetation structure, as most spectral information is retrieved from the topmost parts of the canopy [2].Radar sensors also have critical limitations in assessing AGB over tropical regions, given that their sensitivity to AGB variation decreases in the old-growth forests [7].On the other hand, light detection and ranging (LiDAR) systems offer three-dimensional information on canopy structure that is better suited for AGB modelling than optical or radar imagery [8].In particular, airborne laser scanning (ALS) has been found to be useful for estimating forest biophysical parameters, such as leaf area index (LAI), canopy height, canopy cover and AGB (e.g., [8,9]).
Although ALS has been used increasingly for AGB estimation and mapping in Africa [10][11][12], most studies to date have focused on tropical forests (e.g., [11,13]).However, tropical African landscapes are typically mosaics of multiple land uses, and a large fraction of landscape-level AGB can be located outside the forests [14].For example, [15] determined mean AGB densities for land cover classes using samples of ALS data and a satellite image-based land cover map.Nevertheless, the distribution of AGB is highly variable, even under the same climatic condition and land use types [16,17].Under these circumstances, wall-to-wall ALS coverage provides unique possibilities to recover spatial variations in AGB patterns.
Furthermore, ALS-based AGB maps can be paired with other geospatial data to study the fundamental biotic and abiotic determinants of AGB [6,18,19].Digital elevation models (DEM), climate data, land use and land cover maps and soil databases are commonly used for this purpose [6,16,[18][19][20].Although the underlining mechanisms of the relationships between large-scale drivers (e.g., climate variables) and AGB are currently better understood [21], the landscape scale drivers of the spatial distribution of AGB are still poorly understood.For instance, while the relationships between topography, biophysical factors, climate and AGB have been evaluated in previous studies, inconsistent results were observed in different ecosystems.
The objectives of this study were to (1) generate an AGB map for the Taita Hills using ALS and (2) examine the determinants of AGB using geospatial datasets on topography, soils, climate, land use and statistical modelling at the landscape to the regional scale.The Taita Hills, in southeastern Kenya, provide a unique setting to address our objectives.In a relatively small area, with large topographic variations, the landscape is characterized by a complex land cover mosaic, with large local variations in AGB.

Study Area
The Taita Hills (study area ca.55,000 ha) in southeast Kenya are part of the Eastern Arc Mountains (Figure 1).The highest hilltop (Vuria) reaches 2208 meters above sea level (m a.s.l.).The hills are surrounded by semi-arid scrublands and dry savannas at an elevation of 600-900 m a.s.l.Due to the altitude difference, the hills experience a lower mean annual temperature (18.2 • C) compared with the lowlands (24.6 • C) [22].The region has a bimodal rainfall pattern: long rains from March-May/June and short rains from October-December.Mean annual precipitations between 1986 and 2003 were 1132 mm at 1768 m a.s.l (Mgange) and 587 mm at 560 m a.s.l.(Voi) [23].The higher elevations receive more rain because of the orographic rainfall pattern, while the northwestern slope in the rain-shadow receives less precipitation compared with the southeastern slope [3].
The area belongs to one of the world's most important biodiversity hotspots [1] and hosts some remaining indigenous montane forests [24].The landscape is a highly fragmented mosaic of montane forest patches, exotic plantations and agriculture with agroforestry and horticulture being common [25].
According to [15], the main land use and land cover types in 2011 were cropland (40.5%), thicket (28.9%), shrubland (12.8%), woodland and agroforestry (10.7%), plantation forest (4.6%), indigenous montane forest (1.4%) and bare rock, soil and built-up areas (1.1%).Agricultural expansion has been taking place in the area for decades.Since 1987, cropland area has increased at the expense of thickets and shrublands in the foothills and lowlands, while the hills were taken to agriculture before the 1950s and after [23,24].
The remnant montane forests can be separated into lower and upper montane forests.The lower montane forests are taller and drier compared with the elfin forest like upper montane forests, which receive significant amounts of precipitation through a mist.These mist forests on the top of Vuria, have typically only one layer and are covered by epiphytic mosses and lichens [26,27].
Common tree species in the montane forests include Tabernaemontana stapfiana, Macaranga capensis, Psychotria petitii, Xymalos monospora, Prunus africana, Maesa lanceolata and Albizia gummifera, while Eucalyptus spp., Cupressus lusitanica and Pinus spp.are the main tree species in the plantations.In the lowlands, common tree species include Albizia amara, Sterculia africana, Lannea alata, Acacia tortilis, Commiphora africana and Commiphora incisa.Rocky and sandy areas outside the forests are dominated by Acacia mearnsii, while typical agroforestry species in the hills are Grevillea robusta and Persea americana and Mangifera indica in the lowlands [28,29].elevations receive more rain because of the orographic rainfall pattern, while the northwestern slope in the rain-shadow receives less precipitation compared with the southeastern slope [3].
The area belongs to one of the world's most important biodiversity hotspots [1] and hosts some remaining indigenous montane forests [24].The landscape is a highly fragmented mosaic of montane forest patches, exotic plantations and agriculture with agroforestry and horticulture being common [25].According to [15], the main land use and land cover types in 2011 were cropland (40.5%), thicket (28.9%), shrubland (12.8%), woodland and agroforestry (10.7%), plantation forest (4.6%), indigenous montane forest (1.4%) and bare rock, soil and built-up areas (1.1%).Agricultural expansion has been taking place in the area for decades.Since 1987, cropland area has increased at the expense of thickets and shrublands in the foothills and lowlands, while the hills were taken to agriculture before the 1950s and after [23,24].
The remnant montane forests can be separated into lower and upper montane forests.The lower montane forests are taller and drier compared with the elfin forest like upper montane forests, which receive significant amounts of precipitation through a mist.These mist forests on the top of Vuria, have typically only one layer and are covered by epiphytic mosses and lichens [26,27].
Common tree species in the montane forests include Tabernaemontana stapfiana, Macaranga capensis, Psychotria petitii, Xymalos monospora, Prunus africana, Maesa lanceolata and Albizia gummifera, while Eucalyptus spp., Cupressus lusitanica and Pinus spp.are the main tree species in the plantations.In the lowlands, common tree species include Albizia amara, Sterculia africana, Lannea alata, Acacia tortilis, Commiphora africana and Commiphora incisa.Rocky and sandy areas outside the forests are dominated by Acacia mearnsii, while typical agroforestry species in the hills are Grevillea robusta and Persea americana and Mangifera indica in the lowlands [28,29].

Field Measurements and Aboveground Biomass Computations
The field measurements from a total of 215 circular 0.1-ha sample plots (radius 17.84 m) were used in this study for AGB modelling (Figure 1).The sampling design depended on the field campaign and sampling year.In 2013 and 2014, 150 plots were sampled randomly within 100-ha clusters in the hills and lowlands.In addition, 65 plots were sampled subjectively from high-biomass

Field Measurements and Aboveground Biomass Computations
The field measurements from a total of 215 circular 0.1-ha sample plots (radius 17.84 m) were used in this study for AGB modelling (Figure 1).The sampling design depended on the field campaign and sampling year.In 2013 and 2014, 150 plots were sampled randomly within 100-ha clusters in the hills and lowlands.In addition, 65 plots were sampled subjectively from high-biomass forest areas in 2014 and 2015.Sampling in the forest was guided by the ALS data [31] and visible to near-infrared imaging spectroscopy data (AisaEAGLE) [32] acquired in 2013 in order to cover variation in canopy structure and tree species composition.
Sample plot centers were positioned using a Trimble GeoXH GNSS receiver, and differential correction was made using a GNSS base station setup close to the town of Wundanyi (Figure 1).The slope in field plots varied from 0-50.5 • .Slope correction was done either in the field in order to keep the plot area constant (i.e., 0.1 ha), or sample plot areas were corrected based on mean slope computed from ALS-derived DEM.
All of the tree stems with a diameter (D) at breast height (1.3 m) ≥10 cm were measured for D [18,20], and tree species were identified by a local para-taxonomist.In 2013, tree height (H) was measured for the majority of the trees outside forests by a clinometer, but only random measurements were made in the forests.In 2014 and 2015, H was measured using a laser range finder for at least three trees in each plot (i.e., the tree with minimum, maximum and median D) [32].Furthermore, H was measured for all of the palms.For trees with only D measurement, H was predicted using a two-parameter Curtis height function [33].Non-linear mixed effect modelling and the plot as random effects (e.g., [34]) were used to calibrate the H-D model for each sample plot.The 'nlme' package [35] in the R statistical software [36] was used to fit the model (RMSE 1.47 m, 22.6%).
In total, 6337 trees belonging to 141 species were measured.Tree species-specific equations are poorly available for the species-rich study area (e.g., [37]).However, generic species-specific models were used for Acacia spp., Eucalyptus spp., Pinus spp.and palms [37][38][39]; for the rest of the species, the AGB (kg) was computed based on wood-specific density ( , g cm −3 ), D (cm) and H (m) by the pantropical allometric model of [40].Values of were obtained from the global wood density databases [41,42].Genus level averages were used for the species that lacked species-specific values, and if such could not be derived, a site-averaged value was used.Based on the tree-wise AGB estimates, AGB was calculated for the plots.Descriptive statistics for the field plots are given in Table 1.

Airborne Laser Scanning Data
ALS flights were conducted during several periods in 2014-2015 (Table 2).The data were provided as pre-processed and georeferenced point clouds.LAStools software (rapidlasso GmbH) was used to classify ground, building and vegetation returns and to create 2-m (for normalizing height) and 5-m (for explanatory variables) resolution DEMs.Triangular irregular network (TIN) interpolation was used for interpolating between the ground points [43].
The ALS point cloud elevations were normalized to heights above the ground surface by using the 2-m resolution DEM.Based on the normalized heights, noisy returns (e.g., high points) and returns from the electric lines were identified and removed manually.

Aboveground Biomass Mapping
The vegetation returns were further analyzed using FUSION software [44].For modelling, various ALS metrics representing return height distribution and canopy cover were computed for the plots (Table 3).A 3-m threshold was used to separate understory and ground returns from canopy returns.This threshold was considered to match well with the minimum D used in the field measurements (i.e., 10 cm).For prediction, spatial grids of ALS metrics were generated at a spatial resolution of 30 m × 30 m.The spatial resolution of ALS metrics was determined by the sample plot size (0.1 ha).The correlation coefficients among the ALS metrics are shown in Figure A1.
General models for modelling AGB from ALS metrics are available (e.g., [45,46]).To obtain the highest possible accuracy to map AGB in the study area, the linear regression models were fit using the 'lm' function in R, and the 'regsubsets' function in the package 'leaps' [9,47] was used for exhaustive search of the best models having 1-4 predictors.The coefficient of determination (R 2 ) and the root mean square error (RMSE) were used to evaluate the AGB models.The models were validated using leave-one-out cross-validation.The best model was selected based on the lowest Akaike's information criterion and variance inflation factor (VIF < 4) [48].The topographical and hydrological variables were derived from ALS-based DEM (5 m × 5 m cell).Higher spatial resolution (5 m) was used as the spatial resolution of ALS metrics was considered too coarse for topographic modelling.Slope, aspect, topographic position index (TPI), topographic wetness index (TWI) and river network were computed using System for Automated Geoscientific Analyses (SAGA) GIS (v.2.1.2).The default settings were used if not stated otherwise.Slope ( • ) and aspect ( • ) were created using the 'slope, aspect, curvature' module.TPI is the difference between the elevation of a cell and the mean elevation of the cells in a neighborhood or within a specified radius around it.Positive values represent ridges, negative values valleys, and near zero values are either flat or have a constant slope [49].Six radii between 50 m and 500 m were tested to cover local and landscape-scale variation in the 'topographic position index' module [50].TWI was calculated using the 'SAGA wetness index' module [51].River networks were extracted using the 'channel network' module, which requires elevation and catchment area as input.Catchment area was computed using the 'catchment area (flow tracing)' module [52], and the initiation type greater than 100,000 was used for channel network calculation.Soil type was derived from the Soil Atlas of Africa in the vector format [53].Soil types included calcaric fluvisols (FL), rhodic ferralsols (FR), rendzic leptosols (LP), chromic luvisols (LV) and cambic umbrisols (UM).
Climatic variables included mean annual temperature ( • C) and mean annual precipitation (mm), which were derived from [54] and [23], respectively.Temperatures were collected using small data loggers, which were placed at a height of 1.5 m at 40 sites covering different elevation zones and land cover types from dry, lowlands savannas to humid montane forest [54].A 20-m resolution grid of mean annual precipitation (MAP) data was produced by interpolating monthly rainfall records of 11 ground stations in the Taita Hills from 1987-2005 using the procedures from the ANUSPLINE [55] package [23].
Land use classes (extent of croplands and plantation forest) were extracted from the most recent land use and land cover classification for the area [15,56].A SPOT 4 HRVIR satellite image from October 2011 with a pixel size of 20 m × 20 m has been used for the classification with ten classes.The distribution of buildings was based on the ALS point clouds 'building' class.Misclassified buildings were removed manually with the help of high-resolution satellite imagery available as the base map in ArcGIS [57].Furthermore, the intensity image from ALS and the same imagery were used for digitizing the main roads.The correlation coefficients among the explanatory variables are shown in Figure A2 (Appendix A).
Finally, the AGB map and all of the explanatory variables were aggregated to a 250 m × 250 m cell size for analysis.This was done in order to reduce uncertainties in AGB predictions [6] and explanatory variables and to better capture the landscape-and regional-scale impacts of explanatory variables on AGB.Length (m) of road digitized from the high-resolution imagery -

Statistical Analysis
Boosted regression trees (BRT) [58] were used for examining the relationship between AGB and the explanatory variables.BRT improves the predictive performance by combining many simple models.BRT combines the strength of two algorithms: (1) regression trees and (2) boosting [58].BRT can fit complex nonlinear relationships, and there is no need for prior data transformation.BRT had been used in similar studies for analyzing the effect of topography on biomass (e.g., [59]).
In order to reduce spatial autocorrelation, 25% of the cells were systematically sampled before statistical analysis.Correlations between the explanatory variables were studied, and highly correlated variables were excluded based on the Pearson correlation coefficient (|r| > 0.7) [60].Among the correlated explanatory variables, the variable with the strongest univariate relationship with the response variable was retained [20].The mean annual temperature (MAT) and elevation were correlated with MAP and hence removed (Appendix A).Furthermore, TWI was correlated with the slope.TPIs were correlated with each other, and the TPI500 having the strongest univariate relationship with AGB was retained.
BRT analysis was done using the 'dismo' package [61] in R. The Gaussian error distribution was used.Furthermore, the minimum predictive error was achieved when using a learning rate of 0.008, a tree complexity (interaction depth) of 7, a bag fraction of 0.5 and the tolerance method 'fixed'.

Aboveground Biomass Map
The best AGB model was fitted with an R 2 of 0.88 and RMSE of 52.9 Mg•ha −1 based on the 25th percentile of height values (H.p25) and canopy cover (CC.all.first)(Table 5).The inclusion of additional ALS metrics did not significantly improve the model fit.The scatterplot of the predicted versus field-estimated AGB is presented in Figure 2. A single model performed well across the study area without any signs of systematic over-or under-estimation in any range of the values.The results were back-transformed (squared), and the square of residual standard error (3.78) was added to the predicted values [62].
with the response variable was retained [20].The mean annual temperature (MAT) and elevation were correlated with MAP and hence removed (Appendix A).Furthermore, TWI was correlated with the slope.TPIs were correlated with each other, and the TPI500 having the strongest univariate relationship with AGB was retained.
BRT analysis was done using the 'dismo' package [61] in R. The Gaussian error distribution was used.Furthermore, the minimum predictive error was achieved when using a learning rate of 0.008, a tree complexity (interaction depth) of 7, a bag fraction of 0.5 and the tolerance method 'fixed'.

Aboveground Biomass Map
The best AGB model was fitted with an R 2 of 0.88 and RMSE of 52.9 Mg⋅ha −1 based on the 25th percentile of height values (H.p25) and canopy cover (CC.all.first)(Table 5).The inclusion of additional ALS metrics did not significantly improve the model fit.The scatterplot of the predicted versus field-estimated AGB is presented in Figure 2. A single model performed well across the study area without any signs of systematic over-or under-estimation in any range of the values.The results were back-transformed (squared), and the square of residual standard error (3.78) was added to the predicted values [62].The AGB map predicted at 30 m × 30 m resolution is shown in Figure .3. The largest values were concentrated on montane forest patches in the hills and slopes with northeast and southeast aspects.In general, the AGB density decreases towards the lower elevations.Furthermore, lower densities (less than 31 Mg ha −1 ) are observed in the northwestern part of the area, which is in the rain-shadow of the hills.Some hilltops, for example in Yale and Ngangao, are without tree cover.In the foothills and lowlands, AGB densities are low and less variable compared with the hills, particularly in the northeast and southeast parts of the area.In the lowlands, larger AGB is concentrated in riverine areas, as those provide the consistent wetness necessary for trees.The AGB map predicted at 30 m × 30 m resolution is shown in Figure 3.The largest values were concentrated on montane forest patches in the hills and slopes with northeast and southeast aspects.In general, the AGB density decreases towards the lower elevations.Furthermore, lower densities (less than 31 Mg ha −1 ) are observed in the northwestern part of the area, which is in the rain-shadow of the hills.Some hilltops, for example in Yale and Ngangao, are without tree cover.In the foothills and lowlands, AGB densities are low and less variable compared with the hills, particularly in the northeast and southeast parts of the area.In the lowlands, larger AGB is concentrated in riverine areas, as those provide the consistent wetness necessary for trees.

Determinants of Aboveground Biomass
The AGB map and all of the explanatory variables used to explain the spatial distribution of AGB in the Afromontane landscape are shown in Figure A3.Plantation forests are limited in the northeast and southwest regions, which could be due to the distribution of precipitation in the region.LULC variables, for example buildings, are spatially concentrated, and croplands are really sparsely distributed (Figure A3).
The total explained deviance (cross-validated D 2 ) of the BRT model was 72.7%. Figure 4 shows the relative importance and partial dependence plot for each explanatory variable.The most influential variables were MAP (37.6%), cropland (16.9%) and slope (15.3%).Only 5.9% of the total deviance was explained by rivers, roads and soil type.
According to the partial dependence plots, AGB increases rapidly with increasing precipitation after approximately 1240 mm, which corresponds to elevations above 1350 m.Cropland and slope also play a role in AGB distribution, especially at the lower cropland coverage and the steeper slope (<37 degrees).
The effect of the plantation as predictor variables is observed especially at the lower plantation coverage (i.e., below 65%), and above that, it is constant.Plantation forests were available on small patches, and only a few areas were completely covered.The relative strength of the contribution of TPI500 is high especially at the beginning and end of the range.At the landscape scale (TPI500), high AGB was found both in the upper and lower slope.Aspect clearly shows that AGB is more concentrated in northeast and southeast slopes with decreasing AGB on the leeward side.Furthermore, areas with fewer buildings contain more AGB.
Figure 5 displays two-variable partial dependence plots on some of the most influential variables.Interaction effects of varying degree are presented among these variable pairs.The strongest interactions in the final model were observed between MAP and cropland (Figure 5a), and MAP and slope (Figure 5b).Cropland and slope have important effects on the predicted spatial distribution of AGB when MAP is at the highest and cropland and slope are at the lowest fraction (below 20) and below 37 degrees, respectively.When MAP is less than approximately 1250 mm, cropland and slope have low impacts on AGB.Additionally, the steep slope (up to 37 degrees) in the northeast aspect (Figure 5c) and regions with low cropland fraction (Figure 5d) contain high AGB.However, comparatively less AGB is present in the steepest slope.

Determinants of Aboveground Biomass
The AGB map and all of the explanatory variables used to explain the spatial distribution of AGB in the Afromontane landscape are shown in Figure A3.Plantation forests are limited in the northeast and southwest regions, which could be due to the distribution of precipitation in the region.LULC variables, for example buildings, are spatially concentrated, and croplands are really sparsely distributed (Figure A3).
The total explained deviance (cross-validated D 2 ) of the BRT model was 72.7%. Figure 4 shows the relative importance and partial dependence plot for each explanatory variable.The most influential variables were MAP (37.6%), cropland (16.9%) and slope (15.3%).Only 5.9% of the total deviance was explained by rivers, roads and soil type.
According to the partial dependence plots, AGB increases rapidly with increasing precipitation after approximately 1240 mm, which corresponds to elevations above 1350 m.Cropland and slope also play a role in AGB distribution, especially at the lower cropland coverage and the steeper slope (<37 degrees).
The effect of the plantation as predictor variables is observed especially at the lower plantation coverage (i.e., below 65%), and above that, it is constant.Plantation forests were available on small patches, and only a few areas were completely covered.The relative strength of the contribution of TPI500 is high especially at the beginning and end of the range.At the landscape scale (TPI500), high AGB was found both in the upper and lower slope.Aspect clearly shows that AGB is more concentrated in northeast and southeast slopes with decreasing AGB on the leeward side.Furthermore, areas with fewer buildings contain more AGB.
Figure 5 displays two-variable partial dependence plots on some of the most influential variables.Interaction effects of varying degree are presented among these variable pairs.The strongest interactions in the final model were observed between MAP and cropland (Figure 5a), and MAP and slope (Figure 5b).Cropland and slope have important effects on the predicted spatial distribution of AGB when MAP is at the highest and cropland and slope are at the lowest fraction (below 20) and below 37 degrees, respectively.When MAP is less than approximately 1250 mm, cropland and slope have low impacts on AGB.Additionally, the steep slope (up to 37 degrees) in the northeast aspect (Figure 5c) and regions with low cropland fraction (Figure 5d) contain high AGB.However, comparatively less AGB is present in the steepest slope.4 for a description of the variables.

ALS-Based Aboveground Biomass Map
Using a combination of field plots and ALS metrics, a wall-to-wall high-resolution AGB map was made for the Taita Hills.The AGB model had an R 2 of 0.88 and an RMSE of 52.9 Mg•ha −1 (43%).Some of the previous studies in the tropics have achieved better model fits and accuracies [10,12,16,63].However, a larger plot size has usually been used, which reduces the error [18].Furthermore, the heterogeneous land cover could explain the worse accuracy, but a single model for the whole area was considered appropriate here as systematic over-or under-estimation was not observed.Therefore, averaging AGB estimates to 250 m × 250 m resolution for modelling is likely to improve accuracy considerably [6].The final model was based on percentile height (H.p25) and canopy cover (CC.all.first).The variables are in line with the previous studies in Africa, which have used, for example, mean canopy height, height percentiles and canopy density for modelling [6,[10][11][12].Both modelling accuracy and explanatory variables are similar to the model used for the hills in the study by [15], who developed separate models for hills and lowlands using different ALS data.
Although the general patterns of AGB are very similar, the current map is an improvement over an earlier aboveground carbon map for the Taita Hills [15].The authors in [15] generated their map for analyzing the effect of the land cover change on carbon stocks using a satellite image-based land cover map and land cover type-specific mean carbon density values.However, the continuous AGB estimates are important as AGB depends on many variables, also within land cover classes (see Section 4.2).If converted to carbon, the new map had slightly larger mean carbon densities for montane forests (94.9 Mg•C•ha −1 ), exotic plantation forests (36.Uncertainties in the ALS-based mapping could be decreased by using species-specific allometric equations [37].Due to the absence of local wood density data and species-level allometric equations, beside genera Acacia and Eucalyptus, we used [40] for the rest of the species.Therefore, future studies should attempt to contribute to filling this gap.Furthermore, stems having D < 10 cm were not sampled in this study as they usually contribute only little to the aboveground carbon in mature African tropical forests [20,64].However, in a boreal forest, The author of [65] showed that a minimum DBH of 3 cm produced better AGB modelling results than 10 cm for the young forests, although improvement was small for the mature forests.Therefore, having a smaller DBH limit could improve modelling accuracy slightly in the areas with abundant small trees and shrubs.Additionally, the slope can bias ALS height percentiles [66], which were used in the final AGB model.This could increase AGB predictions in the steep slopes with an effect on the BRT analysis (Section 4.2), and future studies should consider applying a slope correction [67] before computing ALS height metrics.Nonetheless, there was no sign of systematic overestimation in the steep slopes when assessing the model residuals against the slope.

Determinants of Aboveground Biomass Distribution
AGB distribution in the Taita Hills is driven by the quantity and spatial distribution of precipitation.Precipitation is uneven due to the orographic rain pattern causing the eastern slopes and eastern parts of the hills to receive more precipitation than the leeward side in the western region, which is suffering from high water stress.AGB was the greatest at mid-altitude hills where MAP is high in comparison with the lowlands.The authors of [21] observed that precipitation explained 39% of the variation in LAI, and LAI increased significantly with increasing MAP in woody biomes.The bimodal rainfall pattern in the study area provides precipitation for more than half a year.In this study, precipitation explained 37.6% of AGB variation.The authors of [68] observed that AGB in African tropical forests shows a positive relationship with rainfall in the driest nine months and is negative for the wettest three months of the year.Low AGB in the lowlands could be due to higher air temperature, which results in higher respiration, and low precipitation in the lowlands could limit photosynthesis.
Both precipitation and temperature are closely related to elevation.We removed elevation and MAT from the analysis as they strongly correlated with MAP.The authors of [6] observed that elevation and the fraction of photosynthetic vegetation cover together accounted for 27-67% of the spatial variation in aboveground carbon density in areas in the north and south aspects of Madagascar.The authors of [20] revealed that AGB was the greatest in gentle slope and mid-elevation, which explained 63.7% of the variation in aboveground carbon in the Eastern Arc Mountains, especially Tanzania.The authors of [16] highlighted the fact that increasing elevation corresponds to a 53-84% decrease in AGB levels depending on substrate age class.Normally, at higher elevations, growth can be limited by water shortage and reduced temperature [20,69].However, in tropical montane cloud forests, high atmospheric-humidity levels are sustained due to frequent cloud cover and fog [69].
Furthermore, greater AGB was found in areas with smaller cropland fraction and on the relatively steep slope, compared with valleys and ridges.Croplands, if present in the mid-altitude, are concentrated on the northern and southwestern slopes, for example around Vuria and Ngangao.There is an increasing trend in croplands at the expense of thickets and shrublands in the lowlands.The amount of precipitation is low in the lowlands, which is directly linked with productivity.In the hills above 1220 m a.s.l., croplands are converted to shrublands and woodlands, due to the growing number of trees in the farms due to agroforestry and government policies [15].Besides that, cropland in hills has less economic returns from cash crop and high labor demand.This pattern is also observed here as even though the cropland fraction in the area is more than 50%, there is still some AGB (<120 Mg•ha −1 ).Steeper slopes have less cropland as they are not suitable for agriculture due to soil loss during heavy rainfall.The slope was the most important control in tropical landscapes in Costa Rica [19].Convexity, elevation and slope were significantly related to AGB in subtropical broad-leaved forests of Taiwan [70].The results are in line with [19] and [70], as slope also explained more of the variation in AGB in the Afromontane landscape mosaic.In this study, larger AGBs are limited to steep slopes (<37 degrees), compared with shallow and steepest slope in [20] and [71].The steepest slopes and hilltops are covered by bare rocks, and they do not have enough nutrients and soil to hold high biomass forests.For instance, part of the hilltops in Ngangao are covered by Erica mannii and the top of Vuria by short elfin forests.The forest on the hilltops is also reserved for conservation under the Kenya Forest Service [72] and recognized as sacred forests [73].
Plantation of trees on former agriculture lands substantially increases biomass accumulation during the first few years of forest recovery in the restored tropical forest in Costa Rica [74].In addition to the small and fragmented indigenous forest patches, the Taita Hills host plantations of exotic tree species, such as Cupressus lusitanica, Eucalyptus spp., Grevillea robusta and Pinus spp.[24], which explained 8.3% of the AGB variation.Forests found on a steep slope at higher altitudes were less likely to be deforested and less accessible to transport forest materials, while they were more likely to contain the tallest trees harboring the highest AGB.However, the height of trees is low in the upper montane forest, like in Vuria, which could be due to the frequent occurrence of low cloud, the possibility of mechanical damage due to the wind and low radiation (persistence cloud) [69].
Except for plantation and cropland, other land use and land cover variables were found to be relatively insensitive to the spatial distribution of AGB, for example road and rivers.Built-up areas are sparsely distributed and have less impact on AGB.The local people practice agroforestry in the hills in order to preserve the soil, its nutrients and provide shade from trees.There is a low population density in the drier northwestern part of the study area.The authors of [75] showed the dependency of dense population with precipitation in the Taita Hills.In rural Africa, dense populations are linked to croplands with low biomass compared with forests.The authors of [18] observed that for upland tropical rain forest landscapes, AGB was relatively insensitive to soil type and topography.The influence of soil was among the lowest also in this study.
Previous studies have analyzed the relationship between ALS-based AGB and biotic and abiotic explanatory variables from the plot to the landscape level [6,18,19].The relationship between the predictors and response variables differs between biomes and changes over spatial scales and through time [60].The fitted functions are not supposed to be applied to another area as the objective of BRT analysis was to help with understanding the controlling factors of AGB, but not making a predictive model for AGB.Although the regression methods were different, explanatory variables with higher spatial resolution (1 ha) explained 67% of the spatial variation in aboveground carbon density in [6].Additional research is needed to compare the explained deviance using the same regression techniques and explanatory variables at different spatial resolutions.Furthermore, the model could be further improved by including forest stand characteristics, for example stand age or disturbance history.

Conclusions
The Taita Hills comprise a human-dominated landscape with heterogeneous land cover and forest being limited to small patches typically far from the main areas of human settlement.In this study, field estimates of AGB were used for mapping AGB using ALS metrics to construct a wall-to-wall map for the Taita Hills area.ALS metrics representing the height of the canopy and canopy cover (i.e., volume of the canopy) predict AGB with the highest accuracy.AGB is concentrated in native moist evergreen montane forests, which are located in the windward side in the north and southeast aspects of the Taita Hills in higher elevations receiving more rainfall.The high biomass areas are also located on hilltops and steep slopes, as they are too cumbersome to take into agricultural practice.This study revealed the influence of climate, land use and topography in shaping the spatial patterns in AGB in a heterogeneous landscape.Using the coverage of the AGB map and various explanatory variables, it was possible to capture the spatial variation in AGB throughout the highly complex region and a large altitudinal gradient.Hence, the results provide novel insights into the cross-scale nature of determinants of AGB, in a region where multiple factors may interact non-linearly in defining landscape-level AGB.    4 for description of the variables.

Figure 1 .
Figure 1.The location of the study area and sample plots with false color Sentinel-2A Multispectral Instrument (MSI) satellite image from 8 October 2016 [30].

Figure 1 .
Figure 1.The location of the study area and sample plots with false color Sentinel-2A Multispectral Instrument (MSI) satellite image from 8 October 2016 [30].

Figure 2 .
Figure 2. ALS predicted vs. field estimated AGB based on leave-one-out cross-validation.

Figure. 3 .
Figure. 3. The ALS-based AGB map at 30 m × 30 m resolution with hillshade based on DEM in the background.

Figure 3 .
Figure 3.The ALS-based AGB map at 30 m × 30 m resolution with hillshade based on DEM in the background.

Figure 4 .
Figure 4. Single-variable partial dependence plots and smoothed response curves for the explanatory variables where the Y-axis is the fitted function for the response variable (AGB) (scale different).The relative influence of each variable is shown in brackets.The hash marks at the base of each plot represent the deciles of the explanatory variable distribution.Soil types: calcaric fluvisols (FLca), rhodic ferralsols (FRro), rendzic leptosols (LPrz), chromic luvisols (LVcr) and cambic umbrisols (UMcm).See Table4for a description of the variables.

Figure 5 .
Figure 5.The strongest interactions between the explanatory variables in the fitted boosted regression trees (BRT) model: (a) cropland vs. mean annual precipitation (MAP), (b) slope vs. MAP, (c) aspect vs. slope and (d) cropland vs. slope.Note that the plots have different scaling for the Y-axis.

Figure 4 .
Figure 4. Single-variable partial dependence plots and smoothed response curves for the explanatory variables where the Y-axis is the fitted function for the response variable (AGB) (scale different).The relative influence of each variable is shown in brackets.The hash marks at the base of each plot represent the deciles of the explanatory variable distribution.Soil types: calcaric fluvisols (FLca), rhodic ferralsols (FRro), rendzic leptosols (LPrz), chromic luvisols (LVcr) and cambic umbrisols (UMcm).See Table4for a description of the variables.

Figure 4 .
Figure 4. Single-variable partial dependence plots and smoothed response curves for the explanatory variables where the Y-axis is the fitted function for the response variable (AGB) (scale different).The relative influence of each variable is shown in brackets.The hash marks at the base of each plot represent the deciles of the explanatory variable distribution.Soil types: calcaric fluvisols (FLca), rhodic ferralsols (FRro), rendzic leptosols (LPrz), chromic luvisols (LVcr) and cambic umbrisols (UMcm).See Table4for a description of the variables.

Figure 5 .
Figure 5.The strongest interactions between the explanatory variables in the fitted boosted regression trees (BRT) model: (a) cropland vs. mean annual precipitation (MAP), (b) slope vs. MAP, (c) aspect vs. slope and (d) cropland vs. slope.Note that the plots have different scaling for the Y-axis.

Figure 5 .
Figure 5.The strongest interactions between the explanatory variables in the fitted boosted regression trees (BRT) model: (a) cropland vs. mean annual precipitation (MAP), (b) slope vs. MAP, (c) aspect vs. slope and (d) cropland vs. slope.Note that the plots have different scaling for the Y-axis.

Figure A2 .
Figure A2.Pearson correlation between the variables and their scatter plots.

Figure A3 .
Figure A3.ALS-based aboveground biomass (AGB) map (a) and (b-k) explanatory variables at 250m resolution.MAP is the mean annual precipitation, and TPI 500 is the topographic position index with a 500-m radius.Panels (i) and (j) refer to the distance of the river and roads in the given pixel (i.e., 250 m × 250 m).See Table4for description of the variables.

Figure A3 .
Figure A3.ALS-based aboveground biomass (AGB) map (a) and (b-k) explanatory variables at 250-m resolution.MAP is the mean annual precipitation, and TPI 500 is the topographic position index with a 500-m radius.Panels (i) and (j) refer to the distance of the river and roads in the given pixel (i.e., 250 m × 250 m).See Table4for description of the variables.

Table 2 .
Scanner properties and flight parameters of the ALS data.

Table 3 .
Summary of airborne laser scanning (ALS) metrics.Explanatory Variables for Modelling Distribution of Aboveground BiomassA range of explanatory variables was derived from geospatial datasets for modelling AGB.Table4presents the complete list of variables.

Table 4 .
Explanatory variables for aboveground biomass modelling.

Table 5 .
Summary of the final model for AGB mapping.H.p25 = 25th percentile point of all laser return heights above 3 m; CC.all.first= all returns above 3 m/total first returns × 100; *** p < 0.001.