Representing anthropogenic gross land use change , wood harvest , and forest age dynamics in a global vegetation model ORCHIDEE-MICT v 8 . 4 . 2

Land use change (LUC) is among the main anthropogenic disturbances in the global carbon cycle. Here we present the model developments in a global dynamic vegetation model ORCHIDEE-MICT v8.4.2 for a more realistic representation of LUC processes. First, we included gross land use change (primarily shifting cultivation) and forest wood harvest in addition to net land use change. Second, we included sub-grid evenly aged land cohorts to represent secondary forests and to keep track of the transient stage of agricultural lands since LUC. Combination of these two features allows the simulation of shifting cultivation with a rotation length involving mainly secondary forests instead of primary ones. Furthermore, a set of decision rules regarding the land cohorts to be targeted in different LUC processes have been implemented. Idealized site-scale simulation has been performed for miombo woodlands in southern Africa assuming an annual land turnover rate of 5 % grid cell area between forest and cropland. The result shows that the model can correctly represent forest recovery and cohort aging arising from agricultural abandonment. Such a land turnover process, even though without a net change in land cover, yields carbon emissions largely due to the imbalance between the fast release from forest clearing and the slow uptake from agricultural abandonment. The simulation with sub-grid land cohorts gives lower emissions than without, mainly because the cleared secondary forests have a lower biomass carbon stock than the mature forests that are otherwise cleared when sub-grid land cohorts are not considered. Over the region of southern Africa, the model is able to account for changes in different forest cohort areas along with the historical changes in different LUC activities, including regrowth of old forests when LUC area decreases. Our developments provide possibilities to account for continental or global forest demographic change resulting from past anthropogenic and natural disturbances.

Abstract. Land use change (LUC) is among the main anthropogenic disturbances in the global carbon cycle. Here we present the model developments in a global dynamic vegetation model ORCHIDEE-MICT v8.4.2 for a more realistic representation of LUC processes. First, we included gross land use change (primarily shifting cultivation) and forest wood harvest in addition to net land use change. Second, we included sub-grid evenly aged land cohorts to represent secondary forests and to keep track of the transient stage of agricultural lands since LUC. Combination of these two features allows the simulation of shifting cultivation with a rotation length involving mainly secondary forests instead of primary ones. Furthermore, a set of decision rules regarding the land cohorts to be targeted in different LUC processes have been implemented. Idealized site-scale simulation has been performed for miombo woodlands in southern Africa assuming an annual land turnover rate of 5 % grid cell area between forest and cropland. The result shows that the model can correctly represent forest recovery and cohort aging arising from agricultural abandonment. Such a land turnover process, even though without a net change in land cover, yields carbon emissions largely due to the imbalance between the fast release from forest clearing and the slow uptake from agricultural abandonment. The simulation with sub-grid land cohorts gives lower emissions than without, mainly because the cleared secondary forests have a lower biomass carbon stock than the mature forests that are otherwise cleared when sub-grid land cohorts are not considered. Over the region of southern Africa, the model is able to account for changes in different forest cohort areas along with the historical changes in different LUC activities, including regrowth of old forests when LUC area decreases. Our developments provide possibilities to account for continental or global forest demographic change resulting from past anthropogenic and natural disturbances.  Guo and Gifford, 2002;Poeplau et al., 2011;Powers et al., 2011).
Globally, LUC activities have contributed significantly to historical anthropogenic carbon emissions. It is estimated that about 800 Mha (1 Mha = 10 6 ha) of forests were cleared for agricultural purpose and that 2000 Mha of forests were harvested during 1850-1999, giving rise to cumulative emissions of 124 Pg C, or 33 % of the total anthropogenic emissions (Houghton, 1999). Houghton et al. (2012) reviewed LUC emissions from multiple studies and estimated the annual global LUC emissions as 1.1 Pg C yr −1 during 1980-2009, with an uncertainty of 0.5 Pg C yr −1 . Different estimations of historical LUC emissions by dynamic global vegetation models (DGVMs) show a spread as large as 1 Pg C yr −1 (see Fig. 1 in Houghton et al., 2012; see also Hansis et al., 2015, for an even larger range among model estimations). This is partly due to different forcing data used and initial carbon stocks simulated  but also because of different implementations of LUC processes in DGVMs (Prestele et al., 2017). Given the importance of understanding historical LUC emissions in projecting the future land-based mitigation potential, a more realistic representation of LUC processes and land management in DGVMs is desirable.
In most global studies, only net transitions were accounted for in the LUC processes simulated by DGVMs (Le . Changes in land use over each model grid cell are diagnosed as the difference in ground fractions of different land cover types between two consecutive years. At a typical spatial resolution of 0.5 • for global applications (e.g. TRENDY, Sitch et al., 2015;MsTMIP, http://nacp.ornl.gov/ MsTMIP_simulations.shtml), such a scheme has ignored the simultaneous, bidirectional transitions between two vegetation types within the same grid cell (i.e. gross transitions). Such gross transitions can arise from spatial upscaling of LUC data or from certain land use activities. A typical example is shifting cultivation, a form of smallholder subsistence agriculture primarily occurring in tropical regions that involves clearing a forest for a non-permanent agricultural land, which is often abandoned later. Shifting cultivation was historically important in many tropical regions for the subsistence of indigenous people (Hurtt et al., 2011;Lanly, 1985), although more recently it has been in the process of being superseded by more intensified land management . Forest management such as a clear-cut for wood harvest followed by replanting trees is another type of gross transition. Although it does not entail any net change in land cover (forest remaining forest), species choice and forest management can have a significant effect on carbon stocks and fluxes (Erb et al., 2017).
More and more DGVMs started to include gross transitions and we provide an overview of them in Table 1. All models in Table 1 include shifting cultivation and wood harvest except that shifting cultivation is not included in ISAM, and five of them include sub-grid secondary land tiles when accounting for LUC. A recent review by Arneth et al. (2017) found that including processes that have been previously neglected in DGVMs, including gross transitions and other land management processes such as crop harvest and management, can lead to an upward shift of estimated LUC emissions. Their study thus highlights the importance of including these processes. Furthermore, to more robustly account for shifting cultivation and wood harvest, which often have a certain rotation length and mainly involve secondary forests of different ages, it is critical for DGVMs to include subgrid differently aged land cohorts. This feature exists in some DGVMs that combine with a forest gap model (e.g. LPJ-GUESS; Bayer et al., 2017) but it would be difficult to represent forest species change because different tree plant functional types (PFTs) are mixed over a model grid cell. The same also applies for LM3V (Shevliakova et al., 2009). Other so-called area-based DGVMs (Smith et al., 2001) such as ISAM (Jain et al., 2013) and LPX-Bern 1.0 (Stocker et al., 2014) included secondary land tiles in the model but their capability to represent different rotation lengths in land use is limited. In the ORCHIDEE model, sub-grid forest cohorts have been recently included in the ORCHIDEE-CAN branch mainly for forest management purposes , but a combination of both sub-grid land demography and gross land transition is still missing.
Here we present the new model developments in OR-CHIDEE that combine both sub-grid land cohorts and gross LUC. The objectives of this study are (1) to document a new LUC module, including sub-grid vegetation cohorts, forest harvest, and gross LUC in the ORCHIDEE model, that can be run with and without sub-grid age dynamics; (2) to document through an idealized pixel simulation the simulated carbon fluxes from shifting cultivation or land turnover between model set-ups with and without sub-grid age dynamics; and (3) to document the model behaviour and forest age dynamics associated with the historical changes in LUC activities. Whereas the current paper focuses on documenting new model developments and subsequent changes in model behaviour, a companion paper presents a global reanalysis of historical LUC emissions . sphere. The carbon module simulates vegetation carbon cycle processes, including photosynthesis, photosynthate allocation, vegetation mortality and recruitment, phenology, litter fall, and soil carbon decomposition. ORCHIDEE-MICT is a branch initially focusing on improving high-latitude processes (e.g. soil freezing, snow processes, permafrost dynamics, and northern wetlands) but is now under development to include more processes. Of interest for this study is that the grassland management module developed in Chang et al. (2013) is included (r2615). This allows for distinction between natural grassland and pasture that have been mixed together in previous LUC simulations by ORCHIDEE. In ORCHIDEE, land cover types are represented as PFTs, with each PFT being associated with a set of parameters. A typical model simulation consists of two stages: a spin-up stage with stable or constant forcing data until the model reaches an approximately equilibrium state, to mimic an era with no appreciable human perturbation, and a transient stage in which the model is forced with temporally varying forcings (e.g. climate, atmospheric CO 2 , land cover). The LUC module prior to this study accounts for net transitions only (Piao et al., 2009a) and has been used in many applications (e.g. CMIP5, https://portal.enes.org/ models/earthsystem-models/ipsl/ipslesm; TRENDY, . To simulate historical LUC, a spin-up run is initiated with a given initial land cover map (i.e. a PFT map), and then vegetation distribution is updated annually with prescribed PFT map time series during the transient simulation. The LUC module simply compares grid cell fractions of different PFTs between the current simulation year and the next year. Then 12 vegetative PFTs (all standard model PFTs excluding the bare soil PFT) are separated into two groups with expanding versus contracting areas. Carbon stocks and associated carbon fluxes on shrinking PFTs are displaced to expanding PFTs in proportion to their respective surface increments.

Concept of gross transitions in relation to vegetation age structure
The numerical implementation of net transitions is straightforward. However, as explained in the introduction, this scheme omits important sub-grid gross land use transitions. Figure 1 uses an exemplary grid cell to illustrate the dif- ference between the two LUC schemes: one accounting for net transitions only (Fig. 1b), and the other accounting for gross transitions but with no sub-grid cohorts (Fig. 1c, d).
Although the areas of forest and cropland after LUC are identical (Fig. 1b, d), carbon stocks for the same vegetation type (e.g. forest) are different between the two schemes. According to the net transition scheme, the carbon stock of the final forest patch shown in Fig. 1b Stocker et al. (2014). LUC here is similar to in Fig. 1, except that forest is no longer a single ageless patch but consists of two patches of primary and secondary forests, i.e. having an age structure. (a) The same area of forest is converted to cropland as in Fig. 1a but conversion is made from primary forest.
Consequently, a young cropland patch with rich legacy forest soil C is established. Meanwhile, a very young forest patch is established due to the bidirectional gross land flux. Because the model uses multiple sub-grid patches to represent vegetation age structure (or differently aged cohorts), merging of patches with different carbon stocks is no longer necessary. Panel (c) shows an alternative to (a) in which conversion of forest to cropland is made on a secondary forest. Correspondingly, in (d), which shows the post-LUC state of (c), the established young cropland patch will have lower legacy soil C than that in (b).
is an area-weighted mean between the original forest patch not being impacted by LUC and the newly established forest with a low carbon density that results from cropland abandonment. Consequently the carbon stock of the grid cell is expected to be smaller in Fig. 1d than in 1b and LUC carbon emission in Fig. 1d is conversely larger than in Fig. 1b. Figure 1c represents the real land cover state after LUC, while the merging shown in Fig. 1d is only a necessary simplification when no sub-grid cohorts are represented in the model. Ideally, the model capability could be expanded to include cohorts to represent the real world case as in Fig. 1c. In addition, inclusion of sub-grid cohorts would allow not only the distinction between original intact forest and newly established forest but also among different forest cohorts (e.g. primary versus secondary forests) regarding which forest patch to be cleared for cropland. Figure 2 illustrates a case in which gross LUC is combined with sub-grid cohort representation in the model. Here, mul-tiple patches within a grid cell are used to represent cohorts of a single vegetation type but with different ages since establishment. These cohorts often have different carbon stocks either due to different lengths in carbon accumulation time (e.g. for forest) or due to different extents to which legacy soil carbon is present (e.g. for croplands establishing on former forests). The areas subject to gross LUC transition in Fig. 2a and b remain the same as in Fig. 1a (dashed red rectangles), but primary and secondary forests are cleared in Fig. 2a and b, respectively. Thus, LUC emissions from clearing of primary forest are expected to be higher due to its higher biomass stock. Correspondingly, the legacy soil carbon stocks on the cohort of new cropland are also higher (shown in Fig. 2b and d). Figures 1 and 2 have shown the example of LUC transitions between forest and cropland, but other types of LUCs, including forest harvest, can be handled in a similar way. In the case of forest harvest, having cohorts avoids the simplification of merging a young re-established forest after harvest with the original forest, which serves as the exact source of harvest. This can effectively simulate forest management practices that induce rotations of different forest cohorts (e.g. see McGrath et al., 2015, for a forest management history in Europe).

Expansion of ORCHIDEE-MICT capacity to represent sub-grid vegetation cohorts
In order to simulate gross LUC combined with sub-grid vegetation cohorts as illustrated in Fig. 2, we expanded the ORCHIDEE-MICT capability to include sub-grid evenly aged cohorts. This necessitates multiple patches within a grid cell for a single PFT, which inherits most of the parameters from its parent PFT (they still belong to the same PFT and thus are largely physically similar). These patches are named cohort functional types (CFTs) here, to be distinguished from the original plant functional types. In this sense, the original PFTs actually become "meta-PFTs" which were named meta-classes (MTCs). As subsequent LUCs generate differently aged CFTs, the computational demand will be greatly increased. Hence, the number of CFTs within an MTC is limited to a user-defined number. ORCHIDEE-trunk has a feature called "PFT externalization" that allows the creation of a new user-specified PFT by inheriting its parameters from an existing one. A user can then modify specific parameters at their convenience. Based on this feature, the ORCHIDEE-CAN branch (the svn revision number is 2566; Naudts et al., 2015Naudts et al., , p. 2037) has developed representation of sub-grid forest age classes (i.e. equivalent to our CFTs here). Each forest age class is an inheritance of a given forest MTC. There, the transitions from one age class to another were defined by tree diameters. When a forest of a certain age class reaches its diameter limit, it moves into the next age class, and is merged with the existing forest patch of that age class if there is one. All associ-ated biophysical and biogeochemical variables are merged as well following an area-weighted mean approach with a few exceptions for discrete variables such as the applied forest management strategy.
ORCHIDEE-MICT also inherits this externalization feature from ORCHIDEE-trunk. Here we ported the codes of forest age class functionality from ORCHIDEE-CAN to develop the CFT functionality needed for LUC simulation with cohorts in ORCHIDEE-MICT. The code base to include sub-grid forest cohorts was migrated from ORCHIDEE-CAN, with substantial adaptions being made in ORCHIDEE-MICT. Except for this, all other LUC developments have been achieved within the current study. Contrary to ORCHIDEE-CAN (see above), ORCHIDEE-MICT uses woody biomass to delimit different forest cohorts, with older cohorts having a higher woody biomass. Forest grows old by moving from the current cohort to the next one when the woody biomass exceeds the cohort upper boundary. Except for the cohort boundaries, no further cohort-specific parameterizations have been performed, so essentially all cohorts are governed by the same set of biophysical and ecological parameter values. However, in ORCHIDEE-MICT there are indeed some simple aging processes to proximate the key changes when a forest grows old: notably, the net primary production (NPP) allocation to below-ground sapwood decreases with the time since establishment.
In addition, we expanded the concept of CFT to croplands, natural grasslands, and pastures. Cohorts are defined with their soil carbon stocks for these herbaceous vegetation types; this is a definition relevant to LUC emission calculation. Because the directional change of soil carbon largely depends on the vegetation types before and after LUC and on climate conditions Poeplau et al., 2011), ideally agricultural cohorts from different origins should be differentiated. However, to avoid inflating the total number of cohorts and the associated computational demand, as a first attempt, we simply divide each herbaceous MTC into two broad sub-grid cohorts according to their soil carbon stocks and without considering their individual origins. We expect that such a parameterization can accommodate some typical LUC processes, such as the conversion of forest to cropland where soil carbon usually decreases with time, but not all LUC types (for instance, soil carbon stock increases when a forest is converted to a pasture). The biomass or soil carbon thresholds that delineate different CFTs must be properly parameterized in order to have sensible CFT segregation within different contexts of land use change. This will be further detailed in Sect. 2.2.3. In practice, for single-site simulations, the parameterization could be set up via a configuration file enumerating the thresholds for all CFTs. For regional applications, an input file containing spatially explicit thresholds will be used.
The implementation of sub-grid cohort function types as inheritances of meta-classes and the corresponding hierarchy is exhibited in Fig. 3a. Tier 1 of the model parameterization CFT 1,1 CFT 2,1 CFT 2,1 CFT 2,2 … … CFT 9,1 CFT 2,6 CFT 3,1 CFT 3,2 CFT 2,2 … … CFT 3,6 CFT 9,2 … … … CFT 9,1 CFT 9,2 CFT 2,6 … … CFT 9,6 CFT 9,6 CFT 10,1 CFT 10,1 CFT 10,2 CFT 11,1 CFT 11,1 CFT 10,2 CFT 11,2 CFT 11,2 CFT 12,1 Model parameterizaIon hierarchy LUC hierarchy (a) (b) Figure 3. Two parallel hierarchies from the model parameterization and land use change perspective. (a) Sub-grid cohort function types (CFTs) as inheritances of meta-classes (MTCs) and the corresponding parameterization hierarchy. There are in total 14 vegetative MTCs corresponding to four vegetation types. The notation of CFT i,j indicates that it inherits from MTC i and belongs to the j th cohort (Cohort j ). Each forest MTC has six cohorts, with Cohort 1 being the youngest and Cohort 6 the oldest, whereas each herbaceous MTC is set tentatively to have two cohorts. Darker colours indicate older cohorts. (b) Within the gross LUC module hierarchy, Tier 3 remains the level of CFT, but CFTs are reorganized to derive the Tier 2 information based on the level of cohorts under the same Tier 1 as in (a). A cohort baring the notation of Cohort v,i indicates it belongs to vegetation type "v" (where "v" could be forest, natural grassland, pasture, or cropland) and meta-class "i". This reorganization of the hierarchy from the left to the right side is to prepare for properly allocating prescribed LUC transitions first onto the cohort level, then further to different CFTs within each cohort.
hierarchy corresponds to the four basic vegetation types (forest, natural grassland, pasture, and croplands, abbreviated as f, g, p, and c respectively). Tier 2 corresponds to meta-classes in ORCHIDEE-MICT, which contain one bare soil MTC and 14 vegetative MTCs, with each vegetative MTC belonging to one of the four basic vegetation types. Tier 3 corresponds to CFTs. A CFT is noted as CFT i,j to denote that it inherits its parameter values from the MTC i and belongs to the j th cohort. For this study, forest MTCs contain six CFTs and herbaceous MTCs contain two CFTs. The number of CFTs for each MTC is not hard-coded in the model and can be specified by users via a configuration file. With sub-grid cohorts, the model spin-up run is initiated with an input MTC map, essentially the same as in the case without sub-grid cohorts (recall that in Sect. 2.1.1 this MTC map is called a PFT map). But the difference is that the initial prescribed areas (as fractions of grid cell area) of different MTCs are all assigned to their youngest cohorts. During model spin-up forest woody mass will grow to exceed the thresholds of the first cohort, so that forests will move to the second cohort, and so on. At the end of spin-up, all forests thus end up in the oldest cohort of each MTC. The same case applies to herbaceous MTCs, given that cohort thresholds are properly defined (see more details in Sect. 2.2.3).
Natural forest mortality in ORCHIDEE could be either prescribed as a constant rate or dynamically simulated, but in the case of prescribed vegetation cover, mortality takes effect by reducing the amount of existing biomass only, with the coverage of the concerned forest patch being unchanged. Likewise, recruitment increases forest individual density and updates leaf age and other relevant variables, but again, forest coverage remains unchanged. These features are necessary, as the original ORCHIDEE model does not take into account forest demography. As explained in Krinner et al. (2005, p. 8), recruitment sapling biomass is only incorporated when the existing biomass is virtually zero while a larger-than-zero ground coverage is prescribed. These features remain the same when sub-grid cohorts are used, i.e. forest mortality or recruitment does not modify forest cohort ground coverage. In addition, forest mortality and subsequent regeneration due to forest fires are handled in a similar manner. ORCHIDEE-MICT has integrated a prognostic fire module to simulate open grassland and forest fires arising from both natural and anthropogenic ignitions (Yue et al., 2014). Other forest disturbances, such as wind-throw, diseases, and insect outbreaks, are not explicitly considered in ORCHIDEE-MICT. Because of these reasons, after the spin-up, the only way to create secondary cohorts in the model is through LUC.
When entering transient simulations with LUC, younger cohorts will begin to be created. From a modelling perspective, the oldest cohorts in ORCHIDEE-MICT are somewhat equivalent to the primary lands (especially, the oldest forest cohorts are equivalent to primary forests), and other younger cohorts are analogue to secondary lands.

Model developments to include gross land use change and forest harvest, with and without sub-grid cohorts
This section describes the implementation of gross LUC and forest harvest with sub-grid CFTs. We focus on the implementation with sub-grid cohorts because the same LUC process without cohorts could be simply treated as a particular case in which all MTCs have only one single cohort. The module interface is designed to receive forcing information on land area fluxes among four basic land cover types of forest (f), natural grassland (g), pasture (p), and cropland (c), taking into account the current LUC modelling landscape in DGVMs (as briefly reviewed in the introduction) and the availability of LUC reconstructions (e.g. Hurtt et al., 2011). The present developments are intended for the case in which changes in vegetation coverage are only driven by historical LUC activities and so there is no need to use the dynamic vegetation module of ORCHIDEE. This is different from the LUC implementation in JSBACH DGVM in Reick et al. (2013) in which a lot of effort has been devoted to reconciling the vegetation types in the forcing data (primary and secondary natural lands in the Land-Use Harmonization data set version 1 or LUH1 data) and the vegetation distributions simulated by the dynamic vegetation module of JSBACH. We focus on including sub-grid land cohorts in the model and implementing a set of hierarchical rules for which land cohorts are subjected to different LUC processes ( Table 2). The allocation of natural lands into forest versus grasslands in the model, and the reconciliation of LUH1 land cover distribution and model PFT map, are instead handled by independent preparations of reconstructed historical land cover map time series. In order to compare the simulation results from the gross LUC module with the original net-transition-only LUC module, we separate the gross LUC areas into two additive terms: net change equivalent to the original net transition (prescribed by the matrix M net ) and land turnover for the bidirectional equal land fluxes between any pair of land cover types (prescribed by the matrix M turnover ). Similarly, the forest harvest information is prescribed in a third matrix M harvest . For the moment, information for all three LUC types is provided as a fraction of grid cell area. This is a deliberate choice, mainly for the convenience of progressive stage-wise model development. We will come back to the influence of this choice within the land use decision contexts in Sect. 4.
The key processes of the gross LUC module with CFTs are shown in Fig. 4, comprising in total six steps. The LUC module is called at the first day of each year. Input data are the three matrices. M net and M turnover are both square matrices with a size of 4 by 4: , Step 1. Calculate the area for each cohort and each land cover type (i.e., construct the LUC hierarchy in Fig. 3b) Step 2. Construct the final LUC matrix of MnCFT, nMTC Input LUC matrices: Mnet, Mturnover, Mharvest Step 4. Collect legacy carbon pools and fluxes from all source CFTs that collectively comprise the area of the newly established youngest cohort of a given MTC Step 3. Wood product collection Step 5. Establish a proxy patch of the newly created youngest CFT, with its area being the sum of contributions by all source CFTs Step 6  where the element F i>j denotes the land flux from land cover type i to j , with i and j being elements of the vector of [f g p c] T . The diagonal elements correspond to land fractions intact from any land use transitions and are simply ignored in the LUC module. By definition, M turnover is a symmetric square matrix. M harvest is a matrix with only two elements: harvest area from primary and secondary forests. As explained in Sect. 2.1.3, the construction of CFTs within the model follows the model parameterization hierarchy shown in Fig. 3a. The cohort age subjected to LUC is one of the most important considerations in LUC decisions, especially in the context of land turnover and forest harvest. This necessitates a re-organization of the CFTs to derive the LUC hierarchy shown in Fig. 3b, in which Tier 2 information is about areas of different cohorts of the same land cover type, and Tier 3 remains on the level of CFTs. Thus, Step 1 in the LUC module (Fig. 4) is to construct the LUC hierarchy, i.e. to calculate within the model the areas of each cohort for each vegetation type.
When implementing LUC matrices, all information of land transitions between the four basic land cover types must first be downscaled on the cohort tier (i.e. decision on which cohort is subjected to LUC) and then on the CFT tier (i.e. how LUC-affected area is distributed among differ-ent comprising meta-classes within each cohort; refer also to Fig. 3b). This is achieved in Step 2 as shown in Fig. 4. Because all the newly established lands, regardless of their originating LUC process, must belong to the youngest CFT of the MTCs that comprise the target land cover type, the ultimate outcome of Step 2 is a single (large) matrix M nCFT, nMTC (nCFT: no. of CFTs; nMTC: no. of MTCs), which indicates the area transferred from each CFT to the youngest cohort of the MTC concerned. The rules to convert LUC matrices into components of M nCFT, nMTC depend on LUC types and will be explained in detail later. But as long as Step 2 is finished, the remaining steps are rather straightforward.
Step 3 handles forest wood collection (here "collection" rather than "harvest" is used, to avoid the confusion with forest wood harvest, which is a means of forest management), from forest being converted to other land cover types, and forestry harvest (forest remaining forest). We assume that a certain fraction of above-ground woody biomass (i.e. sapwood and heartwood) is lost as instant CO 2 flux into the atmosphere (i.e. due to on-site disturbance), and that the remaining wood is collected as wood product pools.
Step 4 involves the proper displacement of associated carbon stocks and fluxes from the donating CFTs to the newly established (youngest) cohorts of MTCs, after wood collection. Notably, the legacy carbon stocks in litter and soil collected from the donating CFTs are transferred to the newly established youngest CFTs. Then in Step 5, each youngest CFT cohort is established and initialized, with its fraction of grid-cell area being the sum of contributing areas given by each source CFT. Finally, in Step 6, a newly established cohort is merged with the existing youngest CFT cohort if there is one. When merging stocks or fluxes between the newly established and existing CFTs, an area-weighted mean approach is followed: where x is the variable in question (e.g. leaf biomass, soil carbon stock), x new and x existing are the values of the newly established patch and the existing patch before merging respectively, and x merged is the value of the composite patch after merging. The variables area new and area existing are patch areas of the newly established patch and the existing patch respectively. We now return to Step 2, explaining the different rules used to build the M nCFT, nMTC components for different LUC types. We start with M harvest by assuming that it precedes conversion of forest to other land cover types (i.e. land turnover or net LUC). As is explained, the LUC module is designed to receive externally prescribed harvest information, especially from the widely used LUH1 reconstruction (Hurtt et al., 2011), rather than to determine harvest volume internally within the model. The LUH1 distinguishes between harvests from primary and secondary forests and non-forest vegetation but in ORCHIDEE only harvest from forests is considered. The harvest information is provided as both for-est area and harvested biomass in LUH1. Here we used the area information (a deliberate choice that will be discussed in Sect. 4). Because of this, ensuring the consistency between the harvest area in the forcing and that being actually realized in the model is an important consideration. Moreover, as we want to compare simulated LUC impacts between the two model configurations with and without sub-grid cohorts, it is necessary to ensure that exactly the same LUC area is realized in both configurations. This involves a set of decision rules to properly allocate the prescribed harvest area into different forest cohorts ( Table 2).
Implementation of primary forest harvest is straightforward: we always start with the oldest cohort and move sequentially downwards to younger ones if older cohorts are exhausted until the prescribed harvest demand is fulfilled ( Table 2). For secondary forest harvest, we start with intermediately aged cohorts. But if the existing area of intermediately aged cohorts is not sufficient to fulfill the prescribed harvest area, we are left with two options to either search upwards for older cohorts or downwards for younger ones. We decide to first search upward and then search downward, if all cohorts older than the intermediate age still cannot fulfill the prescribed harvest demand (Table 2). This rule allows potential temporal changes in harvested area to be accommodated, as explained in Fig. 5. Under such a scheme, (1) at the very beginning (after spin-up) and before the existence of any secondary forests, harvest will start with the oldest cohort, i.e. corresponding to harvest of primary forests (sometimes, because of the inconsistency between the input harvest information and existing forest cohort structure in the model, secondary forest harvest could be prescribed for pixels in which only primary forests exist in the model).
(2) If harvest area of secondary forests remains stable, then as soon as sufficient intermediately aged cohorts are created via conversion of primary forest to regrowing younger cohorts, a corresponding stable rotation cycle would be maintained in the model as well. (3) If the harvest area increases, the upward searching would allow additional harvest of primary forests (i.e. area subject to the stable rotation is expanded). (4) If the harvest area decreases, moving cohorts from younger to older ones independent of any LUC activities would allow the restoration of older cohorts -e.g. a consequence of abandonment of forest management. (5) Finally, the downward searching for younger cohorts after exhausting all other older cohorts is solely to ensure the consistency between prescribed input harvest area and that actually realized in the model. Hence, this scheme is designed in order to faithfully implement the prescribed harvest areas in the model with an explicit consideration of forest successional states (i.e. primary or secondary). But when this is not possible because of inevitable mismatch between the model and forcing data, harvest areas of primary and secondary forests could mutually compensate for each other in the model to ensure that their prescribed total harvest area remains realized. A number of studies reported that fallow lengths for shifting cultivation could range from a few years to more than 50 years depending on different regions, with the majority being 10-40 years (Bruun et al., 2006;Mertz et al., 2008;Thrupp et al., 1997;van Vliet et al., 2012), and there is a tendency in reduction of fallow lengths possibly because of increased population pressure (van Vliet et al., 2012). Hurtt et al. (2011) assumed a mean residence time of 15 years for shifting cultivation for tropical regions in the LUH1 reconstruction data. Based on these reports, we assume forest clearance for shifting cultivation to occur primarily in secondary forests and treat it similarly as secondary forest harvest when allocating the prescribed LUC area into different cohorts ( Table 2). The only difference is that the destination land cover remains forest in the case of forest harvest but is agricultural land in the case of shifting cultivation. For all other land transfers in shifting cultivation (e.g. pasture to forest), we start exclusively from the oldest cohort and move downwards to younger ones (Table 2). For net LUC, priority is again given to older cohorts followed by younger ones (Table 2).
Finally, we still need to downscale the LUC area in each cohort to its component CFTs. This is done by allocating the LUC area in each cohort to its member CFTs in proportion to the existing area of each CFT.

LUC processes that remain unchanged in the model
ORCHIDEE simulates two wood product pools with a turnover length of 10 years and 100 years. Fractions of above-ground woody biomass as instant on-site losses (F instant ) and entering into the two wood product pools (F 10 yr , F 100 yr ) follow the values in the original nettransition-only LUC scheme (Piao et al., 2009a), as shown in Table 3. Other biomass compartments (i.e. leaves, fine roots, coarse roots, fruits, and reserve pool) are transferred to litter pools during forest harvest or deforestation. Carbon in the two wood product pools is then released into the atmosphere according to the pools' respective turnover time, and this flux contributes to the overall land carbon balance as a source term (see the next section).

Time
Always start with the Cohort 3 . Move sequen@ally to older cohorts when younger ones are exhausted. Further search back to Cohort 2 un@l Cohort 1 to ensure the consistency with prescribed harvest area. Cohorts keep growing and moving toward older ones; old forests will be restored when the harvest area decreases.  (1) first starts with intermediately aged cohort, then moves to older cohorts until the oldest one; (2) if the prescribed harvested area still cannot be satisfied, then the harvest will move back to the even younger cohorts (3) to the youngest one until the prescribed harvested area is fulfilled. Independent of the harvest activity is the movement of forests from younger cohorts to older ones because of growth (grey arrows). (b) Example of cohort dynamics along with temporal changes in the harvest area shown in the black curve: (1) before the onset of any harvest activity (i.e. after the model spin-up), only the oldest cohorts are available so harvest starts with the primary forest; (2) for a stable harvest area, a steady-state cycle involving only secondary forest is established (intermediate secondary cohorts being harvested is represented by the blue arrow, and younger growing cohorts are represented by grey arrows); (3) then with an increase in harvest area, more primary forests are harvested; (4) finally, in this example, the harvest area decreases, and older cohorts are restored.

Forest growth
Other processes relevant to LUC are left unchanged with the original model version. In particular, crop harvest is applied to cropland CFTs with a fraction of 45 % of biomass turnover being harvested in the model and exported outside the ecosystem (Piao et al., 2009a). Pasture CFTs are also harvested in the same fashion. Agricultural harvest and associated fluxes to the atmosphere through food consumption or livestock feeding are assumed to happen locally in the model during the same year of harvest, without considering spatial relocation through international trade. Fires are simulated with a prognostic module, but as explained in Sect. 2.1.3, fire disturbances do not lead to creation of young cohorts, but only their carbon consequences (e.g. emissions, vegetation mortality) are included.

Definition of land-use change emissions (E LUC ) and carbon flux sign convention
The land carbon balance simulated by ORCHIDEE-MICT v8.4.2 (i.e. net biome production or NBP), when land use change is included, is defined as where NPP is the net primary production, and all fluxes with F denoting outward carbon fluxes from the land system (they are assigned a negative sign following the ecosystem convention, indicating that carbon is lost from ecosystems), with F Inst for the instantaneous carbon flux during LUC (e.g. carbon release arising from site preparation, land-clearing burning), F Wood for the delayed carbon release due to wood product degradation, F HR for heterotrophic respiration from litter and soil organic carbon decomposition, F AH for agricultural harvest on both croplands and pastures, and F Pasture for carbon sources from pastures other than harvest, i.e. export of animal production and methane emissions (see Chang et al., 2015, for details). F Inst and F Wood are both fluxes on an annual timescale that depend only on wood mass at the time of forest clearing and the respective wood product degradation rates (see Sect. 2.1.5). F HR is simulated at a time step of 30 min and depends on soil temperature and moisture. F Fire is simulated with a prognostic fire module SPITFIRE (Yue et al., 2015). The LUC emissions (E LUC ) are quantified as the difference in simulated NBP between two paired simulations, with LUC (or a specific LUC process) included in one simulation but not the other one: where NBP LUC and NPB control are NBP simulated with and without LUC. A negative E LUC denotes a carbon source to the atmosphere, i.e. the ecosystem carbon sink is reduced because of LUC. This definition follows Pongratz et al. (2014, p. 178) and is also the same as used in TRENDY  simulations and the existing global carbon budget analysis (Le Quéré et al., 2016). As explained by Pongratz et al. (2014), such a definition quantifies the net LUC flux because it integrates both emissions to the atmosphere (e.g. deforestation) and uptakes by potentially recovering vegetation (e.g. agricultural abandonment). More specifically, this corresponds to the definition "D3" using uncoupled DGVM simulations in Pongratz et al. (2014, Eq. 15c, p. 187), which contains instantaneous fluxes, legacy fluxes, and "loss of additional sink (source) capacity (LOAS)". Instantaneous fluxes refer to the carbon emissions directly arising from LUC, often occurring within the first year since LUC (F Inst in our case). Legacy fluxes arise from the readjustment of carbon stocks to the new type of vegetation and/or the changes in management intensity over time , and LOAS refers to the carbon sinksource difference between the actual land cover after LUC and the otherwise potential one under environmental perturbations. All other flux terms on the right side of Eq. (3) except F Inst contribute to the legacy fluxes and LOAS. Here, as our model development mainly distinguishes the biomass carbon of secondary forests, it is expected that F Inst and F Wood will be the major fluxes to influence the simulated E LUC . To facilitate the demonstration of model behaviour, we refer to F Inst and F Wood collectively as LUC-associated direct fluxes and their variations will be examined in detail by using an idealized grid cell simulation.
The model developments presented here enable us to make two parallel simulations that include LUC: with and without sub-grid age dynamics. Their simulated E LUC can thus be compared to separate the effect of including sub-grid age dynamics. Henceforth, for briefness, we denote the simulation without sub-grid age dynamics as S ageless and the one with age dynamics as S age .

Idealized simulation on a single grid cell
We conducted an idealized grid cell simulation with prescribed land cover and LUC matrices to compare in detail the simulated carbon pools and fluxes between S age and S ageless . The geographical coordinates of the simulation site are 9.25 • S, 18.25 • E on a 0.5 • global grid, in the north of Angola, Africa, where the miombo woodlands are known to be subject to practices of shifting cultivation. The ESA CCI land cover map for the 5-year period of 2003-2007 (https://www. esa-landcover-cci.org/) shows a dominant fraction of tropical deciduous broadleaf forest for this grid cell. Hence, for the idealized experiment, the initial vegetation composition is prescribed as 85 % tropical deciduous broadleaf forests and 15 % C 4 croplands. As we will focus on the LUC impacts, other model forcings (climate, atmospheric CO 2 , etc.) are held constant, with climate input data recycling the year of 1901 (CRUNCEP v5.3.2 climate data) and atmospheric CO 2 concentration being fixed at 350 ppm. The model is tested for a hypothetical scenario of constant annual land turnover with 5 % of grid cell area between forest and C 4 cropland. Forest harvest of the same annual areal fraction is expected to have a largely similar impact. The spin-up was run for 450 years until biomass and soil C stocks reached equilibrium and the mean annual NBP was close to zero without including any LUC. Starting from the spin-up, a transient simulation with the prescribed LUC matrix was performed for 100 years.

Simulation over southern Africa
Subsequently, the model behaviour has been documented for a real-world case over the region of southern Africa (south from the Equator of the African continent). All three LUC types occurred historically in this region, making it ideal to demonstrate model behaviour regarding forest cohort dynamics as presented in Fig. 5. This regional simulation serves a single purpose -to further exemplify model features that cannot be sufficiently demonstrated over a grid cell.
The regional simulation is performed at 2 • resolution for 1501-2005. We used the land use reconstruction from LUH1 covering 1501-2013 (Hurtt et al., 2011, http://luh.umd.edu/ data.shtml#LUH1_Data) re-gridded from the original 0.5 • to a 2 • spatial resolution. From the LUH1 data set we derived the matrices of the three types of LUC: net land use change, land turnover, and wood harvest. Land turnover information is extracted from LUH1 as the minimum land fluxes between two vegetation types. Wood harvest from primary and secondary forests in LUH1 is used, while wood harvest from non-forest is not. Climate forcing data are from CRUNCEP-v5.3.2 at a 2 • resolution. For the spin-up, climate data were cycled from 1901 to 1910, with atmospheric CO 2 concentration fixed at the 1750 level (277 ppm). In the transient simulation, atmospheric CO 2 concentration began to increase in 1750 and climate data were varied starting in 1901. The dynamic vegetation module was turned off in order to apply the prescribed historical LUC. Factorial simulations are conducted to highlight changes in areas of different forest cohorts when different LUC processes are included, as shown in Table 4.
Each forest MTC has six CFTs to represent six cohorts. The woody mass thresholds are set in a way that they correspond roughly to the woody masses at ages of 3, 9, 15, 30, 50 years, and the mature or primary forest (with an age greater than 50 years) during the spin-up simulation for Cohort 1 to Cohort 6 respectively. The Cohort 3 with an age of 15 years is the primary target for secondary forest harvest and land turnover (or shifting cultivation), corresponding to the mean residence time of 15 years of shifting cultivation assumed in LUH1 data (Hurtt et al., 2011). We set two CFTs for each herbaceous MTC with a high and low soil carbon den- Table 4. Factorial simulations to examine forest cohort dynamics when including different LUC processes: net land use change, land turnover, and wood harvest. The plus signs (+) indicate that the corresponding processes (matrices) are included in the simulations. Only simulations with sub-grid age dynamics are carried out, with S0 age having no LUC activities to S3 age including all LUC processes.

Simulations
Net land Land Wood use change turnover harvest S0 age S1 age + S2 age + + S3 age + + + sity. The CFT thresholds of soil carbon stock are the same for all herbaceous MTCs. We first calculate the maximum soil carbon stock of all MTCs (including the forest ones) at the end of spin-up for each grid cell, and cohort thresholds are then taken as this maximum value and its 65 % value. Because the energy balance in ORCHIDEE-MICT is resolved for the average of all CFTs over a grid cell, and the hydrological balance is resolved for three sub-grid water columns (i.e. the water column of bare soil, forest, and herbaceous vegetation), we expect the factors influencing soil carbon decomposition (e.g. soil temperature, soil moisture) to have little variation among CFTs of the same MTC. This justifies the small number of herbaceous CFTs for the sake of computation efficiency. Overall, this feature of separating herbaceous MTCs into multiple cohorts is coded more as a placeholder for the current stage of model development rather than having solid scientific significance. Fully tracking soil carbon stocks of different vegetation types and their transient changes following LUC would require a much larger number of cohorts than used in this study.

Results
3.1 Grid cell simulations with and without sub-grid forest age dynamics 3.1.1 Temporal patterns of biomass carbon stock during the spin-up and transient simulations Figure 6a and b exhibit the evolution of above-and belowground biomass for both S ageless and S age simulations for the spin-up and transient simulation for a test grid cell located in Angola. The results for the S age simulation are shown for individual cohorts (Cohort 1 to Cohort 6 ). For this test an annual forest-cropland turnover of 5 % of the grid cell area was imposed. Figure 6c-h present changes in the ground fractional cover of different forest cohorts during the transient simula-tion. S ageless and S age share the same biomass accretion with time during the spin-up, but S age shows a succession of forest cohorts -with biomass moving from one cohort to the next (Fig. 6a, b). At the end of the spin-up, all biomass is found in Cohort 6 (i.e. the oldest cohort), with an initial forest cover of 85 %. More differences emerge when entering the transient simulation. Above-ground biomass in S ageless shows an initial sharp drop followed by a more gradual decline under constant land turnover because biomass of the single forest patch is constantly diluted by merging with the new forest patch with a low biomass, which is established as a result of land turnover (see also Fig. 1). Below-ground biomass, however, shows a corresponding initial drop but then slightly increases. Eventually, both above-and below-ground biomass stocks in S ageless reach a new equilibrium, which is lower than their values at the end of the spin-up. By contrast, in S age , the fraction of Cohort 6 declines with the start of the transient simulation because of conversion to cropland. This decline continues until the 12th year, after which the remaining Cohort 6 covers only 30 % of the grid cell (Fig. 6h). Younger cohorts are progressively created as forests restore after shifting agriculture abandonment, with the Cohort 1 (i.e. the youngest one) appearing during the initial 6 years after the start of LUC, after which its biomass is moved into Cohort 2 (Fig. 6c, d). Cohort 3 starts to appear at the 12th year when biomass in Cohort 2 moves into it. Then its coverage declines as this cohort, rather than Cohort 6 , is used as the source for shifting cropland, according to the model rule that secondary forest is taken prior to primary forest in land turnover (Fig. 5). After the initial 15 years (the rough age of Cohort 3 ), the fractions of Cohort 1 , Cohort 2 , and Cohort 3 reach a dynamic stable state. As Cohort 3 is being constantly converted to cropland, it has never developed into Cohort 4 or Cohort 5 . This explains the zero fractions of these two latter cohorts in Fig. 6f and g.
While the above-ground biomass continuously grows during the spin-up, the below-ground biomass first increases with time and then slightly declines before reaching the equilibrium value. This is because ORCHIDEE-MICT has a preferential allocation of NPP to below-ground sapwood when forests are young. The small decline in below-ground biomass in the late spin-up stage thus results from an almost stabilized NPP (under a big-leaf approximation), a reduced below-ground allocation, and a constant mortality. Because of this feature, ORCHIDEE-MICT creates a higher belowground biomass in younger forest cohorts (e.g. Cohort 2 and Cohort 3 in Fig. 6a and b) in S age than the single forest patch in S ageless in the transient simulation. However, the aboveground biomass in younger Cohort 2 and Cohort 3 in S age is lower than S ageless . The difference in biomass influences the simulated E LUC between these two simulations, as we will discuss in detail later.

LUC-associated direct carbon fluxes
As shown in Fig. 7a, in S ageless , the instantaneous carbon flux resulting from LUC follows the same temporal pattern as the above-ground biomass, as it is simulated as a fixed fraction of above-ground woody mass (sapwood and heartwood) (see Sect. 2.1.5). In S age , for the initial 12 years, Cohort 6 (undisturbed mature forest) is cleared, so that the instantaneous LUC carbon flux is higher than that in S ageless (where the biomass of the single forest patch is reduced immediately when the land turnover starts). After that, the instantaneous flux shows a stark drop in S age when Cohort 3 enters the land turnover. Since then until the end of the simulation, S age keeps a constantly lower instantaneous flux than S ageless because the LUC-perturbed equilibrium biomass is higher in the latter case (Fig. 6a). As a fixed 10 % of above-ground woody biomass enters the wood product pool with a 10-year turnover time, delayed carbon emissions from wood product degradation in both simulations are smaller than the instantaneous LUC carbon fluxes. They peak around the 12th year after LUC and remain stable afterwards (Fig. 7a). Overall, S age has a higher LUC-associated direct carbon flux than S ageless for the first 12 years and a lower one afterwards (Fig. 7a). The cross point for the cumulative LUC-associated direct fluxes equal in S age and S ageless is around the 20th year (Fig. 7b). When summing over the whole simulation period (100 years), the cumulative fluxes by S ageless are lower in S age by about 11 kg C m −2 , or ∼ 110 g C m −2 yr −1 (Fig. 7b) than S ageless .

LUC emission and its disaggregation into underlying component carbon fluxes
As defined in Eq. (4), the net LUC carbon emission (E LUC ) is diagnosed as the difference in NBP between the LUC simulation and the control one. Since NBP is further a composite flux determined by carbon uptake and releases (Eq. 3), the difference in E LUC age and E LUC ageless can be disaggregated into the effect of each underlying flux, which differs between the LUC simulation and the control simulation. First of all, S ageless (no age dynamics) simulates a larger magnitude (i.e. a larger absolute E LUC value) of mean annual E LUC than S age (with age dynamics) by about 26 g C m −2 yr −1 . Second, for both simulations, the simulated E LUC is an outcome of LUC-associated direct fluxes being compensated for by changes in other fluxes, all of which have an effect to reduce E LUC in this example: NPP, heterotrophic respiration, fire carbon emissions, and agricultural harvest.
NPP is higher in LUC simulations than in the control. This is because young forests are established in the former case (either by merging with existing forest patch or not), leading to a younger leaf age than in the control simulation, which is parameterized to have a higher photosynthetic capacity than older leaves in the model. This suggests the model can some- what integrate the effect of recovering young forests or intermediately aged forests with a higher productivity than the old-growth forests, as reported by Tang et al. (2014) using observation data.
Averaged over the LUC simulation period of 100 years, both S age and S ageless show lower heterotrophic respiration (F HR ) than the control. This is because the biomass stock is lower in the LUC simulations (despite a higher NPP, biomass turnover is accelerated due to site perturbation and wood collection in the process of clearing forest for cropland), causing less litter input and fewer soil carbon stocks (data not shown). The S age simulation shows a much smaller reduction in F HR , mainly because a higher below-ground litter is maintained, which results from a high below-ground litter input out of land turnover, driven by a high below-ground biomass, as explained in Sect. 3.1.1 (Fig. 6a).
Decreases in fire carbon emissions (F Fire , from prognostically simulated natural fires but not land-clearing fires) in the LUC simulations in contrast with the control are because the above-ground litter (dominant fuel for fires) is reduced by land turnover. Reductions in fire emissions, and reductions in heterotrophic respiration, are thus driven by Figure 8. Mean annual carbon flux differences between the LUC and control simulations over 100 years for an annual forestcropland turnover with 5 % of the grid cell area for two model configurations: without (blue) and with sub-grid age dynamics (green). Positive (negative) values indicate contributions to enhanced carbon sink (source) in LUC simulation compared to the control one, either by stronger (weaker) carbon uptake or smaller (stronger) carbon release. E LUC is shown as a negative value here, i.e. the LUC simulation has a lower NBP than the control one, indicating an effect of net carbon source by LUC. the same process, i.e. a reduction in above-ground standing biomass. LUC simulations also result in lower agriculture harvest (F AH , from cropland) although there is no change in the cropland area; this is due to lower biomass in young crop, as the crop harvest is assumed as a constant fraction of the biomass turnover (i.e. routine mortality) at a daily time step. The lower crop biomass in the LUC simulations here is because crop saplings are established on the first day of each calendar year, right before the seasonal biomass peak for the Southern Hemisphere, which artificially reduces the standing biomass.
Overall, the lower E LUC magnitude in S age is a result of the lower LUC-associated direct fluxes having been partly compensated for by a higher heterotrophic respiration. The relative magnitudes between E LUC age and E LUC ageless are dominated by these two fluxes, while other fluxes play a less important role.

Forest cohort area changes as a result of historical land use change over southern Africa
One of the useful features of our model development is to account for sub-grid forest age dynamics as a result of historical LUC, as illustrated in Fig. 9 for southern Africa. When no LUC is included (S0, the control simulation shown in light blue), the areas of all forest cohorts are constant over time. Except that younger cohorts have a very small area (< 0.1 Mkm 2 ) (Cohort 2 and Cohort 3 , probably due to improper cohort thresholds on a very small number of grid cells), almost all forests are found in Cohort 6 , which resembles mature forests. In S1 where only net LUC is considered, the area of Cohort 6 decreases consistently over time due to conversion of forest to other land cover types (Fig. 9a). Occasional increases in areas of other younger cohorts are also present, corresponding to the periods when forest gain happens due to net LUC, for instance, afforestation or reforestation around the 1700s and in the latter half of the 20th century (Fig. 9a). This is consistent with our rule that forest from abandonment of agriculture is established in the youngest cohort ( Fig. 5b -on the right), and progressive movement of forests from younger to older cohorts is also visible as the small waves in the curves of Fig. 9b-f. In the S2 simulation with both net LUC and land turnover, large areas of younger forests, in particular of Cohort 1 and Cohort 2 , begin to appear as a result of continual creation of forests from land turnover and subsequent moving of forests from Cohort 1 to Cohort 2 . Their temporal changes over time follow those of the forest area subject to land turnover, as shown in Fig. 9a (green dashed line). The area of Cohort 3 , however, does not see as much increase as in the two younger cohorts because forests of Cohort 3 are the primary target for clearance in land turnover and thus are incessantly converted back to (shifting) agriculture. As a result, about half of mature forests (Cohort 6 ) are left intact from LUC by 2005 (Fig. 9h). Most interestingly, when there is a decline in the turnover-impacted area around the 1700s (the green arrow in Fig. 9a), a corresponding decline in the area of Cohort 1 is found because these forests move into the next cohort. This pattern of decrease in the current cohort accompanied by the according increase in the next one then propagates into other older cohorts with time, which results in a delayed increase in Cohort 5 around the 1750s (Fig. 9g), and finally in Cohort 6 as well (but less prominent because of its already large area). This demonstrates the model feature of older forest recovery in case of decreased land turnover or wood harvest, as explained in Fig. 5b (right-hand side). Last, when we further include forest harvest in the S3 simulation, because wood harvest area only started to rise in the middle of 20th century, larger areas of Cohort 1 and Cohort 2 are found compared with S2 in the latter half of the last century, and forest area in Cohort 6 is accordingly lower, being converted to younger cohorts as a result of harvest.

Discussion
DGVMs, either used in an offline mode or coupled with climate models, are powerful tools to investigate the role of past and future LUC in the global carbon cycle perturbed by human activities Le Quéré et al., 2016). Therefore, a more realistic representation of LUC processes in these models is a scientific priority. We included two new features in ORCHIDEE-MICT v8.4.2: gross LUC and forest wood harvest, and sub-grid vegetation cohorts. In a recent review (Prestele et al., 2017), proper representation of gross LUC or sub-grid bidirectional land turnover has been identified as one of the three major challenges in implementing LUC in DGVMs for credible climate assessments, despite that these have already been pioneered by some models (Table 1). Large underestimation of LUC emissions would occur when gross LUC is ignored, as is shown by several model results reviewed in Arneth et al. (2017).
Shifting cultivation, or forest wood harvest, or more forest management in general, often involves a stable fallow length or rotation cycle, which involves secondary forests rather than primary ones. In tropical regions, fallow lengths in shifting cultivation range from 10 to 40 years (Bruun et al., 2006;Mertz et al., 2008;Thrupp et al., 1997;van Vliet et al., 2012), with a tendency of reduction in fallow length. In Latin American tropics, agricultural abandonment has already led to prominent growth of secondary forests (Chazdon et al., 2016;Poorter et al., 2016). Forest management, including wood harvest, is more common in temperate and boreal regions. In European forests, rotation lengths depend on tree species, regional climate, and management purposes, ranging from 8 to 20 years in coppicing systems in southern Europe to 80-120 years in northern countries . The prevalence of secondary forests associated with land use and LUC therefore calls for their representation in DGVMs, especially when modelling LUC.
To our knowledge, Shevliakova et al. (2009) performed the first study to include both sub-grid secondary lands and gross transitions in the LM3V model, but the number of PFTs and secondary land tiles are limited in their study (up to in total 12 secondary land tiles compared with 50 in our study). Stocker et al. (2014) included secondary land in LPX-Bern 1.0 but only one tile of secondary land is available. Yang et al. (2010) examined the contribution of secondary forests to terrestrial carbon uptake using the ISAM model by explicitly including secondary forest PFTs, but they did not include the dynamic clearing of secondary forests nor shifting cultivation in LUC. Therefore, none of these studies have included a dynamic decision rule regarding the ages of cohorts to be targeted in different LUC processes or the possibility of targeting different cohort ages in different geographical regions. ORCHIDEE-CAN is especially designed to address forest management and species change. Although certain LUC such as wood harvest and net land cover changes are included, a more comprehensive LUC scheme addressing gross change is missing .
The gross LUC combined with sub-grid cohorts presented here has shown some promising results. We first confirmed that including gross LUC leads to additional carbon emissions. However, these additional emissions tend to be overestimated when secondary forests are not explicitly accounted for. The idealized grid cell simulation explained the mechanism driving such overestimation in S ageless simulations well. The results presented here are closely linked with our model parameterization and in particular the decision rules regarding which forest cohorts to apply for specific LUC processes (Table 2). Land turnover and secondary forest harvest are parameterized to target intermediately aged cohorts as a priority. This is the core mechanism driving the lower LUC emissions when sub-grid forest age structure is accounted for.
As a preliminary effort to demonstrate the model behaviour, the land turnover parameterization is heavily tied with the input LUC forcing data (LUH1), so that the age of Cohort 3 (as the primary target for land turnover) is set as ∼ 15 years, following the assumed mean residence time of shifting cultivation in the LUH1 data set (Hurtt et al., 2011). The model simulations showed that this parameterization is crucial because it largely determines the rotation length in the model and consequently the amount of carbon stocks subjected to LUC and the difference in estimated LUC emissions between the two model configurations (S age and S ageless ). In this regard it should be noted that the information on rotational lengths of shifting cultivation or forest harvest is spatially unbalanced and that at present no systematic global compilation exists. The universal setting used in this study is due to the absence of such a compilation. In fact, because the thresholds in woody mass to distinguish forest cohorts could be configured via a spatial map in the model and such maps could vary among different years, and because the primary cohort target is not hard-coded and can be parameterized as well, it is rather straightforward to apply temporally and spatially different rotation lengths in the model. Such a feature is well considered in the model development design and could be tested when information on spatially and temporally explicit forest rotation lengths or associated biomass thresholds is available.
In the following paragraphs we will discuss the decisions that were marked as deliberate and their potential impacts on modelled LUC stocks and fluxes. First, the LUC module developed is intended for usage within DGVMs and forced with external data sets that provide information on land flows between different land cover types. It is not intended to supersede a LUC model per se, which simulates LUC using other available social and economic information such as population, food demand, wood demand, etc. (Hurtt et al., 2011). In this sense, the LUC module implementation has to inevitably take into account the details of information in forcing data that are available and to reconcile the potential mismatch between the model and forcing data. For example, the LUC module presented here can accommodate forest wood harvest from primary and secondary forests when these two sources are distinguished in the forcing data, but hierarchical decision rules are also made when the model and forcing data disagree (e.g. Fig. 5), for example, when prescribed secondary forest wood harvest can actually harvest a primary forest in the model if all younger cohorts are exhausted.
Second, because of this clearly defined border of the LUC module to use land areas as the input information, model output from ORCHIDEE-MICT can potentially disagree with the socio-economic information used to generate the LUC forcing data. For instance, crop yield simulated by ORCHIDEE may differ with that used to convert food demand and consumption to cropland area, so that simulated crop output or food production may disagree with historical food demand in the real world. The same applies to forestry wood production: simulated harvest wood volume might disagree with the wood volume actually used to generate the harvest area information -the harvested wood biomass information is provided in the LUH1 data set but not used as an input in the current stage of model development. This largely raises the issue of to what extent the information that drives LUC decisions can be internally integrated into DGVMs, for example, to directly use crop production, rather than cropland area, or wood volume, rather than forest harvest area as the model input. One potential obstacle is that statistical information (e.g. on wood volume demand) is often available on a regional basis (FAO global forest resource assessment, http://www.fao.org/ forest-resources-assessment/en/; Eurostat, http://ec.europa. eu/eurostat/data/database), and complex decision-making rules are needed to disintegrate such information on spatial grids that DGVMs are operated on. But in general and over the long term, land use or land management decisions need to be integrated directly into DGVMs. ORCHIDEE-CAN has integrated forest management decisions based on simulated tree diameters and stand density, so that harvested wood biomass is actually a model output that can be validated against historical statistical data (Naudts et al., 2016).
The developments presented here mainly build on a model structure that distinguishes differently aged cohorts. Nonetheless, we have built a better tool to address the impacts of historical LUC on carbon cycle and climate with these developments. Forest demographics, which are shown to have great impact on the current Northern Hemisphere carbon sink (Pan et al., 2011;Piao et al., 2009b), either as a result of active afforestation, agricultural abandonment, or natural regeneration, could then be explicitly investigated. These developments also make it possible to verify modelled global and regional forest age distribution using independent age information from either forest inventory or remote sensing. The model version used here has incorporated the developments in pasture and cropland modules (Chang et al., 2015;Wang et al., 2017). On a regional scale such as Europe, where the comprehensive forcing data are available, it is possible to go beyond the carbon emissions only by LUC activities, but also to include LUC-induced changes in emissions of other greenhouse gases such as methane and nitrogen oxide.

Conclusions
We have presented new developments made in a global vegetation model to include gross LUC and forest wood harvest, in combination with explicit representation of sub-grid forest age dynamics. Furthermore, a set of decision-making rules regarding the land cohorts to be targeted in different LUC processes have been implemented. The presented simulation results are specific of the ORCHIDEE-MICT model, but the methods are generic for other DGVMs. We demonstrated through an idealized pixel simulation that gross LUC leads to additional emissions but accounting for sub-grid land cohorts yields lower emissions than not. Over the region of southern Africa, the model is able to account for changes in different forest cohort areas along with the temporal changes in different LUC processes, including regrowth of old forests when LUC area decreases. Our developments provide the possibility of accounting for forest demography when evaluating LUC impacts on global carbon cycle and climate.
Code availability. The source code for ORCHIDEE-MICT version 8.4.2 is available online (https://forge.ipsl.jussieu.fr/orchidee/ browser/branches/ORCHIDEE-MICT/tags/ORCHIDEE_MICT_ GLUC_8.4.2) but its access is restricted to registered users. Requests can be sent to the corresponding author for a username and password for code access. ORCHIDEE-MICT is governed by the CeCILL license under French law and abides by the rules of distribution of free software. One can use, modify, and/or redistribute the software under the terms of the CeCILL license as circulated by CEA, CNRS, and INRIA at the following URL: http://www.cecill.info.
Data availability. Primary data and scripts used in the analysis and other supplementary information that may be useful in reproducing the authors' work can be obtained by contacting the corresponding author.
Competing interests. The authors declare that they have no conflict of interest.