Abstract

Over the last few decades, remote sensing has emerged as a dependable and cost-effective method for collecting precise data on forest biophysical parameters, aiding in sustainable forest management and global initiatives to combat climate change. This research aimed to develop a model for estimating the above-ground biomass (AGB) of Teak (Tectona grandis L. F.) by combining field measurements with Sentinel-2 earth observation data. The study took place in 36-year-old teak plantation areas within the Sagarnath Forest Development Project in Nepal’s Sarlahi district. Field measurements were conducted using a destructive systematic sampling method, employing 10 × 10 m2 sample plots, and the volume of logs was determined using Newton’s formula. A total of 30 sample plots were used for calibration, while 10 were utilized for validation purposes. The findings revealed that the average AGB per plot was 814 kg (equivalent to 81.4 t ha−1), with a minimum value of 716 kg (71.6 t ha−1) and a maximum value of 1,060 kg (106 t ha−1). The study utilized five independent variables, namely, the Red band, Green band, Blue band, near-infrared (NIR), and normalized difference vegetation index (NDVI) values from Sentinel-2 imagery data, to develop estimation models. Among the 12 models examined, model M10 proved to be the best fit for accurate AGB estimation (adjusted R2 = 0.9809, RMSE = 0.01269, AIC = −170, and p-value = < 8.39e−21). The equation of the best-fitted model was ln (AGB) = A + B × Red + C × Green + D × Blue2 + E × ln (NIR) + F × ln (NDVI), providing an accurate estimate of AGB. Model validation involved a t-test comparing the observed and calculated AGB values for ten sample plots, demonstrating no significant difference (p-value = 0.3662 > 0.05). This model has the potential to facilitate AGB biomass calculations and carbon stock estimates for teak plantations of similar age groups.

1. Introduction

Teak, scientifically known as Tectona grandis L. F., is a highly valuable species of tropical wood that has been cultivated extensively across tropical Asia, far from its original natural range encompassing India, Myanmar, Laos, and Thailand. As of 2010, teak plantations covered a vast area of 6.9 million hectares across 60 countries, with Asia accounting for 88%, Africa for 8%, and the tropical Americas for 4%. Myanmar boasts the world’s largest native teak forest, followed by Thailand [1]. Teak plays a significant role in the Sagarnath Forest Development Project (SFDP), meeting nearly half of the nation’s demand [2]. Due to its rapid growth, teak plantation forests can effectively absorb substantial amounts of carbon dioxide, mitigating the impacts of global warming [3]. Biomass refers to the total mass of organic material, both living and dead, which includes above and below-ground living biomass, leaf litter, and decaying matter [4]. Accurately estimating above-ground forest biomass (AGB) is crucial for assessing carbon stocks, as the above-ground portion generally contains the largest carbon reservoir [5]. The most precise method of calculating biomass data is through the use of allometric equations or destructive sampling in the field [6]. However, these techniques are often expensive, time-consuming, labor-intensive, limited in spatial distribution, and impractical, particularly in remote areas [6, 7]. Remote sensing technologies have now emerged as the primary means of estimating AGB, as they can overcome the limitations of traditional field measurement techniques [7].

Both direct and indirect approaches can be employed to estimate the biomass of an individual tree or a stand. The direct approach, which involves cutting plants into pieces and weighing them to determine biomass, is undoubtedly the most accurate method. However, it is also the most expensive, time-consuming, and labor-intensive [8]. As a result, its application is often limited to small areas and sample sizes [9]. Precise and effective estimation of forest biomass is crucial for implementing carbon emission reduction incentives such as the clean development mechanism [10]. The importance of forests in reducing greenhouse gas emissions has gained increased attention due to growing concerns about climate change and global warming. The international program known as REDD+ (reducing emissions from deforestation and forest degradation and the role of conservation, sustainable management of forests, and enhancement of forest carbon stocks in developing countries) offers developing nations like Nepal a significant opportunity to mitigate climate change, support livelihoods, and protect forest ecosystem biodiversity. The successful implementation of REDD+ relies on accurate assessments of biomass and carbon stocks, as well as the subsequent monitoring of the forest carbon pool [11]. Developing reliable allometric correlations between biomass and other parameters is essential, as they can be applied over large areas to obtain precise estimations [12].

Since the 1980s, several organizations, institutions, and private forests have participated in biomass studies in Nepal. The Forest Research and Training Centre is the primary government agency responsible for conducting biomass studies and developing allometric equations in the country. Biomass investigations and the creation of allometric equations for many tree species in the mid-hills have been carried out by projects such as the Nepal-Australia Forestry Project and the Forest Resource Information System Project [13]. Linear and nonlinear algorithms were developed for biomass estimation using Sentinel-2 products, and the developed models were evaluated using various statistical measures. It is anticipated that this study will serve as a valuable tool for forest managers, forest users, and researchers and contribute to Nepal’s participation in international carbon trade programs. According to Lu [7], the combination of field data and high-resolution remote sensing is crucial for accurate biomass estimation, a view also supported by the research of Gautam et al. [14]. Both authors agree that data integration provides an accurate, precise, and cost-effective monitoring approach for tropical forests. However, collecting field data in remote locations can be time-consuming and labor-intensive [15]. AGB serves as a significant indicator of carbon sequestration in terrestrial forests [16]. Estimating AGB is of primary concern for forest management and policymakers due to the alarming rates of deforestation, fragmentation, and tropical biomass destruction, emphasizing the need to evaluate remaining forest patches and prioritize conservation efforts [17]. Accurate calculations of forest biomass are necessary for greenhouse gas inventories and terrestrial carbon accounting [16].

Biomass mapping across extensive areas provides a detailed understanding of carbon sequestration. Typically, assessments focus on AGB as it constitutes the largest portion of the forest’s total living biomass and is directly affected by deforestation and degradation [18]. Continuous monitoring, assessment, biomass analysis, and mapping of forest regions are essential for biodiversity conservation and sustainable forest resource management [19]. Compared to approaches based on forest inventory plots and extrapolation methods, indirect methods such as the regression equation method offer greater benefits and require less time for assessment. These methods are adaptable and allow for the analysis of AGB fluctuations over time or between different geographic areas, including the effects of disturbances like fire or deforestation [20]. The plantation forest at SFDP is well-established and productive, but there is currently no valid biomass table available to calculate carbon stocks for this species. Determining the biomass of standing trees involves measuring tree components (stem, branches, and leaves) after harvesting, which is time-consuming, expensive, labor-intensive, and tedious. Several universities and research organizations in Nepal have used equations proposed by Sharma and Pukala [21], Tamrakar [22], Chave et al. [10], and other authors. However, these equations have been found to have significant flaws [23]. This study aims to address these shortcomings by developing allometric equations for calculating AGB using Sentinel-2 imagery data. It will also provide a baseline for resource managers at SFDP to estimate the biomass of the teak forest and facilitate proper supervision and monitoring of forest products. The general objective of this study is to model the AGB of teak (T. grandis L. F.) using field measurements and Sentinel-2 imagery. Specific objectives include estimating the AGB of the teak plantation area through field measurements, identifying and correlating remotely sensed tree parameters with calculated AGB, and developing a regression model between AGB and the Sentinel-2 imagery database.

2. Materials and Methods

2.1. Study Area

The research was conducted within the teak plantation area, covering a land area of 7.74 hectares, located in the SFDP in the south-eastern Terai region, specifically the Sarlahi district (Figure 1). The project has undertaken extensive plantations of fast-growing species, including Eucalyptus camaldulensis, Dalbergia sissoo, T. grandis, and other species, encompassing approximately 13,512 hectares [23]. The terrain in the area is characterized by a gentle slope, with an elevation ranging from 60 to 330 m above mean sea level. The study area spans from coordinates 26.9974°N to 85.6749°E. During summer (April/May), the temperature fluctuates between 35 to 45°C, while in the dry winter month of January, it ranges from 10 to 15°C. The average annual rainfall in the region falls between 1,130 and 2,680 mm [24]. The climate of the area is classified as lower tropical, consisting of a mixed hardwood tropical forest, with Shorea robusta G. f. being the dominant species.

2.2. Data Collection

For the purpose of the study (Figure 2), both primary and secondary data were gathered. Secondary data and information were sourced from various internet resources, books, journals, articles, and reports. On the other hand, primary data was obtained through field observations, tree and sample measurements, processing of Sentinel-2 image data, and laboratory analysis.

2.3. Sampling and Measurements

Field measurements were conducted using the systematic sampling method with a sample plot size of 10 × 10 m2. The destructive sampling approach was employed, involving the felling of individual trees across a wide range of diameters at breast height (DBH) and subsequently separating them into boles, branches, and leaves. To obtain measurements, the felled trees were divided into sections with a minimum taper, allowing for the measurement of diameters (at the bottom, middle, and top) as well as lengths of each section. Additionally, three sample discs were collected from each sample plot, representing different parts of the tree (above stump height, top of the tree, and base of branches). These sample discs were obtained from trees located near or at the center of the sample plot. To prevent moisture loss, the sample discs were securely stored inside plastic bags.

The collected sample discs from above the stump, top of the tree, and base of branches were further processed into wood samples measuring 2 × 2 × 2 cm³. The quantities obtained were 10, 10, and 5 samples, respectively. Leaves were also collected from designated plots measuring 0.5 × 0.5 m2, located at the center of each sample plot. The leaves were weighed and then dried in a laboratory. The biomass of the leaves was calculated by multiplying the total number of leaves in a plot by the weight-to-dry ratio. In total, 40 sample plots were established, with 30 of them utilized to calibrate the allometric relationship, while the remaining 10 sample plots were set aside for validation purposes.

2.4. AGB Estimation

To calculate the volume of the main stem, Newton’s formula, as presented in Equation (1) by Mandal et al. [23], was applied. In the volume estimation, the stump height was not considered, and the DBH was measured, including the bark. The formula for volume calculation is as follows:

In Equation (1), S1 represents the basal area of the larger (base) end of the section, Sm represents the basal area of medium (between) portion, and S2 represents the basal area of the smaller (top) end of the section. L denotes the length of the section. Wood density, which is the ratio of the oven-dry mass of a wood sample to its green volume (basic specific gravity), was determined using Equation (2) proposed by Chave et al. [25].

To calculate the AGB of the wood portion, Equation (3) by Moradi et al. [26] was utilized. In Equation (3), AGB represents the above-ground biomass, V represents the volume derived from Newton’s formula, and WD represents the wood density (g cm−3). The total AGB, excluding the stump height, was obtained by summing the biomass of the stem, branches, and leaves.

2.5. Empirical Modeling Approaches

Multiple regression analysis was employed to develop biomass models, with the AGB serving as the dependent variable and the Sentinel-2 image data as the independent variables. The explicit model structures were established using the field inventory data for AGB and various sentinel image data, including spectral bands, vegetation indices, textural images, and subpixel information, as the independent variables.

2.6. Remote Sensing Data

The AGB in the teak plantation area was estimated using Sentinel-2 imagery with a resolution of 10 m. This remote sensing data was of excellent quality and had minimal or no presence of clouds and cloud shadows. The Sentinel-2 image data was obtained from the United States Geological Survey (USGS) Earth Explorer data portal (https://earthexplorer.usgs.gov). Sentinel-2 imagery data is equipped with a multispectral imager that captures 13 spectral bands, with spatial resolutions ranging from 10 to 60 m. For this particular study, only the bands with a spatial resolution of 10 m were utilized for analysis, namely B2 (490 nm), B3 (560 nm), B4 (665 nm), and B8 (842 nm), while the remaining bands were excluded.

2.7. Validation/Accuracy Assessment

The estimated models were evaluated using root mean square error (RMSE), coefficient of determination (R2), and akaike information criterion (AIC).

2.8. Data Analysis

ArcGIS software (Version 10.8) was utilized for creating a study area map and calculating vegetation indices using the Sentinel-2 imagery data. MS Excel was employed for recording field data and performing calculations for tree volume and biomass. The R-Programme and SPSS (version 23) were employed for generating various models.

2.9. Lab Analysis

The weight of wood samples measuring 2 × 2 × 2 cm3 and the quantity of leaves collected from the sample plot were measured both before and after being dried in the laboratory. The leaves were dried at a temperature of 75°C until three consecutive measurements indicated a constant weight. Similarly, the wood samples from the stem and branches were dried at a temperature of 80°C until three consecutive measurements indicated a constant weight, as recommended by Chapagain et al. [27].

3. Results

3.1. Field AGB Estimation

The results revealed that the mean stem volume per plot was 0.855 m3, ranging from a minimum of 0.726 m3 to a maximum of 1.148 m3 (Figure 3). The average oven-dry wood density of the stem was 696 kg m−3, with a minimum of 676 kg m−3 and a maximum of 714 kg m−3 (Figure 4). Analysis of Figure 5 indicated that the average branch volume per plot was 0.313 m3, with sample plot 8 having the lowest volume of 0.250 m3 and sample plot 28 having the highest volume of 0.385 m3. Sample plots 19 and 23 exhibited the lowest branch oven-dry wood density of 627 kg m−3, while sample plot 12 had the highest density of 660 kg m−3. The average oven-dry wood density for branches was 643 kg m−3 (Figure 6). Additionally, the mean weight of fresh leaves per plot was 45 kg, with a minimum of 35.4 kg observed in sample plot 21 and a maximum of 56.2 kg observed in sample plot 24 (Figure 7). Furthermore, Figure 7 illustrates that the average oven-dry weight of leaves per plot was 17.51 kg, ranging from a minimum of 13.81 kg to a maximum of 23.60 kg.

The average AGB measured in the field was 814 kg per plot (equivalent to 81.4 t ha−1), ranging from a minimum of 716 kg per plot (71.6 t ha−1) to a maximum of 1,060 kg per plot (106 t ha−1) (Figure 8). Among the sample plots, 14 (47%) had AGB values ranging from 700.1 to 800 kg, while 11 (37%) and 4 (13%) sample plots had AGB values between 800.1–900 kg and 900.1–1,000 kg, respectively (Figure 9). None of the sample plots had AGB values equal to or less than 700 kg, and only one sample plot (3%) had an AGB value exceeding 1,000.1 kg.

Among the 30 sample plots, the average biomass was 814 kg with a standard deviation of 14.33 kg, and the range of AGB values varied from the lowest value of 716 kg to the highest value of 1,060 kg (Table 1). The median AGB value was 801.7 kg per plot, and the most commonly observed AGB value was 810 kg per plot. The range of AGB values observed in the field was 344 kg per plot, with a standard deviation of 78.5 kg.

3.2. Development of Biomass Model

Multivariate regression analysis was employed to establish a relationship between AGB and independent variables derived from Sentinel-2 earth observation data, including various bands and vegetation indices. The chosen independent variables were the Red band, Green band, Blue band, near-infrared (NIR), and normalized difference vegetation index (NDVI) values obtained from the Sentinel-2 imagery data. These variables were utilized to estimate the AGB. The models exhibited a positive correlation and significance with AGB, as indicated by the Adj. R2 values ranging from 0.9772 to 0.9839. The AIC values varied from −170 to 237, and the RMSE values ranged between 0.0127 and 11.17.

3.3. Selection of Best-Fitted Model

Among the 12 models, model M9 exhibited the highest adjusted R2 value of 0.9839, while models M4, M2, M10, and others followed suit (Table 2). Regarding the RMSE value, model M10 displayed the lowest value of 0.01269, with models M11, M8, M6, and others following suit (Table 2). In terms of AIC values, model M10 had the lowest value of -170, followed by models M11, M6, M8, M12, and others (Table 2). The models were ranked based on the criteria of selecting the lowest value obtained by summing the ranks from the three criteria mentioned above. Consequently, model M10 was deemed the best model for estimating AGB.

The best-fitted model for accurate estimation of AGB among the 12 models was found to be model M10, with an adjusted R2 value of 0.9809, an RMSE value of 0.01269, an AIC value of −170, and a p-value of <8.39e−21 (Tables 3 and 4). Following model M10, models M11, M9, M4, M8, and others showed a good fit for AGB estimation. The R2 value of 0.9809 indicated that the model could explain approximately 98.09% of the variability in AGB. Based on the analysis of R2, RMSE, and AIC values, the best-fitted model for AGB estimation in this study can be represented as follows:where

ABCDEF

−32.07700659−0.002782833−0.0005366756.63504E−095.523833553−0.837778227

The correlation coefficient, multiple R, had a value of 0.99, indicating a strong positive linear relationship between the dependent and independent variables. The coefficient of determination, R2, was 0.98, suggesting that the model fit the data well. The adjusted R2 value of 0.98 indicated that the model adequately accounted for the number of independent variables. The p-value in Table 5 was much smaller than 0.05 (8.39E−21 < 0.05), indicating a significant difference between the means of the groups (Table 4). The linear regression analysis between AGB and vegetation parameters derived from Sentinel-2 data (as shown in Table 3) revealed that most vegetation parameters, including intercept, Red, Green, and ln (NIR), displayed a significant and positive correlation with AGB. Among these parameters, ln (NIR) was identified as the most reliable vegetation index corresponding to AGB.

Upon conducting a graphical examination of the residual plots (Figure 10), model M10 was further analyzed. The analysis revealed that model M10 exhibited a well-distributed pattern. The graphical analysis demonstrated that the majority of residuals were clustered around zero, ranging from 0.03 to −0.03. This observation indicates that the equation was highly suitable and accurately captured the data.

3.4. Validation of Regression Equations

For the purpose of validating model M10, ten sample plots were selected in the same study area using a systematic sampling method. An F-test was employed to compare the variances between the observed AGB and the AGB calculated from the best-fitted model. The assumed hypothesis aimed to test the similarity between the variances of the observed AGB and the calculated AGB, with the null hypothesis (H0) stating that the variance ratio of the observed AGB and the calculated AGB from the model was equal to 1. On the other hand, the alternative hypothesis (H1) posited that the variance ratio of the observed AGB and the calculated AGB from the model was not equal to 1.

F-test to compare two variances:
Data: Observed AGB and calculated AGB
F = 0.47462, num df = 9, denom df = 9, p-value = 0.2822
Alternative hypothesis: true ratio of variances is not equal to 1
95% confidence interval: 0.1178893 1.9108240
The ratio of variances: 0.4746217

The results of the F-test comparing the two variances between the observed AGB and the calculated AGB indicated a p-value of 0.2822, which was greater than the significance level of 0.05 at a 95% confidence level. Therefore, the null hypothesis (H0) was accepted, while the alternative hypothesis (H1) was rejected. This outcome suggested that there was no significant variation between the observed AGB and the calculated AGB, confirming the eligibility of model M10 for further verification using a two-sample t-test. The t-test was conducted at a 95% significance level to assess the significance difference between the observed AGB and the calculated AGB from model M10. The assumed hypothesis for this test was that the null hypothesis (H0) stated there was no significant difference between the observed AGB and the calculated AGB from the model, while the alternative hypothesis (H1) proposed there was a significant difference between the observed AGB and the calculated AGB from the model.

Two sample t-test:
Data: Observed AGB and calculated AGB
t = −0.92699, df = 18, p-value = 0.3662
Alternative hypothesis: true difference in means is not equal to 0
95% confidence interval: −0.09316409 0.03611999
Mean of x = 6.668889
Mean of y = 6.697412

The p-value obtained from the t-test was 0.3662, which was greater than the significance level of 0.05 at a 95% confidence level. As a result, the null hypothesis (H0) was accepted, and the alternative hypothesis (H1) was rejected. This finding indicates that the AGB estimated from field observations and calculated from model M10 were similar, demonstrating the success of model M10 in AGB estimation using the five independent tree parameters. The correlation analysis between the calculated and observed AGB of the ten sample plots revealed a strong coefficient of determination, with an R2 value of 0.87. This suggests that approximately 87% of the observed AGB can be explained by the calculated AGB based on this model (Figure 11).

4. Discussion

This study marks the first attempt to estimate the AGB of teak species in the Sagarnath plantation area using Sentinel-2 imagery data. It represents a significant advancement in the field of forest inventory through the application of remote sensing techniques. The primary objective of this study was to evaluate the usefulness of Sentinel-2 satellite imagery for modeling, forecasting, and mapping AGB. Furthermore, it aimed to explore the potential of utilizing open-source remotely sensed data to support the implementation of the UNFCCC’s REDD+ program in tropical regions. The study employed a destructive sampling method, which is known to be time-consuming and financially demanding when applied to large-scale geographical areas. Similar modeling research studies and biomass estimations have described this method [27, 28]. It is expected that the biomass models derived from these relatively few sampled individuals will sufficiently represent a wide range of sizes, locations, and stand conditions [27, 29].

The findings demonstrate that the oven-dry wood density values ranged from 627 to 714 kg m−3, with an average of 669.5 kg m−3. Comparisons with previous research indicate that teak wood density in a humid tropical site in Costa Rica ranged from 450 to 650 kg m−3 [30], and in East Timor, the mean oven-dry wood density was 607 kg m−3 [31]. Reports from Costa Rica indicate teak wood density ranging from 500 to 650 kg m−3 at 8 years and from 500 to 600 kg m−3 at 4 years [32]. In Bolivia, the average density at 12% constant moisture ranged from 640 to 730 kg m−3 in natural teak forests and plantation teak forests, respectively [33]. In Togo, teak plantations of different ages showed densities of 647, 728, and 779 kg m−3 for ages 11–16, 40−45, and 67–70 years, respectively [34]. During the early stages of plantation development, wood density was found to be more correlated with tree age than with silvicultural management, site, or region [35]. The results indicate that long-rotation teak exhibits higher wood density than short-rotation teak, explaining the slight differences observed in wood density values compared to other literature on teak plantations.

Regarding stem volume, the study revealed an average volume of 85.5 m3 ha−1, with a range of 72.6 to 114.8 m3 ha−1. A study conducted in the dry zone of Sri Lanka reported a volume of 127.8 m3 ha−1 for 29-year-old teak plantations [36]. Thapa and Gautam [37] found volumes of 84.2, 99.9, 125.0, 136.1, and 145 m3 ha−1 for teak stands at ages 7.5, 8.5, 9.5, 10.5, and 11.5 years, respectively, before thinning. The slightly smaller volume estimate in this study can be attributed to differences in tree age and the formula used for volume calculations.

In terms of foliage dry biomass, previous studies reported values ranging from 0.9 kg (225 kg ha−1) to 38.1 kg (952.5 kg ha−1) for teak plantations in Costa Rica [32]. In this study, the foliage dry biomass ranged from 1,381 to 2,360 kg ha−1 at 36 years of age, indicating slightly higher values compared to the earlier report. Biomass estimations in teak plantations in the Dehradun Forest Division, Uttarakhand, and a dry tropical region in India were reported as 69.71 [38] and 77 t ha−1 [39] at 30 years of age, respectively. Rinnamang et al. [40] found total AGB values of 84.

5. Conclusion and Recommendation

This study demonstrated a methodology for estimating AGB in teak forests by integrating Sentinel-2 imagery data and field data. The results revealed that the average AGB in the teak forest area was 814 kg per plot (81.4 t ha−1), with a range of 716 kg per plot (71.6 t ha−1) to 1,060 kg per plot (106 t ha−1). Among the 12 models tested, model M10, which incorporated independent variables such as the Red band, Green band, Blue band, NIR, NDVI, and natural logarithm, yielded the most significant model for AGB estimation in the 36-year-old teak plantation forest of Nepal. Model M10 outperformed compare to other models based on three criteria (Adj. R2 = 0.9809, RMSE = 0.01269, AIC = −170), ranking first. The equation of the best-fitted model (M10) was ln (AGB) = A + B × Red + C × Green + D × Blue2 + E × ln (NIR) + F × ln (NDVI) for accurate AGB estimation, with coefficient values of A = −32.07700659, B = −0.002782833, C = −0.000536675, D = 6.63504E-09, E = 5.523833553, and F = −0.837778227. The parameters used in model M10 are easily obtainable, making it user-friendly for AGB estimation. This model offers a time and cost-effective approach for mapping and monitoring teak forest biomass. The estimated AGB also contributes to national reference scenarios for teak species biomass, carbon stock estimation, sequestration, and the implementation of the REDD+ program. The model is well-suited for use in the lowlands of Nepal and similar eco-regions where teak species thrive. Various opportunities, such as fulfilling demands for veneer, saw logs, firewood, poles, and small timber, exist for promoting teak plantations in the terai and inner terai of Nepal.

Based on the study’s findings, the following recommendations are proposed: to assess the applicability of this model in teak forests with similar stand structures in different physiographical zones, to conduct further research investigating the correlations between field inventory data and Sentinel-2 imagery data in teak plantation forests of various age groups to enhance the utility of the model, and to explore additional research avenues by modifying or expanding the independent variables (bands and vegetation indices) in the model to improve the precision of AGB estimation in forested areas.

Data Availability

Data will be available upon request from the journal.

Conflicts of Interest

The authors declare that they have no conflicts of interest