1 Introduction

Population pressure, agricultural expansion/intensification and development of infrastructure have been suggested as major threats to biodiversity (Bargali et al., 2019; Davidar et al., 2010). Forest vegetation provides a wide range of ecosystem goods and services to the inhabitants (Gosain et al., 2015; Fartyal et al. 2022). Overexploitation of natural resources has created a big gap between the demand and supply of resources. To cope with this situation, plantations of fast-growing species, such as bamboos, are being established in the global scale to meet an increasing demand for goods and services, believing that a large amount of biomass/carbon can be harvested within short period of time (Bargali & Singh, 1991). The afforestation scheme includes monocultures of fast-growing species which are introduced either after clear-felling of the natural forests or within the bunds and boundaries/ mixed with others crops occurring there. It is generally assumed that soil/site-specific characteristics may change as a result of such species introduction (Bargali et al., 1993; Joshi et al., 1997; Kadeba & Aduayi, 1985; Stewart & Kellman, 1982).

Along the Colombian Andes from 900 to 2000 masl, the bamboo species Guadua angustifolia (guadua) is dominant within the so-called bamboo forests (Camargo, 2006). In the coffee region area, despite being highly fragmented (Camargo & Cardona, 2005), remnants of these forests fulfill significant ecological functions (Muñoz-López et al., 2020) and provide raw material for different applications (Arango Arango and Muñoz-López 2021; Correal & Arbelaez, 2010; García & Camargo, 2010; Takeuchi et al., 2009). In fact, guadua represents an important potential within the marketable timber species in Colombia (Cardona, 2012), where enterprises have the possibility to consolidate sustainable production systems (Arango et al., 2017).

Guadua is woody bamboo that has erect habit where the aerial vegetative axes known as culms interconnected by rhizomes in a defined pattern as long necked pachymorph (Judziewicz et al., 1999). This condition implies a density of culms per over 6,000 on average (Kleinn & Morales-Hidalgo, 2006). However, other plant species are also present and shape forest ecosystems with important values of diversity where 172 species belonging to 54 families have been identified (Ramírez-Díaz & Camargo, 2019).

An inventory of bamboo forests performed in the coffee region of Colombia estimated a total of 28,000 ha Kleinn and Morales-Hidalgo (2006). However, during the last decade, these forests have been highly deforested, especially around urban areas, at annual rates of approximately 5% (Muñoz, 2021), which considerably exceeds the national average (0.3%) registered by Ideam (2021).

According to Camargo and Long (2020), bamboo forests provide important ecosystem services (ESs) and therefore are strategic to maintain and guarantee a proper main ecological structure. The ESs provided are related to habitat for plant biodiversity (Ramírez-Díaz & Camargo, 2019), water purification (Chará et al., 2010), landscape restoration (Camargo et al., 2018), bioenergy (Daza Montaño et al., 2013) and carbon storage (Camargo et al., 2010).

Biomass estimation is essential for determining the status and flux of biological materials in an ecosystem and for understanding the dynamics of the ecosystem an (Anderson, 1970; Bargali et al., 1992). The quantity of tree biomass per unit area of land constitutes the primary inventory data needed to understand the flow of materials and water through forest ecosystems (Awasthi et al., 2022; Swank & Schreuder, 1974).

Understanding bamboo ecology is important to know its capacity to store carbon within biomass and is relevant for proper management and conservation (De Campos et al., 2021). Although the proportion of bamboo in the global forest area is small, its contribution to climate change mitigation is significant, especially for characteristics, such as the belowground biomass (Singnar et al., 2021) or the high density of culms per ha. For example, bamboo forests in the coffee region of Colombia include 6284 ± 3128 culms ha−1 on average (Camargo, 2006).

Bamboo may fulfill an important role in climate change mitigation (Yiping et al. 2010). The registered value of total ecosystem carbon (TEC) for bamboo is between 94 and 392 Mg C ha−1, and the aboveground carbon levels from 16 to 128 Mg C ha−1 are lower than that noted for some forest ecosystems (Yuen et al., 2017). However, bamboo is important as a carbon sink, even after only selected culms, where bamboo stores carbon, are harvested (Li et al., 2015).

Values of carbon stock for bamboo forests previously estimated in the coffee region of Colombia come from inventories where culm volume and wood density have been measured. For example, 311 Mg C ha1 was reported by (Kleinn & Morales-Hidalgo, 2006). Additionally, other data are derived from a destructive sampling approach (Arango, 2011), yielding a value of 126 Mg C ha−1. A range of values between 18 Mg C ha−1 and 260 Mg C ha−1 was obtained by Aguirre and Criollo (2020) and calculated from culm volume, wood density and a biomass expansion factor. The estimated carbon values show a high variability associated with features of bamboo forests as well as approaches used for calculations and estimates. Therefore, the challenge is to determine better methods to obtain proper information and estimates.

In addition, studies performed on guadua plantations by Riaño et al. (2002) and Camargo et al., 2010 in Colombia, Aguirre-Cadena et al. (2018) in Mexico and Fonseca and Rojas (2016) in Costa Rica highlighted the importance of these ecosystems for carbon storage and sequestration. The values registered by these studies vary according to the age of plantations and site conditions.

Allometric equations fitted to estimate the biomass of guadua are scarcely reported in the literature. In Bolivia, Quiroga Rojas et al. (2013) fitted an equation for natural bamboo forests located in a national park area at 1000 masl and high annual precipitation (5500 mm). The sample size of this study was 24 plants. Additionally, two specific studies with guadua plantations describe the process to fit allometric equations. One study in Colombia (Riaño et al., 2002) and another study in Costa Rica (Fonseca & Rojas, 2016) included six and twenty-one years of data, respectively. Therefore, equations with the possibility for application in a wider range of ecological conditions are required.

The characteristics of guadua bamboo forests along the coffee region might represent a high variability on the biomass and therefore in their carbon content. Thus, to develop of tools for estimating biomass of guadua as the dominant species within these forests is usefulness for determining their contribution to climate change mitigation.

The aim of this study was to provide a better approximation of the total carbon stored in bamboo forest ecosystems, considering the natural distribution along the elevation gradient in the coffee region of Colombia and including different carbon pools. Furthermore, this study seeks to provide better tools and different approaches for biomass estimation through allometric equations fitted from predictor variables that are easy to measure during inventories.

2 Materials and methods

2.1 Study area

The sites of this study are located along the Colombian Andes on the Central Mountain range in the coffee region zone (Fig. 1). Considering the natural distribution of guadua in this area, six sites from 900 m up to 2000 masl with natural bamboo forests previously included in inventory and monitoring processes were selected for sampling biomass. Two of the sites were located at the lowest elevation (900 m) and associated with flat topography corresponding to alluvial valleys. The remaining four sites were located on hillsides with steeper topography ranging from 1150 m up to 2000 masl.

Fig. 1
figure 1

Zone of study: coffee region of Colombia

Full size image

The climatic conditions of the study area correspond to an average temperature (T) from 17 °C at 2000 masl up to 21 °C at 900 masl. The mean annual precipitation (P) ranges from 2100 at 2500 mm. Predominant soils are Inceptisols at the lower elevation. Along the hillside, most of these soils correspond to Andisols (UTP-GATA, 2020).

3 Data collection and measurements

Within each site, plots of 100 m2 (10 m × 10 m) were selected. These plots were established in the framework of previous projects for inventory and monitoring of bamboo growth and variables to describe the quality of culms. For this study, a total of 23 plots were selected. The number of plots per site varied according to owner approval to permit destructive sampling.

Culms within each plot were accounted for, and data were used to estimate the total density of culms per ha (N). Additionally, the diameter at breast height (dbh) of culms was recorded. Those culms where age was previously registered during monitoring activities were considered for destructive sampling. Thus, a total of 130 of these culms were used for destructive sampling and calculating biomass. The rhizome of 13 culms was also extracted and included for biomass calculations. The age of culms included in destructive sampling varied from 12 to 120 months. Younger and older culms occasionally appeared without leaves and branches or were even broken but were included because they represented the natural conditions of natural bamboo forests.

Each culm selected for destructive sampling was cut and separated into four components, including leaves, branches, culm and rhizome (the last only for 13 culms). Then, a frame of 0.5 m × 0.5 m was randomly thrown thrice within each plot, and litter inside was collected. The total fresh weight (Fw) of each plant component plus litter was obtained in the field using a scale. Then, three samples of each component of approximately 300 g were labeled and sent to the laboratory to determine dry biomass (B). Samples were obtained considering the variability of the components; specifically, those from culms were obtained from the lower, middle and higher parts.

The diameter of culms cut at the first internode height (dint) and breast height (dbh) and the total length (l) were measured. Other variables, such as wall thickness at first internode height (tint) and length of the internode at first internode (lint) and at breast height (libh), were also recorded. Considering the ecological conditions of these bamboo forests, trees found within plots with a diameter at breast height greater than 5 cm (Treedbh) were recorded, and then values were extrapolated to ha (Ntree).

In addition, within each plot, three random soil samples were collected at a depth of 30 cm to determine organic carbon levels. Additionally, at the same points with cylinders of known volume, soil samples were taken at each 5 cm depth up to 30 cm for measuring bulk density (BD). Due to predominant characteristics of soil (depth and high organic matter) along the study area and the volume of cylinders (98 cm3), the portion of sample taken lack of stones or hard material that might influence on increasing soil mass.

4 Laboratory work and estimation of biomass

Samples of all components were brought to the Ecotechnology Laboratory at the Technological University of Pereira and then oven-dried at 105 °C until a constant weight was obtained (Li et al. 2016; Huy et al. 2019). Then, the B of the samples was calculated as follows:

$$B = \frac{{D_{{{\text{ws}}}} }}{{F_{{{\text{ws}}}} }} \times F_{w}$$
(1)

where B is the dry biomass (kg) of the component, Fws and Dw are the fresh and dry weight of the sample, and Fw is the total fresh weight of the bamboo component. Thus, the B values of leaves (Ble), branches (Bbr) and culm (Bcu) were calculated. The aboveground biomass (AGB) in kg was the sum of Ble, Bbr and Bcu:

$${\text{AGB}} = B_{cu} + B_{br} + B_{le}$$
(2)

The biomass of rhizome (Brh) was also calculated in the same manner for the corresponding samples. Using the ratio to AGB, a value of the biomass expansion factor (BEF) was obtained for estimating the total biomass (Bt) in kg:

$$B_{t} = {\text{AGB}} \times {\text{BEF}}$$
(3)

Then, considering N, the average Bt was calculated from the values obtained in a destructive sample, and the values were extrapolated to ha as follows:

$$B_{ba} = Bt \times N$$
(4)

where Bba is the total bamboo biomass (Mg ha−1) and N is the total number of culms per ha. Additionally, litter biomass (Bli) was calculated, and values were extrapolated to ha considering the frame area (0.25 m2).

Then, the proportion of soil organic carbon (SOC) was calculated for samples collected in the laboratory using the wet combustion method of Walkley–Black (Sato et al., 2014). Bulk density (BD) was also calculated after oven drying samples collected in cylinders until a constant weight of soil mass was obtained (Jaramillo, 2002). Thereafter, BD was used to calculate the total soil mass at 30 cm depth, and the value in Mg ha−1 was obtained based on the proportion of SOC.

The value of treedbh was averaged, and tree biomass (Btree) was estimated with the following equation proposed by Segura and Andrade (2008):

$$B_{{{\text{tree}}}} = 10^{{( - 1.11 + 2.64 \times \log \left( {tree_{{dbh}} } \right))}}$$
(5)

For Bba, Bli and Btree, the carbon fraction of 0.47 was applied to obtain the values of carbon pools (Huy and Long 2019). In addition, SOC was added to obtain an approximation of the total ecosystem carbon stock value (Mg ha−1) of bamboo forests (TECaprox):

$${\text{TEC}}_{{{\text{aprox}}}} = \left( {B_{ba } + B_{li} + B_{{{\text{tree}}}} } \right) \times 0.47 + {\text{SOC}}$$
(6)

Descriptive statistics were obtained for all variables calculated and estimated and associated with biomass and carbon pools as well those measured during the sampling. In this sense, correlation analyses between biomass components and variables as N and T were performed. Furthermore, an analysis of variance with elevation (E) as independent factor was also carried out. Analyses were performed using the software Rstudio (R Core Team, 2020).

5 Fitting model equations

The first step was the assessment of the relationships between the dependent variable AGB and the probable predictor variables. Relationships with other components of biomass (Ble, Bbr and Bcu) were also explored. This process included correlation analyses and graphical exploration to identify association patterns between the variables.

The criteria to define the variables included in the models were consistent with the relationships found. Specifically, the criteria included standard variables used as predictors in the literature (e.g., Huy & Long, 2019; Huy et al., 2019b; Picard et al., 2012, 2015) and variables that were cheap and easy to measure.

The form of the assessed allometric equations was consistent with those usually proposed in the literature (e.g., Huy & Long, 2019; Poudel & Temesgen, 2016; Picard et al., 2015; Picard et al., 2012). Therefore, the process started with fitting of the simple linear models and then the power model. The form of the assessed models is as follows:

$${\text{AGB}} = a_{1} + b_{1} \times X_{1} + \varepsilon_{1}$$
(7)
$${\text{log}}\left( {{\text{AGB}}} \right) = a_{2} + b_{2} \times {\text{log}}\left( { X_{2} } \right) + \varepsilon_{2}$$
(8)
$${\text{AGB}} = a_{3} \times X_{3}^{b3} + \varepsilon_{3}$$
(9)

where AGB is the total aboveground biomass in kg, ai and bi are the parameters of the model (i = 1, 2, 3), and Xi is the predictor variable (dbh, l or the combined dbh2 x l), which has also been used in other studies with allometric equations for bamboo (Huy et al., 2019a; Huy & Long, 2019).

Considering heteroscedasticity because the variation in the relationship between AGB and the predictor variable increases with greater values and variance is not constant (Picard et al., 2012), weighted regression models were also fitted. Therefore, weighted variables (1/X2δ), where X = dbh, dbh2 x l and δ the variance function coefficient (Huy et al., 2019a), were included for fitting the models.

Considering that AGB is calculated from the sum of biomass components (Ble, Bbr and Bcu), models independently fitted might be less effective for estimations (Huy et al., 2019a; Sanquetta et al., 2015). Therefore, the seemingly unrelated regression method (SUR) was also performed to simultaneously fit AGB, where the biomass of each component (Ble, Bbr and Bcu) is integrated to the additive systems of biomass equations (Cui et al., 2020). SUR approach may be an alternative to obtain better precision on estimates.

Considering the form of independent models with the better performance and the use of predictor variables easier of measuring during inventories, in SUR approach different combinations of simultaneous systems of equations were assessed and three of them were fitted. These systems had the following general form:

$$B_{cu} = a_{1} \times X_{1}^{b1} + \varepsilon_{1}$$
(10)
$$B_{br} = a_{2} \times X_{2}^{b2} + \varepsilon_{2}$$
(11)
$$B_{le} = a_{3} \times X_{3}^{b3} + \varepsilon_{3}$$
(12)
$${\text{AGB}} = B_{cu } + B_{br} + B_{le } = a_{1} \times X_{1}^{b1} + a_{2} \times X_{2}^{b2} + a_{3} \times X_{3}^{b3} + \varepsilon_{4}$$
(13)
$$B_{cu} = x X_{1}^{b1} + \varepsilon_{1}$$
(14)
$$B_{br} = x X_{2}^{b2} + \varepsilon_{2}$$
(15)
$$B_{le} = a_{3} \times X_{3}^{b3} + \varepsilon_{3}$$
(16)
$$B_{cu} = a_{1} + b_{1} \times X_{1} + \varepsilon_{1}$$
(18)
$${\text{AGB}} = a_{4} \times X_{4}^{b4} + \varepsilon_{4}$$
(17)
$$B_{br} = a_{2} + b_{2} \times X_{2} + \varepsilon_{2}$$
(19)
$$B_{le} = a_{3} + b_{3} \times X_{3} + \varepsilon_{3}$$
(20)
$${\text{AGB}} = a_{4} + b_{4} \times X_{4} + \varepsilon_{4}$$
(21)

where ai and bi are the parameters of the model (i = 1, 2, 3) and Xi is the predictor variable (dbh, dbh2 x l). Equations from 10 to 13 represent the equations system 1, from 14 to 17 the equations system 2 and from 18 to 21, the equations system 3.

6 Model selection

To assess the quality of the equations and select those with higher precision to estimate AGB, statistical tests and graphical analyses were performed. The statistics suggested by (Huy and Long, 2019) and used in this work are as follows:

$${\text{Akaike information criterion}}:{\text{ AIC}} = - 2\ln \left( L \right) + 2p$$
(22)

where L is the likelihood of the model, and p is the number of the parameters of the model.

$${\text{Coefficient of determination}}:{ }R^{2} = 1 - \frac{{\mathop \sum \nolimits_{i = 1}^{N} \left( {Y_{i} - \mathop {Y_{i} }\limits^{ \wedge } } \right)^{2} }}{{\mathop \sum \nolimits_{i = 1}^{N} \left( {Y_{i} - \mathop {Y_{i} }\limits^{\_\_} } \right)^{2} }}$$
(23)

where \(Y_{i}\) and \(\hat{Y}_{i}\) \(\overline{Y}_{i}\) are the observed, predicted and mean values of AGB, respectively.

$${\text{Adjusted coefficient of determination}}:R_{adj}^{2} = 1 - \left( {1 - R^{2} } \right).\frac{{\left( {N - 1} \right)}}{{\left( {N - p} \right)}}$$
(24)

where N is the total number of observations used in fitting the model, p is the number of parameters to estimate, and R2 is the coefficient of determination. For graphical analyses, diagnostic plots of observed versus fitted values and the trend of residuals were used to assess model performance (Huy and Long, 2019). In addition, the error of the models was assessed with the following statistics also suggested by (Huy and Long, 2019):

$${\text{Bias }}\% = \frac{100}{n}\mathop \sum \limits_{i = 1}^{n} \frac{{\left( {y_{i} - \hat{y}_{i} } \right)}}{{(y_{i} )}}$$
(25)
$${\text{RMSE}} = \sqrt {\frac{1}{n}\mathop \sum \limits_{i = 1}^{n} (y_{i} { } - \hat{y}_{i} )^{2} }$$
(26)
$${\text{MAPE}} \% = \frac{100}{n}\mathop \sum \limits_{i = 1}^{n} \frac{{\left| {y_{{\text{i}}} - { }\widehat{{y_{i} }}} \right|}}{{y_{i} }}$$
(27)

where n is the number of samples and yi and \({\widehat{y}}_{i}\) \(\mathrm{and} \overline{y }\) are the observed, predicted and averaged values for the ith AGB sample, respectively.

To compare those models with the transformed response variable (\(log(\mathrm{AGB})={a}_{2}+ {b}_{2} xlog( {X}_{2})+{\varepsilon }_{2}\)), the Furnivall index (Fln) was calculated (da Cunha & Guimarães-Finger, 2012; Picard et al., 2012; Segura & Andrade, 2008) as follows:

$$F_{ln} = {\text{exp}}\left( {\frac{1}{n}\mathop \sum \limits_{i = 1}^{n} \left( {{\text{ln}}({\text{AGB}}_{i} } \right)} \right) \times {\text{RMSE}}$$
(28)

In addition, the root-mean-square percentage error (RMSPE) was used to measure the predictive performance of models (Kralicek et al., 2017; Zell et al., 2014) and as an alternative of comparison with results of other studies on allometric equations for bamboo biomass:

$${\text{RMSEPE}} \% = \sqrt {\frac{1}{n}\mathop \sum \limits_{i = 1}^{n} \left( {\frac{{y_{i} - \hat{y}_{i} }}{{y_{i} }}} \right)^{2 } } \times 100$$
(29)

Fitting analyses were performed with the software Rstudio and R version 4.2.1 (R Core Team, 2020) and the packages lm and nlm for the models independently fitted. For the SUR approach, the simultaneous systems of equations were fitted with the package systemfit (Henningsen & Hamann, 2008).

7 Model validation

For those models with better performance after the fitting phase and considering the possibility of being used according to the predictor variables included, a cross-validation process was conducted. The weighted nonlinear fixed effect models and the three systems of simultaneous equations were performed by the Monte Carlo method (Huy & Long 2019), whereas the K-Folds repeated (James et al., 2013) method was applied to the weighted linear and weighted log-transformed models and permitted us to evaluate their uncertainty and reliability. Therefore, the models were fitted with 80% of the dataset, and cross-validation was performed with the remaining 20%. Then, the model with a smaller value of cross-validation error bias (%), RMSE, MAPE (%) and RMSEPE (%), is considered better.

8 Results and discussion

8.1 Values of the variables measured during the sampling

Tables 1, 2 and 3 summarize the characteristics assessed during the sampling, including site conditions and dendrometric and stand variables of bamboo forests as well as those related to biomass and carbon. The sites sampled are included within the conditions of the natural distribution of guadua bamboo forests, which correspond to humid, warm and predominantly lowlands of Andean ecosystems (Table 1).

Table 1 Climatic variables and elevation associated with sampling sites
Full size table
Table 2 Dendrometric and stand variables of the assessed natural guadua bamboo forests
Full size table
Table 3 Biomass and carbon content in the different pools associated with guadua bamboo forests
Full size table

Dendrometric and stand variables showed a wide range of values. The predominant culm age was four years, and the older was up to 10 years. The presence of trees species within plots, as mentioned above, certainly demonstrates the heterogeneity associated with natural bamboo forests where there are presence of other plant species, whit different growth patterns and habits shaping a high diversity forest ecosystem. Conditions that different among plantations are usually more homogeneous in terms of composition and structure (Table 2).

Culm represents 67% of Bba and is the more important component for the allocation of biomass. Some values of 0 for leaves and branches correspond to the youngest plants included, which were completely elongated and recently had lost culm sheath but without leaf or branch development. A few of the samples were old culms that were dry and broken in the upper part (Table 3).

The AGB values range from 9.6 kg to 92.7 kg, which confirms the high variability associated with and described above for dendrometric variables. However, differences when values are expressed per ha are not as notable. In addition, as other pools are integrated to obtain the value of TECaprox, the range of differences decreases, demonstrating less variability among stands than within them.

The analysis of variance showed a significant (p < 0.05) effect of elevation on AGB where lower values correspond to elevations around 2000 m and the higher were found along the range between 1200 and 1450 m.

The SOC accounts for an important proportion of the TECaprox, which is approximately 45%. Considering the contribution of Ctree and Cli, this value would reach approximately 50%. The distribution of C among the different pools shows the particularities of natural guadua bamboo forests in terms of a more complex structure and composition, where bamboo is an important component but not the only component (Table 3).

The dendrometric variables assessed showed values within the range of other studies performed for natural guadua bamboo stands (e.g., Camargo, 2006; Kleinn & Morales-Hidalgo, 2006). Most of these forests are eventually harvested, especially for domestic purposes (Camargo et al., 2007). Therefore, the conditions of structures with culms at different ages are maintained, including those over four years old, to reach the proper conditions for commercial uses (Maya et al., 2017; Rodríguez et al., 2010).

The average AGB and Cba carbon values found (131.4 ± 47 Mg Cbaha−1) correspond to the highest values registered in the literature for 70 species (Yuen et al., 2017). The SOC estimate of up to 30 cm depth with values over 100 Mg C ha−1 might also be considered high according to the values usually estimated for soil depths less than 50 cm (80–160 Mg C ha−1) (Yuen et al., 2017).

AGB and Cba values are comparable with values previously registered for guadua (Aguirre & Criollo, 2020; Arango, 2011; Quiroga Rojas et al., 2013); however, a remarkably wide range between 55.2 171.2 Mg C ha−1 and 171.2 Mg C ha−1 is noted and is related to the range of N between 4186 culms ha−1 and 7700 culms ha−1.

In this study, a BEF was estimated that represents the ratio between AGB and Crh. If this value is used as an approximation to the ratio between belowground biomass (BGB) and AGB, the average value obtained (0.71) shows an important contribution of Crh for allocation biomass. This value is comparable with values for other species, such those reported by Singnar et al. (2021). This value is greater than that reported by Yuen et al. (2017) for this bamboo species, and this difference is likely due to the fact that the samples were obtained from young plantations, such as those assessed by Riaño et al. (2002).

The TECaprox for the guadua bamboo forests assessed shows an interesting value of 259.7 Mg C ha1 ± 43.3, which is comparable to values estimated for tropical forests (Saatchi et al., 2011) and greater than that noted for South America (Sullivan et al., 2017). Regarding bamboo forests, ecosystems might also be considered within the higher values registered by Yuen et al. (2017). This condition represents the important relevance of guadua bamboo forests in the face of climate change.

Compared with Moso bamboo species, which represent the largest covered area and are important for commercial use, these species account for approximately 14 kg C culm and 54.6 Mg Cba ha−1 (Zhuang et al., 2015) with total carbon ranges from 87.83 Mg C ha−1 to 119.5 Mg C ha−1 (Xu et al. 2018a). In addition, the higher registered values here enhance the importance of guadua bamboo forests.

9 Model fitting and validation

Correlation analyses between dendrometric variables and AGB showed the strongest and most significant relationship (p < 0.05) with dbh, dint, l and dbh2l (p < 0.05). These variables were also plotted against AGB, Bcu, Br and Ble to explore the type of relation between variables (Fig. 2). AGB and Bcu tended to increase as greater dbh, dint, l and dbh2l. However, the relationship of dendrometric variables with Bbr and Ble did not show a clear pattern. Therefore, for the independent models only AGB was the response variable, which is in fact the variable of interest.

Fig. 2
figure 2figure 2

Relationship between AGB, Bcu, Bbr and Ble with dbh a, b, c and d with dint, e, f, g and h, with l i, j, k and l and with dbh2l m, n, o and p

Full size image

In addition, variables, such as N and temperature, showed a significant correlation (p < 0.05) with AGB. This condition was considered when including these variables within multiple linear regression; however, their contribution for fitting was very low.

The fitting of simple linear models shows better results for those where weighted residuals were included. Table 4 presents the parameters and statistics used to assess the performance of the models. Weighted linear Models 2 and 4 explain more variability in AGB and show less error. Model 2 is fitted with dbh, whereas Model 4 is fitted with dbh2l. Although the error of Model 4 is slightly less, the fact of including l as a predictor might represent a limitation given that the measurement of l for guadua is highly time-consuming. In addition, residuals and comparisons between the fitted and observed values are shown in Fig. 3.

Table 4 Fitted linear models (\(AGB={a}_{1}+ {b}_{1} x {X}_{1}+{\varepsilon }_{1}\)), parameters and error
Full size table
Fig. 3
figure 3

Plot of linear weighted model; a and c observed vs. fitted AGB values; b and d residuals vs. fitted AGB values; a and b the predictor variables are dbh and c and b dbh2 l for the linear model and the weighted linear model

Full size image

The weighted log-transformed linear models fitted with dbh and dbh2l showed similar precision and error to the weighted simple linear models described above. Once again, the model with the predictor variable dbh2l showed slightly better performance. The calculated Fln values for Models 5 and 6 were 9.97 and 8.77, respectively. This value is useful for comparisons with other models where Fln corresponds to RMSE. Table 5 shows the statistics of the models assessed.

Table 5 Fitted weighted log-transformed models (\(log(AGB)={a}_{2}+ {b}_{2} xlog( {X}_{2})+{\varepsilon }_{2}\)), parameters and error
Full size table

In addition, the back transformation model (\(AGB=CF x (a x {X}^{b})\)) fitted for Models 5 and 6 with dbh and dbh2l as predictor variables (X) with the parameters and statistics shown in Table 6 and Fig. 4 depicts the variation and distribution of weighted residuals at predicted biomass.

Table 6 Back transformation model (\(AGB=CF x (a x {X}^{b})\)) for the fitted weighted log-transformed models (5 and 6) (\(log(AGB)={a}_{2}+ {b}_{2} xlog( {X}_{2})+{\varepsilon }_{2}\)), parameters and error
Full size table
Fig. 4
figure 4

Plot of back transformation model (\(AGB=CF x (a x {X}^{b})\)) where a) and c are observed vs. fitted AGB values and b and d are residuals vs. fitted AGB values. The predictor X variable is dbh for a and b) and dbh2l for c and d

Full size image

After fitting the weighted power models using the same predictor variables, the precision and error of the models were also sufficient according to the statistics and the variation observed in the predicted values and residuals (Table 7 and Fig. 5).

Table 7 Weighted power models fitted (\(AGB={a}_{3} x {X}_{3}^{b3}+{\varepsilon }_{3}\)), parameters and error
Full size table
Fig. 5
figure 5

Plot of weighted power model (\(AGB={a}_{3} x {X}_{3}^{b3}+{\varepsilon }_{3})\) where a and c are observed vs. fitted AGB values and b and d are residuals vs. fitted AGB values. The predictor variable is dbh for a and b and dbh2l for c and d

Full size image

For the SUR approach, values of error for prediction of AGB resulted with a better performance that de models independently fitted, even using only dbh as a predictor variable (systems 1, 2 and 3). However, the poor fitting effect of Bbr and Ble had a high influence on the performance and the parameter ai of both equations resulted not significant (p > 0.05) in equations systems 1 and 2. To improve this shortcoming, in equations system 3, as predictor variable was included the variable dbh2l in equations for Bbr and Ble. Thus, the performance was better and error decreases; however, the contribution of Bbr and Ble was also poor, although the parameters were significant (p < 0.05). In Table 8, 9 and 10 are shown the parameters and statistics for assessing the models. Furthermore, in Fig. 6 are depicted the observed versus the fitted AGB values and as the residuals versus the fitted AGB values as well.

Table 8 Fitted simultaneous equations system 1, parameters and error
Full size table
Table 9 Fitted simultaneous equations system 2, parameters and error
Full size table
Table 10 Fitted simultaneous equations system 3, parameters and error
Full size table
Fig. 6
figure 6

Plot of the observed vs. fitted AGB values and residuals vs. fitted AGB values from the simultaneous systems of equations fitted by the SUR approach, a and b system 1, c are and d system 2 and e and f system 3. The predictor variable is dbh

Full size image

Cross-validation performed for Models 2, 4, 5, 6, 7, 8 and the simultaneous equations system 3, shows proper precision for estimating AGB and consistent values of error regarding those estimated in the fitting phase (Table 11). For the log-transformed Models 5 and 6, the Fln values calculated after cross-validation were 9.86 and 8.63 (28.9% and 23.9% of RMSEPE), respectively. These values show the lower values of error in comparison with other models. As evidenced in the fitting phase, lower values of error were found when models were fitted with the predictor variable dbh2l, although in the simultaneous equations system 3, the equation for estimating AGB fitted from dbh had the lower RMSE.

Table 11 Cross-validation statistics of the assessed AGB models
Full size table

Models fitted with higher performance from the predicted variables dbh and dbh2 l are consistent with previous results for assessing allometric equations for bamboo (Huy et al., 2019a; Yuen et al., 2017; Zhuang et al., 2015). The inclusion of dbh2 l is also promoted by Picard et al., (2015) because collinearity between both predictors is avoided under this form. However, for practical reasons, dbh is preferable given that it is easier to measure in inventories (Camargo and Arango Arango, 2012) and with other bamboo species (Huy et al., 2019a; Li et al., 2016). Other variables, such as dint and l, showed a lower performance in models and were not suggested for use. In addition, l is a complicated parameter to measure because culms are bent and intermingled. In fact, allometric equations are typically used for estimating culm length (Camargo and Kleinn, 2010).

The models fitted with weighted linear (2 and 4), weighted log-transformed linear (5 and 6) and weighted power models (7 and 8) exhibit similar prediction performance with lower error when the predictor variable is dbh2l (4, 6 and 8). Considering the literature (pe. Huy and Long, 2019), the weighted power model might be suggested. However, conditions and simplicity may justify the application of the other two equations as an alternative. In fact, Picard et al. (2015) suggest seeing the data before suggesting a power model as the best alternative. In addition, the error of the weighted log-transformed linear model was slightly lower. For practical reasons, Models 2, 5 and 7 fitted from dbh should be preferably considered for use given that culm length is not easy to measure.

The environmental factors included in the linear models (e.g., elevation, precipitation, N) did not show a significant effect on the model with the exception of N. However, the model did not substantially improve. Although the coefficient of determination was greater, the error slightly increased. These results are similar to some studies where N was included in models (Devi & Singh, 2021; Xu et al., 2018) and with culm age (Zhang et al. 2014). The relationship between age and biomass was restricted to the dimension of the plant. Therefore, given that diameter does not increase and its length is reached during the first six months (Camargo, 2006), biomass probably has no particular effect. Additionally, Aguirre and Criollo (2020) found that intrinsic variables influenced by management, such as N, have a greater effect on biomass than environmental factors.

The simultaneous equations system 3 showed possibilities for being used for prediction of AGB and Bcu; however, the poor contribution of Bbr and Ble to the performance of the models is a disadvantage. Although the use of simultaneous equations is widely used for estimates biomass (Sanquetta et al., 2015; Kralicek et al., 2017; Huy et al., 2019a; Cui et al. 2020) as a reliable approach, in this study the relationship between Bbr and Ble and the predictor variables used did not contribute to a better performing of models. The characteristics of guadua forest as above described, have incidence on the variability of culms where the proportion of leaves and branches may vary depending on the location (distance to the stand border), maturity and N, independently on the predictor variables used.

10 Conclusions

Given their capacity to store carbon, guadua bamboo forests in the coffee region of Colombia represent a key resource for climate change mitigation. Under natural conditions, these forest ecosystems enhance this capacity, and other pools are usually not included, such as litter, SOC and trees, which additionally represent part of the vegetal biodiversity associated with these forests.

Allometric fitted equations independently and simultaneous system (3) provide proper precision and a good explanation of variability for the prediction of biomass (AGB) according to the conditions assessed and considering the intrinsic high variability associated with these ecosystems. Those models fitted from dbh were easier and more practical to apply in inventories and management plans typically performed for these forests. Therefore, it might be considered as a complementary measure to integrate AGB estimation in inventories and plan management. An interesting advantage was also that the site variability of the natural distribution affected the fitting models; therefore, the lack of precision depended on those intrinsic factors related to the structure of natural stands.