A sub-canopy structure for simulating oil palm in the Community Land Model

Introduction Conclusions References


Introduction
Land-use changes in South-East Asia's tropical regions have been accelerated by economy-driven expansion of oil palm (Elaeis guineensis) plantations since the 1990s 5 (Miettinen et al., 2011). Oil palm is currently one of the most rapidly expanding crops in the world (Carrasco et al., 2014) and Indonesia as the largest global palm-oil producer is planning to double its oil-palm area from 9.7 million ha in 2009 to 18 million ha by 2020 (Koh and Ghazoul, 2010). Since oil palms favor a tropical-humid climate with consistently high temperatures and humidity, the plantation expansion has converted 10 large areas of rainforest in Indonesia in the past two decades replacing mainly intact forests (47 %), logged forests (22 %), and agroforests (21 %) and also happening on carbon-rich peat soils (∼ 10 %) (Carlson et al., 2012).
Tropical deforestation caused by the expansion of oil palm plantations has significant implications on above-and below-ground carbon stocks and on land-atmosphere bio-15 physical and biogeochemical processes (Nogueira et al., 2014). Undisturbed forests have a significantly higher capacity to store carbon than disturbed or managed vegetation (Luyssaert et al., 2008). The greenhouse gas balance of oil palms is still uncertain due to incomplete monitoring of the development of young oil palm plantations, and lack of understanding of the carbon and nitrogen exchange between oil palm, soil and the at-Introduction Oil palm is a perennial evergreen crop following the Corner's architectural model (Hallé et al., 1978), each phytomer carrying a large leaf (or frond) axillating a fruit bunch. A number of phytomers emerge successively (nearly two per month) from a single meristem (the bud) at the top of a solitary stem. These phytomers develop in series throughout their life cycles, being on top of the crown at emergence, then, progressively 5 being covered by new phytomers, until senescence. Each has its own phenological stage and yield, according to respective position in the crown. Carbon (C) and nitrogen (N) allocation also has to be determined for each phytomer specifically according to its phenology. The corresponding gross primary production (GPP) by leaves and carbon output by fruit harvests both exhibit seasonal dynamics in response to environmental 10 drivers and agricultural managements. In order to capture the inter-and intra-annual dynamics of land-atmosphere fluxes in the oil palm system, a new structure and dimension detailing the sub-canopy phenological, C and N allocation and agricultural management processes has to be added to the current integrated plant-level biogeochemical parameterization in the land models. This specific refinement needs to remain 15 compliant with the current model structure though, and be simple to parameterize.
In this study, we develop a sub-canopy phytomer-based structure for simulating an "oil palm plantation" land use within the framework of the Community Land Model (CLM4.5). The new phytomer-based structure is designed for a special oil palm PFT, or even for oil-palm-like plantations (e.g. coconut, date palm etc.). It introduces a sub-20 canopy phenological and physiological parameterization, so that multiple leaf and fruit components operate in parallel but at delayed steps, separated by a thermal period. More specifically, at the PFT level vegetative growth is represented by leaf development of different age cohorts (40 ranks) of phytomers as well as the development of shared stem and root. Reproductive growth of the PFT consists of simultaneous fruit-25 filling at different age cohorts of bunches, and yield includes multiple harvests per year of the mature bunches at the lower end of the canopy. The phenology is driven by thermal time and available C and N are allocated accordingly to the vegetative and reproductive components for the whole plant and then partitioned between phytomers. The functions implemented for oil palm combine the characteristics of both trees and crops, such as the woody-like stem growth and turnover but the crop-like vegetative and reproductive allocations which enable fruit C and N output. Agricultural practices such as transplanting, fertilization, and leaf pruning are also represented.
Overall, this new parameterization on oil palm aims to allow future investigation of the 5 effects of forest (forest PFTs already available in CLM) to oil palm conversion on carbon, water and energy fluxes between land and atmosphere across regional to global scales.
The main objectives of this paper are to: (i) describe the development of the new subcanopy structure for oil palm in the CLM including its phenology, carbon and nitrogen 10 allocation, and yield output, (ii) optimize model parameters using field-measured leaf area index (LAI) and observed long-term monthly yield data from a mature oil palm plantation in Sumatra, Indonesia; and (iii) validate the model against independent data from eight oil palm plantations of different age in Sumatra, Indonesia. 15 The crop model component of CLM4.5 implements the interactive crop management parameterizations from the land surface model AgroIBIS (Kucharik and Brye, 2003;Levis et al., 2012;Drewniak et al., 2013). It has explicit representations of crop type, planting, harvesting, fertilization, and irrigation. The managed crop types currently simulated and validated in CLM4. 5  Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | emergence and declines continuously until zero when grain-fill begins. Since grain-fill, the LAI declines drastically because of background litter-fall controlled by leaf longevity. Furthermore, the entire yield represented by grain C and N is currently routed into the litter pool in CLM4.5.

Model development
For the adequate description of oil palm functioning, these simplifications are not 5 sufficient and the structure of the CLM crop model should be adapted. First, oil palm is a perennial evergreen crop and thus the annual phenological cycle has to be extended to perennial. Second, allocation functions have to be modified accordingly for an evergreen pattern. Third and most important, a multiple growth and yield scheme that can respond to seasonal climatic dynamics has to be coined for oil palm's continuous 10 fruit-filling and harvests throughout a year. Lastly, it would be more logical to output the harvested fruit C into a separate export C pool instead of routing it to litter. This is necessary to enhance the decomposition computation and allows the validation against observed yield data and soil carbon data. A sub-canopy structure can meet most of the above needs by dividing the canopy 15 into multiple layers or multiple "big leaves", each having its own phenology, growth and yield (Fig. 1a). An oil palm produces successively 20-30 phytomers every year according to phyllochron (the period in heat unit between initiations of two subsequent phytomers) (Fig. 1b). Each phytomer carries a large leaf (frond) and an axillated inflorescence bunch. A maximum of 40 expanded leaves, each growing up to 7 m long, 20 are usually maintained in plantations by pruning management, eliminating senescent leaves at the bottom of the crown. There are also around 60 initiated phytomers developing slowly inside the bud, all being heterotrophic. The eldest ones, emerging at the top of the crown and unexpanded yet, are named "spear" leaves. Each phytomer can be considered a sub-PFT component that has its own prognostic leaf growth and 25 fruit yield capacity like a PFT but having (1) the stem and root components that are shared by all phytomers, (2) the soil water content, nitrogen resources, and resulting photosynthetic assimilates that are also shared and partitioned among all phytomers, and (3) a vertical structure of the foliage, with the youngest at the top and the oldest at the bottom of the canopy. Within a phytomer the fruit and leaf components do not compete for growth allocation because leaf growth usually finishes well before fruitfill starts. However one phytomer could impact the other ones through competition for assimilates, which is controlled by the C and N allocation module according to their respective phenological stages.

5
The phytomer configuration is similar to the one already implemented in other oil palm growth and yield models such as the APSIM-Oil Palm model (Huth et al., 2014) or the ECOPALM yield prediction model (Combres et al., 2013). However, the implementation of this sub-canopy structure within the CLM modeling framework is the first attempt among land surface models. It incorporates the ability of yield prediction, like 10 an agricultural model, beside that it allows the modeling of both biophysical and biogeochemical processes as a land model should do, e.g. what is the whole fate of carbon in plant, soil and atmosphere if land surface composition changes from a natural system to the managed oil palm system?
In this model description, we focus on the new phenology, allocation and agricultural 15 management functions developed specifically for the oil palm PFT. Photosynthesis, respiration, water and nitrogen cycles and other biophysical processes are the important functions already implemented in CLM4.5 and they are not modified (except N retranslocation scheme) for the current study. Details about these functions can be found in Oleson et al. (2013). The following diagram shows the major model developments 20 in this paper and their coupling with existing modules within the CLM4.5 framework (Fig. 2).

Phenology
In CLM, three types of phenologies (evergreen, seasonal-deciduous and stressdeciduous) are available for natural PFTs. These phenology algorithms provide prog-25 nostic controls on the seasonal timing of new leaf growth (onset) as well as leaf litter-fall (offset) and root and stem turnovers. sufficient to consider the 3-phase phenological cycle for the plant whole life. Leaves are assumed to grow until grain-fill and harvest implies removal of the reproductive part as well as of all the vegetative components incl. leaf, stem and root. The oil palm, as a perennial evergreen crop, is fundamentally distinct from the natural vegetation and other types of crops. It has to sustain a high level of LAI while producing fruits continuously (nearly two fruit bunches per month) for many years, including a juvenile stage. In order to simulate such phenology, the vegetative growth must be partially decoupled from the reproductive growth. For example, leaf growth and maturity could be finished well before fruit-fill and leaf senescence could happen after a harvest event. Moreover, each phytomer has its own phenology that controls its leaf and fruit development (Fig. 1).
We designed six post-planting phases for the phytomer phenology: (1) from leaf initiation to leaf expansion, (2) from leaf expansion to leaf maturity, (3) from leaf maturity to fruit-fill, (4) from fruit-fill to harvest, (5) from harvest to leaf senescence; and (6) from leaf senescence to pruning (Fig. 1b). The modified phenology module controls the life 15 cycle of each phytomer as well as the stem and root turnover for the whole plant. Different phytomers are assumed to follow the same phenological steps and the same thermal length for each phase. Air temperature is the key variable or "clock" for the phytomer phenology as measured by GDD. For tropical crops such as oil palm, a new GDD variable with 15 • C base temperature and 25 • -days daily maximum (Corley and 20 Tinker, 2003;Goh, 2000;Hormaza et al., 2012) is accumulated since planting (abbr. GDD 15 ). The six phenological phases are signaled by respective GDD requirements, except that pruning is controlled by the maximum number of live phytomers according to plantation management (

Planting and leaf initiation
Planting is implemented in the similar way as in the CLM4.5 crop phenology except that GDD 15 is tracked since planting and an option of transplanting is enabled. An initial phytomer emergence threshold (GDD init ) is prescribed for attaining the first leaf initiation after planting (Table 1). When GDD init is zero, it implies transplanting from nursery in-5 stead of seed sowing in the field. Oil palm seedlings usually grow in nursery for 1-2 year before being transplanted into the field. Therefore, in this study GDD init is set to zero and the first new phytomer is assumed to initiate immediately after transplanting in the field. An initial total LAI of 0.15 is assigned to the existing expanded phytomers, whose leaf sizes are restricted to be within 10 % of PLAI max (Table 2). 10 The oil palm phytomers initiate as leaf primordia in the apical bud and then appear as leaves on the stem successively according to relatively stable intervening periods, termed plastochron (the duration in terms of heat unit (GDD) between successive leaf initiation events) and phyllochron (the rate of leaf emergence from the apical bud). Here for simplicity, the phyllochron is assumed equal to the plastochron. As the apical buds 15 in palms usually do not start to accumulate dry mass immediately after physiological initiation but wait until several phyllochrons before expansion (Navarro et al., 2008), we define leaf initiation as the start of active accumulation of leaf C in this model, so that the phenological steps and C and N allocation process can be at the same pace.
A parameter phyllochron is prescribed with an initial value of 130 • -days at planting 20 with reference to GDD 15 and it increases linearly to 1.5 times at 10 year old (Huth et al., 2014;von Uexküll et al., 2003). Given GDD init and phyllochron, a heat unit index H init p for triggering leaf initiation can be calculated for each new phytomer when a preceding phytomer initiates: where subscripts p and p+1 refer to successive phytomers and 1 refers to the first new phytomer initiated after planting. As the GDD accumulates since planting, new phytomers will be turned on in sequence when GDD 15 > H init p , and will enter the 6-phase life cycle one by one. The timing of later phenological steps for each new phytomer is determined at the time of 5 initiation by adding the length of a corresponding phase period (Table 1). Each newly initiated phytomer is assigned a negative rank of −N and remains packed in the bud until the next phase of leaf expansion is triggered. The oldest unexpanded phytomer (spear leaf), right before expansion, has a rank of −1. The GDD period between leaf initiation and expansion is used to calculate the number of bud phytomers that have 10 already initiated before transplanting, i.e. N = GDD exp phyllochron .

Leaf expansion
During the phase from initiation to leaf expansion, leaf C already starts to build-up in the bud or spear leaf but it remains photosynthetically inactive. The thermal threshold for leaf expansion is calculated by H a phytomer ranked −1, the leaf starts to expand and becomes photosynthetically active. Its rank changes to a positive value of 1, while the ranks of other phytomers all increase by 1 at the same time. The expansion phase lasts for roughly 5-6 phyllochrons until leaf maturity (Legros et al., 2009).
Hereafter, the pre-expansion and post-expansion growth periods, distinguished 20 by negative and positive ranks, are treated separately so as to differentiate nonphotosynthetic and photosynthetic increases in leaf C. The following post-expansion phases and their thresholds are determined with reference to H exp p .

Leaf maturity
Another phenological step is added for the timing of leaf maturing so as to control Introduction reaches maturity well before fruit-fill starts on the same phytomer. Therefore, we set the parameter GDD L.mat to be smaller than GDD F.fill ( + GDD F.fill . At this point, the phytomer enters reproductive growth. Growth allocation increases gradually for the fruit component while leaf C and LAI remain constant on the mature phytomer until senescence. Due to the fact that most inflorescences on the initial phytomers within 2 years after planting are male (Corley and Tinker, 2003), another threshold GDD min is used to control the beginning of first fruiting on the palm. Only when GDD 15 > GDD min , the mature phytomers are allowed to start fruit-filling.

Fruit harvest and output
Fruit harvest occurs at one time step when a phytomer reaches fruit maturity, measured 15 by a heat unit index H F.mat p = H exp p +GDD F.mat . Since GDD build-up is weather dependent and phyllochron increases through aging, the harvest interval is not constant. New variables track the flow of fruit C and N harvested from each phytomer to PFT-level crop yield output pools. The fruit C and N outputs are isolated and are not involved in any further processes such as respiration and decomposition, although their fate is 20 largely uncertain.

Litter fall
For oil palm, leaf litter-fall is performed in two phases: senescence and pruning. Senescence is simulated as a gradual reduction in photosynthetic leaf C and N on the bottom exp p + GDD L.sen . These phytomers are allowed to stay on the palm until pruning is triggered. Their senescence rates are calculated as the inverse of the remaining time until the end of a phytomer's life cycle (GDD end ). Leaf C removed during this phase is not put into the litter pool immediately but saved in a temporary pool C senescent leaf until pruning, while the photosynthetic LAI of 5 senescent phytomers are updated at every time step. The reason to do this is that each oil palm frond is a big leaf attached tightly to the stem and its leaflets do not fall to the ground during senescence unless the whole frond is pruned. Thus, the dynamics of soil litter pool and decomposition process could be represented better with this function. Nitrogen from senescent phytomers is remobilized to a separate N retranslocation pool that contributes to photosynthetic N demand of other phytomers and avoids supplying excessive amount of N to the litter. The proportion of N remobilized from senescent leaves before pruning is adjusted by the length of senescent period (GDD end − GDD L.sen ) with a given pruning frequency, and the rest N goes to the litter pool.
Pruning is conducted at one time step if the number of expanded phytomers (incl. senescent ones) exceeds the maximum number allowed (i.e. mxlivenp). All senescent phytomers are subject to pruning at the time of harvest and their remaining C and N together with the temporary C senescent leaf pool are moved to the litter pool immediately. The frequency and intensity of pruning is determined through the combination of 20 mxlivenp, GDD L.sen and phyllochron. A larger mxlivenp gives lower pruning frequency and a smaller GDD L.sen results in more senescent leaves being pruned at one time. Besides, since phyllochron increases by age, the rate of phytomer emergence decreases and thus pruning frequency also decreases when the plantation becomes older. 25 Unlike other crops, the oil palm stem is represented by two separate pools for live and dead stem tissues (Fig. 1a). Although the stem of oil palm is not truly woody, field observations have found that the stem section below the lowest phytomer only contains less than 6 % of live tissues in the core of trunk for transporting assimilates to the roots ( van Kraalingen et al., 1989). This is similar to the stem of most woody trees that largely consists of functionally dead lignified xylem. Therefore, conversion from live to dead stem for oil palm follows the CLM stem turnover function for trees, except that 5 the turnover rate is slightly adjusted to be the inverse of leaf longevity (in seconds), such that when a leaf is dead the stem section below it will mostly become dead. Leaf longevity is around 1.6 years measured from leaf expansion to the end of senescence. The oil palm fine-root turnover follows the CLM scheme for trees and crops which also uses a turnover rate as the inverse of leaf longevity. When the maximum plantation age (usually 25 years) of oil palm is reached and a new rotation cycle starts, the whole PFT is turned off and all C and N of the leaves, stem and roots go to litter. Existing fruit C and N of mature phytomers go to the fruit output pools. The PFT is then replanted in the next year and enters new phenological cycles.

15
In the CLM model plant growth, yield, stem turnover and litter-fall are simulated as the expansion and shrinking of C and N pools of the leaf, stem, root and grain/fruit components, and mediated by the developmental stage (phenology), atmospheric forcing, and N supply along with other environmental factors. The fate of newly assimilated carbon from photosynthesis is determined by a coupled C and N allocation routine.

20
Potential allocation for new growth of various plant tissues is calculated based on allocation coefficients and C : N ratios for each tissue type ( Table 2). The C : N ratios link C demand and N demand so that a N down-regulation mechanism is enabled for scaling down potential GPP depending on N availability from soil mineral N or retranslocated N pool.

25
Here allocation refers to the partitioning of remaining photosynthetic carbon for the new growth of various tissues, after accounting for their respiration costs and N downregulation. Photosynthesis calculation is based on the model of Farquhar et al. (1980) Stomatal conductance is linked to net leaf photosynthesis and scaled by the relative humidity and CO 2 concentration at the leaf surface as well as soil water stress (Collatz et al., 1991;Sellers et al., 1996). Only expanded leaves and the remaining photosynthetically active part of senescent leaves are involved in the photosynthesis process. Maintenance respiration is calculated as a function of temperature and the C and N content of each live vegetation component (Ryan, 1991). The major fraction of stem tissue is dead and thus excluded from respiration processes. CLM assumes a common base rate of maintenance respiration per unit N content for all live tissue types. The carbon cost of maintenance respiration at daytime is subtracted from available photosynthetic C pool before its allocation to new growth. Growth respiration cost is set 10 as a constant percentage (25 %) of allocated C to new growth (Oleson et al., 2013). The above photosynthesis and respiration mechanisms, already existing in CLM4.5, are well coupled with the new growth allocation mechanism and the sub-canopy structure for oil palm (Fig. 2). A two-step allocation scheme is designed for the multilayer sub-canopy phytomer 15 structure and according to the new phenology described in Sect. 2.1. First, available C (after subtracting respiration costs) is partitioned to the root, stem, overall leaf, and overall fruit pools at the PFT level. Available N is allocated accordingly by the C : N ratios. Then, the C and N allocated to the overall leaf or fruit is partitioned between different phytomers at the sub-PFT level. Details are described below. 20

PFT level allocation
Corley and Tinker (2003) described an "overflow" model that oil palm produces vegetative dry matter at a relatively constant rate after maturity and only the excess assimilates go to the reproductive pool, which can guide the partitioning of available assimilates between vegetative and reproductive pools. In our module, C and N allocation at 25 the PFT level is treated at two distinct phenological phases. Before fruiting starts (i.e. GDD 15 < GDD min ) all the allocation goes to the vegetative components. The initial allocation coefficients ( at a certain age determined by d mat . Real-time allocation to stem is partitioned to its live and dead stem tissues by the ratio F live stem . The live portion is then involved in stem turnover (Sect. 2.1.7). The equations reflect changed vegetative allocation strategy that reduces stem and root growth rates and shifts resources to leaf for maintaining LAI and increasing photosynthetic productivity when fruiting starts. The three vegetative alloca-10 tion ratios A leaf , A stem and A root always sum to 1.
At the second phase when fruit-fill begins a fruit allocation ratio A fruit is introduced, which is relative to the total vegetative allocation unity. To represent the dynamics of reproductive allocation effort of oil palm during its development, we adapt the stem allocation scheme for woody PFTs in CLM, in which increasing net primary production 15 (NPP) results in increased allocation ratio for the stem wood (Oleson et al., 2013). We use the similar formula for reproductive allocation of oil palm so that it increases with increasing NPP: where NPP mon is the monthly sum of NPP from the previous month and 100 20 (g C m −2 mon −1 ) is the base monthly NPP when the palm starts to yield (Kotowska et al., 2015). Parameters a and b adjust the base allocation rate and the slope of curve, respectively (Table 2). This function generates a dynamic curve of coefficient A fruit increasing from the beginning of fruiting to full vegetative maturity ( Fig. A1 in Appendix). The average value of A fruit is around 1 (relative to total vegetative allocation), 25 resulting nearly half of available assimilates allocated to the fruit pool and the rest to the vegetative pools. 4561 Total C and N allocation to the overall leaf pool (Eqs. A2 and A5) is divided to the displayed and storage pools by a fraction lf disp (Table 2) according to the following equation.
The leaf storage here is only intended to differentiate the pre-expansion (negative ranks, photosynthetically inactive) and post-expansion (positive ranks, photosynthetically active) growth phases of leaves and is not used for supporting other tissues. These two phases and two types of pools coexist throughout oil palm's life at the PFT level, although at the sub-PFT level each phytomer only transits from one phase to 15 another through time.

Sub-PFT level allocation
Total leaf and fruit allocations are partitioned to different phytomers, respectively, according to sub-PFT level allometry. Fruit allocation per phytomer is distributed to each phytomer using a sink size index: Real-time leaf allocations to the displayed and storage pools are distributed evenly to expanded and unexpanded phytomers, respectively, at each time step. When a phytomer enters the expansion phase, C and N from its leaf storage pools transfer gradually to the displayed pools during the expansion period. Therefore, a transfer flux is added to the real-time allocation flux and they together contribute to the post-expansion leaf growth. LAI is calculated only for each expanded phytomer according to a constant specific leaf area (SLA, Table 2) and prognostic amount of leaf C accumulated by phytomer n. In case it reaches the prescribed maximum (PLAI max ), partitioning of leaf C and N allocation to this phytomer becomes zero.

Other parameterizations
Nitrogen retranslocation is performed exclusively during leaf senescence and stem turnover. A part of N from senescent leaves and from the portion of live stem that turns to dead is remobilized to a separate N pool that feeds plant growth or reproductive demand (Sect. 2.1.6). Nitrogen of fine roots is all moved to the litter pool during root 20 turnover. We do not consider N retranslocation specifically during grain-fill that is designed for annual crops, which in effect turns the N contents in all vegetative parts to levels similar to crop residues in order to meet organ demand (Drewniak et al., 2013). For oil palm, once fruit-fill starts it continues year around at different phytomers. All the vegetative components have to be sustained at the same time. Therefore, it is not Introduction The fertilization scheme for oil palm is also adapted to the plantation management in our study area, which applies fertilizer biannually, only since 6 year old after planting, assuming each fertilization event lasts one day. Currently the CLM4.5 model uses an unrealistically high denitrification rate, which results in a 50 % loss of the unused available soil mineral nitrogen each day under conditions of excess N (Oleson et al., 5 2013). This caused our simple biannual regular fertilization nearly useless because peak N demand by oil palm is hard to predict given its continuous fruiting and vegetative growth and most fertilized N is thus lost in several days. The high denitrification parameterization has been recognized as an artifact (Drewniak et al., 2013;. According to a study on a banana plantation in the tropics (Veldkamp and Keller,10 1997), around 8.5 % of fertilized N is lost as nitrogen oxide (N 2 O and NO). Accounting additionally for a larger amount of denitrification loss to gaseous N 2 , we modified the daily denitrification rate from 0.5 to 0.001, which gives a 30 % annual loss of N due to denitrification that matches global observations (Galloway et al., 2004).
The irrigation option is turned off because oil palm plantations in the study area are 15 usually not irrigated. The vegetative structure is updated for oil palm including a special calculation of stem area index (SAI) and canopy top and bottom height allometry (Table B1 in Appendix), which are involved in biophysical processes in CLM model (Fig. 2). Other input parameters for the new oil palm PFT such as its optical, morphological, and physiological characteristics are estimated based on a literature review and 20 field observations ( and dry weight and SLA of whole leaf samples were measured for different sites. They were converted to LAI as above. Oil palm fresh bunch harvest data were collected for a whole year from July 2013 to July 2014. Harvest records from both PTPN-VI and the 8 validation sites were converted to harvested carbon (g C m −2 ) with mean wet/dry weight ratio of 58.65 % and C content 60.13 % per dry weight according to 20 C : N analysis (Kotowska et al., 2015). clay Acrisol. Soil texture such as sand/silt/clay ratios and soil organic matter C content were measured at multiply soil layers (down to 2.5 m) (Allen et al., 2015). They were used to create two sets of surface input data for the Harapan (H) and Bukit Duabelas (B) regions separately. ated spin-up procedure (Koven et al., 2013) was performed till 1850 with constant CO 2 and nitrogen deposition to get a steady-state estimate of soil C and N pools. Simulation continued on this equilibrium condition but was forced with dynamic CO 2 and NCEP meteorological data until 1990. Restart files from 1990 were used as soil initial condition for the following simulations.

20
The surface was then replaced with 100 % oil palm PFT, while the soil C and N pools were maintained from the initial condition. The new sub-canopy phytomer structure, phenology and allocation schemes were turned on. Recent simulations were all forced with dynamic CO 2 and 3 hourly ERA Interim reanalysis atmospheric data until 2014 (Dee et al., 2011). Input meteorological variables include surface incident solar ra- A simulation from 2002 to 2014 at the PTPN-VI site was used for model calibration. Additional eight simulations were run for the sites HO1, HO2, HO3, HO4, BO2, BO3, BO4, BO5 with different surface input files (for soil texture) and different atmospheric forcing files for the H and B plots, respectively. The simulations started from different years (1996,1997,1999,2000,2001,2002,2003,2004) when the palms were planted 5 at the individual sites. Outputs from these simulations were used to validate the model in terms of LAI and yield.

Calibration of key parameters
Both the PFT level and phytomer level LAI development were calibrated with field observations in 2014 from a chronosequence approach (space for time substitution) using oil palm samples of three different age and multiple phytomers of different rank (Sect. 3.1). Simulated yield outputs (around twice per month) were calibrated with monthly harvest records of PTPN-VI plantation from 2005 to 2014. Cumulative yields were compared because the timing of harvest in the plantations was largely uncertain and varied depending on weather and other conditions.

15
To simplify model calibration, we focused on parameters related to the new phenology and allocation processes. Phenological parameters listed in Table 1 were determined according to field observations and existing knowledge about oil palm growth phenology (Combres et al., 2013;Corley and Tinker, 2003) as well as plantation management in Sumatra, Indonesia. Allocation coefficients in Table 2 were more uncertain 20 and they were the key parameters to optimize in order to match observed LAI and yield dynamics. Measurements of oil palm NPP and its partitioning between fruit, canopy, stem, and root from the eight sites (Kotowska et al., 2015) were used as a general reference when calibrating the allocation coefficients.
Parameters related to photosynthesis, stomatal conductance and respirations were 25 set at similar levels as those of other crops, except that leaf traits such as PLAI max and SLA were determined by field measurements. Other parameters such as C : N ratios of

Sensitivity analysis
Performing a full sensitivity analysis of all parameters used in simulating oil palm (more than 100 parameters, though a majority are shared with natural vegetation and other 5 crops) would be a challenging work. As with calibration, we limited the sensitivity analysis to a set of parameters introduced for the specific PFT and model structure designed for oil palm (Tables 1 and 2). Among the phenological parameters, mxlivenp and phyllochron are closely related to pruning frequency (Sect. 2.1.6) but they should not vary widely for a given oil palm breed and plantation condition. Therefore, they were fixed at the average level for the study sites in Jambi, Sumatra (Table 1). GDD init was kept zero because only the transplanting scenario was considered for seedling establishment. We tested two hypotheses of phytomer level leaf development based on the other phenological parameters: (1) considering the leaf storage growth period, that is, the bud & spear leaf phase is explicitly simulated with the GDD parameter values in Table 1 and 15 lf disp = 0.3 in Table 2, (2) excluding the storage growth period by setting GDD exp = 0 and lf disp = 1 so that leaf expands immediately after initiation and leaf C and N allocation all goes to the photosynthetic active pools. The sensitivity of allocation and photosynthesis parameters in Table 2 were tested by adding or subtracting 10 or 30 % to the baseline values (calibrated) one-by-one and calculating their effect on final cumulative yield at the end of simulation (December 2014). In fact, all the allocation parameters are interconnected because they co-determine photosynthesis capacity and respiration costs as partitioning to the different vegetative and reproductive components varies. This simple approach provides a starting point to identify sensitive parameters, although a more sophisticated sensitivity analysis is 25 needed in the future effort.
Parameter PLAI max is only meant for error controlling, although in our simulations phytomer-level LAI never reached PLAI max (see Fig. 4 in results) because environmen- tal constraints and nitrogen down-regulation already limited phytomer leaf growth well within the range. Therefore, this parameter was not analyzed for its sensitivity. The C : N ratios and some photosynthesis and respiration parameters (Table A1) were evaluated thoroughly in Bilionis et al. (2015). Since we do not consider specific N retranslocation during fruit-fill, some C : N parameters are not used for oil palm and the aspect of N 5 content in different plant tissues is not prioritized for this sensitivity analysis.

Validation
In this study, we only validated the model structure and model behavior on simulating aboveground C partitioning and flux as represented by LAI and fruit yield. Independent LAI and yield data in 2013-2014 from the eight different sites were compared with 10 model predictions from eight simulations using the above model setting and calibrated parameters.

Calibration with LAI and yield
In model calibration with the PTPN-VI plantation, the PFT-level LAI dynamics simulated 15 by the model incorporating the pre-expansion phase matches well with the LAI measurements for three different ages (Fig. 3). Simulated LAI for the PFT increases with age in a sigmoid relationship. The dynamics of LAI is also impacted by pruning and harvest events because oil palms invest around half of their assimilates into fruit yield. Oil palms are routinely pruned by farmers to maintain the maximum number of expanded 20 leaves around 40. Hence, when yield begins 2-3 years after planting, LAI recurrently shows an immediate drop after pruning and then quickly recovers. Simulations without the pre-expansion storage growth phase show an unrealistic fast increase of LAI before 3 years old, much higher than observed in the field. At older age after yield begins, LAI drops drastically and recovers afterwards. Although the final LAI could stabilize at a similar level, the initial jump and drop of LAI at young stage do not match field observations and cannot be solved by adjusting parameters other than GDD exp . Hereafter, the following simulations were all run using the pre-expansion phase. The phytomer level LAI development also matches well with field measurements (Fig. 4). This sub-PFT level LAI dynamics clearly marks the 6-phase phenology cycle, 5 including gradual increment during phytomer development and gradual declining during senescence. The figure only shows photosynthetic LAI that updates at every time step. Leaf C accumulation before leaf expansion is not shown here. Only when the palm becomes mature, phytomer LAI could come closer to the prescribed PLAI max (0.165). However, during the whole growth period from 2002 to 2014 none of the phytomers 10 have reached PLAI max , which is the prognostic result of the carbon balance simulated by the model.
The cumulative yield of baseline simulation has overall high consistency with harvest records (Fig. 5). The mean percentage error (MPE) is only 4 %. The slope of simulated curve increases slightly after 2008 when the LAI continues to increase and 15 NPP reaches a high level (Fig. A1). The harvest records also show the same pattern after 2008 when heavy fertilization began (456 kg N ha −1 yr −1 ). If a high denitrification rate is used, fertilization becomes useless even given the already high dosage and the plant cannot catch up with observation after 2009.
Comparing the monthly trend of yield (Fig. 6), the per-month harvest records exhibit 20 strong zig-zag pattern. One reason is that oil palms are harvested every 15-20 days and summarizing harvest events by calendar month would result in uneven harvest times per month, e.g. two harvests fall in a previous month and only one in the next month. It would be more appropriate to use yield data per harvest event to compare with model output. Yet it still shows that harvests at PTPN-VI plantation dominated from 25 October to December whereas in the earlier months of each year harvest amounts were significantly lower. The simulated amount of yield per harvest event has less seasonal fluctuation, but it responds to the drought periods (Fig. 6) correlation exists between simulated yield and the mean precipitation of a 60 day period before each harvest event (Pearson's r = 0.15).

Sensitivity analysis
The leaf nitrogen fraction in Rubisco (F LNR ) is shown to be the most sensitive parameter (Fig. 7), because it determines the maximum rate of carboxylation at 25 • C (V cmax25 ) to-5 gether with SLA (also sensitive), foliage nitrogen concentration (CN leaf , Table B1) and other constants. Given the fact that F LNR should not vary widely in nature for a specific plant, we constrained this parameter within narrow boundaries to give a V cmax25 around 100, which is similar to that shared by all other crop PFTs (100.7) and higher than forests (around 60) in CLM. SLA (0.013) was fixed by field measurements. The 10 value is only representative of the photosynthetic leaflets. The heavy rachis that supports leaflets is about 75 % of the weight of a whole phytomer and contains a large proportion of lignified tissue (Corley and Tinker, 2003). Thus, the rachises are considered as a part of the stem like tree branches and included in stem turnover, which reduces respiration cost compared to when it is maintained as a part of the metaboli-15 cally active leaf. However, bias could exist as the physiology of rachis esp. during the juvenile phase of bud & spear leaf is likely very different from stem. We use this parameterization because it fits the structure of the CLM model. The initial root allocation ratio (a i root ) has signification influence on yield because it modifies the overall respiration cost along the gradual declining trend of fine root growth across 25 years (Eqs. A1, Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | stem turnover rate (Sect. 2.1.7), and therefore it brings higher respiration cost and lower yield. Decreasing the fruit allocation coefficient a results in a higher base rate of A fruit according to Eq. (2), whereas increasing coefficient b brings up the rate of change and final magnitude of A fruit if NPP rises continuously. Their relative influence on yield is lower than the leaf allocation coefficients because of the restriction by NPP dynam-5 ics (Eq. 2, Fig. A1). Parameters lf disp and transplant have negligible effects. lf disp has to work together with the phenological parameter GDD exp to give a reasonable size of spear leaf before expansion according to field observation (Sect. 5.2). Varying the size of seedlings at transplanting by 10 or 30 % does not alter the final yield, likely because the resulting initial LAI is still within a limited range (0.1-0.2) given the baseline value 10 0.15.

Model validation with independent dataset
The LAI development curves for the eight oil palm sites follow similar patterns since field transplanting in different years (Fig. 8a). The final LAIs at 2014 match poorly with observations at each individual site, but the average LAI values between simulations 15 and field measurements are comparable (MPE = 10 %) (Fig. 8b). There is large uncertainty in field LAI estimates because we did not directly measure LAI at the plot level but only sampled leaf area and dry weight of individual phytomers. Plantation management is also very different between these 8 smallholder plantations. For example, H plots applied significantly higher fertilization than B plots, while the model prescribes 20 the same fertilization for all plots (Sect. 2.3). The simulated annual yields match closely with the average yield of the eight sites measured in 2013-2014 (MPE = −4 %) but the model predicted variability across the sites is much lower than field records (Fig. 9). Modelled yield generally increases with plantation age, which can be explained by the increasing fruit allocation rate A fruit with 25 increasing LAI and NPP (Fig. A1). We do not have data to test an aging decline function of growth and yield and assume the oil palm plantations remain productive for 25 years (Age max ) before replanting. Site-to-site variation in the field harvest data is still difficult 4572 Introduction

Tables Figures
Back Close

Full Screen / Esc
Printer-friendly Version

Interactive Discussion
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | to interpret, because the measured yields do not correlate with age and LAI measurements, possibly reflecting unknown differences in fertilization or other management practices.
Overall, the simulated LAI and yield seem reasonable on an average, but has smaller site-to-site variability than observation as fertilization and other management impacts 5 were set constant in the model.

Discussion
Calibration and validation with multiple site data demonstrate the functionality of the new PFT and its sub-canopy structure for simulating growth and yield of the unique oil palm plantation system within a land surface modeling context. Results showed 10 a limitation of the model, that is, the difficulty to capture the large site-to-site variation in yield and LAI. The discrepancy was very likely due to insufficient representation of management (e.g. fertilization, harvest cycle) in the model, which has been shown to be crucial for determining oil palm growth and yield (Euler et al., 2015). Field observations in the 8 smallholder plots (B and H plots) show that increased harvest cycle 15 decreases yield because of fruit respiration costs. Other factors such as pruning, soil conditions (only two categories for H and B plots) and possibly different oil palm progenies could also result in difference in the average size and number of leaves and fruits per palm, but the model uses uniform parameterization of these aspects across sites. Especially the amount and timing of fertilization vary largely from plantation to 20 plantation and from year to year. Thus a more complex dynamic fertilization scheme needs to be devised and evaluated thoroughly (together with the new denitrification rate) with additional field data, which we lack at the moment. The seasonal variability in yield simulated by the model well corresponds to precipitation data but it is difficult to interpret the difference with harvest records due to the artificial zig-zag pattern. Moreover, the harvest records from plantations do not necessarily correspond to the amount of mature fruits along a phenological time scale due to varying harvest arrangements,

Tables Figures
Back Close

Full Screen / Esc
Printer-friendly Version

Interactive Discussion
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | e.g. fruits are not necessarily harvested when they are ideal for harvest, but when it is convenient. Observations of mature fruits on a tree basis would be more suitable to compare with modeled yield, but such data are not available. Some studies have also demonstrated important physiological mechanisms on oil palm yield including inflorescence gender determination and abortion rates that both respond to seasonal climatic 5 dynamics although with a time lag (Combres et al., 2013;Legros et al., 2009). The lack of representation of such physiological traits might affect the seasonal dynamics of yield simulated by our model, although these mechanisms are rarely considered in a land surface modelling context. Nevertheless, the multilayer phytomer-based structure, the extended phenological 10 phases for a perennial crop PFT and the two-step allocation scheme are distinct from existing functions in land models. Their significance and implications as shown in the results are discussed below. In the future, a fuller picture of the carbon, water and energy fluxes over the oil palm landscape will be examined with the oil palm module presented here and be evaluated with Eddy Covariance flux observation data.

The unique phenological phase for leaf storage growth
The pre-expansion growth phase is proved necessary for simulating phytomer leaf growth and PFT LAI in a prognostic manner within the CLM framework. Real-time allocation of assimilates to the leaf C pool is depending on incoming solar energy, N availability and other environmental forcing variables. In our model implementation, the 20 leaf C storage pool accumulated by each phytomer during its bud & spear leaf stage is released gradually, which comes together with the real-time leaf C allocation, to support the increment of photosynthetic active leaf C pool within the leaf expansion period. Therefore, the storage growth period and abundant storage C pool provides an efficient buffer to support phytomer development and maintain overall LAI when resources are 25 limited during fruit-fill or environmental constraints. It also avoids an abnormally fast increase of LAI in the young palms before fruiting when C and N allocation is dedicated to the vegetative components. ops unrealistically fast at young age and then enters an emergent drop once fruit-fill begins (Fig. 3). This is because the plant becomes unable to sustain leaf growth just from its current photosynthetic assimilates when a large portion is allocated to fruit, together with removal of leaf C by pruning. The pre-expansion and post-expansion growth phases are usually not represented 5 together in other studies. For example, the APSIM-Oil Palm model (Huth et al., 2014) considers only post-expansion leaf growth for a period of 5 phyllochrons (∼ 2 months) and it prescribes a logistic leaf area expansion pattern that is not constrained by assimilates supply. The ECOPalm model (Combres et al., 2013) assumes a 7 month leaf growth period with constant growth rate solely within the pre-expansion phase, that is, 10 at negative ranks in bud & spear leaf stage. According to field observations, it is more realistic to include both pre-expansion and post-expansion growth for a phytomer to reach full vegetative maturity. Furthermore, differentiating the two contrasting phases of leaf growth could avoid abrupt increase in photosynthesis if a phytomer with full dry mass shifts from photosynthetically inactive to active status at one step. We imple-15 ment the leaf senescence function for the similar purpose (Sect. 2.1.6). The realistic parameterization of senescence differentiates the photosynthetic capacity of a leaf at the bottom layer of canopy from those on the top so as to avoid unwanted drastic reduction in photosynthetic LAI that would happen if the bottom leaves were turned off immediately.

Allocation strategy
Allocation patterns for common annual crops are easier to simulate because their LAI is often assumed to decline during grain-fill (Levis et al., 2012) cannot sustain vegetative maintenance and growth. The dynamics of A fruit as a function of monthly NPP (Eq. 2) is meant to capture the increasing yield capacity of oil palms during maturing at favorable conditions (often the case in oil pam plantations). The average near 1 : 1 ratio to partition available assimilates to the reproductive and vegetative pools matched closely with field observations (Kotowska et al., 2015). In the simulations A fruit increased rapidly from the base rate within one year after fruiting began, without impairing allocation for LAI development. This is because within the same period NPP also increased nearly 50 % (Fig. A1). Correlating A fruit with NPP thus matches the "overflow" model and follows ecophysiological evidence of source limitation of oil palm yield (Corley and Tinker, 2003). It avoids unrealistic increase of A fruit 10 under stressed conditions if a simple time-dependent or vegetative-size-dependent fruit allocation function is used. When severe stress conditions do happen, this NPP-related function is expected to decrease fruit allocation and shifts resources to the vegetative components.
The pre-expansion phase is proved to be a crucial strategy that enables oil palm to 15 maintain a stable LAI when a large ratio of reproductive allocation and leaf pruning are co-occurring. We set the value of lf disp to be 0.3 for the post-expansion phase which means a majority of available C and N fluxes are allocated for storage growth. Together with the length of each growth phase set by GDD exp and GDD L.mat , they determine that roughly 60-70 % of leaf C in a phytomer is accumulated before its expansion with the 20 remaining part developed until leaf maturity. This is comparable to observations on coconut palm that dry mass of the oldest unexpanded leaf accounts for 60 % of that of a mature leaf (Navarro et al., 2008).

The sub-canopy structure
We devised the sub-canopy structure for the phenology and allocation mechanisms 25 in order to capture the growth and yield dynamics of the oil palm plantation system. Phytomer-to-phytomer relationship is considered more "social" rather than independent (Corley and Tinker, 2003 ranks of the canopy while most photosynthesis happens on the top layers of young expanded phytomers which are exposed to more sunshine but not bearing fruits yet. This implies the young phytomers are more than self-sufficient and they share excessive photosynthetic products with older phytomers to support their fruit-fill demand. All the photosynthetic phytomers also support the storage growth of unexpanded phy-5 tomers that are in the bud & spear leaf stage, in addition to maintaining stem and root growths. Sharing a common pool of photosynthetic assimilates for C and N allocation also implies a competition among phytomers by their different demand levels at different phenological stages. Thus, the allocation ratios are demand-driven as controlled by phenology, whereas the actual allocation rates (in terms of g C s −1 ) are constrained by 10 assimilates supply which is the result of climatic forcing, N down-regulation, and other environmental variables. The sub-canopy structure could be also useful for considering light competition among expanded phytomers. However, in the current paper the CLM default onelayered sun/shade big-leaf canopy model is used for radiative transfer and photosyn-15 thesis, which is based on the two-stream approximation for an integrated canopy (Sellers, 1985). It also derives the sunlit and shaded fractions for the whole canopy using an analytical solution by Dai et al. (2004), but the light profile is not resolved through the canopy depth. Thus, specific light competition among different phytomers and influence on their photosynthetic rates is not considered. But the sun-shade canopy is 20 a relevant and accurate simplification according to De Pury and Farquhar (1997) and Ryu et al. (2011). It has been cross-validated with refined 3-D models on coconut palm by Roupsard et al. (2008). Calculating only canopy-level photosynthesis is sufficient for the two-step social-sharing allocation strategy in which a common assimilates pool is partitioned to different phytomers regardless of the light profile.

25
Nevertheless, resolving light profile at sub-canopy layers could possibly improve the accuracy of deriving canopy-level photosynthesis and energy fluxes. A classical multilayer radiative transfer scheme by Norman (1979)  lations are currently incorporated into the CLM model as an option to work with the oil palm modules, allowing simulating sub-canopy light profile.

Summary
The development of the perennial oil palm PFT including its structure, phenology, and carbon and nitrogen allocation functions was proposed for modeling an important agri- 5 cultural system and land use transformation in Indonesia. This paper demonstrates the ability of the new oil palm module to simulate the inter-annual dynamics of vegetative growth and fruit yield from field planting to full maturity of the plantation. The multilayer canopy structure and associated sub-PFT level phenology and allocation strategy are necessary for this perennial evergreen crop which yields continuously on multiple phy-10 tomers. They are potentially applicable to similar plantations (coconut palm, date palm, etc.). The pre-expansion leaf storage growth phase is proved essential for buffering and balancing overall vegetative and reproductive growth. Average LAI and yield were satisfactorily simulated for multiple sites, and the model, designed initially to be run for large areas also proved to be attractive for point simulations, thanks to the overarch-15 ing framework of the CLM model and the underlying ecophysiological processes that were introduced for oil palm. The point simulations here provide a starting point for calibration and validation at large scales. The oil palm module is aimed to capture the C and N dynamics within the plantation system from the beginning of its establishment till full vegetative maturity and final 20 rotation (replanting). To be run in a CLM regional or global grid, the age class structure of plantations needs to be taken into account. This can be achieved by setting multiple replicates of the oil palm PFT, each planted at a point of time at a certain grid. As a result, a series of oil palm cohorts developing at different grids could be configured with a transient PFT distribution dataset, which allows for a quantitative analysis of the effects of land-use changes, specifically rainforest to oil palm conversion, on carbon, water and energy fluxes. This will contribute to the land surface modeling community for simulating this structurally unique, economically and ecologically sensitive, and fast expanding oil palm land cover.

Appendix A
Allocation equations for the vegetative components before and after fruiting starts: 1. From leaf emergence to begining of fruiting where DPP Age max ≤ 1, and DPP is the days past planting.
where a 2 leaf equals the last value of A leaf calculated right before fruit-fill starts and DPP 2 is the days past planting right before fruit-fill starts. 0 ≤ DPP−DPP 2 Age max ×d mat −DPP 2 ≤ 1.  (Combres et al., 2013;Corley and Tinker, 2003;Hormaza et al., 2012;Legros et al., 2009   to P n indicate expanded phytomers and P −1 to P −n at the top indicate unexpanded phytomers packed in the bud. Each phytomer has its own phenology, represented by different colors corresponding to: (b) the 6-phase phytomer phenology: from initiation to leaf expansion, to leaf maturity, to fruit-fill, to harvest, to senescence and to pruning. Phytomers initiate successively according to the phyllochron (the period in heat unit between initiations of two subsequent phytomers). Introduction

Tables Figures
Back Close

Full Screen / Esc
Printer-friendly Version

Interactive Discussion
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Figure A1. Reproductive allocation rate (A fruit ) as a function of monthly sum of NPP from the previous month (NPP mon ) according to Eq. (2) in Sect. 2.2.1. A fruit is relative to the vegetative unity (A leaf + A stem + A root = 1).