Journal 
Introduction
Bhutan has constitutionally declared to maintain at least 60% of territory under forest cover. Currently, Bhutanese forests cover more than 70.5% of the land surface and are net carbon (C) sinks (NSSC 2011). Bhutanese government has declared to remain C neutral ultimately and has also recently participated in the joint declaration for mitigating global climate change with EU. For this reason, estimating forest C dynamics in Bhutan can contribute to understanding global forest C sink (IPCC 2014a).
In order to estimate the forest C dynamics in Bhutan, accumulation of empirical data is essential. However, Bhutan currently lacks adequate forest C data. Under this circumstance, there are two options for quantifying forest C dynamics: using default values and modeling (IPCC 2014b). Using default values is a simple option, but it seems highly uncertain. Modeling can estimate forest C dynamics at diverse scales while most forest C models require much amount of data (Lee et al. 2010).
The Forest Biomass and Dead organic matter Carbon (FBDC) model, which is a developed forest C model in Korea, is able to simulate forest C dynamics with small amount of required data (Lee et al. 2014). The FBDC model estimated the C dynamics in Korean forests at diverse spatial and temporal scales (Park et al. 2013a, 2013b; Yi et al. 2013; Lee et al. 2014, 2015, 2016). In spite of the lack of empirical data, it is possible to estimate forest C dynamics in Bhutan with the FBDC model. In this study, we attempted to provide a possibility of using the FBDC model to Bhutan and some empirical knowledge for annual forest C dynamics in Bhutan. The FBDC model was newly adjusted for Bhutanese forests and changes in C stocks of biomass, litter, dead wood, and mineral soil were separately provided.
Materials and methods
1.Study site
This study was conducted in Bumthang and Paro regions of Bhutan. According to World Bank Group, the mean air temperature and annual precipitation were 11.93°C and 1861 mm, respectively, during 1990-2012 (http://sdwebx.worldbank.org). Tsuga dumosa and Juniperus spp. were selected for estimating forest C dynamics in Bhutan because they are major tree species of temperate zone in Bhutan. Tree data (species, diameter at breast height (DBH), tree height, and tree ring) were collected from six circular plots with 12.62 m radius (0.05 ha). The most dominant tree species were selected from each plot, on the basis of the basal area. Tsuga dumosa and Juniperus spp. dominated 2 and 4 plots, respectively. The age of individual tree was determined by number of tree ring and it ranged from 25 to 121 years. The age of plot was calculated by averaging age of each individual tree in the plot. The DBHs and tree heights were also measured and they were 11-77 cm and 0.5-15.4 m, respectively.
2.The FBDC model
The FBDC model can simulate forest C dynamics: growth of biomass, input of litter and dead woody debris, and decomposition of dead organic matter (Yi et al. 2013; Lee et al. 2014). The model comprises 5 biomass pools and 8 dead organic matter pools (Lee et al. 2014). Growth of biomass is formulated with yield tables and allometric functions, and C dynamics in dead organic matter are simulated with model parameters (Lee et al. 2014). Detailed formulations and explanations are shown in Lee et al. (2014) and Yi et al. (2013).
3.Parameterization of the FBDC model and simulation
Parameterization of the FBDC model is essential to simulate forest C dynamics in Bhutan. For the parameterization, growth functions and model parameters should be substituted. Due to the lack of allometry data in Bhutan, most of the growth functions and model parameters in the original FBDC model were used (Table 1, Table 2, Table 3). The growth functions and parameters for Pinus densiflora and Chamaecyparis obtusa in Korea were used as those for Tsuga dumosa and Juniperus spp. in Bhutan, respectively, because they are included in the same genus and family. The growth functions of Juniperus spp. were newly developed because the FBDC model does not include those functions of Chamaecyparis obtuse (Table 1). In order to calibrate heterogeneity in productivity among the plots, the growth modifier was used (Lee et al. 2014). The growth functions, stocking volume, and mean forest age at each plot were used to estimate the growth modifiers.
After the parameterization for Bhutanese forests, the forest C dynamics during 80 years were simulated. In order to initialize C stocks in each pool at forest age of 0, the spin-up process, which is a common method for initialization, was conducted with 80-year-interval harvest as Lee et al. (2014) did. Then, the actual forest C dynamics were simulated for 80 years.
Results and discussion
1.Forest C dynamics
The total C stocks (Mg C ha-1) ranged from 118.35 to 200.04 with an average of 168.41 (Fig. 1). Especially, the biomass and mineral soil were dominant components of total forest C stocks, explaining 36.0 and 56.0% of total C stocks, on average. The average biomass and mineral soil C stocks (Mg C ha-1) were 60.80 and 94.25. The biomass showed high variation (3.40-88.13 Mg C ha-1) with respect to stand age. In contrast, the C stocks (Mg C ha-1) in mineral soil showed low variation, ranging from 91.45 to 97.90. The C stocks (Mg C ha-1) in litter and dead wood varied from 4.24 to 24.95 and from 1.99 to 20.31, which were generally less than those in biomass and mineral soil. The average C stocks (Mg C ha-1) in litter and dead wood were 9.19 and 4.17, accounting for 5.5 and 2.5% of the average total C stocks, respectively.
In order to verify the model estimates, the estimated values were compared to values from previous studies. The estimated total C stock (168.41 Mg C ha-1) was within the previous reports of total forest C stock for temperate forests (150-217 Mg C ha-1) (Dixon et al. 1994; IPCC 2001). However, according to Global Forest Resources Assessment for Bhutan (FAO 2011), the average C stocks (Mg C ha-1) in biomass, litter, and mineral soil were 103.42, 18.47, and 69.56, respectively. Compared to those values, the model seemed to underestimate the C stocks in biomass and litter. In contrast, the C stock in mineral soil was overestimated. Forests in Bhutan actually grow under various climatic and environmental conditions (Ohsawa 1987). The uncertainty might result from the fact that this study only dealt with the C dynamics of a single temperate forest type. Meanwhile, the measured C stocks in Bhutan also showed high variation.
At the study plots, the measured biomass C stocks (Mg C ha-1) ranged from 23.91 to 101.87. Also, mineral soil C stocks in Bhutan showed large variation (56.50-196.76 Mg C ha-1; Tshering, unpublished). It was found that C stocks were highly dependent on forest age, as described in previous studies (Zhou et al. 2009; Lee et al. 2016). Thus, systemic forest C inventory system is urgently required to verify the model estimates.
The simulation exhibited that the C balance was closely related to forest age (Fig. 2). During the early development stage (1-20 years), the forests were C sources (0.82 Mg C ha-1 yr-1). During the forest age of 20-80 years, the forests sequestered the atmospheric C at a rate of 1.12 Mg C ha-1 yr-1. The annual C sequestration rate converged to zero as the forests highly matured. Our estimate was less than an annual C sequestration in Nepal (1.88 Mg C ha-1 yr-1; Banskota et al. 2007). It might be attributed that tree sequesters more C during the intermediate growing stage. The annual C sequestration was around 1.91 Mg C ha-1 yr-1 during the forest age of 20-50 years in our study site (Fig. 2).
The annual changes in C stocks highly varied among the C pools (Fig. 1). The biomass C pool especially dominated the annual changes in C stocks in this study and it is consistent to previous studies (Lee et al. 2014, 2016). The biomass C stock generally increased with increasing forest age at the rate of 1.07 Mg C ha-1 yr-1. However, the increment decreased with increasing forest age as forests mature. The C stocks in litter and dead wood were decreasing at a rate of 0.95 and 1.07 Mg C ha-1 yr-1 until the forest age of 20 years with increasing forest age at the early development stage while these C stocks increased at a rate of 0.10 and 0.02 Mg C ha-1 yr-1 with increasing forest age after that stage. It might be attributed to the changes in C balance between decomposition and accumulation of dead organic matter (Liski et al. 2005; Kurz et al. 2009; Lee et al. 2014). In contrast, the mineral soil C stock seemed relatively constant with forest age.
2.Implications
We piloted the FBDC model for estimating forest C dynamics in Bhutan with limited data. The simulation results showed empirical knowledge on forest C dynamics in Bhutan. It showed a possibility that annual Bhutanese forest C inventory can be reported with modeling (Tier 3), which is one of the most reliable reporting system in IPCC Good Practice Guidance (IPCC 2014b).
In order to report Bhutanese forest C inventory at a national scale, field measurements should be conducted in parallel for major tree species and forest types. Even though this study estimated annual changes in C stocks of two major tree species at a plot scale, the estimates have limited applicability because Bhutanese forests consist of various tree species and forest types. Especially, those compositions of forests change with altitudinal gradient, from tropical zone to alpine zone (Ohsawa 1987; Wangda and Ohsawa 2006). Thus, basic forest data for major tree species and forest types are required to improve the FBDC model estimates for establishing Bhutanese forest C inventory. Furthermore, the model-estimated C stocks can be verified well with measured C stocks.
Conclusion
We provided the C stocks and their changes in biomass, litter, dead wood, and mineral soil in Bhutanese forest although available data was very limited. The forest C dynamics highly varied by forest age and C pools. Especially, the biomass C stocks increased with forest age while the mineral soil C stocks seemed insensitive to forest age. Our results showed a possibility for reporting annual forest C inventory in Bhutan with reliable modeling method. Field work for major tree species and forest types are required to estimate Bhutanese forest C dynamics at a national scale.
