Allometric Models for Estimation of Forest Biomass in North East India

Abstract

In tropical and sub-tropical regions, biomass carbon (C) losses through forest degradation are recognized as central to global terrestrial carbon cycles. Accurate estimation of forest biomass C is needed to provide information on C fluxes and balances in such systems. The objective of this study was to develop generalized biomass models using harvest data covering tropical semi-evergreen, tropical wet evergreen, sub-tropical broad leaved, and sub-tropical pine forest in North East India (NEI). Among the four biomass estimation models (BEMs) tested AGB_est = 0.32(D²Hδ)^0.75 × 1.34 and AGB_est = 0.18D^2.16 × 1.32 were found to be the first and second best models for the different forest types in NEI. The study also revealed that four commonly used generic models developed by Chambers (2001), Brown (1989), Chave (2005) and Chave (2014) overestimated biomass stocks by 300–591 kg tree⁻¹, while our highest rated model overestimated biomass by 197 kg tree⁻¹. We believe the BEMs we developed will be useful for practitioners involved in remote sensing, biomass estimation and in projects on climate change mitigation, and payment for ecosystem services. We recommend future studies to address country scale estimation of forest biomass covering different forest types.

1. Introduction

Land based climate change mitigation strategies have received much global attention in the recent past, due to the large sink capacity and economic viability [1]. Among the terrestrial ecosystems, tropical and sub-tropical forests are considered central to global terrestrial carbon (C) stocks [1,2]. Tropical and sub-tropical regions are well-recognized for losing forests due to agricultural expansion and bio-energy production [3]. Information on tropical forest biomass and C fluxes is gaining both economic and political currency in renewable energy development, C credit markets and research in global environmental change. Since the approval of the REDD⁺ (reducing emissions from deforestation and forest degradation) during the Conference of the Parties (COP 19) in November 2013, the Warsaw Framework has become a formal mechanism for creating performance-based C financing. The Warsaw Framework requires consistency in methods, definitions, and comprehensiveness in measurement, reporting and verifying emissions by sources and removals by sinks, and changes in C storages [4].

In spite of the significant progress in biomass estimation methods, quantifying C stocks in tropical and sub-tropical forests is still challenging. A large degree of uncertainty exists in measured C stocks and fluxes in tropical and sub-tropical forests [5,6,7,8,9]. Some of the uncertainty results from the lack of consistencies in methods, thus leading to widely varying results even among similar studies. Based on remotely sensed data, two independent studies, namely Harris et al. [5] and Baccini et al. [8], published maps of tropical forest C storages, which are widely used in REDD⁺ projects [6]. For 2000–2005, Harris et al. [5] reported gross emission of 0.81 Pg (Petagram, 1 Pg = 10¹⁵ g) C year⁻¹ compared to 2.22 Pg C year⁻¹ by Baccini et al. [8]. Although the maps developed in both the studies used the same LiDAR data from the Geoscience Laser Altimeter System (GLAS) across the tropics, the maps reveal substantial differences in total biomass stocks, with little consistency in the direction of differences [6,9]. Differences of this magnitude are reason for concern, not only in policy formulation, but also global climate change science [9]. According to Mitchard [6] one of the causes of the differences between the two maps is the difference in the biomass estimation models (BEMs) used to estimate biomass from the ground plots. The choice of BEMs can significantly influence local, regional and global biomass estimates. In addition, the choice of BEMs poses a practical limit to the accuracy with which remote sensing methods can predict regional biomass [6]. This highlights the need for the development of region-specific BEMs.

The introduction of C credits under REDD has financial implications to C stock in tropical and sub-tropical forest ecosystems [10]. Therefore, accurate estimation of C loss and sequestration is fundamental to the initiatives of managing forested ecosystems for reducing CO₂ emissions [11] and model applications to climate change studies [12]. Accuracy in estimates of forest C stocks are limited by the challenge of developing robust models to estimate tree biomass [13]. Different direct and indirect methods are used for biomass estimation. In the direct methods, a sample of trees in a given area is harvested and measured for estimation of dry weight in different tree components (e.g., trunk, branches, and leaves). Direct methods can be expensive, especially when dealing with large sample areas and several species [14]. In the indirect method, biomass is usually estimated using BEMs, which relate measurable variables such as total tree height, diameter at breast height, and woody density to total tree biomass [15,16,17]. Therefore, the use of indirect methods is often preferred over direct methods [18].

Multi-species pan tropical models have been developed for estimation of above ground biomass (AGB) for major forest types [16,19]. However, such models may not accurately predict biomass of forests in different ecological regions of the world [20,21]. Species-specific models [16] and LiDAR technology can also be unreliable for application to mixed species stands.

In the tropical and sub-tropical parts of North East India (NEI), forests cover 66% of the total geographical area [22]. These forests have significant influence on regional and national C balance. In NEI, different species-specific models have been developed for Pinus kesiya [23], Hevea brasiliensis [24], and Barringtonia acutangula [25]. Given the uncertainty in biomass estimation, these species-specific models may not have wider application in mixed species tropical forests. Biomass and C stock for diverse forest ecosystems in NEI have also been estimated [26,27,28,29] using various generic models including those developed by Brown et al. [30], Chambers et al. [31] and Chave et al. [16,19]. However, the accuracy of biomass estimates using these models has rarely been tested. Therefore, accurate estimation of biomass and C stock in NEI will form the baseline for regional/national C balance datasets.

The objective of this study was to develop generalized biomass models using harvest data covering tropical semi-evergreen, tropical wet evergreen, sub-tropical broad leaved, and sub-tropical pine forest in NEI. The aim of this study is to present generalized BEMs for NEI for use by practitioners involved in projects on climate change mitigation, and payment for ecosystem services.

2. Methods

2.1. Descriptions of the Study Region

This study covered NEI, which consists of eight different states: Assam, Arunachal Pradesh, Mizoram, Meghalaya, Manipur, Nagaland, Sikkim, and Tripura. NEI covers 26.3 million hectares (M ha) equivalent to 8% of total geographical area of India [32]. NEI is situated at the confluence of the Indo-Chinese, Indo-Malayan, and Indian bio-geographical realms. Due to this unique geographical location, NEI represents numerous forest types falling within one of the biodiversity hotspots of the world, the Indo-Burma biodiversity hotspot [33]. In NEI, 17.2 M ha of land is covered with forests, which constitutes ~25% of India’s total forest area [32] and is represented by five broad forest types based on the elevation, forest structure and composition: (i) tropical semi-evergreen (up to 600 m), (ii) tropical wet evergreen (up to 900 m), (iii) sub-tropical broad leaved (900–1900 m), (iv) sub-tropical pine (1000–3500 m), and (v) alpine temperate (above 3500 m) [34,35].

2.2. Sampling Strategies

Sample tree data (a total of 303 trees) were collected from four major forest types (

Table 1. Dominant and co-dominant tree species sampled for the present study.

Forest Types	Altitudinal Range (Meters)	Species
Alpine Temperate	>3500	Data not available for this study
Sub-Tropical Pine	1000–3500	Pinus kesiya Royle ex Gordon, Pinus roxburghii Sarg.
Sub-Tropical Broad Leaved	900–1900	Schima wallichii Reinw. ex Blume, Quercus oblongata D. Don, Ficus benghalensis L., Machilus gamblei King ex Hook.f., Mallotus philippensis (Lam.) Muller.-Arg., Myrica sapinda Wall., Terminalia myriocarpa Van Heurck & Mull., Terminalia chebula (Gaertn) Retz, Toona ciliata M.J. Roem, Juglans regia L., Alnus nepalensis D. Don
Tropical Wet Evergreen	Up to 900	Tectona grandis L.f., Macaranga denticulata (Blume) Muller. -Arg., Mesua ferrea L., Dipterocarpus turbinatus C.F.Gaertn
Tropical Semi-Evergreen	Up to 600	Albizia procera (Roxb.) Benth., Syzygium cumini (L.) Skeels, Macaranga peltata (Roxb.) Muller. -Arg., Bauhinia variegata L., Artocarpus chama Buch. –Ham.

, Figure 1) accounting for over 90% of the forest cover in NEI. These four forest types span over tropical to sub-tropical climates. Alpine zone that represents temperate forests of NEI accounts for only 5%–7% of total forest cover of NEI. Therefore, this forest type was not included in this study.

To cover diverse tree sizes for harvest, seven diameter classes were formed: 10.1–20 cm, 20.1–30 cm, 30.1–40 cm, 40.1–50 cm, 50.1–60 cm, 60.1–70 cm, and 70.1–90 cm. Tree sizes ≤10 diameter at breast height (D) were not considered for harvest. Sample trees were harvested from a minimum of four different diameter classes for each of the four different forest types. This sampling strategy was adopted because for certain forest types, trees were distributed up to maximum of four diameter classes. Then selected trees were cut at ground level and total height (m) was measured. After felling, tree components were separated in to leaf, branch, and bole, and fresh weight of each component was measured in the field with a digital balance. Sub-samples (1 kg) of each component were collected, taken to the laboratory and dried at 65 ℃ until a constant weight was reached. Then, the fresh weight to dry weight ratio was used to compute the dry mass of the total tree.

2.3. Model Development

Development of an empirical biomass model is subject to the use of appropriate independent variables and covariates that are likely to influence biomass and the chosen variables. Although the selection of independent variables still remains a matter of debate [17], D and total tree height (H) are the most commonly used variables when developing BEMs [16,17]. The compound forms of D and H with or without wood density (δ) are also widely used in BEMs [13,19]. The conventional approach in the development of BEMs has been to use models with only fixed effects, i.e., without considering covariates. As biomass estimates may vary not only with the fixed effects, but also covariates such as forest type, stand density, site quality, competition, etc., it is important to include such covariates as sources of additional variation in the model. In this study, we used both conventional and a linear mixed modeling (LMM) framework where covariates are considered in model development. The general form of the LMM model is given as:

$${\mathrm{Y}}_i=\beta {\mathrm{X}}_i+u_i{Z}_j+{\epsilon }_i.$$

(1)

where Y_i is the n-dimensional response vector, β is the p-dimensional parameter vector for fixed-effects, X_i is the nxp design matrix for fixed-effects, u_i is the q-dimensional vector of parameters for random effects; Z_i is the nxq design matrix for random-effects, ε_i is the n-dimensional error vector [36]. It is assumed that the expectation (E) is E(ε_i) = E(u_i) = 0, and the variances (V) and co-variances (Cov) are: Var(u_i) = R_i, Var(u_i) = D_i, Cov(ε_i, u_i) = 0, respectively. Furthermore, ε_i and u_i are assumed to be normally distributed as: ε_i ~ N(0,R_i) and u_i ~ N(0,D_i) [36].The rationale for the LMM is that both the fixed and random parameters can be entered in the model simultaneously, thus providing consistent estimates of parameters and their standard errors than the conventional method [37].

In this study, we had only information on forest type, but information was lacking on stand density, site quality and competition. Initially, we entered a random intercept and forest type as random effects in the model, but this created problem with convergence of algorithms and testing the significance of the random intercept. Although convergence criteria were met, the final hessian was not positive definite when the intercept was entered as a random effect. Therefore, we used only forest type as the random effect. Here, four forest types, namely, sub-tropical pine, sub-tropical broad leaved, tropical wet evergreen and tropical semi-evergreen forest were considered in the model.

Both the conventional and LMM versions of four commonly used BEMs [17] were compared here. In the conventional framework, the four models are formulated as follows:

Model 1:

$\ln \left(AGB\right)=\ln \left(a\right)+bln\left(D\right)+\epsilon$;

Model 2:

$\ln \left(AGB\right)=\ln \left(a\right)+bln\left(D^{2}H\right)+\epsilon$;

Model 3:

$\ln \left(AGB\right)=\ln \left(a\right)+bln\left(D^{2}H\delta \right)+\epsilon$;

Model 4:

$\ln \left(AGB\right)=\ln \left(a\right)+bln\left(D\right)+cln\left(H\right)+dln\left(\delta \right)+\epsilon$

where AGB is total above-ground biomass in kg dry matter tree⁻¹, a is the intercept, b, c and d are slope parameters, D is diameter at breast height in centimeters, H is total height in meters, δ is wood density (or specific gravity) in kg m⁻³ and ε is the random error. In all models we had H measurements, but δ was obtained from Global Wood Density Database [38]. Upon back transformation, each model was multiplied with a correction factor (CF) estimated from the residual variance or mean square error (MSE) using the formula CF = exp(MSE/2) [39]. The CF accounts for the back transformed of the error and a step in log-transformed data in BEMs [17].

In the LMM framework, the models were formulated as follows:

Model 1:

$\ln \left(AGB\right)=\ln \left(a\right)+bln\left(D\right)+uF+\epsilon$;

Model 2:

$\ln \left(AGB\right)=\ln \left(a\right)+bln\left(D^{2}H\right)+uF+\epsilon$;

Model 3:

$\ln \left(AGB\right)=\ln \left(a\right)+bln\left(D^{2}H\delta \right)+uF+\epsilon$;

Model 4:

$\ln \left(AGB\right)=\ln \left(a\right)+bln\left(D\right)+cln\left(H\right)+dln\left(\delta \right)+uF+\epsilon$

where AGB, D, H, δ, a, b, c, d and ε are defined as in the conventional. The added quantity is the parameter u, which is the random effect of forest type (F). In all LMM models, it is assumed that u and ε are uncorrelated, random and normally distributed variates with mean 0 and variance of 1. In these models, there are two variance components; the variance associated with the random effect (${\sigma }_u^2$) and the residual variance (${\sigma }_{\epsilon }^2$).

To ensure comparability of parameters, we estimated parameters of both the conventional and LMM models using the PROC MIXED procedure in the SAS system. The SAS codes used for estimation of parameters of the conventional and LMM models are given in

Table A1. PROC MIXED codes for estimating parameters of the conventional models.

/*Code for fitting Model 1*/ Proc mixed data = Biomass method = REMLcovtest ratio ic; Class Forest; Model lnAGB = lnD/solution outp = check cl; Run;

/*Code for fitting Model 2*/ Proc mixed data = Biomass method = REMLcovtest ratio ic; Class Forest; Model lnAGB = lnDDH/solution outp = check cl; Run;

/*Code for fitting Model 3*/ Proc mixed data = Biomass method = REMLcovtest ratio ic; Class Forest; Model lnAGB = lnWDDH/solution outp = check cl; Run;

/*Code for fitting Model 4*/ Proc mixed data = Biomass method = REMLcovtest ratio ic; Class Forest; Model lnAGB = lnDlnHlnW/solution outp = check cl; Run;

and

Table A2. PROC MIXED codes for estimating parameters of the LMM models.

/*Code for fitting Model 1*/ Proc mixed data = Biomass method = REMLcovtest ratio ic; Class Forest; Model lnAGB = lnD/solution outp = check cl; Random Forest; ESTIMATE “1” intercept 1| Forest 1; ESTIMATE “2” intercept 1| Forest 0 1; ESTIMATE “3” intercept 1| Forest 0 0 1; ESTIMATE “4” intercept 1| Forest 0 0 0 1; Run;

/*Code for fitting Model 2*/ Proc mixed data = Biomass method = REMLcovtest ratio ic; Class Forest; Model lnAGB = lnDDH/solution outp = check cl; Random Forest; ESTIMATE “1” intercept 1| Forest 1; ESTIMATE “2” intercept 1| Forest 0 1; ESTIMATE “3” intercept 1| Forest 0 0 1; ESTIMATE “4” intercept 1| Forest 0 0 0 1; Run;

/*Code for fitting Model 3*/ Proc mixed data = Biomass method = REMLcovtest ratio ic; Class Forest; Model lnAGB = lnWDDH/solution outp = check cl; Random Forest; ESTIMATE “1” intercept 1| Forest 1; ESTIMATE “2” intercept 1| Forest 0 1; ESTIMATE “3” intercept 1| Forest 0 0 1; ESTIMATE “4” intercept 1| Forest 0 0 0 1; Run;

/*Code for fitting Model 4*/ Proc mixed data = Biomass method = REMLcovtest ratio ic; Class Forest; Model lnAGB = lnDlnHlnW/solution outp = check cl; Random Forest; ESTIMATE “1” intercept 1| Forest 1; ESTIMATE “2” intercept 1| Forest 0 1; ESTIMATE “3” intercept 1| Forest 0 0 1; ESTIMATE “4” intercept 1| Forest 0 0 0 1; Run;

, respectively. PROC MIXED uses the restricted maximum likelihood (REML) estimation method. The REML method is known to produce consistent estimates of the variance-covariance matrix [37] and asymptotic standard errors of the covariance parameter estimates.

The performance of the different models was assessed using the bias-corrected Akaike information criterion (AICc) and the coefficient of determination (R²). In the LMM framework, the R² is not routinely estimated due to theoretical problems in its definition and practical difficulties in implementation [40]. As such, PROC MIXED does not normally report it. Therefore, R² was estimated as: R² = 1 − RV_m/RV_i, where RV_m is the residual variance of the full model, and RV_i is the residual variance of empty (intercept-only) model [40].

LMMs are based on the theory of empirical best linear unbiased predictors (EBLUPs) of the random effects and the best linear unbiased estimates (BLUEs) of the fixed effects [37]. The advantage of EBLUPs estimation is that it does not require normality of the random effects [41]. As such, LMMs provide an efficient approach to small area (or domain) estimation by incorporating random effects that account for dissimilarities between domains (forest type in our analysis). Taking advantage of this property, localized AGB predictions that are specific to each forest type were generated using the intercepts of the different forest types for the two highly ranked LMM model. Finally, the standardized residuals generated using the linear mixed effects model were plotted against the independent values to check for heteroscedasticity in residuals. Outliers were detected by checking the studentized residuals, and values below −2 or +2 were considered as outliers [17]. Heteroscedasticity will typically be manifested as residuals whose magnitude is correlated with that of the response variable and a plot of the residuals against the predicted values will reveal a megaphone pattern if the errors are not homogenous [17].

2.4. Model Validation

Cross validation is usually recommended to determine how accurately BEMs will perform when applied to an independent dataset. Usually 5-fold or 10-fold cross validation provides a good balance between bias and variance [17]. Therefore, a 10-fold cross-validation was employed to evaluate the predictive performance of the four selected models. The goodness of fit criteria was calculated for the validation dataset using the lava and forecast packages of the R package (Appendix A

Table A3. R scripts for the four models.

Model1<-lm(lnAGB ~ lnD, data = Biomass_data) Model2<-lm(lnAGB ~ lnDDH, data = Biomass_data) Model3<-lm(lnAGB ~ lnWDDH, data = Biomass_data) Model4<-lm(lnAGB ~ lnD + lnH + lnW, data = Biomass_data)

The R script for cross validation using 10-fold using the lava packages: cv(list(Model1, Model2, Model3, Model4), k = 10, data = Biomass_data)

Various goodness of fit indices were calculated using the following codes using the “forecast” package of R: Model1_results < -t(data.frame(CV(Model1))) Model2_results < -t(data.frame(CV(Model2))) Model3_results < -t(data.frame(CV(Model3))) Model4_results < -t(data.frame(CV(Model4))) Model_results < -rbind(Model1_results, Model2_results, Model3_results, Model4_results)

). Specifically, the adjusted coefficient of determination (Adj R²), root mean square of error (RMSE), the AICc and Bayesian information criterion (BIC) were estimated from the 10-fold cross validation using the forecast packages of R (see Appendix A).

2.5. Comparison with Generic Models

We also compared our highest ranked LMM model with the following commonly used generic models developed by Brown et al. (1989) [30], Chamber et al. (2001) [31], Chave et al. (2005) [16] and Chave et al. (2014) [19] for broadly defined forest types. We chose these models for comparison with our own models because they are commonly used for biomass estimation in NEI by previous researchers [26,27,28,29].

Chamber’s model is given as

$$AGB=exp\lfloor -0.370+0.333\mathrm{lnD}+0.9333ln{D}^2-0.122ln{D}^3\rfloor$$

(2)

The model proposed by Brown et al. (1989) for tropical regions is:

$$AGB=13.2579-\left(4.8945D\right)+0.6713D^2$$

(3)

Chave’s model 1 was based on his equation for tropical wet forests (Chave et al. 2005) given as:

$$AGB=0.0776{\left(\rho D^{2}H\right)}^{0.940}$$

(4)

Chave’s model 2 is Equation (4) of Chave et al. (2014) said to be the best-fit pan tropical model. It is given as:

$$AGB=0.0673{\left(\rho D^{2}H\right)}^{0.976}$$

(5)

We determined the appropriateness of these generic models by comparing with our highest ranked model using the R², RMSE, relative prediction error (Error), mean absolute percentage error (MAPE) and AICc. We also compared the 95% CI (confidence interval) of the slopes (b) of the regression of measured against fitted values to assess whether or not significant prediction errors exist. If significant prediction errors exist b ≠ 1, and the 95% CI of b will not cover 1 [17].

3. Results

Parameter estimates of models 1–4 generated using the conventional method and the linear mixed effects modeling (LMM) framework are given in

Table 2. Comparison of models fitted using the conventional (CONV) method with models fitted within a linear mixed effects (LME) framework where forest type was used as a random effect. The Akaike information criterion (AIC), adjusted R² (Adj R²) and RMSE were used for comparing CONV with LMM.

		Model Parameters					Adj
Model		ln(a) (SE) †	b (SE)	c (SE)	d (SE)	AIC ‡	R2	RMSE
1	CONV	−2.12 (0.34)	2.32 (0.10)			769	0.62	0.851
	LME	−1.73 (0.45)	2.16 (0.10)	-		705	0.705	0.749
2	CONV	−2.30 (0.35)	0.82 (0.04)			771.7	0.619	0.852
	LME	−1.64 (0.44)	0.74 (0.04)	-		734.5	0.675	0.785
3	CONV	−1.83 (0.32)	0.82 (0.04)			770.3	0.621	0.850
	LME	−1.25 (0.43)	0.75 (0.04)	-		725.7	0.685	0.773
4	CONV	−2.12 (0.39)	2.00 (0.16)	0.43 (0.16)	0.37 (0.34)	763.2	0.627	0.842
	LME	−1.21 (0.49) ns	2.22 (0.16)	−0.08 (0.16) ns	0.87 (0.30)	698.8	0.710	0.740

^† Figures in parentheses are asymptotic standard errors (SE) of parameter estimates; ^‡ AIC values in bold face indicate the better model when the conventional was compared with LMM; ns: not significant: Note: Model 1: ln(AGB) = ln(a) + b*ln(D); Model 2: ln(AGB) = ln(a) + b*ln(D²H); Model 3: ln(AGB) = ln(a) + b*ln(D²Hδ); Model 4: ln(AGB) = ln(a) + b*ln(D) + c*ln(H) + d*ln(δ).

. The AICc and adjusted R² show that models 1–4 fitted using LMM are superior to those fitted using the conventional method (Table 2). Therefore, all inferences hereafter will be based on the models fitted using the mixed effects framework. Estimates of the covariance parameters for the different LMM models are given in

Table 3. Estimates of covariance parameters and their significance for the different models.

Model §	Covariance Parameter	Variance Estimate	Z Value †	p Value ‡
1	Forest	0.38 (0.32)	1.18	0.1187
	Residual	0.56 (0.05)	12.21	<0.0001
2	Forest	0.30 (0.26)	1.16	0.1236
	Residual	0.62 (0.05)	12.20	<0.0001
3	Forest	0.33 (0.28)	1.17	0.1218
	Residual	0.60 (0.05)	12.20	<0.0001
4	Forest	0.45 (0.38)	1.18	0.1182
	Residual	0.55 (0.05)	12.16	<0.0001

^† The Z value is the Wald Z-test for covariance parameter estimates; ^‡ p value represents the significance of the Z test of the hypothesis that the variance of the effect is 0; ^§ Specification of models is as in Table 2.

. In all models the variance of random effect was not significantly different from zero, while the residual variance was significantly larger than zero (Table 3). The largest residual variance was recorded in Model 2.

Table 4. Comparison models using cross validation goodness of fit criteria including the adjusted coefficient of determination (Adj R²), root mean square of error (RMSE), the bias-corrected Akaike information criterion (AICc) and Bayesian information criterion (BIC).

Model	Adj R2	RMSE	AICc	BIC
1	0.62	0.728	−93.7	−82.7
2	0.62	0.729	−93.1	−82.0
3	0.62	0.726	−94.4	−83.4
4	0.63	0.713	−98.5	−83.8

Note: Model 1: ln(AGB) = ln(a) + b*ln(D); Model 2: ln(AGB) = ln(a) + b*ln(D²H); Model 3: ln(AGB) = ln(a) + b*ln(D²Hδ); Model 4: ln(AGB) = ln(a) + b*ln(D) + c*ln(H) + d*ln(δ).

provides goodness-of-fit statistics of the cross validation of models 1–4. Models 1–3 were indistinguishable in terms of RMSE, AICc, BIC and R² (Table 4). Although the cross-validation goodness of fit criteria indicates that model 4 is the best, the mixed model analysis revealed that parameters a and d in Model 4 are not significantly different from zero. Therefore, Model 3, the one with the next smaller RMSE, AICc and BIC, was chosen as the best model. Model 3 is given as $AG{B}_{est}=0.32{({\mathrm{D}}^2\mathrm{H}\mathsf{\delta })}^{0.75}\times 1.34$ in the arithmetic domain where 1.34 is the correction factor. The second best option would be the power law model given as $AG{B}_{est}=0.18D^{2.16}\times 1.32$ in the arithmetic domain. The AGB values predicted using the different models are shown in Figure A1. In all models, the residuals did not reveal heterogeneity of variance (Appendix A Figure A1).

The LMM estimates of the intercept (and SE) of Model 3 were −1.01 (0.34) for sub-tropical pine, −0.83 (0.36) for sub-tropical broad leaved, −1.10 (0.31) for tropical wet evergreen, and −2.08 (0.31) for tropical semi evergreen forest. The predictions produced using these estimates for each forest type are presented in Figure 2. For Model 1, the forest type-specific intercept were −1.34 (0.34) for sub-tropical pine, −1.46 (0.37) for sub-tropical broad leaved, −1.49 (0.31) for tropical wet evergreen, and −2.64 (0.32) for tropical semi-evergreen forest, and corresponding the predictions for each forest type are presented in Figure 3. The discrepancies in predictions between the conventional and the LMM models were largest in tropical wet evergreen forest (Figure 2d), and the discrepancies increased with tree diameter (Figure 3a). On the other hand, the smallest discrepancy was observed in tropical semi-evergreen forest subtropical pine forest (Figure 2d and Figure 3d).

With high R² and low AICc, RMSE and MAPE values, our highest rated model was also superior to the generic model developed by Brown and the two pan-tropical models developed by Chave (

Table 5. Comparison of the predictions of our highest rated model with the generic models using the various goodness of fit criteria and the slope (b) of the measured above ground biomass (AGB) with the fitted values.

Model	R2	RMSE	Error	MAPE	AICc	b (95% CL)
Our highest rated model	0.869	141	197.4	235.7	3001.8	1.06 (1.01–1.11)
Chamber’s model	0.870	140	591.4	595.4	2997.9	0.33 (0.32–0.35)
Brown’s model	0.848	151	303.6	314.8	3045.7	0.59 (0.56–0.61)
Chave’s model 1	0.829	161	299.6	313.2	3081.8	0.42 (0.40–0.45)
Chave’s model 2	0.821	165	372.8	382.7	3096.5	0.32 (0.31–0.34)

). Although our highest rated model was comparable with Chamber’s model in terms of R², AICc, and RMSE, it was superior in terms of MAPE and prediction error (Table 5). With b < 1 all the generic models also have significant prediction errors (Table 5). Compared to our highest rated model, Chamber’s model and Chaves model 2 severely overestimated AGB especially for trees with DBH exceeding 40 cm (Figure 4a). The largest average deviation from our model prediction (595.5 kg tree⁻¹) was recorded with Chamber’s model, while the lowest (191.4 kg tree⁻¹) was recorded with Brown’s model (Figure 4b). However, the deviations increased with increasing tree diameter (Figure 3b). Irrespective of tree diameter, when errors in estimation were evaluated, Chamber’s, Brown’s, Chave Model 1 and Chave Model 2 overestimated biomass stock by 591.4, 303.6, 300 and 372.8 kg tree⁻¹, respectively, while our highest rated model overestimated biomass stock by 197.4 kg tree⁻¹.

4. Discussion

Our analysis indicated that models 1–4 fitted using the LMM approach is superior to those fitted using the conventional method. This is probably because mixed effects models can account for the clustered or nested structure of data in forestry. Interestingly, the discrepancies in predictions between the conventional and the LMM models were large in tropical wet evergreen forest.

Our analysis also indicated that the models with wood specific gravity (Model 3 and 4) fitted the data better than models without it although wood specific gravity was obtained from the literature. This is in disagreement with analyses by Stegen [42] that shows inconsistent relationship between forest biomass and wood specific gravity. Although recent analysis has also shown that wood specific gravity has a very weak contribution to AGB [43,44], we believe its inclusion can improve prediction. The use of fewer explanatory variables has been recommended for ease in model application and validation [17]. In the present study D and H were directly measured and could be included in BEMs. In the study area model 3 was found to be more appropriate than Model 1, 2 and 4. Although the cross validation shows that Model 4 is slightly better than Model 3, parameters a and c of Model 4 were not significantly different from zero. Therefore, Model 4 cannot be reliably used for predictive purposes. The limitation of our highest rated model is that it requires H and δ data, which are often not available in many situations. In situations where measured height (H) and wood density (δ) are lacking, second highest rated model, i.e., $AG{B}_{est}=0.18D^{2.16}\times 1.32$, may be used for biomass estimation in NEI.

Although species-specific BEMs have often been applied across multiple sites, they are not necessarily applicable to other species, especially those with differing wood densities [16] and plant functional types [45]. The generic models [16,19,30,31] currently in use in NEI severely overestimated the biomass stock when applied to our dataset. This suggests that our model is better suited for AGB estimation for diverse forest types in NEI. The differences in biomass stock estimation between our highest rated model and the generic models may be attributed to differences in the forest types. The generic models were widely used because of the high goodness of fit to the sample data used to develop them. Models usually have a grossly inflated performance in-sample compared to their performance in follow-up studies. This is called the winner’s curse [46] in statistics literature. It must be noted that a good model fit does not necessarily translate into good predictions of AGB at the landscape level or outside the study area [47]. The ability of a model to describe the data at hand (in sample fit) is sometimes confused with predictive power (out-of-sample fit). The value of cross validation is to avoid bias in predictions of biomass [45]. Cross validation also provides a better method for model assessment as it estimates how accurately a model will perform when applied to an independent dataset. It will also curtail problems such as over-fitting [17]. Therefore, we caution against overdependence on generic models without validating their suitability in new areas.

5. Conclusions

Availability of appropriate biomass models for species-diverse forest ecosystems was a major constraint in the study of C balance in NEI. The models developed in the present study can be applied for multi-site and multi-species evaluation of stand biomass and C assessment at local and regional scale. Such models can also find application in studies of plant community and forest structure, and remote sensing methods. Future studies may incorporate harvest data from Alpine zone with larger data set for country scale estimation and validation. We also recommend future studies to address country scale estimation of forest biomass covering different forest types.