Synthesising Existing Forest Inventory Datasets for Estimating Historical 1 Aboveground Biomass Stocks, Growth and Mortality in Logged-over Tropical 2 Dipterocarp Forests of Kalimantan, Indonesia

Background : Extensive forest inventory data is available from commercial timber companies. For 14 this study, over 20,000 plots were compiled for North, East and West Kalimantan provinces, with 15 more than 17,000 of these exceeding our quality assurance tests. This study aimed to: (1) explore 16 the potential use of existing permanent sample plots and forest inventory data established and 17 measured by timber concessions; (2) assess uncertainties of aboveground biomass (AGB) estimates 18 using various allometric models; (3) analyse the dynamics of AGB in logged-over dipterocarp 19 forests; (4) analyse AGB stocks and emission factors in tropical dipterocarp ecosystems. 20 Methods : Two types of forest monitoring datasets measured by timber companies in Indonesia were 21 compiled and assessed in this study: permanent sample plots (PSPs) (24 1-ha plots), and the overall 22 periodic timber inventory (OPTI) (17,301 plots). We compared various allometric equations for 23 estimating AGB of the plots and developed a simple AGB equation using basal area (BA) as 24 predictor. We further evaluated the AGB growth and mortality using the PSP plots.


Results:
We found that the model using only tree diameter (D) as a predictor variable tended to be 26 unbiased when aggregating the estimates at larger plots. We also found that BA per hectare could 27 explain the variation of AGB at plot level (adjusted r 2 = 0.911; root mean square error [RMSE]: 28 27.8). We overlaid the OPTI plot with the land cover map and estimated the mean AGB of the 29 associated land cover classes. The mean AGB of primary dryland forest, secondary dryland forest 30 and bush classes were 281.1 + 4.0 Mg/ha, 231.5 + 1.7 Mg/ha and 179.0 + 5.0 Mg/ha, respectively. 31 Nine years after logging, the mean AGB is still lower than the mean AGB two years after logging. 32 The growth rate (2.5%) was still lower than the mortality rate (3.1%), and recruitment (0.2%) did 33 not occur until seven years after logging. 34

Conclusions:
The results of this study suggest that the existing forest monitoring data should be 35 incorporated into the carbon accounting system at district, province and national level to improve 36 the estimation of forest biomass and emission factors related to forest degradation and deforestation. 37 However, there is a need for data quality assessment prior the analysis and a standardised platform 38 for nation-wide forest inventory database is therefore required. 39

40
Tropical lowland dipterocarp forests (LDF) play a crucial role in economic development, ecological 41 service and carbon balance. Over the last several decades, the tropical LDFs have reduced in size 42 and quality because of rapid deforestation and forest degradation (Stibig et al., 2014). The causes 43 are mostly economic, such as excessive timber extraction and conversion to agricultural land and 44 plantation (Gibbs et al., 2010;Geist and Lambin, 2001). To reduce deforestation and forest 45 degradation in tropical regions, the United Nations Framework Convention on Climate Change 46 (UNFCCC) developed an incentive mechanism under the climate change mitigation framework to 47 reduce emissions from land use and forestry (UNFCCC, 2010). Accurate estimates of carbon 48 dioxide (CO 2 ) emissions from deforestation and forest degradation, as well as sequestration from 49 forest conservation, sustainable forest management and forest carbon stock enhancement, are 50 required to assess the impact of Reducing Emissions from Deforestation and Forest Degradation 51 (REDD+) activities and contribution to global climate change mitigation efforts (Herold and 52 Skutsch, 2009). 53 A major source of emissions in tropical forests is not only deforestation but forest degradation 54 (Huang and Asner, 2010). Forest degradation is associated with a reduction in quality of forest 55 structures or biophysics, such as timber stock, biomass or biodiversity (Lund, 2009). The main 56 causes of forest degradation are selective cutting and small-scale illegal logging (Hosonuma et al., 57 2012). Compared to tropical deforestation, the impact of forest degradation in a unit area is smaller 58 and more difficult to monitor (Lambin, 1999), but occur over substantial number of hectares 59 (Gaveau et al., 2014) and contribute to a larger greenhouse gas emission than from deforestation 60 (Pearson et al., 2017). A better understanding of forest dynamics in selectively-logged forests is 61 crucial for policy development related to sustainable forest management and climate change 62 mitigation. 63 Monitoring of REDD+ related activities is still problematic and uncertain because of limited 64 research on the tropical region (Herold et al., 2011). Although the monitoring of deforestation and 65 forest carbon stock are more advanced because of the availability of relevant satellite imageries, the 66 scope of existing studies is mostly global, with medium to low-resolution satellite imageries used 67 (Baccini et al., 2012;Saatchi et al., 2011). Therefore, to quantify the carbon impact of forest 68 degradation and carbon stock enhancement activities on a regional scale using global scale studies 69 is common (Penman et al., 2016). 70 IPCC suggested two methods for estimating emissions from forests: stock-difference and gain-loss 71 methods. The first method requires aboveground biomass (AGB) stocks at two different 72 measurements to generate emission factors. The later method requires AGB growth and loss 73 measurement to estimate annual net increment. In the forest reference emissions level document, 74 Indonesia acknowledges the need to utilise the existing data on forest inventory and forest 75 permanent sample plots (PSPs) established and measured by timber companies (MoEF, 2015). 76 Almost half of forest estates in Indonesia were designated as production forests for timber 77 where OPTI and PSP datasets were compiled 132

Data Collection 133
The data collection was conducted between February and May 2012 under the auspices of the 134 FORCLIME-GIZ, a German-Indonesian cooperation in forestry and climate change. The datasets 135 used in this study were compiled from timber concessions operating in project pilot provinces (West 136 Kalimantan, East Kalimantan and North Kalimantan). Our interest was in collecting raw and digital 137 data on PSP and OPTI datasets from the existing, active timber concessions operating in the study 138 area. 139 First, we contacted all responsible government institutions, such as Forestry Services at province 140 and district levels. Officially, all timber concessions must report all forest inventory and permanent 141 sample data as a requirement for getting approval for their cutting plan proposal. Second, we sent 142 official letters to the existing timber concessions to request related datasets. As some concessions 143 are no longer operational, only active concessions were targeted for data compilation. 144

PSP Datasets 145
In each 1-ha plot of 100 m × 100 m, trees with diameter at breast height (DBH) > 10 cm were 146 measured, their scientific names, marked and labelled for long-term monitoring. The plot 147 boundaries were clearly marked with poles in each corner. The circumference at breast height 148 (CBH) or 20 cm above buttresses and total tree height (H) were measured for each tree in the plot. 149 The DBH is CBH divided by phi value (3.14). Within the 1-ha plots, the trees were recorded in each 150 subplot of 10 m × 10 m, allowing us to regroup the trees into smaller plots of 20 m × 20 m; 30 m × 151 30 m; 40 m × 40 m or 50 m × 50 m (Fig. 2). The smaller plots were used to evaluate the accuracy 152 of AGB estimation at various plot size. However, we ignored the potential spatial autocorrelation 153 between subplots, which may inflate the degree of freedom and slightly bias the estimates. We compiled digital files of PSP data from five timber concessions in East and North Kalimantan 158 (Table 1). For this study we selected only datasets that had complete DBH and scientific name 159 records. Only datasets from two companies out of five met the requirements totalling 24 hectare 160 plots, which were measured independently. For the first dataset, a total of 18 1-ha plots were 161

OPTI Datasets 169
For the purpose of allowable annual cut planning, each timber concession in Indonesia is required 170 to conduct OPTI for the whole concession area (MoF, 2007). The OPTI collect only the DBH and 171 commercial tree names of each recorded tree, in a nested plot of circular subplot, 10 m × 10 m 172 subplot, 20 m × 20 m subplot and 20 m × 125 m for saplings, poles, trees d < 35 and large trees d > 173 35 cm, respectively. Additionally, the physical appearance and quality of trees were also recorded. 174 The plots were systematically distributed with a distance of about 900-1000 metres, depending on 175 the size of the concession. 176 We successfully compiled the OPTI datasets measured between 2009 and 2011 from 33 timber 177 concessions in East, North and West Kalimantan provinces. The dataset consisted of 20,133 plots 178 (Table 2). Duplicate plots from the same timber concessions and plots from a non-native timber 179 plantation were excluded (3.7%) (Filter 1). Because each dataset was measured independently by 180 the timber concession, data checking for quality assurance was essential. We removed plots with 181 data inconsistency, such as large trees recorded in small subplots, redundant tree numbers within

201
Due to the occlusion of tree canopies, measuring tree height accurately in tropical forests is difficult. 202 Figure SM1 showed a low precision of tree height measurement in the plots. In the case of canopy 203 occlusion, the tree height was often estimated using a local D-H model. The overlay between the 204 tree height data with the regional D-H model for lowland forests in Indonesia (Manuri et al. in press) 205 showed strong agreement. We performed outlier analysis based on the studentised residuals of the 206 regression between the measured and modelled tree heights. Residual values larger than 2 and 207 smaller than -2 were outliers (Sileshi, 2014). We estimated the tree height of the outliers using the 208

Development of AGB-BA Models 218
Assuming the AGB estimates using complete variables are the most accurate, we developed an 219 AGB model based on BAs, which was found to be simple but with relatively high precision 220 (Burrows et al., 2002;Slik et al., 2010). We assessed the performance of the AGB-BA models at 221 individual tree and plot levels (0.04 ha, 0.09 ha, 0.16 ha, 0.25 ha and 1 ha plots) using the first 222 measurement of 24-ha PSP data sets. 223

Assessing AGB Dynamics in Logged-Over Dipterocarp Forests 224
We used the selected PSP dataset, which has complete measurement of DBH, tree height, species 225 identification and long measurement (more than five years). Net annual AGB increment (IAGB) was 226 calculated as the net annual AGB change due to growth (GAGB), recruitment (RAGB), mortality 227 (MAGB) and shrinkage (SAGB) of trees that annually averaged across the monitoring period, assessed 228 using the following equations Alder (1995)

Estimating AGB Stocks of Logged-over and Primary Dipterocarp Forests 237
For logged-over forests we used an OPTI and PSP dataset compiled from existing timber 238 concessions in the study area. We used the D1 and DGH AGB model (Table 3)

Testing Allometric and Biometric Models for Estimating AGB of Logged-over Forests 246
We analysed the regression between the best-predicted AGB (using the DGH model) and the 247 predicted AGB values (using a model with fewer predictor variables) ( Figure 3). As expected, the 248 DG model performed similarly with the DGH model at tree level. The DG model was still relatively 249 unbiased in small plots, but tended to be larger in larger plots. A similar trend occurred with D2 250 models, which used diameter and wood density class as predictor variables. The D1 model, which 251 used only diameter as a predictor variable, tended to be unbiased, aggregating the estimates at larger 252 plots (e.g., 2500 m 2 ). 253 254 Figure 3. Regression between best-predicted total AGB and the predicted total AGB using less 255 predictor models at tree individuals and plot levels.  1 ha.). We found that the non-linear models were better than the linear models, in terms of the 262 normality of residuals distribution (result not shown). Models with additional WD as a predictor 263 variable were only slightly better than the AGB model using only BA as a predictor variable (Figure  264 4). We decided to use the model with BA alone as a predictor variable, as it is also not practical to 265 estimate average WD at plot level. We found that the best AGB-BA model is the model for 0.25-ha 266 plot (adjusted r 2 =0.911; RMSE: 27.8). The 0.25-ha plot is coincidently the same size as the OPTI 267 plot. 268 16 269 Figure 4. Regression between the best-predicted AGB (best AGBp) with predicted AGB using non-270 linear BA models(AGB-BA) 271

Estimating AGB Stock of Logged-over Forests 272
We computed AGB of OPTI dataset using the D1 equation, which was tested to be less biased when 273 applied to the 0.25 and 1 ha PSP plots (Figure 3). The result was compared to the AGB calculated 274 in previous studies in primary forests using pan-tropical equations, and our study in PSP using the 275 DGH equation. The AGB estimates of OPTI dataset using the D1 equation were still in agreement 276 with the AGB estimates in logged-over forests and primary forests, where the plots have a BA less 277 than 40 m 2 /ha ( Figure 5). However, the estimates using D1 equation tended to be lower than that of 278 previous studies in primary forests, especially where the plots have a BA of more than 40 m 2 /ha. In 279 contrast, our predicted AGB using a AGB-BA model (AGB = 6.37×BA 1.206 ) developed from PSP 280 dataset was in better agreement with the estimates from logged-over forests with BA less than 25 281 m 2 /ha and primary forests that have BA more than 40 m 2 /ha. The AGB estimates of AGB-BA model 282 were lower than the studies in primary forests, especially in the plots that have BA between 25 -40 283 m2. 284 285 Figure 5. Regression  We compared our BA estimates using our datasets (OPTI and PSP) with BA estimates from 291 previous studies in lowland tropical dipterocarp forests. Our estimates were lower than the estimates 292 of previous studies in primary forests, but in accordance with the estimates from logged-over forests 293 ( Figure 6). The mean BA of OPTI and PSP dataset is comparable with the BA estimates of medium 294 and high-impact logged-over forests, respectively (Sist and Nguyen-Thé, 2002). The results suggest 295 that the OPTI dataset is generally consistent and reliable. 296 297

Figure 6. BA distribution (mean and standard deviation) from PSP (logged-over forests), OPTI 298 (mixed forest cover) and literature (primary and logged-over forests) 299
We used forest cover type information recorded during the field measurement for estimating mean 300 AGB. Primary dense forest had the highest mean AGB (371.2 Mg/ha) and non-forest had the lowest 301 mean of AGB (148.2 Mg/ha). The mean AGB of other forest types ranged from 206.9 to 289.7 302 Mg/ha. Unexpectedly, the bush had a mean AGB higher than the secondary dense forests (Table 4). 303 About 45% of the plots did not have information on forest cover type. 304 Based on the Tukey test on least square mean differences, only mean AGB from primary dense 305 forest and non-forest were significantly different to other forest cover types. Other forest cover 306 types, ranging from primary medium forest to bush, were not significantly different from each other. 307 We further overlaid the plots with geographical references with land cover map 2009 derived from 311 satellite imagery classification (MoF, 2012). Only three land cover classes were represented by 312 more than 100 plots, meaning they had high confidence level of the mean AGB estimates (i.e., 313 primary dryland forest, secondary dryland forest and bush [ The estimate of each land cover class was lower than the estimate based on similar forest cover 322 classes derived from the field. This is because the land cover classification used a visual delineation 323 method, which classifies several pixels close to each other as one entity, based on the majority of 324 pixels. If the resolution is low, small portions of different pixels will likely be classified as different 325 classes of majority pixels. Therefore, the mean AGB of forest classes will be lower than the forest 326 classes based on plot identification, potentially due to the inclusion of small patches of low-density 327 forest or logging-effected areas. 328

AGB Dynamics in Selectively-logged Dipterocarp Forests 329
Six 1-ha PSP datasets were used to analyse AGB dynamics after logging. The datasets have the 330 longest measurement period: six time measurements from two to nine years after logging. At first 331 measurement (two years after logging) the mean total AGB was 258.3 Mg/ha. At the second 332 measurement (three years after logging), this was reduced to 240.72 Mg/ha due to the mortality of 333 large and medium trees (Figure 7). The mean AGB continued to decline towards 229.45 Mg/ha at 334 the fourth measurement (five years after logging) before it finally increased at the fifth (seven years 335 after logging) and the sixth (nine years after logging) measurements (Table 6). 336 annual mortality was high in the earlier measurement period, and decreased towards the end of 343 measurement. In contrast, the annual growth in the earlier measurement was very low, then 344 continued to increase to 10.9 Mg/ha/year at seven years after logging. 345 At nine years after logging, the mean AGB was still lower than the first measurement two years 346 after logging, because the growth rate was still lower than the mortality rate. Also, recruitment did 347 not occur until seven years after logging, which contributed only 0.49 Mg/ha annually. Surprisingly, 348 the shrinkage was higher than the recruitment (0.74 Mg/ha/year). 349

350
There have been few attempts to evaluate the potential use of existing forest inventory for assessing 351 historical forest carbon stocks and biomass growth in Indonesia. This study explores the potential 352 of existing forest inventory plots for quantifying carbon stocks and growths in logged-over forests. 353 We successfully compiled data from 20,133 OPTI plots. After data filtering for consistency and 354 outlier checks, we found that about 17,301 plots (85.9%) were reliable. Most of the removed plots 355 had unrealistically large total BAs or AGB, because of frequent recording errors or unexpected 356 existence of very large trees in the plots. A validation process must be conducted to address this 357 issue. The current process for validating the OPTI result is performed only if there is a discrepancy 358 between the plan and the implementation, mostly due to administrative matters (MoF, 2009). The 359 selection of plots for field validation should also be based on the outlier plots. 360 Our analysis of the distribution of AGB based on the forest cover information from the field found 361 that they were not consistent, so the mean AGB values among forest cover types were mostly not 362 significantly different. It seems the definition of each forest cover type overlapped and were hard 363 to distinguish in the field, thus confusing field crews when defining forest cover. We also found that 364 57.3 % of the OPTI plots did not record geographical position. Because of Indonesia's large 365 geographical size, it is suggested that a geographic coordinate system be used for easy compilation 366 and comparison among OPTI databases. Also, it seems that the coordinate positions of the plots do 367 not represent the actual plot position, as can be observed from the fully systematic distribution of 368 the plots. In relatively difficult terrain with limited accessibility, reaching a plot or placing a plot as 369 the plan is often problematic. Therefore, 100% similarity between the planned and actual plot 370 position seems to be unrealistic. The actual location of the plot using GPS in the field is not only 371 useful for documentation and revisiting for field validation, but also for validating estimation of 372 forest metrics based on remote sensing imageries. A relationship model between AGB and BA could be used for estimating AGB stock from historical 380 forest inventory plot summaries. Our finding suggests that the AGB-BA non-linear model is better 381 than the linear model for estimating a wide range of BA values representing areas with scarce trees 382 to dense forests, with very large trees in tropical dipterocarp forests. This is different from a study 383 on eucalyptus forests, where a linear model of AGB-BA achieved similar accuracy (Burrows et al., 384 2002). The main reason is that tree diameter and height ranges in tropical dipterocarp forests are 385 much higher than in eucalyptus forests. 386 Similar to OPTI, the PSP datasets were utilised in a limited manner within the timber concession 387 for yield and annual allowable cut regulation. Twelve out of the compiled 19 PSP series were not 388 in digital format. Only four series were in digital format, with botanical identification at least at 389 genus level. The mortality rate is still higher than the accumulation of growth and recruitment rate 390 nine years after logging. A study in dipterocarp forests in East Kalimantan found that the highest 391 mortality rate occurred one to three years after logging (Susanty et al., 2015), which is in agreement 392 with our study. Mortalities in selectively-logged forests occurred even after eight to 17 years after 393 logging (Cannon et al., 1994), due to damage from logging (Nguyen-The et al., 1998) and wind 394 disturbance after fragmentation (Laurance et al., 1998). This implies at the need for managing long-395 term forest plot monitoring database in logged over forests, which established and measured by 396 timber concessions in Indonesia since 1995 (Tata et al., 2010). Unfortunately, since 2009 397 development and measurement of PSP is no longer a requirement for timber concessions when 398 applying for cutting permits in Indonesia. 399

400
Using existing OPTI datasets, we are able to estimate mean AGB stock with high confidence, as 401 well as using them for estimating AGB based on land cover map. In this study we developed 402 models using BA per plot as a predictor variable for estimating AGB with high precision and low 403 bias. Our estimates of AGB in primary forests were calculated based on forest cover information 404 from field plots, and land cover maps derived from satellite imagery classification, which were 405 useful for estimating emission factors from the land and forest cover types. 406