Nutrient Balance as a Tool for Maintaining Yield and Mitigating Environmental Impacts of Acacia Plantation in Drained Tropical Peatland—Description of Plantation Simulator

: Responsible management of Acacia plantations requires an improved understanding of trade-offs between maintaining stand production whilst reducing environmental impacts. Intensive drainage and the resulting low water tables ( WT ) increase carbon emissions, peat subsidence, ﬁre risk and nutrient export to water courses, whilst increasing nutrient availability for plant uptake from peat mineralization. In the plantations, hydrology, stand growth, carbon and nutrient balance, and peat subsidence are connected forming a complex dynamic system, which can be thoroughly understood by dynamic process models. We developed the Plantation Simulator to describe the effect of drainage, silviculture, fertilization, and weed control on the above-mentioned processes and to ﬁnd production schemes that are environmentally and economically viable. The model successfully predicted measured peat subsidence, which was used as a proxy for stand total mass balance. Computed nutrient balances indicated that the main growth-limiting factor was phosphorus (P) supply, and the P balance was affected by site index, mortality rate and WT . In a scenario assessment, where WT was raised from − 0.80 m to − 0.40 m the subsidence rate decreased from 4.4 to 3.3 cm yr − 1 , and carbon loss from 17 to 9 Mg ha − 1 yr − 1 . P balance shifted from marginally positive to negative suggesting that additional P fertilization is needed to maintain stand productivity as a trade-off for reducing C emissions.


Introduction
Peatlands extend from the polar regions to the humid tropics, and may hold over 1000 Pg of carbon (C) [1,2], more than the total pool of C in the atmosphere. In pristine peatlands the input of dead organic matter exceeds the rate of biological decomposition of the organic residues and therefore pristine peatlands tend to store C [3]. However, the cultivation of peatlands, which requires intensive water management and control of plant nutrition, may turn peatlands into a C source [4][5][6][7][8] and cause other environmental problems, such as peat subsidence [5,9], increased nutrient export to water courses [10], and increased fire risk [11].
A responsible approach to the bio-economy of cultivated tropical peatlands calls for an improved understanding of the trade-offs between maintaining crop production whilst reducing negative environmental impacts. In Southeast Asia, ombrogenous peatlands cover approximately 25 Mha [12]. Half of this area is under some form of agricultural use-either for large-scale plantations of oil palm or Acacia for pulp and paper production, or for smallholder farming [13]. Mitigating C losses and associated peat subsidence in these landscapes is of increasing importance in the context of global climate change agenda and regional land use planning [14].
When a peatland is drained, the same biogeochemical processes that enhance crop production are also responsible for increased C emissions, subsidence and nutrient exports. Lowering of the water table (WT) increases the rate of organic matter decomposition and consequently nutrient release. It is likely that this is the main mechanism through which drainage improves crop production; for example, on infertile ombrogenous peatlands, decomposition of organic matter provides the main source of nutrients for tree stands [15], with several field experiments demonstrating that drainage improves both nutrient availability in peat [16][17][18] and the growth rate of trees. On the other hand, there will also be a direct growth response attributable to amelioration of root oxygen (O 2 ) stress following the removal of excess water. Since most deep peats in South East Asia are composed of woody, fibric materials that have a low bulk density and high macroporosity, there is typically favourable aeration in the rooting zone above the WT and particularly in peatlands that have been drained. These arguments highlight the role of nutrition rather than O 2 stress in the drainage-induced growth response of forest stands on ombrogenous peat, with previous studies demonstrating that tree growth is frequently limited by phosphorus (P), potassium (K) and, at infertile sites, by nitrogen (N) availability [19][20][21].
If an enhanced nutrient supply is the driver behind the drainage-induced growth response, it provides us with tools to maintain crop production and to mitigate harmful environmental effects with a designed combination of WT control and careful nutrient management of a stand. In tropical peatlands, it has been estimated that raising plantation WT from −80 cm to −40 cm, would decrease the CO 2 emissions by 50 Mg ha −1 yr −1 [7] and the rate of peat subsidence by 1.7 cm yr −1 [9]. At the same time, however, nutrient release from decomposition of organic matter will decrease. This is shown for example by Marwanto et al. [22], who observed a sharp decrease in soil water cation concentrations after rewetting a peat profile following a dry period under oil palm plantation in Kalimantan, Indonesia. In order to maintain productivity at higher WT, it is therefore likely to be necessary to replace nutrients that are no longer being generated through peat decomposition via artificial fertilization. Fertilization can be designed to meet the plant nutrient needs, thereby fulfilling the production demand, whilst potentially decreasing the leaching loss of other nutrients. On the other hand, the use of artificial fertilizers incurs an economic cost, and risks greenhouse gas emissions displacement to fertilizer production or increased N 2 O emissions from the plantation. In the case of pulpwood plantations, however, there is considerable potential to recover nutrients from harvested biomass during the pulp production process, and to return these to the plantation as wood ash, presenting the opportunity for a productive, low-emission, 'circular' nutrient management system if implemented successfully.
Nutrient management requires a thorough understanding of stand growth and nutrient demand, litter input, and litter and peat decomposition and consequent nutrient release; and finally how these processes depend on the dynamics of the WT. In tree plantations on peatland, stand biomass growth and yield, nutrient supply, C balance, stand total mass balance, peat subsidence and peat hydrology are connected to each other in a tight and recursive manner resulting in a complex dynamic system. Attempts to optimize crop production targets within responsible environmental constraints therefore require a thorough understanding of this system, which can be provided by dynamic process models. Many process models accounting for plantation growth and yield, photosynthesis, and water and nutrient use in mineral soils have been developed [23][24][25][26], but these models are not particularly suitable for peatland sites, where hydrological conditions, and environmental and plant nutrition problems are different from mineral soils.
In this paper, we describe the Plantation Simulator, a model that describes biomass production, carbon and nutrient balance with respect to different plantation management operations such as drainage, silviculture, fertilization, and weed control. It allows us to find production schemes that are socially and environmentally bearable; socially and economically equitable; and environmentally and economically viable. Using a Monte Carlo -approach we ran the model with a large variety of plausible input data and tested the model performance against published data from similar conditions. After demonstrating plausibility, we studied the system behavior under different WT in order to optimize plantation growth by manipulating nutrient availability, whilst minimizing peat C loss and subsidence.
We apply the Plantation Simulator to outline management scenarios for an Acacia plantation, where requirements for WT management have recently changed. After the severe fire-related haze (air pollution) season in 2015, the Indonesian government released a regulation obliging all managers of plantations on peat to maintain an average WT within −0.40 m of the peat surface with the primary aim of reducing fire risk whilst also seeking to reduce emissions of greenhouse gases and other pollutants (SK.22/PPKL/PKG/PKL.0/7/2017). Previously, standard Acacia plantation operating procedures have generally aimed at keeping the WT in the region of −0.70 m. Given the rapid implementation of the requirement for higher WT in plantations, there has been little opportunity to gather datasets on the implications that implementing higher WT management would have for plantation growth and yield. This lack of experimental evidence can, however, be overcome through the application of a modelling study that merges theoretical understanding of biogeochemical processes. A higher WT will reduce the rate of organic matter decomposition and the subsequent release of nutrients; and is therefore likely to reduce stand growth. However, it will introduce environmental benefits by decreasing rates of peat subsidence and greenhouse gas emissions. In this study we ask: 1.
how the nutrient balance of an Acacia stand will change under a higher WT regime; 2.
what kind of nutrient management should be applied to maintain the yield under higher WT management; and 3.
how much the rate of peat subsidence and CO 2 efflux will be decreased under this management regime?
Uncertainties of the model are discussed and the needs for future field measurements to decrease these uncertainties are suggested.

Acacia Plantation in Drained Tropical Peatlands
Acacia crassicarpa (Leguminosae) is a fast-growing, nitrogen-fixing tree that is the principal pulpwood plantation species grown on peat soils in SE Asia. It is tolerant of a wide range of soil types and suitable for low fertility sites [23]. The typical plantation rotation period between planting of tree seedlings to harvest, is four to five years, and a closed canopy develops at around 12 to 18 months. Tree height at harvest is in the range from 19 to 24 m. In order to provide suitable rooting conditions, peatland plantations are drained by a network of canals, typically 5 to 8 m wide and 500 to 800 m apart, and by smaller field drains. One dose of artificial P fertilizer is applied at the time of planting [24]. Weeding occurs several times before canopy closure and at two years of age any self-seeded or invasive tree species are removed from the plantation.

Study Outline
The Plantation Simulator is based on a mass balance (biomass, C, N, P, K) computation of drained tropical peatland. We explicitly account for organic matter produced by the plantation, returned to soil as litter, and decomposed and lost as CO 2 efflux to the atmosphere. The processes driving the mass balance are complex, and comprehensive experimental data for model validation are not available due to confidentiality of data owned by plantation companies. Therefore, the model is tested against published data and by filling the data gaps with required assumptions. First, the overall model structure is described (detailed description in Appendix A), then the model performance in mass balance calculation is shown against published subsidence data. Because the subsidence reports contain little information about the stand properties and dynamics, we used a Monte Carlo approach, deriving unknown site properties from a probability distribution and computing a large simulation set. The results of the simulation set were compared to available independent, published subsidence data sets from similar land use, climatic conditions and peat formations [5,9]. Subsidence is a suitable variable for model testing because it integrates input and output of organic matter through a long time period, whereas direct measurements of the C balance of tropical peatlands are limited to a very small number of shorter-term studies. After having demonstrated plausibility from the mass balance viewpoint, the nutrition simulation results were evaluated against known fertilization schemes. Thereafter, the model was applied for two WT schemes, and the nutritional consequences of the different WT were identified ( Figure 1).

Overall Model Description
The Plantation Simulator describes an Acacia plantation as a system composed of tree stand, soil, nutrient, and hydrology modules ( Figure 2). Basic plantation inventory data such as tree species, planting density, site index describing the growth potential of the site [25], and literature-derived parameters are required as inputs for the model. Aspects of the plantation management regime including weeding, and fertilization are given as input parameters to the simulation. The Plantation Simulator computes growth and yield of the Acacia stand, organic matter decomposition and stand nutrition and allows quantification of the impact of plantation management and drainage on these processes. A one-month time step is used. A recursive loop returns information concerning decision criteria, i.e., production, carbon and nutrient balance and peat subsidence, to the model user, thereby enabling adjustment of management procedures in support of improved production and environmental management practices ( Figure 2). A detailed model description is provided in the Appendix A. Figure 2. Structure of the Plantation Simulator. The plantation system (dark grey panel) accounts for growth and yield of the stand, water table and biogeochemical processes within the site, and returns target variable information to the management system. Target variables, i.e., current yield, nutrient status, carbon balance and subsidence rate are evaluated against decision criteria and the management practices can be adjusted to reach better production and environmental targets. Codes A1, A11, A12, etc. refer to chapters in the detailed model description in Appendixes A.1, A.1.1, A.1.2, etc.

Deduction with Plantation Simulator
The key output of the Plantation Simulator is a detailed dynamic nutrient balance (Appendix A.5.3) of the Acacia stand. The balance contains N, P and K demand (Appendix A.5.1) and supply (Appendix A.5.2) as monthly values. Nutrient demand and supply are calculated independently; therefore, the results of growth and yield are conditional on nutrient supply: Following the Liebig's law of the minimum, growth can occur only if enough resources are available. Using the nutrient balance, the plantation managers can identify which nutrient is needed and when and can construct a fertilization regime to sustain the maximal growth.
Nutrient demand (Appendix A.5.1) is calculated from stand growth (Appendix A.1.1) and litter production (Appendix A.2). Net nutrient demand is derived from nutrient concentrations and growth of Acacia biomass components: foliage, bark, branch, fis, coarse roots, and stem wood ( Figure 3). Thereafter, the gross total nutrient demand is calculated as the sum of the net demand and the nutrients lost through monthly litterfall. The competing nutrient demand of weeds is calculated accordingly. Example of a phosphorus balance of an Acacia stand. Stand P storage (lower left) is computed from the biomass components (bark, branch, foliage, stem, fis, coarse roots) and their nutrient concentrations. Net demand (upper left) is the monthly change in nutrient storage. Gross demand (upper middle) accounts for net uptake and the uptake that replaces the nutrients lost in litterfall. Cumulative gross demand (lower middle) shows that the stand needs P at an approximate rate of 280 kg ha −1 during three rotations (250 kg for above-ground components and 30 kg for the below-ground components shown below the X-axis). The instantaneous nutrient supply (upper right) shows how much P becomes available from internal retranslocation, atmospheric deposition, fertilizer dissolution and decomposition of litter, coarse woody debris and peat. The cumulative nutrient supply is shown in the lower right panel. The dotted red line in the rightmost panels indicates the total nutrient demand of the stand.
The Plantation Simulator accounts for the main nutrient sources: Retranslocation from senescing litter, atmospheric deposition, nutrient release from decomposition of stand and weed litter, coarse woody debris and peat, and nutrient dissolution from fertilizers. Leguminous atmospheric N-fixation is also considered for Acacia. Each component in the demand and supply can be followed separately, thereby facilitating distinction of the most important processes and detailed component-wise consideration of uncertainty. Nutrient balance is calculated as both cumulative and instantaneous monthly values. The cumulative balance is used if we can assume that all released nutrients are stored in the root zone and are available for plant uptake; whereas the instantaneous balance is valid under the assumption that only nutrients released in the current time step are available for uptake. Cumulative and instantaneous balances represent extreme ends of the uptake scenarios; therefore, interpretation of the results should rely on both balances. The cumulative and instantaneous nutrient balances (Appendix A.5.3) are calculated on an area-based nutrient supply implying that each square meter of the compartment is equally relevant as a nutrient source for a tree. This holds well for a mature plantation with a closed canopy but is not valid for the initial plantation stage when tree seedlings have small root systems. Therefore, we also calculate a tree-wise nutrient balance, which takes account of the area from which a tree can extract nutrients. This is dependent on the canopy size ( Figure 4). After the canopy closure the tree-wise nutrient balance equals to the stand nutrient balance, but before the closure the tree-wise balance leads to nutrient deficiency more easily than the stand-wise nutrient balance. . Tree-wise nutrient balance for N, P and K. The tree-wise nutrient balance allows us to consider available nutrient reservoirs within the reach of the root system before canopy closure. Therefore, even though the area-based nutrient balance would suggest an adequate nutrient supply, the individual tree might suffer deficiency. The red area between the supply and demand curves indicates deficiency.

Model Testing
We used the data presented by Hoojier et al. [5] and Evans et al. [9] to test the performance of the Plantation Simulator. Table 3 in Hoojier et al. [5] presents the mean and standard deviation of WT and subsidence from 21 Acacia crassicarpa sites located on peatlands in Riau, Indonesia. The paper contains, without a specified location, the mean and standard deviation of the peat bulk density, and a record that the bulk density increases by 45 kg m −3 in 15 years. Evans et al. [9] present a uniquely large dataset of peat subsidence with respect to mean WT (see Figure 2 in [9]). However, this paper does not contain information about the peat bulk density; we therefore derived the value for initial bulk density from older plantation reported in [5].
Neither of the studies give explicit information about the tree stand nor its development. We simulated the Acacia stand dynamics using the growth and yield model for Acacia mangium [26,27], assuming a constant rate of mortality, and adjusted the dominant height development. Classically the rotation length is set to the age when the mean annual increment (MAI) and the current annual increment (CAI) of the stand intersect. We adjusted the dominant height model parameter values so that they give a height of 21 m (typical height of harvestable stand) at the age of 5 years and produce an intersection of MAI and CAI between the ages of 4 and 5 years. The growth and yield model and its parameterization are presented in detail in the Appendixes A.1.1 and A.1.2.
Nutrient demand and litter production depend on how trees allocate biomass to bark, branches, foliage, roots and stems as all these components have distinctively different nutrient contents. A full set of allometric equations for all these biomass components are not, to our knowledge, publicly available for Acacia crassicarpa. Therefore, we composed the allometric functions calculating the aboveground biomass function specially made for A. crassicarpa [28] and then assuming that biomass is distributed among the components as described by Krisnawati et al. [29] for A. mangium. A more detailed description of the procedure, the biomass equations and parameters are presented in Appendix A.2.
A Monte Carlo simulation set was conducted to fill the data gaps with plausible site and stand characteristics. We lifted parameter values from normal distribution with given mean and standard deviation (Table 1). For mortality rate we assumed a uniform distribution. We ran the model for 15 yrs, i.e., three full rotations. A monthly WT record was generated from the reported mean WT and assuming 20% standard deviation from the mean with annual dry/wet season cycle [5]. Because each simulation starts from planting, all the stands at a given time point are of equal age, unlike in the measured subsidence dataset. Therefore, a random 10 yrs subset of the simulation results was taken to ensure that the results are independent of the stand age, and thus comparable to the reported field subsidence data. The simulated subsidence accounts for the peat surface elevation change resulting from the simulated decrease in soil organic matter stock, long-term compaction, and short-term shrinking and swelling caused by WT fluctuation. A detailed description of the subsidence calculation is given in the Appendix A.3.4. The overall ecosystem C balance was calculated as the C storage change in soil and the change in stand above-ground biomass C, excluding the commercially usable stems (Appendix A.3.5). After demonstrated plausibility from the mass balance viewpoint we constructed the N, P and K balance for each simulation and identified the adequacy of nutrient supply with respect to the stand nutrient demand. The timing and quantity of nutrient deficiency was recorded, and the need for fertilization was identified. The effects of each of the input variables in the Monte Carlo simulation (site index, mortality rate, peat bulk density and mean WT) on the subsidence rate, P balance and C balance were assessed using frequency distribution plots.

Scenario Assessment
The parameterized and tested model was applied for two WT schemes which describe the effects of a change in WT from a mean depth of −0.80 m observed by Jauhiainen et al. [7] for Acacia plantation to a mean of −0.40 m, based on the Indonesian government regulation (SK.22/PPKL/PKG/PKL.0/7/2017). The nutritional consequences of raising the WT were identified. To cover different site and stand characteristics, we ran 200 simulations with varying site index, mortality rate and peat bulk density for low (mean WT −0.8 m) and high WT (mean WT −0.4 m) scenarios. Consequences for stand subsidence, C balance and fertilizer need were quantified.

Testing the Mass Balance Using Published Subsidence Data
The Plantation Simulator successfully predicted both the mean and the range of the peat subsidence in an Acacia plantation, based on standard plantation water and stand management practices ( Figure 5). The simulation set shows that the model satisfactorily covers the variation in the Hoojier et al. [5] dataset, which included detailed site-specific input data. However, the model missed the extremes in the Evans et al. [9] dataset, where considerably less input information was provided. In both cases, the mean subsidence was remarkably well predicted. In the simulation set the unknown site variables (site index, mortality rate and peat bulk density) were lifted from probability distributions, which propagates variation into the model target variables as described in Figure 6. Subsidence rate was correlated with WT and with initial peat bulk density. The good match between the measured and modelled subsidence provides a plausible basis for the nutrient balance calculation, because the nutrient dynamics in deep-peat sites depend on the same biogeochemical processes as the mass balance dynamics reflected by the subsidence.  . Effect of site index (m at index age of 5 years) (a,e,i), mortality rate (stems ha −1 month −1 ) (b,f,j), peat bulk density (kg m −1 ) (c,g,k) and mean WT (cm, positive down) (d,h,l) on annual subsidence rate (cm yr −1 ) (row 1), tree-wise P balance (kg ha −1 ) (row 2) and site C balance (kg ha −1 yr −1 , negative value indicating C source) (row 3). The figure was constructed from the Monte Carlo simulations using the Evans et al. [9] dataset.

Nutrient and C Balances
Stand and tree-wise N, P and K balances were extracted from the Monte Carlo simulations based on the Evans et al. [9] dataset ( Figure 7). The results indicate a large variation in the nutrient balances. With standard plantation management practices, N and K supply were mainly adequate for production of plantation biomass. In contrast, P supply is likely to be the growth limiting factor. On average, P demand and supply were equal but in half of the cases P balance was negative (Figure 7b,e). Site index, mortality rate and WT affected the P balance ( Figure 6 e,f,h). Increasing nutrient demand with increasing site index led to a more negative P balance. The positive correlation between P balance and mortality rate indicates decreasing competition for nutrient resources with increasing mortality. However high mortality opens gaps in the stand and allows nutrient uptake by ground vegetation, which can depress P balance in the highest mortality rates (Figure 6f). Higher WT tended to have a negative impact on P balance but a positive impact on C balance (Figure 6h,l). . N, P and K balance calculated based on standard plantation management practices, using data from the Evans et al. [9] dataset. The simulation period covers 15 years including three full rotations from planting to harvesting. The upper panels describe the instantaneous stand-wise N (a), P (b) and K (c) balance (difference between nutrient supply and demand), lower panels describe the cumulative tree-wise N (d), P (e) and K (f) balance. Negative values in all cases indicate that nutrient supply is smaller than the demand. The red area shows the range, and the solid black line shows the mean.

Scenario Assessment: Raising the Water Table
We applied the Plantation Simulator to identify the nutrient balance change caused by raising WT from a mean of −0.80 m to −0.40 m. Here we concentrate on P, because it was identified to be the main growth limiting nutrient. The cumulative stand-wise P balance remained positive, but decreased by 50 kg ha −1 . The tree-wise cumulative P balance changed from being marginally positive (5 kg ha −1 ) to negative (−17 kg ha −1 ), implying that the internal supply of P from peat mineralization and nutrient recycling would not be sufficient to achieve the rates of tree growth assumed in the simulation (Figure 8a,b). Under the −0.40 m scenario, it would therefore be necessary to close this nutritional gap, of an estimated10 to 17 kg P ha −1 rotation −1 , with increased fertilization. On the other hand, raising the WT reduced subsidence from 4.4 cm yr −1 to 3.3 cm yr −1 (Figure 8c), and reduced carbon loss from 17 to 9 Mg ha −1 yr −1 (Figure 8d).

Need for a Simulation Model
Rates of nutrient consumption are high in tropical short rotation tree plantations, and there is a risk that this could lead to the long-term depletion of soil nutrient resources [24]. This risk is exacerbated for pulpwood plantations that are typically located on infertile soils, and better understanding of the nutrient balance of these systems is therefore needed [30]. Complete nutrient balance studies for tropical plantations are scarce, with most studies dealing only with single nutrients, and focused on mineral soil sites [24,31]. Nutritional issues in plantations growing on peatlands are therefore poorly understood. Subsidence rates in drained tropical peatlands are high, typically 3-5 cm yr −1 [9]. Considering that plantation drainage depths are adjusted over time to maintain a constant WT relative to the peat surface, and that root layer thickness is constant, this continuously exposes "new" organic matter reservoirs to decomposition and nutrient release, but with a high environmental cost. Raising WT decreases peat oxidation and subsequently limits the organic matter decomposition and subsidence, but at the same time reduces the supply of nutrients from decomposing peat [22] which may lead to stand nutrient limitation.
Managing plantations on peat to balance this trade-off between stand productivity and environmental sustainability requires a holistic approach to nutrient management that considers both anthropogenic and natural processes, and the effects of plantation management on these. Therefore, nutrient management of a plantation should incorporate water management, fertilization, weeding, litter and slash management including debarking of harvested stems, and evaluation of nutrient export with harvested biomass (e.g., [30,32]). The Plantation Simulator incorporates these activities within an integrated biogeochemical model that makes visible the different nutrient sources, such as internal translocation of nutrients, circulation of nutrients from litter back to the stand, competition by weeds, and fertilization (see Figure 3). By linking multiple nutrient cycles the model can identify limiting nutrients for growth, helping to focus on nutrient management practices, as argued by Mendham and White [24], and to optimize fertilizer use in order to minimize the associated environmental and economic costs.
A further strength of the Plantation Simulator is that it provides a causal explanation for subsidence, extending empirical studies in which subsidence is predicted solely as a function of WT [5,7,9]. Therefore, it enables a search for optional management schemes connected to combinations of WT control, fertilization, slash management, planting density, and weeding that all affect the nutrient and mass balance and thus help to maximise stand productivity whilst minimizing subsidence and associated CO 2 emissions. Compared to other plantation models [23][24][25][26], Plantation Simulator succeeds to introduce critical characteristics of peatland plantation production into the simulation, such as P and K balance, the role of peat and other organic matter decomposition in the nutrient supply, and the specific the environmental problems, including subsidence and elevated C emissions.

Model Performance
In the test runs, the Plantation Simulator was able to replicate both the mean and the range of peat subsidence in a tropical pulpwood plantation. With the given range of inputs, the model produced subsidence rates ranging from 2.0 to 8.2 cm yr −1 . This covers the previously reported peat subsidence rates for tropical plantations [5,9,33,34]. According to experimental studies, the aboveground biomass in a mature A. crassicarpa stand can vary from 51 to 81 Mg ha −1 [35] or 63 ± 8.3 Mg ha −1 [28]. Using the Plantation Simulator and applying the mean site index (21 m) and mortality rate (9 stems ha −1 month −1 ), we obtained a mean above ground biomass of 66 Mg ha −1 at maturity, very close to the reported biomass values. The good fit between the observed and modelled subsidence and biomass accumulation suggest that the stand mass balance is plausibly described and allows us to proceed to an assessment of nutrient balances (Figure 1).
The Evans et al. [9] subsidence dataset represents stand mass balance under normal plantation management, where the mean WT was managed at −0.7 m during the monitoring period that extended up to 10 years. In the long run, site productivity reflects a balance between growth resource supply and demand; thus, the growth of biomass is scaled according to the supply of the minimum growth factor following Liebig's law of the minimum [25]. Figure 7 clearly suggests that N and K supply are adequate, while P is growth limiting. P shortage was evident during the first three years of the rotation. This is likely due to the competition with weeds, small root systems in the early development of the stand, and a delay between litter fall and nutrient release from the stand litter [24,36]. This finding fits well with the current plantation fertilization practice, which has demonstrated the benefits of P fertilization at time of planting [24,30]. Bich et al. [37] emphasise the importance of P fertilization in Acacia plantations as it enhances microbial activity, decomposition rate, and microbial N fixation.

Scenario Assessment
When WT is raised, the model shows that nutrient supply from the decomposing peat is considerably reduced, and the stand growth finds, without compensating fertilization, a new equilibrium at lower growth level. In our scenarios, the mean WT was raised from −0.80 m (mean in [7]) to −0.4 m, reflecting recent Indonesian government regulations. The gap between the scenarios in P balance varied from 10 kg ha −1 rotation −1 to 17 kg ha −1 rotation −1 (Figure 8a,b). According to Mendham et al. [32] a universal P dose of 10 kg ha −1 is recommended in the planting of A. mangium in South Sumatra. Assuming that a similar dose is applied to A. crassicarpa, this suggests that (at most) compensation of the raised WT would require a doubling of the current fertilization rate.

Uncertainty
Each module of the Plantation Simulator contains different levels of uncertainty. The growth and yield module (Appendix A.1.1) represents a traditional whole-stand model, where the stand growth is driven by H dom and diameter distribution change. Even though currently tree-level growth and yield models are preferred due to their flexibility, stand level models are still suitable for modelling singe-species and even-aged stands [38], such as pulp wood plantations. In this application we could not use the experimentally derived growth and yield parameters for A. crassicarpa because of the confidentiality of plantation company data. Instead we modified parameter values using those reported for A. mangium [26,27] with the exception that the H dom parameters (b 11 and b 12 , Table A1) were adjusted so that the production scheme became realistic in terms of rotation time. While we are unable to rigorously evaluate how well these parameters fit for A. crassicarpa in the absence of direct observations, we did find that the simulation as a whole produced plausible results. Omitting the parameterization problem, the growth and yield module per se contains rather low uncertainty, as the growth is easily verifiable using standard plantation inventory data. Thereafter, the relevant parameters affecting mortality rate, height, basal area and diameter development, can be easily modified accordingly. This is typically the first step in the application of the Plantation Simulator.
Biomass models (Appendix A.2) are usually based on much smaller datasets than growth and yield models, because measurement of biomass components is laborious [38]. Therefore, uncertainty of the biomass module may be slightly higher than that in the growth and yield module, although the allometric relationships within two stands of the same species tend to be rather stable if the stand density is similar [25]. This is particularly true for plantations with controlled genetic material and standardized stand management practices. It is clear, however, that a species-specific additive biomass model that includes bark, branch, foliage, roots and stems for A. crassicarpa would improve the accuracy of the nutrient balance estimates.
Decomposition is an essential part of nutrient cycling [30], and a major source of nutrients in drained deep peat sites. Decomposition of organic matter is a complex phenomenon, where not only original substances decompose, but also new substances are generated [39]; therefore, process-based decomposition models typically have rather a large uncertainty if site-specific input data are limited [40]. The Plantation Simulator considers separately decomposition and storage for CWD (Appendix A.3.1), stand litter (Appendix A.3.2) and peat (Appendix A.3.3). CWD model [41] applies empirical approach where the decomposition rate depends on mean air temperature, stem diameter and initial wood density. Stand litter is decomposed using the mechanistic Romul model [42]. An empirical carbon emission function determined for Acacia plantations in tropical peatlands was used to estimate the total decomposition including all organic matter components [7]. To allocate decomposition to CWD, stand litter and peat, we subtracted the calculated C emissions of CWD model and Romul model from the total C emission. The remaining emission represented the peat decomposition. This structure decreases uncertainty of the decomposition module as a whole, because it frames the total decomposition closely to the range observed by Jauhiainen et al. [7] for A. crassicarpa plantation on tropical peat and in the same area.
Nutrient demand module (Appendix A.5.1) contains less uncertainty than the nutrient supply module (Appendix A.5.2), because nutrient demand is closely related to stand development. A simulation with site index (21 m) and mortality rate (9 stems ha −1 month −1 ), shows that two thirds of the P uptake ends up in foliage and branches. This emphasizes the importance of leaving logging residues on site after harvesting [37]. In our application the litterfall rate is estimated based on the longevity of each biomass component, with parameters based on expert knowledge. The strength of this approach is that the litter production follows the growth and biomass component models and all components will eventually be transferred to litter. This narrows down the uncertainty range. Nutrient supply from the decomposing stand litter and from decomposing peat were the main sources of P supply contributing three quarters of the cumulative P supply during one rotation. The uncertainties in the decomposition models propagate to the nutrient supply estimates.
In this application C balance was calculated as a C storage change in soil and biomass, but omitting the contribution of stems, which are used for pulp production. Thus, the C balance includes biomass of weeds and unmerchantable parts of trees, and subsequently provides a C balance estimate coherent to our mass balance approach. The IPCC emission factors for tropical Acacia plantation on drained peatland are based on soil C storage change (−25 to −15 Mg CO 2 -C ha −1 yr −1 ), equivalent to 15 to 25 Mg CO 2 -C ha −1 yr −1 of emission [5,7,43,44]. In our estimates that account for the unmerchantable parts of trees, the balance varied between −20 and −10 Mg CO 2 -C ha −1 yr −1 (Figure 6i) indicating a slightly smaller C source than the IPCC values. Using the Plantation Simulator and applying the mean site index (21 m), mortality rate (9 stems ha −1 month −1 ), and mean WT of −0.8 m, the C balance accounting only for soil Equation (A32) was −19 Mg CO 2 -C ha −1 yr −1 , omitting the stems Equation (A33)-14Mg CO 2 -C ha −1 yr −1 and including the stems Equation (A34) −11 Mg CO 2 -C ha −1 yr −1 . The C balance estimates can be considered realistic, and therefore the Plantation Simulator can be a useful tool in producing more accurate net C emission estimates, because it takes into account the effect of stand properties, as well as stand and WT management. Future development of Plantation Simulator includes incorporation of hydrology, photosynthesis and net primary production modules to the model, and changing from stand-based to grid-based solution to facilitate simultaneous simulation of large areas.

Informed Consent Statement: Not applicable.
Data Availability Statement: Plantation simulator has been programmed using Python 3.7. The source code is openly available in https://github.com/annamarilauren (accessed on 15 February 2021). The repository also contains a parameter file and readme-file for instructions to running the model.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Appendix A. Plantation Simulator Description
Appendix A.1. Stand Appendix A.1.1. Stand Growth and Yield Growth and yield model, running in monthly time step, is the core of the Plantation Simulator. Growth and yield are computed using a classical dominant height-diameter distribution approach presented by Forss et al. [26,27] originally for upland Acacia mangium. The stand level variables including dominant height (H dom , m) and basal area (G, m 2 ha −1 ) are computed first. In this application we apply constant rate of mortality (M r , stems yr −1 ), and thus the stocking (N s , stem number) is given as: where N ini is the initial stocking (stems ha −1 ) and A is the stand age in years. The H dom (A) is given so that it produces anamorphic growth curves for different site index (S I ): where S I is the site index, which is determined as H dom at index age (A ind ), here 5 years; b 21 and b 22 are shape parameters for the height development (Table A1). The parameters are taken from Forss et al. [26,27] and shown in Table A1. Basal area (G, m 2 ha −1 ) is calculated using H dom and N s as: where b 30 , b 31 , b 32 , b 33 are parameters. The basal area weighted mean diameter (d wmean , cm) was computed from G and N s as follows: In the following step, a two-parameter Weibull distribution was used to describe the breast height diameter (d bh , cm) distribution for all time steps. The distribution parameters b W and c W were derived for each time step from the mean diameter and age as: where b 50 , b 51 , b 60 , b 61 are parameters. Number of stems in diameter class d (n(d), diameter in cm) is then given as: Thereafter, the individual tree height (h), volume (v) and merchantable volume (v m ) are computed for each diameter class in the distribution. The parameter values for the individual tree height model (b 100 and b 101 ) depend on d wmean , H dom and N as:

. Parameterization of H dom Model
A species specific H dom model for Acacia crassicarpa is not publicly available. Therefore, we deduced the shape parameters b 21 and b 22 for Equation (A2) using the following nonconfidential facts: H dom at harvesting time is 17. . . 25 m and rotation time from planting to harvesting is 5 yrs. Optimal rotation time is at age when the mean annual increment (MAI) and current annual increment (CAI) intersect. First, we set the site index to 21 m and M r to 9 stems ha −1 month −1 , computed the growth and yield model and plotted MAI and CAI. Then we iteratively adjusted b 21 and b 22 until MAI and CAI intersected between age 4 and 5 yrs ( Figure A1). Here, MAI in the intersection point equals to 33 m 3 ha −1 yr −1 , which fits to the range reported for A. mangium [24]. Figure A1. Growth and yield model parameterized with Acacia mangium as described by Forss et al. [26,27], and adjusting the H dom model Equation (A2) to produce the optimal rotation time to 5 years to estimate the Acacia crassicarpa production.

. Biomass and Litter
Biomass, its development and distribution among different tree components (bark, branch, foliage, and stem wood) is essential information when nutrient balance is calculated. Full set of allometric equations for biomass components were not available for Acacia crassicarpa. Bi et al. [45] show additive biomass functions, where biomass depends on tree diameter, for several Eucalyptus and Acacia species: where M i is dry biomass (kg) in biomass component i (foliage, branch, and bark), β i0 , β i1 are biomass parameters (see Table A2). The total stand biomass in the component i is computed by integrating over the diameter distribution Equation (A7). Chen et al. [28] present an equation for aboveground biomass for A. crassicarpa, but the parameter values are not directly applicable, because the biomass is dependent on both tree diameter and height.
To obtain β i0 , β i1 , we first calculated the total aboveground biomass using the equation presented by Chen et al. [28]; and plugging in the tree diameter and tree height calculated using Equation (A10). The relative proportion of biomass component i from the total aboveground biomass were obtained from allometric equations for Acacia mangium [29]. Then we calculated all biomass components i for a diameter range 1 cm to 20 cm; and then fitted parameters β i0 and β i1 using minimum sum of squares -method. Parameter values are presented in Table A2. Stem biomass was obtained by multiplying the stand volume with the wood density (500 kg m −3 ). Root biomass was assumed to be 20% of the above-ground biomass. Fine biomass was estimated to be 5% of the total root mass. Weeds and weeding play an important role in stand carbon and nutrient balance. Plantation simulator contains Weed module, which in a simple way describes development of weed biomass, litter production and interaction with stand as a function of time. The above ground weed biomass follows S-shape curve Equation (A2) and is parameterized according to approximate field observations. Below-ground weed biomass is set proportional to the above ground biomass (below/above ground ratio was set to 0.5).
where M weeds is the dry above ground biomass of weeds (kg ha −1 ), b 130 and b 131 are parameters (Table A2), and A weeds is the weed ages, i.e., time elapsed from the previous weeding. Here, the parameters values for b 130 and b 130 were 6000 and 1, respectively. This indicates that the weed biomass would approach asymptotically 6000 kg ha −1 , and that in A weeds of 1 year, M weeds would be 2200 kg ha −1 . The weed below-ground biomass was assumed to be 50% of the above-ground biomass. Shading of the canopy reduces weed mass. This relationship is described with Maximum green mass -parameter (gm max , value here 8000 kg ha −1 ), which sets the upper limit for the sum of the stand leaf mass and the above-ground weed mass. When the total green mass would exceed the gm max , the above ground weed mass is reduced accordingly. Weeding kills all the above and below ground weed biomass and transfers it to litter. Thereafter the development of weeds starts from the beginning. Equation (A12) describes the net biomass in each biomass component i (M i ) at given time t. To obtain the amount of new biomass constructed in each time step we used the M i and the longevity (life span, l i , years) information. If, say foliage biomass, in time step t is 5000 kg ha −1 , and the longevity of foliage is 6 months (here 6 time steps), the current M f oliage must be a sum of the foliage growth in current time step and the five time steps before. Arranging all the time steps into a linear system of equations we can solve the biomass growth (g i , kg ha −1 time step −1 ) as: 0 0 0 · · · 1 0 0 0 0 0 · · · 1 1 0 0 0 0 · · · 1 1 1 where the right hand side of the equation is the time series of biomass component i (vector M with length equal to number of time steps (tn), Equation (A14)), the second term in the left hand side of the equation is the unknown biomass growth time series (denote vector g, length tn), and the first term is a matrix (A, dimensions tn by tn) describing the time steps contributing to the M i (t). The diagonal denotes the current time step, and the left-hand elements from the diagonal are filled with ones equal to number of time steps fitting inside the longevity value. Following the foliage example, each row contains 6 ones. Now, the biomass growth can be solved as: Biomass exceeding the longevity age is located to litter. Litter production in biomass component i from the living biomass (L i l , kg ha −1 dt −1 ) is obtained from the biomass growth (g i , kg ha −1 ) and the longevity (life span, l i years) for each time step (t) as follows: Before the litterfall, the senescing tissues transfer a proportion of the nutrients into the rest of the plant in retranslocation. Only the litter from the living biomass is subject to restranslocation. Therefore, the nutrient J (J = N, P, K) content in the living litter component i l in plant p (Acacia, weeds) is: where J i l is nutrient J mass in biomass component i from living biomass (kg ha −1 ), L i l is litter mass in component i (kg ha −1 ), R J p is a proportion of retranslocation for nutrient J for plant p, and J i is nutrient J concentration in the litter component i (%). The parameter values are available in Table A2.
where J i d is the nutrient mass in biomass component i from dead biomass (kg ha −1 ). Stems are transferred to coarse woody debris pool (L CWD , kg ha −1 ). When weeding is done, both above and below-ground weed mass is transferred to the respective litter pools without retranslocation.

Appendix A.3. Soil and Peat
Organic matter is decomposed in three different procedures. A simple exponential decay model is applied for coarse woody debris as described by MacKenzen et al. [41]. The decomposition rate depends on temperature, stem diameter and initial wood density. The litter originating from living or dead biomass and produced by the stand and weeds is assumed to decompose under aerobic conditions. Widely used Romul-decomposition model [42] was applied to compute for the carbon and nutrient dynamics of the litter. An empirical carbon emission function, determined for Acacia plantations in tropical peatlands, was used to estimate the total heterotrophic C emissions from soil [7]. Peat decomposion was estimated by subtracting the C emissions produced by the litter decomposition from the total emissions.
General inputs to the decomposition module are soil temperature (T soil ) and water table (WT). In tropics the fluctuation on air temperature is small, and T soil is more dependent on shading than meteorological conditions. Therefore, T soil in open area (T open ) and under a closed canopy (T close ) were given as parameters. Monthly T soil was obtained from the above-ground green biomass as: where M f oliage is Acacia stand foliar mass, M weed above is the above ground weed biomass and gm max is the maximum green mass in the stand (see Table A2). In this study, we did not apply hydrological simulation (because lack of data), but instead we gave the monthly WT as an exogenic input. Monthly WT values were lifted from a normal distribution with given mean and standard deviation.

. Coarse Woody Debris (CWD) Decomposition Module
CWD decomposition is accounted by using single-exponential model presented by MacKenzen et al. [41]. The rate of decomposition (k CWD ) is dependent on mean annual temperature, diameter and density of stem: where T a mean is the mean air temperature (°C), d is the stem breast height diameter (cm) and ρ wood is the density of wood (kg m −3 ) (Table A3). The monthly input cohort of CWD (L CWD , kg ha −1 ) comes from the mortality computed in growth and yield module. The CWD module accounts for the current mass, and the storages and release of C, N, P and K from the CWD. The mass of the remaining wood on the peat (P wood cohort , kg ha −1 ) in a CWD cohort produced in time step t for the remaining time steps is computed as: where [t] is a time vector (in years) from the current time step until the end of simulation. The total amount of CWD is obtained as a sum of all P wood cohort . The decayed wood in the cohort (D wood cohort ) was computed as: Conversion factor from CO 2 to C 0.2727 ctomass Conversion farcor from C to dry biomass 2 The time derivative of D wood cohort represents the decomposition in the time step, and the C emission for the time step (C e f f lux CWD kg ha −1 ) is obtained by summing over all existing CWD cohorts: where c C C WD is the carbon content of CWD, n cohorts is the number of CWD cohorts present in the stand. Above ground CWD and coarse roots are computed using the same algorithm. N, P and K release from CWD were assessed as a function of the decomposed CWD following the graphs presented by Palviainen and Finér [46]. Interpolation functions between the proportion of the remaining mass and nutrients were constructed from the following arrays: where f N(M), f P(M), f K(M) are interpolation functions. The remaining mass of N, P and K in each CWD cohort were obtained by applying Equation (A20) together with the interpolation functions, and the total nutrient storage in CWD (N CWD , P CWD , K CWD kg ha −1 ) was obtained by summing over all the cohorts. Nutrient release from each CWD cohort was computed as a time derivative of the remaining nutrient variable. The total nutrient release from CWD was gain by summing the cohort-wise nutrient release (N rel CWD , P rel CWD , K rel CWD , kg ha −1 time step −1 ).

Appendix A.3.2. Stand Litter Decomposition, Romul Module
Acacia stand and weeds produce litter on and under the soil surface. Because the plantations locate on drained peatlands it is reasonable to assume that the stand litter is decomposed mainly under oxic conditions. Widely used decomposition model Romul was applied to compute the organic matter, carbon and nutrient storage of these plant remnants, and the carbon and nutrient release from them [42,47]. Romul describes the decomposition of mass and release of nutrients as succession from fresh litter (JL) through partly decomposed JF storage into stabile humus (JH) ( Figure A2). The rate of succession and mineralization (k1. . . k6) depends on ambient temperature and moisture, nitrogen, ash and lignin content of the remnant and soil pH. Variables, computation of the rate parameters and nitrogen processes are described in details by Chertov et al. [42]. Temperature and soil moisture are given as a time-step-wise array input to the model. Soil pH is given as a parameter (Table A3).
Input to Romul is given separately for litter originated from living and dead Acacia and weeds and in foliage, branch, bark and fis components ( Figure A2). The input variables include litter mass (L i l , L i d , see Table A2), and amounts of N, P and K (J i l , J i d , see Table A2) in the litter. The nutrient contents of litter originating from living biomass are subject to nutrient retranslocation before the litter fall, and therefore they have smaller nutrient content than the litter originating from dead biomass (mortality or weeding). The JL decomposition is computed separately for each litter cohort, but all above ground litter are lumped together for the computation of JF and JH materials. Similarly, all the below ground litter are lumped for the computation of the below ground JF and JH materials. This module accounts for the current mass of organic material above and below ground, and the current carbon efflux originating from the Acacia and weed litter remnants.
The original version of Romul accounts for decomposition of mass and release of N. Here we assume that the release of P and K follow the dynamics of mass loss. Soil organic matter mass components, nutrient storage, CO 2 efflux and nutrient release are the main output of the Romul module (Table A3). Figure A2. Flow chart of Romul-decomposition model as described by Chertov et al. [42]. The aboveand below-ground litter fractions and their nutrient contents are imported from the biomass and litter module. Romul describes the decomposition of mass and release of nutrients as a succession from fresh litter (JL) through partly decomposed material (JF) to stabile humus (JH). The transfer rates between the pools (k2, k4 and k5) and carbon release to the atmosphere (k1, k2 and k6) are controlled by ambient temperature nitrogen, ash and lignin content of the remnant and soil pH.

Appendix A.3.3. Peat Decomposition Module
Total heterotrophic CO 2 efflux is computed as a function of water table depth using empirical equation presented by Jauhiainen et al. [7]. The measurements behind the study were conducted in April Asia peatland plantation under Acacia crassicarpa stand of varying age. Chamber technique was applied in the study. The total CO 2 efflux (CO 2 e f f lux peat , kg ha −1 time step −1 ) through a unit area of peat is calculated as [7]: CO 2 e f f lux tot = (71.1 * WT + 23.15) * 1000 * pdt * aQ 10 , where WT is water table depth (m, positive downwards), pdt is time step in years (30/365), and aQ 10 is temperature adjustment: where Q 10 is a temperature sensitivity factor, T soil is soil temperature and T re f is reference temperature (Table A3). The model in Jauhiainen et al. [7] represents the total heterotrophic CO 2 efflux through the soil surface, and does not differentiate the decomposition of peat from the decomposition of the Acacia remnants or coarse woody debris. In Plantation Simulator these components are kept apart because the CWD, Acacia and weed litters are the main C inputs to soil and therefore slow down the total C loss from soil. Thus, the peat decomposition is obtained by subtracting the heterotrophic CO 2 efflux computed for the Acacia remnants and CWD from the total CO 2 efflux.
CO 2 e f f lux peat = CO 2 e f f lux tot − CO 2 e f f lux stand litter − CO 2 e f f lux CWD , (A29) CO 2 e f f lux peat was converted to C release and then to organic matter decomposition. Thereafter, the release of nutrient J (J = N,P,K) was calculated stoichiometrically using the organic matter decomposition and the peat nutrient concentrations (Table A4).
J release peat = CO 2 e f f lux peat * co2toc * ctomass * J peat 100 , where J release peat is the N,P,K release in peat decomposition (kg ha −1 time step −1 ), co2toc is conversion factor from CO 2 to carbon (12/44) and ctomass is conversion factor from carbon to mass (Table A4).
from fertilizers. Atmospheric deposition is accounted as a constant input to the system according to parameter values N deposition , P deposition , and K deposition (kg ha −1 year −1 , see Table A4). When simulating plantations situated on thick peat deposits, the mineral weathering can be omitted. In Acacia stand the microbial fixing of N can be important component of nutrition, because Acacia belongs to N fixing species. In optimal conditions, N fixing can provide up to 60 % of the stand N demand. However, N fixation is inhibited by low pH in the root layer, and pH in tropical peat can be as low as 3.5. The share of microbial fixing is given as a parameter in the simulation (Table A4).
where J microbial f ix is the nutrient supply through microbial fixing (kg ha −1 time step −1 ), p J f ix is the share of microbial fixing from the nutrient demand (for N given as parameter, Table A4, and 0 for P and K). Nutrient retranslocation (J retransi , kg ha −1 time step −1 ) for nutrient J (J = N, P, K) from the litter component il in plant p (Acacia, weeds) is: where J i is nutrient J mass in biomass component i (kg ha −1 ), L i l is litter mass in component i (kg ha −1 ), R J p is a proportion of retranslocation for nutrient J for plant p, and J i is nutrient J concentration in the litter component i (%). The parameter values are available in Table A2. Nutrient release from decomposition of organic matter was derived as a sum of nutrient J release from CWD (J release C WD ), stand and weed litter (J release l itter ), and peat (J release p eat ). In Plantation simulator the information about fertilization is given as mass of fertilizers given per tree at any time step. The cumulative amount of nutrient J released from fertilizer (J release f ert , kg ha −1 ) in soil is calculated using the following exponential decay function: where F J nut cont is the fertilizer nutrient content (in N, P 2 O 5 and K 2 O), dose is the amount of fertilizer given to a single tree (g tree −1 ), conv is a conversion factor from P 2 O 5 to P (0.42), and K 2 O to K (0.83), N s is standing stocking (stems ha −1 ), k f ert is the decay rate of nutrient J in the fertilizer F, and [t] is a time array. Release of nutrients in time step t was obtained as a time derivative of Equation (A39). Total nutrient supply (J tot s upply kg ha −1 time step −1 ) in time step t was obtained as a sum of all the nutrient supply components.
J tot supply = J deposition + J microbial f ix + ∑ i J retrans i + J rel CWD + J release litter + J release peat d Appendix A.5.3. Nutrient Balance Nutrient balance constructed from it's components enables us to study the relative importance of each nutrient source and to identify where the nutrients are most needed. Area-based nutrient balance assumes that all nutrients are equally available for the stand no matter where the resources are located. This is most likely valid when the canopy is closed. Time-step-wise nutrient balance (J balance , kg ha −1 time step −1 ) compares the nutrient supply and demand from the same time step as: J balance = J tot supply − J gross demand t .
This balance implicitly assumes that only the nutrients released at the same time step are available as a growth resource, and no nutrient storage in soil takes place. Cumulative nutrient balance, achieved by integrating Equation (A41) over the time, assumes that all nutrients that were not used by the stand are available as a growth resource for the future time steps, and thus no nutrient loss takes place. The reality is somewhere between these approaches, therefore it is useful to study both of the balances at the same time.
In the early development of the stand tr systems are small and only the resources close to the tree are available. Tree-wise balances balance was constructed by weighing the nutrient supply with tree coverage. The tree coverage was obtained as: Thus, the tree-wise nutrient supply was: J tree supply = p tree coverage J deposition + J rel CWD + J release litter + J release peat + J microbial f ix + ∑ i J retrans i + d dt J release f ert , (A43) Tree-wise nutrient demand (J tree demand , kg ha −1 time step −1 ) was constructed using Equation (A43) but omitting the weed components. The tree-wise nutrient balance was obtained by subtracting J tree s upply from J tree d emand , and thereafter integrated over the time to obtain the cumulative tree-wise nutrient balance. Figure A3. Left upper panel: Cumulative litter input to soil from living biomass, from mortality and CWD, left lower panel: C input to soil from living biomass (green area), from mortality (light grey area) and CWD (dark grey area). Upper middle: Organic matter storage in soil including stand residue L, F, LH, H, SH component computed by Romul-submodel (colored areas), CWD (light grey area) and peat (dark grey area); and subsidence as affected by soil organic matter change refer to right y-axis, compaction and shrinking/swelling. Lower middle: cumulative C emissions from stand litter, CWD and peat. Lower right: C storage in stand and weed biomass, in stand total biomass and in stand biomass excluding stems described above the x-axis. Lowermost line below the x-axis describes the change in soil C storage, and the next dotted line up denotes C balance.