Recent and projected impacts of land use and land cover changes on carbon stocks and biodiversity in East Kalimantan, Indonesia

: Large-scale land use and land cover (LULC) changes can have strong impacts on natural ecosystems, such as losses of biodiversity and carbon. Future impacts, under one or multiple future scenarios, can be estimated with the use of LULC projections from land use change models. Our aim is to quantify LULC change impacts on carbon stocks and biodiversity in the West Kutai and Mahakam Ulu districts in East Kalimantan, Indonesia. Hereto, we used LULC data from 1990 to 2009 and land use change model projections up to 2030 under four contrasting LULC change scenarios differing along two axes: land development (limited vs. unlimited) and zoning (restricted vs. unrestricted), explicitly considering the uncertainties in the land use change model. For the LULC change impact calculations, three quantitative indicators were evaluated: aboveground biomass (AGB) (for carbon stocks), closed-canopy forest patch size distribution and plant species richness (for biodiversity). Subsequently, we statistically assessed whether the motivation to opt for a specific scenario was conclusive given the uncertainty in the indicator values. We found that under the limited development scenarios the projected AGB decrease towards 2030 was insignificant, plant species richness was projected to decrease significantly by ฀3%, and closed-canopy forest patches mainly of 100–1000 ha were projected to become fragmented. The effect of zoning was insignificant under these scenarios. The difference between the limited and unlimited development scenarios was significant, with the projected impacts under the unlimited development scenarios being much higher: AGB was projected to decrease 4–30%, plant species richness 10–40%, and the closed-canopy forest was projected to completely loose its typical patch size distribution. The effect of zoning on these scenarios was positive and significant. These results suggest that the most sustainable pathway for East Kalimantan, given our indicators, would be Large-scale land use and land cover (LULC) changes can have strong impacts on natural ecosystems, such as losses of biodiversity and carbon. Future impacts, under one or multiple future scenarios, can be estimated with the use of LULC projections from land use change models. Our aim is to quantify LULC change impacts on carbon stocks and biodiversity in the West Kutai and Mahakam Ulu districts in East Kalimantan, Indonesia. Hereto, we used LULC data from 1990 to 2009 and land use change model projections up to 2030 under four contrasting LULC change scenarios di ﬀ ering along two axes: land development (limited vs. unlimited) and zoning (restricted vs. unrestricted), explicitly considering the uncertainties in the land use change model. For the LULC change impact calculations, three quantitative indicators were evaluated: aboveground biomass (AGB) (for carbon stocks), closed-canopy forest patch size distribution and plant species richness (for biodiversity). Subsequently, we statistically assessed whether the motivation to opt for a speci ﬁ c scenario was conclusive given the uncertainty in the indicator values. We found that under the limited development scenarios the projected AGB decrease towards 2030 was insigni ﬁ cant, plant species richness was projected to decrease signi ﬁ cantly by ∼ 3%, and closed-canopy forest patches mainly of 100 – 1000ha were projected to become fragmented. The e ﬀ ect of zoning was insigni ﬁ cant under these scenarios. The di ﬀ erence between the limited and unlimited development scenarios was signi ﬁ cant, with the projected impacts under the unlimited development scenarios being much higher: AGB was projected to decrease 4 – 30%, plant species richness 10 – 40%, and the closed-canopy forest was projected to completely loose its typical patch size distribution. The e ﬀ ect of zoning on these scenarios was positive and signi ﬁ cant. These results suggest that the most sustainable pathway for East Kalimantan, given our indicators, would be to limit land development, mainly large-scale cash-crop cultivation. If land development cannot be limited, the implementation of restricted development zones is advised. The methodologic novelty of our approach is that we propagate uncertainties from a land use change model to the impact assessment and test the signi ﬁ cance of di ﬀ erences between future scenarios, in other words we test if a potential policy instrument has a signi ﬁ cant (positive) e ﬀ ect on the studied indicators and may thus be worth implementing.


Introduction
Land use and land cover (LULC) change throughout the world is driven by multiple natural and anthropogenic factors. The demand for agricultural land for the production of food, feed, fibre and fuel is the most prevalent LULC driver (Lambin and Meyfroidt, 2011). Although globally the agricultural area has stabilized around 37% of the total land area in the past decade, it is still growing with about 1 or 2 percent per year in many Latin-American, Asian and African countries (The World Bank Group, 2018). This growth is expected to continue in the coming decade (OECD/FAO, 2015AO 2015. Expansion of the total agricultural area or a change between agricultural systems (e.g. a change in management level or crop type) impacts natural ecosystems and agroecosystems. Depending on, for example, the type of transition, https://doi.org/10.1016/j.ecolind.2019.04.053 Received 20 July 2018; Received in revised form 7 February 2019; Accepted 15 April 2019 its point in space and in time, and the viewpoint and interest of the beholder, such impacts may be positive, negative or neutral.
Whereas past impacts can be measured, if sufficient data before and after the change are available, future impacts cannot be measured by definition. They can, however, be estimated with the use of LULC projections from land use change models. This is especially useful when evaluating the effect of a set of potential future scenarios with different policy instruments seeking to mitigate or minimize negative impacts (Verstegen et al., 2017). A number of studies has quantified the impacts of multiple scenarios of projected agricultural expansion on carbon stocks (van der Hilst et al., , 2018Lawler et al., 2014), biodiversity (Lawler et al., 2014;Visconti et al., 2016), and water (Lin et al., 2007;Rajib et al., 2016) for different case study areas and time frames.
It is recognized that the LULC projections generated by land use change models contain uncertainties stemming from multiple sources, such as errors in input data, over-simplified or missing model transition rules (model structure uncertainty), and inaccurately estimated parameter values within these rules. In the past few years, efforts have been made to try to quantify these uncertainties by creating probabilistic/ stochastic land use change models and performing error propagation (e.g. Verstegen et al., 2016a, Krüger andLakes, 2016) or by doing multi-model ensemble runs (e.g. Alexander et al., 2016).
Obviously, the uncertainties in the LULC projections cause uncertainties in the derived impacts. We deem it important that uncertainties in the land use change models are propagated to the indicators, because this results in a better description of the magnitude and robustness of impact projection and allows for tests on the significance of differences between the impacts of different future scenarios; important for reducing risks and enhancing resilience (Dale et al., 2018). Although the landscape ecology community is strong in spatial statistics, and uncertainty quantification in past impacts is common practise (e.g. Pfeifer et al., 2017), uncertainty quantification in future impacts is not. The opportunity is there: previous work in which future impacts are quantified actually use stochastic land use change models (e.g. van der Hilst et al., , 2018Lawler et al., 2014), but the uncertainties in the land use projections are not propagated to the impacts. This leads to unrealistic impact estimates. Furthermore, as significance tests need probability distributions to give meaningful results, deterministic impact analyses leave the policy maker with the question whether impacts are significantly different under different policy instruments, i.e. to what extent the instruments are effective.
Our aim is to quantify LULC change impacts in West Kutai and Mahakam Ulu districts in Kalimantan, Indonesia under four contrasting future land use change scenarios, differing in the amount of agricultural land development (limited vs. unlimited) and zoning (restricted vs. unrestricted). The relevance of studying LULC change impacts in this study area is discussed in the next section. In line with what is predominantly studied in impact analyses (see e.g. Danielsen et al., 2009;Thomas et al., 2013) and with the criteria advised by Dale and Beyeler (2001), two types of impacts are considered: changes in carbon stocks and changes in biodiversity. As a quantitative indicator for carbon stocks we use aboveground biomass (AGB) per unit area. As an indicator for biodiversity, plant species richness per unit area is computed. Species richness is the most commonly used biodiversity indicator (Immerzeel et al., 2014), and we focus on plants because modelled plant distribution maps were available for Kalimantan from a study of Raes et al. (2013). In addition, patch size distribution of closedcanopy forest is computed to have a quantitative indicator of ecological connectedness in general (Jaeger, 2000). The indicators are calculated based on data from 1990 to 2009, and on the outputs of the PCRaster Land Use Change model, PLUC, (Verstegen et al., 2012) for the four scenarios from 2009 to 2030. To quantify the uncertainty in these indicators, we use error propagation to propagate the model structure uncertainties from PLUC and combine these with error information from inputs for the impact assessments. We also apply a sensitivity analysis (Convertino et al., 2014) to compute the individual contribution of these two components (PLUC and the impact assessment) to the total variance in each indicator. Our specific research questions are: 1) How do the projected impacts vary over the four scenarios for the three indicators?, and 2) Are the differences between the four scenarios for 2030 significant, in other words, is the motivation to opt for a specific scenario conclusive given the uncertainty in the indicator values?

Case study area
West Kutai and Mahakam Ulu districts (which were a single district before 2012) are located in the western part of the province of East Kalimantan, Indonesia, and have a total land area of 3.3 Mha (Fig. 1). From the 1950s onward, the Indonesian government tried to incentivise migration of people from Java and Bali into the less densely populated Kalimantan (and other provinces) through a transmigration program. One feature of the program was granting the people with logging concessions and land development licenses. Kutai Barat and Mahakam Ulu received 60% of the people that migrated to East Kalimantan, which was the largest share (Clauss et al., 1988).
Because of these processes, LULC changed in about one-third of the study area between 1990 and 2009, leading to a 9% decrease in forest area (van der Laan et al., 2018). Whereas the original forests in the region are dominated by closed canopy forests mostly from the family Dipterocarpaceae (Slik et al., 2002;Toma et al., 2005), the current situation has multiple forest types with different canopy covers, ranging from primary forests with a closed canopy cover to secondary forests with a medium open and very open canopy (Fig. 1). Changes have taken place particularly in the lowlands in the south-east. The higher altitudes up to ∼2200 m in the north-west of the area are still covered with primary and secondary forest. To mitigate carbon emissions and  (2014)).
biodiversity losses, Indonesia has extended its nationwide Moratorium of 2011 on new licenses for concessions on primary forests and peatlands (Presidential Instruction 10/2011; 6/2013; 8/2015). The still growing local population, and the international and domestic demands for cash crops (mainly oil palm and rubber), food crops and timber, could lead to further LULC changes in East Kalimantan in the future, resulting in impacts on carbon stocks and biodiversity. Projecting the impacts quantitatively under different potential policy landscapes can help to inform decision makers about the effectiveness of policy instruments with respect to sustainability.

LULC change scenarios and land use change model
In previous work, we developed four LULC scenarios for West Kutai and Mahakam Ulu districts (van der Laan, 2016) on the four extremes of two axes: land development: limited (L) vs. unlimited (uL), and zoning: restricted (R) vs unrestricted (uR). The narrative behind these scenarios is that amount of land development can be influenced for example through the transmigration program as well as through yield improvements, whereas zoning can be installed for example through a Moratorium (see Section 2), thereby controlling deforestation in particular locations. The four scenarios combine land development and land zoning restrictions, as follows: 1. limited developmentrestricted zoning (LR), 2. limited developmentunrestricted zoning (LuR), 3. unlimited developmentrestricted zoning (uLR), and 4. unlimited developmentunrestricted zoning (uLuR).
For this previous study, LULC maps were acquired for 1990, 2000 and 2009 at a 250 × 250 m resolution (van der Laan, 2016(van der Laan, , 2018Budiman et al., 2014). Ten LULC types were distinguished in these maps, namely closed canopy forest, medium open canopy forest, very open canopy forest, grass and shrubs, cleared land, mixed agriculture (including rice and other food crops), rubber (mostly smallholder), pulpwood plantation, monoculture oil palm, mining (mostly coal) and settlements (Fig. 1).
We used the land use change model PLUC (Verstegen et al., 2012), an open source model implemented in Python with the PCRaster Python framework (Karssenberg et al., 2010) for future LULC change projections. PLUC simulates the locations of LULC change for a selected number of LULC types. The projected change is a function of an exogenous time series of demands for land or products and a weighted set of drivers of location (such as slope, distances to primary roads, secondary roads, rivers and towns, etc.) and no-go areas per active LULC type. For the West Kutai and Mahakam Ulu district case study, we used six active LULC types: mixed agriculture, rubber, pulpwood plantation, oil palm, mining, and settlements. The weights of the individual drivers of location for these active LULC types were identified using a forward (Wald) stepwise binary logistic regression (van der Laan, 2016). PLUC is a stochastic model, able to take into account uncertainty in inputs, model structure and parameters (Verstegen et et al., 2014). For the West Kutai and Mahakam Ulu district case study, only uncertainty in the model structure was accounted for by adding for each LULC type uniformly distributed random noise as a driving factor of location, with a weight of 1-R 2 , where R 2 is the proportion of the variance in the past LULC change that is explained by the other drivers of location in the logistic regression model of that LULC type. This resulted in random noise weights varying from 0.31 (for pulpwood plantations) to 0.53 (for mining). Using PLUC, annual maps of LULC were projected for 2010 until 2030 under the four scenarios in a Monte Carlo assessment with an ensemble of 100 members. In our analysis in the current paper, we use the map ensembles for 2020 and 2030 to derive the indicators, see next sections. Fig. 2 gives an overview of the steps to calculate the aboveground biomass (AGB). As pre-processing, we estimate the mean AGB (t/ha) for each LULC type on the maps. To do so, we use a 'strata approach' (Gibbs et al., 2007) in which a mean and standard deviation of AGB are attributed to each LULC type.

Carbon stock indicator: aboveground biomass
For the estimation of the AGB of oil palm plantations and the three forest types we collect data in the field in East Kalimantan in 2011 and 2012. It is important to include both diameter at breast height (DBH) and tree height (H) in the estimate of AGB, as tree height varies strongly for any given tree diameter in tropical trees (Chave et al., 2014;Feldpausch et al., 2012). So, DBH and H are measured for each dead and alive tree, palm and liana in 42 plots of 0.2 ha over a range of undisturbed and degraded (by logging and fire) forest types and in 6 plots of 0.2 ha in oil palm plantations, (see for further details Appendix 1). Next, the plot AGB is estimated using allometric equations of Chave et al. (2005) for moist tropical forest (Eq. (1)) and Frangi and Lugo (1985) for oil palm (Eq. (2)). (1) With: AGB est = estimated AGB (kg) per tree ρ = wood density (0.66 t m −3 , regional value based on Brown (1997)) D = tree DBH (cm) H = tree height (m) With: AGB est,oil palm = estimated AGB (kg) per oil palm tree H stem = stem height (m) Subsequently, we sum the AGB values (kg) of all growth forms for each 0.2 ha plot, convert from kg to tonnes, and from 0.2 ha to 1 ha, to come to total AGB in t/ha. Variation in the mean AGB can exist due to, for example, the variations of growth forms within one LULC type (monoculture, mixed), the type of land clearance or disturbance (including logging and fire) or the management level. We use our field data in combination with AGB values from literature to calculate a frequency distribution around the mean AGB assuming a normal distribution (see Appendix 1 and 2). For the remaining LULC types, we do the same but with literature values only (see Appendix 2). As Step 1 of the carbon stock impact analysis (Fig. 2), a random AGB value is drawn from its normal distribution for each LULC type and assigned to all cells of that LULC type in the LULC map, because the assumed distribution represents uncertainty about the mean and not variation over space. This step is executed for each scenario and for each year in the analysis : 1990, 2000, 2009, 2020 and 2030. For the past data, a single LULC map is available per year, whereas for the projections the ensemble of 100 maps exists, from which one is selected each time. This step is repeated 500 times (5 times per LULC map in the ensemble for the projections) per scenario per year. Subsequently (step 2a), the local mean AGB and standard deviation (sd) are calculated over all 500 realizations. We also calculate the local coefficient of variation (sd/mean), as a relative measure of the AGB uncertainty, such that uncertainty is comparable between indicators. Local, in this context, means per cell, as opposed to the mean AGB aggregated over the whole study area, which is calculated per realization in the next step (2b in Fig. 2). All these averages together are then ordered and represented as a boxplot.

Biodiversity indicator: dense forest patch size distribution
Dense closed canopy forest is the 'undisturbed' state of the natural ecosystem in Kalimantan. The change in dense forest patch-size distribution, which is the frequency of patches in each patch-size category, is used as an indicator of disturbance for both plants and animals. To detect a dense forest 'patch', we used the Moore neighborhood (8 surrounding cells) to define a contiguous group of cells of the dense forest LULC type (Fig. 2, steps 3 and 4). The sizes (m 2 ) of all patches are calculated for each LULC map in the Monte Carlo ensemble, using PCRaster map algebra functions (Schmitz et al., 2013). Mean and standard deviation in number of patches per size category are plotted on a log-log scale as is common for patch size distributions (Meloni et al., 2017). Based on the output patch sizes, we decided to use fifteen size categories between 10 4 and 10 11 m 2 , equally spaced on a logarithmic scale, so the class boundaries are (in m 2 ) 1.0•10 4 , 3.2•10 4 (=1.0•10 4.5 ), 1.0•10 5 , 3.2•10 5 (=1.0•10 5.5 ), etc. There is no need to calculate patch sizes multiple times per LULC map, as is done for the other indicators, because the calculation contains no stochastic variables besides the LULC map itself.

Biodiversity indicator: plant species richness
As an input to calculate the plant species richness, we use existing plant distribution maps derived from species distribution models (SDMs), produced by Raes et al. (2013). SDMs are statistical models that associate locations where species are present with environmental predictors (Elith and Leathwick, 2009). SDMs result in a spatial prediction of habitat suitability for a given species, which can then be converted to species distribution maps. The plant families included in this study are: Dipterocarpaceae (includes Meranti hardwood, 21% of total), Leguminosae (legume family, 19% of total), Lauraceae (laurel family, 16% of total), Moraceae (mulberry-or fig family, 15% of total), Fagaceae (beech family, 10% of total), Ericaceae (heather family, 9% of total), Myristicaceae (nutmeg family, 7% of total), Sapindaceae (soapberry family, 4% of total) (Raes et al., 2013). These plant families include trees, shrubs and grasses, see Table 1 and Appendix 3.
To estimate how LULC change impacts plant species richness per raster cell, we estimate the number of plant species that remained in each raster cell after LULC change has occurred (Table 1). Hereto, we assumed a species remaining range (in %) for each LULC type, based on the method of Lucey et al. (2015). The species remaining range is defined by the minimum and maximum possible percent of the original presence of each plant family that is expected to remain after the conversion to the LULC type. The species remaining range is a function of the type of land cover (i.e. natural, mixed agriculture, monoculture, cleared), the management and disturbance level of the land cover (none, low, medium, high), the growth forms that may occur in the particular LULC type (trees, shrubs or herbs), and the growth forms that occur in a particular plant family. By defining a species remaining range, instead of a single value, we account for the fact that species richness is not only impacted by LULC change, but also by e.g. overexploitation, invasive species, pollution and climate change (MEA, 2005). The assumptions behind the species remaining ranges are explained in more detail in Appendix 4.
The eight 'undisturbed' (pre-1990) species distribution maps (Raes et al., 2013) have a resolution of 9.8 × 9.8 km at the equator, which is much coarser than the resolution of the LULC maps, which are nominal maps at 250 × 250 m. Therefore, (step 5a in Fig. 2)w efirst calculate the fractions of LULC types within each 9.8 × 9.8 km raster cell. In the same step, we draw a random value from a standard uniform distribution, U(0,1), for each cell. For each LULC typeplant family combination, we take the value in the species remaining range according to this random draw. For example, when a random value of 0.4 is drawn for a cell that has 10% shrubs: 3 + 0.4 * (10-3) = 5.8% of all Myristicaceae remain, 25 + 0.4 * (100-25) = 55% of all Ericaceae and Leguminosae remain, and 0% of all other families remain in the shrub land within that cell (see Table 1 for the ranges). Next, we multiply the percentages of these species with the fraction of the cell that is covered by shrubs (10% in the example), and finally with the cell values in the 'undisturbed' species distribution maps for the corresponding families. This is done for all LULC types and all plant families, which are finally summed to obtain the total number of plant species remaining per cell.
In the same way as for carbon stocks, this is repeated 500 times (5 times per LULC map in the ensemble for the projections) per scenario per year (step 5b in Fig. 2). Subsequently (step 6a), the mean, standard deviation (sd), and coefficient of variation are calculated per cell over all 500 realizations. And (step 6b), the mean number of species remaining aggregated over the whole study area is calculated per realization and represented as a boxplot.

Significance test and sensitivity analysis
We test whether there is a significant difference in the spatially aggregated AGB and plant species richness between the five points in time, 1990, 2000, 2009, 2020, and 2030, and the four scenarios, to answer the research question whether the motivation to opt for a specific scenario is conclusive given the uncertainty in the indicator values. Hereto, the Kruskal-Wallis test (one-way ANOVA on ranks) (Kruskal and Wallis, 1952) is applied on the areal AGB and plant species richness means. The Kruskal-Wallis test is chosen because it is non-parametric (normality of the resulting AGB and plant species richness distributions cannot be ensured) and tests over independent samples, which our scenarios are. A significant Kruskal-Wallis test indicates that at least one scenario is significantly different from the others. If this is the case, we apply post-hoc tests after Nemenyi (1963) with Chi-squared approximation for independent samples to identify which scenario is different from which other scenario(s).
The Kruskal-Wallis test is not applied to the dense forest patch size distribution, because it is not a single value indicator. To test for one point in time, each patch size category would have to be compared for each scenario, leading to six scenario combinations times twelve size categories, that is 72 p-values. Including the five points in time would lead to 360 p-values.
The significance analyses are done in R (R core team, 2018) with the 'stats' and 'PMCMR' (Pohlert, 2014) packages. The significance of pairwise comparisons between AGB and plant species richness in different scenarios and at the five points in time are plotted with the compact letter display method (Piepho, 2004), using the 'multcompView' package (Graves et al., 2015). To assign the letters in this plot, we use a significance level (p-value) of 0.01, to be conservative about differences. Yet, the p-values of post-hoc tests are provided in Appendix 5,so the reader has the option to see the result for any other significance level.
Finally, a global sensitivity analysis is performed to compute the contributions of the land use change model and the impact assessment models to the total uncertainty in each indicator in 2030: spatially aggregated AGB, local AGB, dense forest patch size distribution, spatially aggregated plant species richness, and local plant species richness. Hereto, we use the Sobol' method (Sobol', 1993in: Convertino et al., 2014. The Sobol' method represents the contributions of all model components (here the land use change model and the impact assessment model) as a fraction of the total variance in the output(s) (here each indicator). To calculate these fractions, we run all our analyses twice more: once with a single output of PLUC instead of all 100 MC samples, such that only uncertainty in the impact assessment model is included, and once with the mean of the plant species richness range and the mean of the AGB distribution per land use type, such that only uncertainty in the land use change model is included. For each run the indicator variances are computed for the year 2030; the 'local' variance is the sum of the variances in all raster cells. The indicator variances in both runs are divided by the variance of the same indicator in the original run (with both components uncertain) to get the contribution.
The sum of the fractions of all components is one (100%) for additive models and less than one for non-additive models, where the difference (one minus the sum) designates the influence of interactions between the components on the output variance (Convertino et al., 2014). In our case study, we do not expect interactions, because the land use change model and the impact assessment models are not tightly but loosely coupled (Fig. 2), so there are no feedbacks between them.

Impacts of LULC change on carbon stocks
We found a large variation in AGB for every LULC type in the literature (see Appendix 2 for an overview). The highest mean AGB (dry matter) was reported for closed canopy forest (411.0 ± 123.4 t/ha) and the lowest for mining (14.0 ± 7.6 t/ha). Furthermore, we found a mean AGB of 205.5 ± 69.7 t/ha for medium open canopy forest and of 102.2 ± 88.0 t/ha for very open canopy forest. For rubber, we found a mean AGB (129.7 ± 128.4 t/ha) that was higher than the mean AGB of  very open canopy forest and of pulpwood (88.1 ± 64.7 t/ha) and oil palm plantations (54.2 ± 34.5 t/ha). Lower mean AGB values are reported for shrubs, grassland and settlements (27.7 ± 28.1 t/ha). Between 1990 and 2009, LULC change resulted in a mean AGB decline in our study area from about 280 to 230 t/ha (Fig. 4), which is about 18%. The decline took place mainly in the South-East, North-East and North-West of Mahakam Ulu, and in the South-East and South-West of West Kutai, where the expansion of agriculture and mining was projected to occur (green colours in Fig. 3). For every raster cell the coefficient of variation (cv) shows the variation in local AGB estimates (size of circles in Fig. 3). The cv of the AGB in 1990 ranges from ≤30% in the undisturbed North, to up to 80% in the slightly disturbed South of Mahakam Ulu region, to > 80% in West Kutai. Although the local cv (relative measure) are the lowest (mean of 42%) in 1990 compared to the other years, the AGB interquartile ranges over the whole area (absolute measure, length of the box plot in Fig. 4, left panel) are the highest for 1990.
Between 2009 and 2030, the projected LULC change follows different trajectories under the four scenarios, with different impacts on AGB for at least one scenario as indicated by a strongly significant overall Kruskal-Wallis p-value < 2.2 10 −16 . The two scenarios with limited development, LR and LuR, show mean AGB patterns similar to The two scenarios with unlimited development, uLR and uLuR, have mean AGB patterns distinct from each other and from 2009, both over space (Fig. 3) and over time (Fig. 4). The strongest decline in mean AGB under the uLR scenario is projected to occur in the centre of West Kutai and in the western valley in Mahakam Ulu (Fig. 3). Under the uLuR scenario, the decline occurs in the entire West Kutai, and in Mahakam Ulu everywhere but in the northern highlands (Fig. 3). The mean AGB over the whole area decreases from around 230 t/ha in 2009 to 220 t/ ha in 2030 (∼4%) under the uLR scenario and to 160 t/ha (∼30%) under the uLuR scenario, with a steep decline in the last 10 years. Local cvs in projected AGB are high everywhere under the uLuR scenario, except where cultivation is made impossible by the landscape characteristics such as slope (Fig. 3). Under the uLR scenario, the cvs are low in the zones where expansion is prohibited by the zoning. The Nemenyi post-hoc test strongly rejects the hypothesis of similarity of mean AGB between uLR and uLuR, and between them and the two limited development scenarios (LR and LuR) (different letters Fig. 4, see also Appendix 5). Only between LuR and uLR, the difference is not significant in 2030 (both have letter e in 2030 in Fig. 4).
The global sensitivity analysis for 2030 shows that the variance in the mean AGB over the whole area is sensitive to the uncertainty in the impact assessment model only (100%), i.e. the AGB estimates per LULC type. On the other hand, the variance in the local AGB does also depend on the land use change model (Fig. 5). The contribution of the uncertainty in the land use change model to the variance in the local AGB varies between the scenarios, from 4.0% in LR to 18% in uLuR.

Impacts of LULC change on biodiversity
Recent patch size distributions of dense closed canopy forest (1990,2000,2009, right-hand side of Fig. 6) have the largest number of patches, around 50, in the patch size categories of 100-300 and 300-1000 ha. Over the subsequent size categories, the number of patches decreases exponentially (straight line in log-log plot Fig. 6). The three largest size categories, up to 3 Mha, have one patch each, of which the largest one (just over 1 Mha) is split in two halves in 2009, showing that the remaining part of continuous forest started to fragment between 2000 and 2009. In the three smallest size categories of 3-100 ha, approximately 25 patches were present between 1990 and 2009.
Between 2009 and 2030, dense forest fragmentation follows different trajectories under the four scenarios (four left panels in Fig. 6). Under the two scenarios with limited development, LR and LuR, patches in 100-300 and 300-1000 ha categories are split up, leading to more small patches of up to 3 ha. The patches in largest size categories that were present in 2009 mostly stay intact, i.e., are projected not to be fragmented under the limited development scenarios, although there is some uncertainty about this in the LuR scenario for the category 10000-30000 ha. This uncertainty in patch size is controlled by the land use change model only (Fig. 5).
Under the uLR scenario, patch fragmentation in the intermediate size categories is more pronounced. In this scenario in 2030, the small size categories even have a higher number of patches than the most prominent 100-300 and 300-1000 ha categories between 1990 and 2009. The large size patches remain intact in 2030 under this scenario. In the uLuR scenario, the largest patch remains in 2030, but the single largest patch is split up after 2020. The number of very small patches increases to about 1000 in 2030 under the uLuR scenario.
According to the SDMs in combination with the 1990 LULC map, the highest plant species richness in our region can be found in the Mahakam Ulu region in the recent past, with 300 up to 600 plant species per raster cell. In West Kutai, the local variation is larger with 35 (in the East) up to 600 (in the West) species (Fig. 7). The mean over the whole region is around 415 plant species per cell. Between 1990 and 2009, the mean species richness over the whole region was estimated to decrease from 415 to 385 per cell (∼7%) (Fig. 4, right panel). The losses mainly occur in the South-West of West Kutai, and the local cvs are relatively low with a mean of 0.11 in 1990 (Fig. 7). The relative decline in plant species richness in this period is the largest for the family Sapindaceae (10%), followed by Dipterocarpaceae (9.7%), Myristicaceae (9.2%), Lauraceae (8.4%), Moraceae (7.6%), Leguminosae (7.2%), Fagaceae (6.7%) and Ericaceae (3.4%) (Fig. 8).
Between 2009 and 2030, LULC change follows different trajectories under the four scenarios, with different impacts on species richness (overall Kruskal-Wallis test with a p-value < 2.2 · 10 −16 ). Under the two scenarios with limited development, LR and LuR, a significant decrease in mean species richness from 385 to 375 (3%) is projected between 2009 and 2030. The two scenarios show mean plant species richness patterns similar to each other, both over space (Fig. 7) and over time (Fig. 4). Accordingly, the Nemenyi post-hoc test does not reject the hypothesis of similarity in species richness between LR and LuR in 2030 (both have letter d in Fig. 4).
The two scenarios with unlimited development, uLR and uLuR, have mean plant species richness patterns different from each other and from the species richness in 2009, both over space (Fig. 7) and over time (Fig. 4). The mean number of species over the whole area decreases from around 385 in 2009 to 345 in 2030 (∼10%) under the uLR scenario and to 230 (∼40%) under the uLuR scenario. The latter scenario has a steep decline in the last 10 years, similar to the AGB. The Nemenyi post-hoc test strongly rejects the hypothesis of similarity of species richness between these scenarios, and between them and the two limited development scenarios (LR and LuR) in 2030 (Fig. 4, see also Appendix 5).
The global sensitivity analysis for 2030 shows that the variance in the mean plant species richness over the whole area and the variance in the spatially explicit species richness are similar and both mainly controlled by the impact assessment model, i.e. by the estimated ranges of species remaining after a LULC change (Fig. 5). The contribution of the uncertainty in the land use change model to both variances varies over the scenarios from 0.18% (aggr) and 0.22% (local) in uLR to 4.3% (aggr) and 4.4% (local) in uLuR (Fig. 5).

Impacts of LULC change on carbon stocks
Literature reports that on average native forests in the region have an AGB of about 411+/− 123 t/ha, while secondary forests have slightly lower AGB values (see Appendix 2) and might be more resilient to losses in AGB (Poorter et al., 2016). Between 1990 and 2009, we found an 18% decrease in AGB. Our projections point to a loss of AGB in West Kutai and Mahakam Ulu, ranging from 0% (LR and LuR) to 4% (uLR) to 30% (uLuR) between 2009 and 2030. The uLuR scenario will likely lead to an extreme unprecedented deforestation, similar to other areas in Kalimantan (Carlson et al., 2012b, e.g. Carlson et al., 2012a, in Papua New Guinea (e.g. Bryan et al., 2010 found losses of 25% over 30 years) and in Sumatra (e.g. Kotowska et al., 2015). The significant differences in policy scenarios suggest that limiting the development, of mainly large-scale cash-crop cultivation, is the most effective measure to reduce negative impacts on AGB, and that restricted zoning is effective under high development scenarios. On the contrary, under the limited development scenarios, restricted zoning does not significantly increase (or reduce) AGB in 2030. This is good news from an economic point of view, as it means that enforcement of protected areas does not interfere with agricultural profits, up to a certain amount of land development.  Fig. 4) and 'local' is spatially explicit (cf. Figs. 3 and 7).

Impacts of LULC change on biodiversity
The patch size distributions between 1990 and 2009 suggest that the size categories of 100-300 and 300-1000 ha are the typical patch sizes of our study area, characterizing the ecosystem at the chosen scale (Kéfi et al., 2014). Under all scenarios, some of the patches in these previously most prominent categories are fragmented into very small patches. In the uLuR scenario this occurs to such an extent that the patch size distribution approaches the shape of a power law (linear on the log-log plot in Fig. 6). This means that the dense closed canopy forest achieves a 'scale-free' patch-size distribution, without any typical patch size (Kéfi et al., 2014). The higher number of small patches also implies a relatively large edge length, associated with additional disturbance (Harper et al., 2005). This is in line with recent studies showing that deforestation in the tropics is resulting in large scale fragmentation and loss of typical patch size distributions (Taubert et al., 2018), with long lasting impacts (Haddad et al., 2015), even in forests that remain intact (Betts et al., 2017).
A recent review report from the IUCN Oil Palm task force found that about 99% of tree and plant species were affected by oil palm conversion (Foster et al., 2011), and 90% of mammal species have lost habitat (Meijaard et al., 2018). Our projections show a decrease in plant species richness between 2009 and 2030 ranging between the four scenarios from 3% (LR and LuR) to 10% (uLR) to 40% (uLuR). Although our simulations do not allow to model extinctions of species, the projected impacts are already substantial in terms of habitat loss, species richness decrease and fragmentation. The positive effects of limiting development on species richness were clearly larger than the effects of zoning; however, we also find an important role of zoning on species losses under unlimited development, particularly because of the overlap between (remote but) suitable areas for oil palm and high biodiversity areas (Fitzherbert et al., 2008). This finding is in line with IUCN's suggestion that zonal restrictions could prevent major biodiversity losses and that species extinctions could be avoided if plantations were developed in secondary forest or degraded areas (Meijaard et al., 2018).
The list of plant families analysed in this study includes the top five families in Kalimantan in terms of number of endemic species, namely the: Dipterocarpaceae (160 endemic species/272 total number of species), Ericaceae (118/132), Myristicaceae (69/115), Moraceae (62/175), and Fagaceae (44/89) (van Welzen and Slik, 2009). Substantial decreases in species richness amongst these plant families, as were found in this study, can result in a serious threat to endemic species if land development continues at a high pace. The largest loss of species toward 2030 was projected in the Dipterocarpaceae family. This was expected since the lowland forests, usually dominated by species of the family Dipterocarpaceae (Ashton, 2004;Whitmore, 1984), are relatively easy to reach and their timber is commercially valuable (Slik et al., 2003). A continued decrease in Dipterocarpaceae-dominated forests can be disastrous, since such forests are among the richest worldwide in terms of flora and fauna (Whitmore, 1984). Deforestation of such forests is particularly alarming when occurring in Kalimantan, because this is where the greatest diversity of Dipterocarpaceae occurs (Ashton, 1982).

Comparison between indicators
The highest decreases in AGB and plant species richness, and the highest fragmentation rates all occur under the same scenario, uLuR. Similarly, LR and LuR are the best scenarios for all three indicators, i.e. with the least negative impacts. These findings suggest that there are no trade-offs between conserving carbon stocks and biodiversity in our study area over our four scenarios, given the applied indicators. Contrary to this, Murray et al. (2015), found low carbon stocks along with higher species richness. Yet, our conclusion is only valid for our study area under the given scenarios; a more in-depth approach is required to study whether this remains valid under other conditions (see e.g. Verstegen et al., 2017).
The model projections are less variable for plant species richness than for AGB (much lower cv for plant species richness than AGB; e.g. mean of 0.11 in 1990 for species vs. 0.42 in 1990 for AGB) (Figs. 3 and  6). The first reason for this is that we have pre-1990 species distribution maps, which give us a precise undisturbed number of species value to start with, whereas for AGB the initial values had to be estimated first, causing additional uncertainty. The second source of the lower cv in species richness, compared to AGB, is the lower resolution (250 m vs. 9.8 km, see Fig. 2). The lower resolution cancels out local variation in allocated agricultural expansion stemming from the model structure uncertainty of the land use change model. For example, consider the situation that in some of the 100 LULC maps for 2030 cell x is selected for expansion, and in others one of his neighbours. When scaling the LULC maps up to the 9.8 km cells of the biodiversity impact analysis (step 5a in Fig. 2), this variation is likely to be 'lost', when both cell x and his neighbour fall into the same 9.8 × 9.8 km cell. The third reason is that the variation in species impact (% of species remaining after LULC change, Table 1) is assumed to be local variation, whereas for AGB variation in the mean is assumed (see Fig. 2). The local variation is cancelled out at the spatially aggregated level of the boxplots. The latter also causes the difference in the contribution of the uncertainty of the land use change model to the total variance: the variance in the mean AGB over the whole area is only sensitive to the uncertainty in the impact assessment model (Fig. 5), because, in our method, AGB only depends on the area of the land use types and not on their location. The initial plant species distribution maps are spatially explicit, so the location where the land use is being changed determines the impact. The contribution of the uncertainty in the land use change model to the total variance in the indicators is higher for scenarios with a higher amount of land use change, as more change means more 'actions' by the land use change model. No interaction effects between the two models were found, as expected in our case of loose model coupling.

Limitations and future work
We are able to show that propagating uncertainties from a land use change model to the impact assessment is straightforward and can help to answer questions about the significance of differences between future scenarios. However, we are aware that our estimation of local impacts, average impacts and uncertainty ranges have room for improvement in many respects, which will be discussed here with respect to the land use modelling, the carbon stock impacts, and the biodiversity impacts. Regarding the land use change modelling, future work  Raes et al. (2013). Colour shows the mean over all realizations per cell. The size of the yellow points shows the coefficient of variation (cv) per cell. In the class ranges, a square bracket indicates inclusion, while a round bracket indicates exclusion.
could take into account land use change model input and parameter uncertainty besides model structure uncertainty, or even try to quantify the probability of more abrupt changes in system and thus in model structure (Müller et al., 2014;Verstegen et al., 2016b). Also, modelling a larger area or multi-scale modelling with integrated land use change models at different scales could lead to a generalization of insights (see e.g. Verstegen et al., 2016a). In this respect, we stress that it is crucial to make models and/or model outputs available to a wider audience, in order to allow others than the land use change modellers themselves to conduct model-based analyses. Outputs of stochastic models can be sizeable, so making the models themselves open, e.g. with executable research compendia (Nüst et al., 2017), is a promising option to stimulate such openness.
Impact assessments of LULC change on carbon stocks should also include estimates of impacts in below ground carbon and soil organic carbon (SOC). Below ground carbon was not considered in our study because estimating initial below ground carbon requires good quality soil maps, which were not available for the area. Estimates on emissions from the below ground carbon stock component correspond to ∼32% of the global carbon stock, but are less robust (Doetterl et al., 2015), and we know that the type of LULC change might drive carbon dynamics in different directions van der Putten et al., 2009). For example, Don et al. (2011) showed increases and decreases in SOC with different land-use types and their transitions in tropical areas, and estimated that the effects of land use changes might be underestimated by 28% if SOC is not accounted for.
With respect to biodiversity, future work may focus on quantifying the accuracies of the initial species richness, including other biological species, and collect more data to be able to apply other indicators that more explicitly calculate ecosystem functionality. Moreover, as more research develops, we could have a sounder theoretical and empirical support for the species reduction ranges (Table 1) and/or have a better estimation of these ranges as the indicator is very sensitive to them (Fig. 5). Further it would be important to account for time lags of impacts as well as to make an explicit connection between patches and species which is now not possible due to the mismatch in resolutions.
The latter would allow for the calculation of connectivity metrics (Calabrese and Fagan, 2004).

Conclusions
In this study, we aimed to estimate the impacts of four contrasting land use and land cover (LULC) change scenarios on carbon stocks and biodiversity, explicitly taking into account the propagated uncertainties from the land use change model used to run the scenarios. We studied the impacts of agricultural expansion in West Kutai and Mahakam Ulu districts in Kalimantan, Indonesia. Using LULC data from 1990 to 2009 and LULC projections towards 2030 from the land use change model PLUC (Verstegen et al., 2012), we study changes in carbon stocks using aboveground biomass (AGB) as indicator, and changes in biodiversity using dense closed-canopy forest patch size distribution and plant species richness as indicators. The LULC projections were done for four scenarios on the combined extremes of two axes: land developmentlimited (L) vs. unlimited (uL) -and zoningrestricted (R) vs unrestricted (uR).
The answers to our two research questions led to the conclusion that the most sustainable pathway for East-Kalimantan in terms of carbon stocks and biodiversity is to limit land development as this leads to the lowest negative impacts. Limiting development can be achieved by, for example, limiting transmigration to East-Kalimantan, bridging yield gaps and/or better land use planning with respect to potential yields of the soils (e.g. Meijaard et al., 2018). We thereby acknowledge that these strategies might have other negative (or positive) impacts, locally as well as elsewhere (spill-over effects). Under unlimited development, zoning becomes important, particularly because of the overlap between (remote but) suitable areas for oil palm and high biodiversity areas (Fitzherbert et al., 2008).
With this paper, we have demonstrated a general methodology characterized by error propagation from a land use change model to the impact assessment. The modelled variances in patch size were fully controlled by the land use change model. The variances in the other indicators were found to be controlled mainly by the impact assessment models, and up to 18% by the land use change model, underlining that land use change model uncertainties should not be ignored in projected impact assessments. Furthermore, we have presented a new method to test the significance of differences between future scenarios, that is, to test if a potential policy instrument has a significant (positive) effect on the studied indicators, and is thus advisable to be implemented.