Articles | Volume 13, issue 1
Geosci. Model Dev., 13, 185–200, 2020
Geosci. Model Dev., 13, 185–200, 2020

Development and technical paper 27 Jan 2020

Development and technical paper | 27 Jan 2020

Accounting for forest age in the tile-based dynamic global vegetation model JSBACH4 (4.20p7; git feature/forests) – a land surface model for the ICON-ESM

Accounting for forest age in the tile-based dynamic global vegetation model JSBACH4 (4.20p7; git feature/forests) – a land surface model for the ICON-ESM
Julia E. M. S. Nabel1, Kim Naudts1,a, and Julia Pongratz1,b Julia E. M. S. Nabel et al.
  • 1Max Planck Institute for Meteorology, 20146 Hamburg, Germany
  • anow at: Department of Earth Sciences, VU University Amsterdam, Amsterdam, the Netherlands
  • bnow at: Ludwig-Maximilians-Universität München, Munich, Germany

Correspondence: Julia E. M. S. Nabel (,


Natural and anthropogenic disturbances, in particular forest management, affect forest age structures all around the globe. Forest age structures in turn influence key land surface processes, such as photosynthesis and thus the carbon cycle. Yet, many dynamic global vegetation models (DGVMs), including those used as land surface models (LSMs) in Earth system models (ESMs), do not account for subgrid forest age structures, despite being used to investigate land-use effects on the global carbon budget or simulating biogeochemical responses to climate change. In this paper we present a new scheme to introduce forest age classes in hierarchical tile-based DGVMs combining benefits of recently applied approaches the first being a computationally efficient age-dependent simulation of all relevant processes, such as photosynthesis and respiration, using a restricted number of age classes and the second being the tracking of the exact forest age, which is a prerequisite for any implementation of age-based forest management. This combination is achieved by using the tile hierarchy to track the area fraction for each age on an aggregated plant functional type level, whilst simulating the relevant processes for a set of age classes. We describe how we implemented this scheme in JSBACH4, the LSM of the ICOsahedral Non-hydrostatic Earth system model (ICON-ESM). Subsequently, we compare simulation output to global observation-based products for gross primary production, leaf area index, and above-ground biomass to assess the ability of simulations with and without age classes to reproduce the annual cycle and large-scale spatial patterns of these variables. The comparisons show decreasing differences and increasing computation costs with an increasing number of distinguished age classes. The results demonstrate the benefit of the introduction of age classes, with the optimal number of age classes being a compromise between computation costs and error reduction.

1 Introduction

Land use, particularly forest management, substantially influences the age structure of global forests (Pan et al.2011; Erb et al.2017). More than 19 Mkm2 of forest area, i.e. about 15 % of global ice-free land, are under some kind of management (Luyssaert et al.2014), with 65 % being under regular harvest schemes and another 7 % being intensive plantations (Erb et al.2017). Often, management practices make use of rotation cycles, as is common in shifting cultivation (Boserup1966) or even-aged forest management strategies that historically were common in temperate forests and are still the dominant management type in boreal forests (Kuusela1994; Puettmann et al.2015; Kuuluvainen and Gauthier2018). Forest age structures are also influenced by other natural or anthropogenically caused disturbances such as fires, windthrows, droughts, pests, and insect outbreaks (e.g. Soja et al.2006; van Mantgem et al.2009; Dore et al.2010; Pan et al.2011, 2013).

Changes in forest age structure in turn influence biophysical and biogeochemical interactions with the atmosphere, through changes in land surface properties such as albedo and carbon uptake (e.g. Juang et al.2007; Dore et al.2010; Sun et al.2010; Kirschbaum et al.2011; Pan et al.2011, 2013; Poorter et al.2016; Erb et al.2017). Forest age structure changes can influence the susceptibility to the environment and to environmental changes. It is, for example, hypothesised that the response of forests to increasing atmospheric CO2 ceases as the forest matures because other resources than CO2, such as water and nutrients, become growth limiting (Körner2006).

It is crucially important to include forest age structures in estimates of the effects of land use on the global scale, primarily (1) due to the aforementioned large extent and substantial effects of forest undergoing changes in age structure; (2) due to the aim of global studies to include forest management effects in addition to anthropogenic land cover changes; and (3) because global studies usually only have a very coarse resolution. One example of such global studies is the estimate of global land-use emissions for the annual global carbon budget (Le Quéré et al.2018), which is conducted with dynamic global vegetation models (DGVMs). Here, 10 out of the 16 participating DGVMs account for wood harvesting. Other examples are studies estimating biogeochemical and/or biophysical interactions between the land surface and the climate system. Such studies are conducted with Earth system models (ESMs) including their land surface models (LSMs), many of which taking part in the coupled model intercomparison projects (CMIPs). Here, considerations of forest age structure might in particular be important for future scenarios that often include strong land-based mitigation measures, such as forest management and afforestation (e.g. in CMIP6's land use intercomparison project LUMIP; see Lawrence et al.2016). Global studies, in particular the computationally expensive ESM simulations, inevitably need to be conducted on coarse horizontal resolutions, typically only about 0.5 to 2. Land use will thus usually only happen on fractions of the grid cells, creating the need to represent subgrid forest age structures. The importance of subgrid forest age structures is also underlined by recent studies stressing the role of forest (re-)growth for the historical and future terrestrial carbon uptake (e.g. Kondo et al.2018; Krause et al.2018; Yao et al.2018) and by studies simulating smaller land-use emissions when accounting for secondary forests (e.g. Shevliakova et al.2009; Yue et al.2018a). Despite all this evidence, many DGVMs, and particularly also those used as LSMs in ESMs, do not account for subgrid forest age structures (Pongratz et al.2017).

There are categorically different approaches of how subgrid forest age structures can be implemented in DGVMs, depending on whether these models are individual- and cohort-based or tile-based. In the class of individual- and cohort-based models (referred to as vegetation demographic models in Fisher et al.2018), subgrid structures are inherently provided. Examples are ED-derivatives (Fisher et al.2015, 2018), LPJ-Guess (Smith et al.2001; Bayer et al.2017), and the SEIB-DGVM (Sato et al.2007). In the (larger) class of tile-based models (also referred to as area-based in Smith et al.2001), subgrid structures are not inherently provided. In these models each tile describes an average individual per plant functional type (PFT). Examples of this class of DGVMs are CABLE (Haverd et al.2018), Class-CTEM (Melton and Arora2016), ISAM (Yang et al.2010), JSBACH (Reick et al.2013; Mauritsen et al.2019), LM3 (Shevliakova et al.2009), LPX-Bern (Stocker et al.2014b), different versions of ORCHIDEE (Naudts et al.2015; Yue et al.2018b), and others. In our study we focus on the second class of DGVMs, as they are more commonly used as LSMs in ESMs. One reason that tile-based models are more commonly used is simply that they have lower computational costs. Another reason is their often historically conditioned top–down development, which facilitates a fully coupled execution within the corresponding ESM.

To expand tile-based DGVMs to represent subgrid forest age structures, two approaches have recently been developed. The more frequently applied approach has been to increase the number of tiles in such a way that a certain number of age classes or structurally similar stands can be distinguished. A pioneer study is found in the paper by Shevliakova et al. (2009), using the model LM3 with a fixed number of in total 12 secondary land tiles for all PFTs and a similarity-based merging of tiles in order to maintain the number of tiles despite further land use or disturbances. Comparably, ORCHIDEE-MICT (Yue et al.2018b) introduced a fixed number of six tiles per woody PFT, with tile merging upon exceeding predefined woody biomass boundaries. In ORCHIDEE-CAN three tiles per woody PFT with tile merging upon exceeding diameter thresholds have been used, while further within-stand structuring has been applied in each tile by accounting for a user-defined diameter distribution (Naudts et al.2015). An increase in tiles has also been chosen in ISAM (Yang et al.2010) and LPX-Bern (Stocker et al.2014a, b); in these models, however, only one additional tile per PFT has been introduced in order to distinguish primary and secondary vegetation. A common drawback of the hitherto existing implementations is the missing traceability of the actual age of the forests as soon as tiles are merged. Merging of tiles, however, is a necessity when the number of age classes is constrained by computation costs.

The alternative approach for extending the number of tiles to represent subgrid forest age structure in tile-based DGVMs is to keep the information about the forest structure in a separate module. For ORCHIDEE-FM (Bellassen et al.2010), for example, ORCHIDEE has been coupled to a forest management module (FMM). FMM takes the tile wood increment calculated in ORCHIDEE as input, allocates the increment to tracked individual trees, conducts self-thinning and forest management, and feeds back the leaf area index (LAI), biomass, and litter to the tile. A comparable coupling is described in Haverd et al. (2018), where the DGVM CABLE is coupled to the Population Orders Physiology (POP) module for woody demography and disturbance-mediated landscape heterogeneity (Haverd et al.2014). POP has a detailed description of the forest structure and simulates the growth of age or size classes of trees competing for soil resources and light. For each forest tile, POP gets the stem net primary production (NPP) from CABLE and returns woody vegetation height, mortality, and sapwood mass. Whilst such a use of a separate module principally enables tracking the exact age of the forest in a grid cell, it has the drawback of using average tile information to compute simulated processes, such as photosynthesis. Age dependencies of these processes can thus not be represented.

In this paper we bridge the two approaches for extending tile-based DGVMs to represent subgrid forest age in the sense that we present a way to trace the actual age of the forests in a grid cell despite following the first approach using a restricted number of additional tiles and thus required merges. The suggested approach is applicable for any tile-based DGVM, provided the tiles are structured in a hierarchical way. We describe the implementation of our approach in the DGVM JSBACH4 and use the new model version to conduct test simulations with different numbers of age classes and age distribution schemes. Subsequently, we compare the different simulation results against observation-based data to investigate the compromise between computation costs and error reduction.

2 Methods


The DGVM JSBACH4 is used as the LSM in the ICON-ESM (Giorgetta et al.2018). In addition, JSBACH4 is developed with a flexible interface, such that it is also usable within MPI-ESM1.2 (Mauritsen et al.2019) and as a standalone model driven by climate data. JSBACH4 is a re-implementation of JSBACH3, the original LSM used in MPI-ESM1.2 (Mauritsen et al.2019), but with a more flexible and extendable structure via a hierarchical representation of tiles (Fig. 1). The implementations described in this paper are based on the current version (4.20p7), which includes most of the processes implemented in JSBACH3, such as land physics, photosynthesis, carbon allocation, and natural disturbances, but is still lacking JSBACH3's representation of anthropogenic land cover change (Reick et al.2013). Furthermore, the current version does not provide an infrastructure for the horizontal exchange of properties among tiles, such as, for example, the movement of area fractions from one PFT to another.

As an important amendment to the current version (4.20p7), we ported a new JSBACH3 development, which we implemented in a recent independent study (Naudts et al.2020): while previous JSBACH3 versions assumed a PFT-dependent but constant maximum leaf area index (LAI), that is the LAI value that can maximally be reached at the peak of a season, Naudts et al. (2020) introduced a dependency of the maximum LAI on the available leaf biomass. Such a dependency is a prerequisite for simulating forest regrowth and thus for the introduction of age classes.

Figure 1The hierarchical tile structure of JSBACH4. In our study, the default tile structure of JSBACH4 (in black) has been extended by a variable number N of forest age classes (ACs) below each of the K forest plant functional types (PFTs; in blue).


Figure 2Unified Modelling Language diagram showing the relation between tiles, forest plant functional types (PFTs), and forest age classes (ACs) in JSBACH4–FF, together with a selection of example state variables and functions (in italics). Forest PFTs and ACs are distinct types of tiles, with each PFT having N associated ACs. Each tile hosts certain state variables, for example the grid-cell fraction that it covers, as well as functions, for example to navigate the tile hierarchy. Different types of tiles add further variables and functions. Tiles representing ACs host variables and functions required for processes calculated on the lowest level of the tile hierarchy, such as photosynthesis or carbon allocation. Tiles representing PFTs host variables to maintain meta-information, for example the fractPerAge vector, which contains the fraction covered by each age, i.e. one entry per year up to the maximum tracked age (maxAge). Furthermore, PFTs host functions altering more than one associated AC, for example ageing or harvest.



As outlined in the introduction, we aimed for a scheme that allows subgrid forest age structures to be introduced in hierarchical tile-based models in a computationally efficient way, i.e. using a restricted set of age classes, but nevertheless with exact tracing of the age of the forest. For our new JSBACH4 git feature “forests” (JSBACH4–FF) we took advantage of the existing hierarchical tile organisation of JSBACH4, which allows us to introduce different processes and associated state variables on different levels of the tile hierarchy (Fig. 2). Already existing processes which are specific to the development of an age class, such as photosynthesis and carbon allocation, are still executed on the leaves, i.e. now per age class. However, processes related to several age classes, such as natural disturbances and the newly introduced ageing and harvesting processes, are implemented on the PFT level, which we use to manage associated age classes and to maintain meta-information about the forest age structure.

In JSBACH4–FF we introduced a fixed user-defined number N of age classes preset in the configuration file for all forest PFTs (Fig. 1). In addition, an upper age bound per age class ACK (maxAK) as well as a total maximum age (maxAge) were introduced, which also have to be predefined. maxAge determines the oldest age up to which the age of an area is tracked, i.e. the length of the fractPerAge vector (Fig. 2). Area fractions with ages exceeding maxAge are not further distinguished and are refereed to as old-grown forest. For a maximum age of 150 years, for example, each forest PFT would contain a vector with a length of 150 to track the exact age of the entire forest area up to 150-year-old forest. Each ACK covers a certain interval of years [maxAK−1, maxAK) (Fig. 3), with the youngest AC (AC1) always covering the range of year 0 to 1 and the oldest ACN covering all forest older than maxAN−1, i.e. [maxAN−1, ), with maxAN-1 maxAge.

In JSBACH4–FF we implemented three processes on the PFT level which can cause shifts of area fractions from one AC to another (Fig. 3), each tracking changes in age fractions in the fractPerAge vector.


The newly implemented process of forest “ageing” happens annually: upon ageing each tracked forest fraction gets 1 year older. Yet, a shift from one age class to the next age class only happens for the area of the oldest age (maxAK-1-1) of an age class ACK−1; i.e. only the forest area which upon getting 1 year older exceeds the upper age bound maxAK−1 of the ACK−1 needs to be shifted into ACK. Thanks to the tracking of the age in the fractPerAge vector, the exact area fraction with age maxAK-1-1 is known.


In the current version, we implemented harvest as a clear-cutting of a certain fraction of an AC, which can happen annually. Harvest causes a shift of the harvested fraction of the affected AC to the youngest AC. Since the exact age of each forest fraction is tracked in the fractPerAge vector, age-based harvest rules can be specified.


Following JSBACH3 (Brovkin et al.2009) wind and fire disturbances can happen daily in JSBACH4 and are assumed to clear certain area fractions of vegetation. In JSBACH4–FF disturbances were partly re-implemented. While the calculation of the disturbed area is still conducted on each leaf, i.e. on each AC (calc_disturbed_area in Fig. 2), the movement of area between ACs is managed on the PFT level (manage_disturbed_area in Fig. 2). This separation was required since the state variables used to determine the disturbed area, such as available fuel, need to be derived on the lowest layer of the hierarchy. Disturbances are realised as shifts of fractions of affected ACs to the youngest AC (AC1).

Figure 3Visualisation of forest age-class (AC) boundaries, the fractPerAge vector, and processes causing shifts from one AC to another in JSBACH–FF. Each AC covers a certain interval of ages. The first AC contains all forest younger than 1 year; an arbitrary ACK covers [maxAK−1, maxAK), i.e. the ages in the right-open interval of the upper age boundary of the previous AC (maxAK−1) and its own upper age boundary (maxAK); finally the last class ACN covers all forest older than and including the upper age boundary of the next younger class (maxAN−1), with maxAN−1 being smaller or equal to the total maximum age (maxAge). The information on the area covered by the different ages (indicated in red) is tracked in the fractPerAge vector of the associated forest PFT. Three processes can lead to the movement of area fractions among ACs: ageing leads to the movement of the fraction exceeding the maximum age of an AC; harvest and disturbances lead to the movement of fractions to the first AC.


Since JSBACH4 so far does not provide the infrastructure for the horizontal exchange of properties among tiles (Sect. 2.1), we had to implement such an infrastructure in JSBACH–FF, redistributing area fractions and associated state variables such as the carbon pools or the maximum LAI. This redistribution has been realised in two steps: in the first step, each of the three processes described above determines required redistributions on the PFT level (schedule_state_changes in Fig. 2) according to Eqs. (1) and (2). Here, fa is the area fraction moved from one AC to another, e.g. upon ageing; VS is the value of the state variable of the source ACS and delta_VS and delta_VT are the scheduled changes of the state variable on the source ACS and the target ACT, respectively.


In the second step, the redistributions are applied to the ACs (apply_state_changes in Fig. 2) according to Eq. (3). Here, VAC represents the affected state variable of the age class AC, delta_VAC the scheduled change, fa is the incoming and fc the current area of the age class AC.

(3) V AC = V AC f c + delta _ V AC ( f c + f a )

2.3 Simulation set-up and measures of model performance

The main purpose of JSBACH is to conduct global applications, often in a mode coupled to an ESM. Therefore, we assess the ability of different set-ups without and with different numbers of age classes to reproduce the annual cycle and large-scale spatial patterns of gross primary production (GPP), forest LAI, and forest above-ground biomass (AGB) by comparing simulated variables against global observation-based products for different seasons and regions. We conducted simulations following a protocol described below (Sect. 2.3.2), with the aim to simulate actual 2010 forest age distributions and forest states. Forest GPP, forest LAI, and forest AGB simulated for 2001 to 2010 were compared to GPP, LAI, and AGB data based on observations using a normalised root mean squared error (NRMSE; Sect. 2.3.5). In addition, we created Taylor diagrams for each variable, season, and region (see Figs S5.5–S5.11 in the Supplement).

2.3.1 Observation-based data

We used 2010 MODIS LAI (Myneni et al.2002) and GPP data obtained from machine learning methods trained on flux-tower measurements (Tramontana et al.2016). These LAI and GPP datasets had already been mapped to JSBACH's forest PFTs in a previous study (see Supplement Nyawira et al.2016). From these datasets seasonal means were calculated and expressed per forest area by dividing them by the sum of the forest cover fractions used for the JSBACH4–FF simulations (Sect. 2.3.2). The AGB per forest area (Avitabile et al.2016) was downloaded from the GEOCARBON data portal (, last access: 19 October 2017) and remapped to T63 using the conservative remapping operator of the climate data operators (CDOs, version 1.9.5). Figures S4.2–S4.4 in the Supplement show maps of the preprocessed observation-based data.

Table 1Conducted simulations with number of age classes and applied age distribution scheme. The “+1” in the number of age classes refers to the youngest age class, which always covers the years 0 to 1 in JSBACH4–FF.

a EAS: equal age spacing. b IAS: increasing age spacing.

Download Print Version | Download XLSX

2.3.2 General simulation set-up

We conducted simulations with JSBACH4 (4.20p7) feature/forest in standalone mode hosted within the MPI-ESM environment (see Supplement Sect. S1 for further information). We used JSBACH's default set-up comprising 12 PFTs, of which 4 are of a forest type: tropical evergreen and deciduous forest (TE and TD) and extratropical evergreen and deciduous forest (ETE and ETD). Simulations started in 1860 from scratch, i.e. with empty vegetation carbon stocks, and were run up to 2010. Empty carbon stocks are a simplification used in the absence of global knowledge on the state of the forest in 1860 but have no influence on our results since in simulations with JSBACH4 (4.20p7) LAI, GPP, and AGB only depend on the age since the last clearing event, not on the history before that. The starting date of 1860 was chosen such that it covers at least one full cycle of regrowth, as the oldest age resolved in the simulations matches that of the observation-based data (Poulter et al.2018). We used T63 resolution (192 longitudes × 96 latitudes; 1.9×1.9), a climate forcing based on GSWP3 (Kim et al.2012), and CO2 from the collection of greenhouse gas concentrations for CMIP6 (Meinshausen et al.2017). To obtain forest age distributions comparable to those observed for 2010 (Poulter et al.2018), harvest was conducted following prescribed maps (Sect. 2.3.4) and natural disturbances were switched off in order to not additionally alter forest age. The simulations were conducted with a static land-use map for 2010, based on TRENDYv4 JSBACH3 output (Le Quéré et al.2015). The TRENDYv4 JSBACH3 simulation started from a potential vegetation map extrapolated from remote sensing (Pongratz et al.2008) and was forced by the Land-use Harmonization dataset LUH1 (Hurtt et al.2011). We conducted simulations with different numbers and distributions of age classes (Sect. 2.3.3). All simulations were conducted on Mistral, the high-performance computing system of the German Climate Computing Center (DKRZ), using an identical number of CPUs.

2.3.3 Simulated number of age classes and age distribution schemes

Table 1 lists the conducted simulations. We used different numbers of age classes and two different age distribution schemes described below. The selected numbers of age classes are arbitrary; however, the finest resolution into 15+1 age classes was motivated by the age map discretisation in Poulter et al. (2018) using 15 age classes that cover 10 years each. In addition to simulations with age classes, we performed one simulation only using PFTs, i.e. without age classes.

The two applied age distribution schemes are defined as follows:


The equal age-spacing (EAS) distribution scheme spreads the age classes evenly over the age space. For example, a maximum traced age of 150 years distributed evenly over 10+1 age classes (EAS11 in Table 1) results in age classes covering 15 years with the following upper age bounds: 1, 16, 31, 46, …, 136, . This distribution scheme was motivated by the equal spacing applied in the forest age map by Poulter et al. (2018).


The increasing age-spacing (IAS) distribution scheme uses an increasing age range, i.e. younger age classes cover smaller intervals in the age space than older age classes. The upper age bound of a forest age class K (uLimK) is defined following Eq. 5.


with maxAge being the maximum age and N being the number of age classes. A maximum age of 150 years distributed with IAS over 10+1 age classes (IAS11 in Table 1) results in age classes with the following upper age bounds: 1, 3, 8, 16, 26, 39, 55, 74, 95, 119, 2. This second distribution scheme was motivated by the fact that young forests usually have larger incremental changes in most variables than old ones (see, e.g., Amiro et al.2006; Martínez-Vilalta et al.2007; Leslie et al.2018).

Both distribution schemes are applied in a static way, i.e. the age-class boundaries do not change during runtime. Figure 4 shows the division into age classes resulting for the different simulation set-ups for an example grid cell in Canada.

Figure 4Division into age classes (ACs) for the different simulations listed in Table 1 (EAS: equal age spacing; IAS: increasing age spacing). Purple lines mark the upper age boundary of each age class. The blue bars show the relative fraction for each year resulting for an example cell in Canada (lat 47.5639, long 286.875) in the simulation year 2010. Note that no harvest was conducted in the final simulation year 2010; therefore the smallest age class is empty.


2.3.4 Harvest maps

Harvest maps were derived such that the observed 2010 forest age distribution given by Poulter et al. (2018) is reached in the final simulation year 2010 for simulations using age classes. The observed forest age map of course not only reflects forest harvest, but all processes influencing the age of a forest, i.e. also natural disturbances and anthropogenic land cover change. Because assigning the observed age structure to forest harvesting vs. natural disturbances vs. anthropogenic land cover change would come with uncertainties and is not relevant for our study (as only the affected fraction of an age class matters, independently of the underlying causes), we apply only forest harvest to obtain the observed age distribution.

The map by Poulter et al. (2018) provides a grid with a 0.5 resolution of the global forest age distribution of four forest PFTs: needleleaf evergreen (NE) and needleleaf deciduous (ND), as well as broadleaf evergreen (BE) and broadleaf deciduous (BD). The map uses a discretisation into 15 age classes, covering 10 years each, with the last class containing all area with an age >140 years. In a preprocessing step, the map was remapped to T63 using the conservative remapping operator of the CDOs. Subsequently, the area sums of the two evergreen and the area sums of the two deciduous PFTs from Poulter et al. (2018) were used to derive the age-class maps for JSBACH's evergreen and deciduous PFTs, respectively, following Eq. (6):

(6) cf _ age i _ PFT k = ( cf _ age i _ N + cf _ age i _ B ) cf _ PFT k i = 1 15 ( cf _ age i _ N + cf _ age i _ B ) ,

where i is one of the 15 age classes from Poulter et al. (2018) and k refers to one of the four forest PFTs of JSBACH. N (needleleaf) and B (broadleaf) refer to the PFTs in Poulter et al. (2018), where either both are evergreen (in case PFTk is evergreen) or both are deciduous PFTs (in case PFTk is deciduous).

From these age maps, we derived harvest maps for each simulation year such that the simulated age distribution in simulations with age classes conforms with the observed one in 2010, assuming that the fractions given by Poulter et al. (2018) are equally distributed over the 10 years covered by each age class (see Supplement Sect. S2 for more details). In the first (1860) and in the last simulation year (2010) no harvest was conducted.

In different simulation types – with or without age classes – the same harvest maps were used, but different forest management schemes were applied. In simulations with age classes, clear-cutting according to the fractions in the harvest map was taken from the oldest age class. In the simulation without age classes, the PFT simulation, we used the same harvest fractions as in the simulations with age classes, but harvest was applied as done in JSBACH3 (Reick et al.2013), i.e. by diluting the wood carbon of the harvested PFT tile.

Table 2Selected regions used for the comparison of simulation results and observation-based data and their latitudinal boundaries and indices.

Download Print Version | Download XLSX

2.3.5 NRMSE

We calculated the area-weighted root mean squared error (RMSE) according to Eq. (7) based on difference maps between “OBS”, the observation-based data (Sect. 2.3.1), and “SIM”, the results of each simulation (see Table 1). The RMSE was calculated for 2001–2010 simulation output averages, separately for each variable “V” and three selected regions “R”. Each selected region defines a latitudinal band, including all forested land on a subset of the 96 latitudes and the entire 192 longitudes (see Table 2 for the regions, their latitudinal boundaries and the latitude indices b1 and b2). For GPP and LAI the four seasons “S” (DJF, MAM, JJA, SON) were calculated separately. The RMSE for each variable, region, and season was subsequently normalised with the range (max–min) observed for that variable, region, and season (Eq. 8).

(7) RMSE V , S , R = k = 1 192 m = b 1 b 2 ( OBS V , S , lon ( k ) , lat ( m ) - SIM V , S , lon ( k ) , lat ( m ) ) 2 AREA lon ( k ) , lat ( m ) AREA R 192 ( b 2 - b 1 + 1 )
(8) NRMSE Max - Min , V , S , R = RMSE V , S , R Max ( OBS V , S , R ) - Min ( OBS V , S , R )

To more easily assess changes in performance when increasing the number of age classes the different NRMSEMax−Min values were subsequently aggregated per variable by averaging over the regions (for AGB) and in addition over the seasons (for GPP and LAI), using equal weights.

2.3.6 Computational costs

In addition to determining the NRMSE for different variables, we also determined the computation costs of the different set-ups. We calculated the average CPU time recorded for the simulation years 2001–2010. Whilst absolute computation times are of less interest here, particularly since JSBACH4 is still highly under development and currently does not reach the targeted performance, relative differences among the set-ups depict the costs of the introduction of subgrid forest age structures.

3 Results and discussions

Having forest age classes in JSBACH4–FF facilitates a finer discretisation in each grid cell and is a precondition for any implementation of age-based forest management. The number of age classes in JSBACH4–FF is flexible, and in the following we describe the evaluation of simulation results using different numbers of age classes and age distribution schemes and discuss the compromise between computation costs and error reduction, when selecting a certain number of age classes (Sect. 3.1). Subsequently, we more closely examine differences between an example simulation with a selected number of age classes and a simulation only using PFTs, i.e. without age classes, to investigate the benefits of having age classes in JSBACH4–FF (Sect. 3.2). Finally, we discuss assets and drawbacks of alternative schemes introducing age classes in tile-based DGVMs (Sect. 3.3).

3.1 Evaluation

In this section we use the NRMSEMax−Min for different regions/seasons as an aggregated measure to compare the different simulation set-ups. A closer examination between a simulation with and without age classes including a spatially explicit comparison follows in Sect. 3.2.

Introducing age classes improves the comparison to observation-based data for nearly all compared variables, regions, and seasons (Fig. 5), with the only notable exception of the AGB in the boreal region, where the PFT simulation was more similar to the observation-based data than the simulations with age classes (Fig. 5c). For most comparisons, the NRMSEMax−Min indicates a small but distinct improvement over not representing a forest age structure for all simulated numbers of age classes and both age distribution schemes.

Figure 5Evaluation of the conducted simulations (Table 1) with observation-based data by means of the normalised root mean squared error (NRMSEMax−Min, Sect. 2.3.5). Depicted are calculated NRMSEMax−Min (Sect. 2.3.5) values for each simulation for the gross primary production (GPP; panel a), the leaf area index (LAI; panel b), and the above-ground biomass (AGB; panel c). The NRMSEMax−Min is calculated as the root mean squared error of observation-based data and simulation results, normalised with the range (max–min) of the according variable for each of the selected regions (Table 2) and for LAI and GPP also for each of the four seasons.


Averaging the NRMSEMax−Min, giving each region and each season the same weight, results in an NRMSEMax−Min decreasing with the number of age classes for GPP and LAI (Fig. 6a and b) but saturating for larger numbers of age classes. This shape holds for all regions, with a faster decrease and an earlier saturation for the Northern Hemisphere temperate and tropical regions than for the boreal region (Fig. S3.1a–f). The NRMSEMax−Min for AGB shows a slowly saturating increase with the number of age classes for the boreal region (Fig. S3.1g) and only small differences among the different numbers of age classes in the Northern Hemisphere temperate and the tropical regions (Fig. S3.1h and i). The observed increase in NRMSEMax−Min for the boreal AGB is due to an increased underestimation when accounting for more young forest, as is also discussed below (Sect. 3.2). Apart from the boreal AGB comparison, all comparisons show a smaller NRMSEMax−Min for simulations using the IAS distribution scheme (Figs. 6 and S3.1), i.e. a distribution applying an increasing age space (visualised in Fig. 4). This decrease in NRMSEMax−Min is due to the finer discretisation of younger age classes which have fast-changing LAI and GPP, which saturates for older age classes (see, e.g., Fig. 7 for LAI). In summary, a finer discretisation, particularly of the younger age classes, is leading to values closer to the observation-based data, even though the benefit of increasing the number of age classes is slowly saturating towards larger numbers of age classes (Fig. 6).

Figure 6Change in the normalised root mean squared error (NRMSEMax−Min, Sect. 2.3.5) and in the CPU time when increasing the number of age classes. Panels (a) to (c) show the averaged NRMSEMax−Min (see also Fig. 5). Averaging has been conducted giving equal weights to all selected regions (Sect. 2.3.5) for AGB (c) and in addition to all four seasons for GPP and LAI (a, b). Figure S3.1 in the Supplement shows the same data separately for each region. Panel (d) shows the computation time required per simulation year averaged over the years 2001–2010.


We performed the averaging of the NRMSEMax−Min to more easily assess differences among the simulations performed (Table 1). For this, we equally weighted the selected regions because we wanted to equally account for these regions, which strongly differ in simulated PFTs and land–atmosphere interactions. Alternatively, we could have weighted the regions by area, which would have led to an increasing weight of the tropical region and thus to an earlier saturation of the NRMSEMax−Min with increasing age classes.

Comparisons of required CPU times show a near-linear increase with an increased number of age classes (Fig. 6d) and neither a difference between the two age distribution schemes nor a striking offset as compared to the PFT simulation. A near-linear increase with an increasing number of age classes was expected, since the processes requiring most of the computing time, such as the calculation of photosynthesis, carbon allocation, and respiration, are conducted on the age classes. The absence of a striking offset comparing the PFT simulation with the age-class simulations indicates that the introduced organisational overhead on the PFT level in simulations with age classes, i.e. tracing of the exact forest age and redistributions of area fractions and other state variables among tiles, does not dominate the computation times.

As expected, the optimal number of age classes is a compromise between computation costs and reduction of the error, which is a logical and commonly observed aspect when dealing with discretisation in models (see, e.g., Nabel2015; Fisher et al.2018). In the end, the choice of the number of age classes to be used in a JSBACH4–FF simulation will depend on the application. Simulations comparing different forest management regimes in detail might, for example, aim for a fine discretisation, while more general simulations covering long time spans might tend to aim for fewer age classes. For the remaining parts of this paper, one set-up has been selected as an illustrative example: IAS11 (see Table 1), i.e. the simulation with 10+1 age classes and the age distribution scheme with increasing age space. This set-up is a compromise between the error reduction for GPP and LAI comparisons on the one hand and CPU time on the other. However, the main findings will not depend on the exact number of age classes selected, particularly not as long as they are in the saturation part of the decreasing function regarding GPP and LAI comparisons.

3.2 On the benefit of having age classes

The evaluation with observation-based GPP, LAI, and AGB data showed that simulations with age classes were closer to the observation-based data for nearly all comparisons (Figs. 5 and 6). Spatially explicit comparisons of the results from the PFT simulation and observation-based data (“OBS-PFT” in Figs. S4.2–S4.4, column 2) indicate several areas of underestimation (red) and of overestimation (blue) for all variables. In Fig. 7 we compare results of a JSBACH4–FF simulation with age classes, a simulation only using PFTs, as representative of a DGVM without forest age structures, and the observation-based data of spring (MAM) LAI. The comparison is done for illustrative grid points that were selected to cover areas of both over- and underestimation and to represent different typical land-use histories or forest management regimes, resulting in different age distributions: a grid point with an age distribution matching historically continuous clear-cuttings and some more recent changes in land-use intensity in Germany (Fig. 7a); a grid point with uniform age distribution resulting from a continuous, steady clear-cutting in Finland (Fig. 7b); a grid point with untouched old-grown forest on the one hand and young managed forest on the other hand in India (Fig. 7c); a grid point with intensive harvest or disturbances in the south-east of the US (Fig. 7d); a heavily deforested example in east South America resulting nearly exclusively in young forest (Fig. 7f); a grid cell with recent afforestation in China (Fig. 7e); and a grid cell with predominantly old-grown forests in central Africa (Fig. 7g). In general, the simulation with age classes results in smaller GPP, LAI, and AGB values (Figs. 7 and S4.2–S4.4, column 3), which is expected, since GPP, LAI, and AGB are non-linearly increasing and saturating with age (see, e.g., Fig. 7). Therefore, a harvested age-less forest in the PFT simulation has higher values for these variables than a fraction-weighted average of an age-structured forest in the same grid cell in the simulation with age classes (Fig. 7). Since the simulation with age classes generally results in smaller GPP, LAI, and AGB values, overestimations can be alleviated, causing a decrease in the NRMSEMax−Min, while underestimations can become more severe, causing an increase in the NRMSEMax−Min. The comparison of the differences between observation-based data (OBS) and the PFT simulation results on the one hand and OBS and the IAS11 simulation results on the other hand accordingly shows a higher similarity in several areas where the PFT simulation indicated overestimation (areas which are blue in columns 2 and 4 in Figs. S4.2–S4.4) and less similarity in some areas with underestimation (areas which are red in columns 2 and 4 in Figs. S4.2–S4.4). Figure 7 shows several grid-point examples with increased underestimations of spring LAI (Fig. 7a and e), reduced overestimations (panel b and c), and grid points where the previous overestimation is now replaced by a slight underestimation (Fig. 7d and f). Globally, reduced overestimations become particularly visible for LAI in the east of South America, and for several seasons also for example over China, North America, and Europe (Fig. S4.3d, h, i, p). For GPP (Fig. S4.2d, h, i, p) and AGB (Fig. S4d) the pattern is more mixed, with reduced overestimations particularly in the east of North America and China and partly for the east of South America. In addition, there are several areas of under- and overestimation which are very similar in the two simulations (areas coloured in column 2 and white in column 4 in Figs. S4.2–S4.4). These are particularly areas with predominantly old-grown forests, i.e. without a distinct age structure, such as central Africa, the central Amazon, and Siberia, where the PFT and the age-class simulation led to similar results (see, e.g., Fig. 7g). In summary, simulations using age classes led to a decrease in the simulated GPP, LAI, and AGB values due to their non-linear increase with a saturation for older ages. This caused a decrease in the NRMSEMax−Min in areas where the PFT simulation was biased high and an increase in the NRMSEMax−Min in areas where the PFT simulation was biased low. Thus, if such a forest age structure would be implemented in a DGVM being predominately biased low, the difference to the observation-based data could increase. In this context, caveats regarding the observation-based data themselves also need to be raised. A known caveat regarding MODIS LAI data is the problem of reflectance saturation in dense canopies making the reflectance insensitive to changes in LAI (Myneni et al.2002). This problem which is particularly relevant to the tropical region could lead to a general high bias of the model compared to the observation-based data. However, since this problem is more typical for denser, old-grown forests, this high bias would also occur in simulations with age classes. Regarding the GPP data from Tramontana et al. (2016), a recent study by Besnard et al. (2018) criticised that the applied empirical upscaling techniques do not directly consider forest age, making it unclear how well they can capture age-related dynamics. In their study, Besnard et al. (2018) advocate the development of alternative global datasets considering forest age as a predictor.

Importantly, besides the error reduction observed for JSBACH4–FF simulations, the newly implemented forest age structure adds further functionality to JSBACH4–FF. It facilitates keeping the coarse resolution required in ESM simulations while nevertheless capturing some of the sub-grid-scale heterogeneity that is important to better resolve several of the simulated processes. Furthermore, the forest age structure is a precondition for any implementation of forest management regimes while simultaneously accounting for differences in the productivity and the standing stocks. The grid-point examples shown in Fig. 7 highlight the relevance of a distinction of age classes, since they demonstrate the non-linear relationship between LAI and forest age. A similar relationship can be found for AGB and GPP. Consequently, the ability to distinguish age classes enables a more accurate simulation of the biogeochemical consequences of land use and particularly prescribed harvest regimes. For example, harvesting of younger age classes will lead to lower land-use emissions, as also described in other studies (e.g. Shevliakova et al.2009; Yue et al.2018a). Similarly, being able to distinguish forest age classes will also affect biophysical land–atmosphere interactions, since younger forests, for example, have lower LAI and higher albedo (e.g. Bright et al.2013). A constantly thinned age-less forest will therefore always lead to a lower albedo than a young forest regrowing after clear-cutting.

Figure 7Example grid points comparing 2001–2010 mean spring leaf area index (MAM LAI) from the simulations without (PFT) and with age classes (IAS11) to observation-based data (OBS). The map in the centre shows the difference of differences between the observation-based data and the simulations (abs(OBS-PFT)-abs(OBS-IAS11)); i.e. it shows where the results from the simulation with age classes (IAS11) deviate less (blue) or more (red) from the observation-based data than the PFT simulation results (see also Figs. S4.2–S4.4, column 4). Dashed lines in the map mark the three selected regions (see Table 2). The plots (a)(g) show the LAI of selected PFTs (ETD: extratropical deciduous; ETE: extratropical evergreen; TD: tropical deciduous; TE: tropical evergreen) as well as their according area fractions per age class and per year at the labelled grid points. Centre latitude, longitude, and grid-cell cover fraction (cf) of the depicted PFT are indicated. The x axis reflects the age from 0–151 (purple) with the age classes (black) indicated at the age centres. The two right y axes represent the bars: depicted are the 2010 area fractions relative to the area of the depicted PFT. Blue bars are per age class (black y axes) and depict the fraction of each age class (i.e. one bar per age class); the yellow framed purple bars depict the fraction of each age (i.e. one bar per year). The left y axis depicts the LAI. Stars mark the simulated LAI per age class and the lines the LAI of the depicted PFT – blue dashed line: IAS11 simulation; black line: PFT simulation; green line: 2010 value from the observation-based data. Note: (1) the age-class LAI is only depicted for age classes having non-zero fractional cover over the whole time span 2001–2010 (this is not the case for the age classes 9 and 10 in panels c, f, and g); (2) age and age-class fractions of classes 2–8 in panel (g) are very small and therefore not visible above the x axis; (3) since we did not apply any harvest in the final simulation year 2010, the first year and accordingly the youngest age class are always empty.

3.3 Limitations and alternative schemes

In the previous sections we compared a JSBACH4–FF simulation when only using PFTs with simulations including forest age classes and discussed associated trade-offs and benefits. In this section we discuss limitations and advantages of the applied and of alternative schemes.

Since JSBACH is a tile-based DGVM, the introduction of an individual- or cohort-based approach as used in some other DGVMs (e.g. Sato et al.2007; Fisher et al.2015; Bayer et al.2017) would be very complex. Regarding forest age structures these models have the essential advantage of naturally providing forest demography (Fisher et al.2018). Due to their complexity, however, they are less commonly used as fully coupled LSMs for ESMs. Being fully coupled with an ESM, however, is one major aspect and purpose of JSBACH, which historically has been part of the MPI-ESM (Mauritsen et al.2019) and is now also part of the ICON-ESM (Giorgetta et al.2018).

For tile-based DGVMs, there is at least one option mentioned in the literature that provides an alternative to simply increasing the number of tiles: the coupling of a separated module dealing with the woody demography (see, e.g., Bellassen et al.2010; Haverd et al.2018). On the one hand, this approach shares the advantage with individual- or cohort-based DGVMs that it provides a forest demography and thus principally enables the tracking of forest age. On the other hand, this approach has the important limitation of still calculating key land surface processes at the aggregated tile level; i.e. in this approach, processes such as photosynthesis and respiration, are not computed for separate age classes. This restriction impairs the calculation of biogeochemical and biophysical interactions, due to the non-linearity of forest growth and the associated non-linear relationships of those key processes with forest age (e.g. as depicted for spring LAI in Fig. 7). This limitation can only be avoided by increasing the number of tiles.

Building on the approach of increasing the number of tiles, the scheme suggested in this paper adds an important benefit of the alternative schemes by explicitly tracking forest age. It thereby enables the implementation of age-based forest management schemes that historically were common in temperate forests and are still the dominant management type in boreal forests (Kuusela1994; Pan et al.2011; Puettmann et al.2015; Kuuluvainen and Gauthier2018). Another advantage of the explicit tracking concerns the discretisation error. While the presented approach does require frequent area-weighted merges in order to maintain a limited number of age classes, it only requires shifting the actually affected parts of an age class and not entire age classes or “cohorts”, as was common in previous applications (e.g. Shevliakova et al.2009; Yue et al.2018b). Upon ageing, for example, in our approach only those fractions of an age class will be shifted that are actually at the age limits of an age class. An important restriction of the approach presented in this paper is that it is only applicable for a DGVM with a tile hierarchy and would not be applicable for a DGVM with a flat tile organisation, such as JSBACH3, since the different layers of the tile hierarchy are used to introduce different processes. In a DGVM with a flat tile organisation, the PFT level associated with the age classes on the leaf level of the tile hierarchy would be missing, which we use for the management of the forest age classes and for the tracking of the exact forest age.

With regard to previous studies that increased the number of tiles in order to introduce a more detailed representation of the forest state, our evaluation indicates that the number of additional tiles used in previous applications might have been too few. Solely separating primary and secondary forests (e.g. Yang et al.2010; Stocker et al.2014a, b) or introducing only a few age classes or cohorts (e.g. Shevliakova et al.2009; Yue et al.2018b) might not be sufficient to discretise non-linear relationships with forest age (see, e.g., Fig. 7, and also Fig. 6), at least not on the coarse resolutions that are common in global model studies dealing with human land use (e.g. Le Quéré et al.2018).

In this paper, we presented two different approaches to distribute the age space onto the available age classes: the equal age distribution scheme EAS, which spreads the age classes evenly, and IAS, a scheme that increases the age space with increasing age. The evaluation indicated the second approach to be superior to the first (Fig. 6), which can be explained by the finer discretisation of younger age classes that more accurately resolves the steep part of the non-linear age-dependent relationship of GPP, LAI, and AGB (see, e.g., Fig. 7). There are, however, other possible age distribution schemes. One could, for example, use smaller age classes for old ages in addition to the smaller age classes used for young ages in the IAS scheme. With such a scheme, one could better cover age-related declines as described in Zaehle et al. (2006) and Bellassen et al. (2010). Another possibility would be to replace the static distribution schemes that are equally applied to all grid cells with a dynamic scheme creating individual distributions for each grid cell. In such a dynamic scheme, age classes could be defined depending on the demand for each grid cell, with merging based on similarity criteria (see, e.g., Shevliakova et al.2009); i.e. age classes sharing similar values for a selected variable (e.g. GPP) could be merged creating space for new age classes covering an age space with less similar values. Such an approach could potentially further decrease the discretisation error, particularly for cells with only infrequent disturbances or harvest events. A drawback could be an increase in the organisational overhead caused by the similarity tests required for each merging step. However, the additional computational effort is not expected to be very large, considering that currently the organisational overhead seems to be very small (no striking offset as shown in Fig. 6d) and particularly since in the current set-up dynamic merges would only be required once a year.

4 Summary and outlook

In this paper we described a new scheme to introduce forest age structure in a hierarchical tile-based DGVM and presented its implementation in JSBACH4. JSBACH4–FF allows key land surface processes to be simulated in dependence of forest age and, simultaneously, to trace the exact forest age, which is a precondition for any implementation of age-based forest management schemes in JSBACH4–FF.

JSBACH4 itself is still highly under development regarding infrastructure and processes integrated from JSBACH3. In the version used for this paper (4.20p7), the representation of natural and anthropogenic land cover change, in particular, has not yet been ported from JSBACH3. Upon implementation, new processes will have to be integrated with the age structure. In addition, other developments would be desirable: harvest, for example, has so far only been implemented as area clear-cutting, following the implementation of other disturbances in JSBACH3 (see Brovkin et al.2009). For a representation of different forest management strategies including intermediate thinning before a final felling, an implementation of forest thinning would be required (Otto et al.2014; Naudts et al.2015). Anthropogenic thinning could be implemented in JSBACH4–FF by keeping the number of individuals as a state variable for each age class that is manipulated upon thinning, with anthropogenic thinning overruling the already implemented self-thinning.

Despite planned and potential extensions, together with the newly implemented age classes, JSBACH4–FF already now provides a valuable tool to study forest management effects, particularly due to its integration with the ICON-ESM.

Code and data availability

The hosting MPI-ESM model version (MPI-ESM 1.2.01p1) is made available under a version of the MPI-M Software License Agreement and can be obtained after registration from (last access: 20 December 2019, Max Planck Institute for Meteorology2019). Data and scripts used in the analysis, the JSBACH4 (4.20p7; git feature/forests) code, a patch to the hosting MPI-ESM required to run JSBACH4–FF, and other supplementary information are archived by the Max Planck Institute for Meteorology (, last access: 20 December 2019, Max Planck Society2019) and can be obtained by contacting


The supplement related to this article is available online at:

Author contributions

JN, KN, and JP initiated the implementation of age classes in JSBACH. JN planned and conducted the implementation, designed and performed the simulations, and wrote the first draft of the paper. All co-authors contributed to the analysis and edited the paper.

Competing interests

The authors declare that they have no conflict of interest.


All simulations have been conducted at the German Climate Computing Center (DKRZ; allocation bm0891). The authors would like to thank the two anonymous reviewers for their helpful comments on the discussion paper, Stiig Wilkenskjeld for reviewing the paper prior to submission, Veronika Gayler for support regarding JSBACH3, Reiner Schnur for support regarding JSBACH4, and Thomas Raddatz for discussions regarding the implementation of the age classes.

Financial support

This research has been supported by the DFG (grant no. PO1751/1-1).

The article processing charges for this open-access
publication were covered by the Max Planck Society.

Review statement

This paper was edited by Carlos Sierra and reviewed by two anonymous referees.


Amiro, B., Orchansky, A., Barr, A., Black, T., Chambers, S., III, F. C., Goulden, M., Litvak, M., Liu, H., McCaughey, J., McMillan, A., and Randerson, J.: The effect of post-fire stand age on the boreal forest energy balance, Agr. Forest Meteorol., 140, 41–50,, 2006. a

Avitabile, V., Herold, M., Heuvelink, G. B. M., Lewis, S. L., Phillips, O. L., Asner, G. P., Armston, J., Ashton, P. S., Banin, L., Bayol, N., Berry, N. J., Boeckx, P., de Jong, B. H. J., DeVries, B., Girardin, C. A. J., Kearsley, E., Lindsell, J. A., Lopez-Gonzalez, G., Lucas, R., Malhi, Y., Morel, A., Mitchard, E. T. A., Nagy, L., Qie, L., Quinones, M. J., Ryan, C. M., Ferry, S. J. W., Sunderland, T., Laurin, G. V., Gatti, R. C., Valentini, R., Verbeeck, H., Wijaya, A., and Willcock, S.: An integrated pan-tropical biomass map using multiple reference datasets, Glob. Change Biol., 22, 1406–1420,, 2016. a

Bayer, A. D., Lindeskog, M., Pugh, T. A. M., Anthoni, P. M., Fuchs, R., and Arneth, A.: Uncertainties in the land-use flux resulting from land-use change reconstructions and gross land transitions, Earth Syst. Dynam., 8, 91–111,, 2017. a, b

Bellassen, V., Maire, G. L., Dhôte, J., Ciais, P., and Viovy, N.: Modelling forest management within a global vegetation model – Part 1: Model structure and general behaviour, Ecol. Model., 221, 2458–2474,, 2010. a, b, c

Besnard, S., Carvalhais, N., Arain, M. A., Black, A., de Bruin, S., Buchmann, N., Cescatti, A., Chen, J., Clevers, J. G. P. W., Desai, A. R., Gough, C. M., Havrankova, K., Herold, M., Hörtnagl, L., Jung, M., Knohl, A., Kruijt, B., Krupkova, L., Law, B. E., Lindroth, A., Noormets, A., Roupsard, O., Steinbrecher, R., Varlagin, A., Vincke, C., and Reichstein, M.: Quantifying the effect of forest age in annual net forest carbon balance, Environ. Res. Lett., 13, 124018,, 2018. a, b

Boserup, E.: The conditions of agricultural growth: the economics of agrarian change under population pressure, vol. 4, Earthscan, London, 1966. a

Bright, R. M., Astrup, R., and Strømman, A. H.: Empirical models of monthly and annual albedo in managed boreal forests of interior Norway, Climatic Change, 120, 183–196,, 2013. a

Brovkin, V., Raddatz, T., Reick, C. H., Claussen, M., and Gayler, V.: Global biogeophysical interactions between forest and climate, Geophys. Res. Lett., 36, L07405,, 2009. a, b

Dore, S., Kolb, T. E., Montes-Helu, M., Eckert, S. E., Sullivan, B. W., Hungate, B. A., Kaye, J. P., Hart, S. C., Koch, G. W., and Finkral, A.: Carbon and water fluxes from ponderosa pine forests disturbed by wildfire and thinning, Ecol. Appl., 20, 663–683,, 2010. a, b

Erb, K.-H., Luyssaert, S., Meyfroidt, P., Pongratz, J., Don, A., Kloster, S., Kuemmerle, T., Fetzel, T., Fuchs, R., Herold, M., Haberl, H., Jones, C. D., Marín-Spiotta, E., McCallum, I., Robertson, E., Seufert, V., Fritz, S., Valade, A., Wiltshire, A., and Dolman, A. J.: Land management: data availability and process understanding for global change studies, Glob. Change Biol., 23, 512–533,, 2017. a, b, c

Fisher, R. A., Muszala, S., Verteinstein, M., Lawrence, P., Xu, C., McDowell, N. G., Knox, R. G., Koven, C., Holm, J., Rogers, B. M., Spessa, A., Lawrence, D., and Bonan, G.: Taking off the training wheels: the properties of a dynamic vegetation model without climate envelopes, CLM4.5(ED), Geosci. Model Dev., 8, 3593–3619,, 2015. a, b

Fisher, R. A., Koven, C. D., Anderegg, W. R. L., Christoffersen, B. O., Dietze, M. C., Farrior, C. E., Holm, J. A., Hurtt, G. C., Knox, R. G., Lawrence, P. J., Lichstein, J. W., Longo, M., Matheny, A. M., Medvigy, D., Muller-Landau, H. C., Powell, T. L., Serbin, S. P., Sato, H., Shuman, J. K., Smith, B., Trugman, A. T., Viskari, T., Verbeeck, H., Weng, E., Xu, C., Xu, X., Zhang, T., and Moorcroft, P. R.: Vegetation demographics in Earth System Models: A review of progress and priorities, Glob. Change Biol., 24, 35–54,, 2018. a, b, c, d

Giorgetta, M. A., Brokopf, R., Crueger, T., Esch, M., Fiedler, S., Helmert, J., Hohenegger, C., Kornblueh, L., Köhler, M., Manzini, E., Mauritsen, T., Nam, C., Raddatz, T., Rast, S., Reinert, D., Sakradzija, M., Schmidt, H., Schneck, R., Schnur, R., Silvers, L., Wan, H., Zängl, G., and Stevens, B.: ICON-A, the Atmosphere Component of the ICON Earth System Model: I. Model Description, J. Adv. Model. Earth Sy., 10, 1613–1637,, 2018. a, b

Haverd, V., Smith, B., Nieradzik, L. P., and Briggs, P. R.: A stand-alone tree demography and landscape structure module for Earth system models: integration with inventory data from temperate and boreal forests, Biogeosciences, 11, 4039–4055,, 2014. a

Haverd, V., Smith, B., Nieradzik, L., Briggs, P. R., Woodgate, W., Trudinger, C. M., Canadell, J. G., and Cuntz, M.: A new version of the CABLE land surface model (Subversion revision r4601) incorporating land use and land cover change, woody vegetation demography, and a novel optimisation-based approach to plant coordination of photosynthesis, Geosci. Model Dev., 11, 2995–3026,, 2018. a, b, c

Hurtt, G. C., Chini, L. P., Frolking, S., Betts, R. A., Feddema, J., Fischer, G., Fisk, J. P., Hibbard, K., Houghton, R. A., Janetos, A., Jones, C. D., Kindermann, G., Kinoshita, T., Klein Goldewijk, K., Riahi, K., Shevliakova, E., Smith, S., Stehfest, E., Thomson, A., Thornton, P., van Vuuren, D. P., and Wang, Y. P.: Harmonization of land-use scenarios for the period 1500–2100: 600 years of global gridded annual land-use transitions, wood harvest, and resulting secondary lands, Climatic Change, 109, 117,, 2011. a

Juang, J.-Y., Katul, G., Siqueira, M., Stoy, P., and Novick, K.: Separating the effects of albedo from eco-physiological changes on surface temperature along a successional chronosequence in the southeastern United States, Geophys. Res. Lett., 34, L21408,, 2007. a

Kim, H., Yoshimura, K., Chang, E., Famiglietti, J. S., and Oki, T.: Century long observation constrained global dynamic downscaling and hydrologic implication, AGU Fall Meeting Abstracts, GC31D-02, 2012. a

Kirschbaum, M. U. F., Whitehead, D., Dean, S. M., Beets, P. N., Shepherd, J. D., and Ausseil, A.-G. E.: Implications of albedo changes following afforestation on the benefits of forests as carbon sinks, Biogeosciences, 8, 3687–3696,, 2011. a

Kondo, M., Ichii, K., Patra, P. K., Poulter, B., Calle, L., Koven, C., Pugh, T. A. M., Kato, E., Harper, A., Zaehle, S., and Wiltshire, A.: Plant Regrowth as a Driver of Recent Enhancement of Terrestrial CO2 Uptake, Geophys. Res. Lett., 45, 4820–4830,, 2018. a

Körner, C.: Plant CO2 responses: an issue of definition, time and resource supply, New Phytol., 172, 393–411,, 2006. a

Krause, A., Pugh, T. A. M., Bayer, A. D., Li, W., Leung, F., Bondeau, A., Doelman, J. C., Humpenöder, F., Anthoni, P., Bodirsky, B. L., Ciais, P., Müller, C., Murray-Tortarolo, G., Olin, S., Popp, A., Sitch, S., Stehfest, E., and Arneth, A.: Large uncertainty in carbon uptake potential of land-based climate-change mitigation efforts, Glob. Change Biol., 24, 3025–3038,, 2018. a

Kuuluvainen, T. and Gauthier, S.: Young and old forest in the boreal: critical stages of ecosystem dynamics and management under global change, Forest Ecosystems, 5, 26,, 2018. a, b

Kuusela, K.: Forest resources in Europe 1950–1990, vol. 1, Cambridge University Press, 1994. a, b

Lawrence, D. M., Hurtt, G. C., Arneth, A., Brovkin, V., Calvin, K. V., Jones, A. D., Jones, C. D., Lawrence, P. J., de Noblet-Ducoudré, N., Pongratz, J., Seneviratne, S. I., and Shevliakova, E.: The Land Use Model Intercomparison Project (LUMIP) contribution to CMIP6: rationale and experimental design, Geosci. Model Dev., 9, 2973–2998,, 2016. a

Le Quéré, C., Moriarty, R., Andrew, R. M., Canadell, J. G., Sitch, S., Korsbakken, J. I., Friedlingstein, P., Peters, G. P., Andres, R. J., Boden, T. A., Houghton, R. A., House, J. I., Keeling, R. F., Tans, P., Arneth, A., Bakker, D. C. E., Barbero, L., Bopp, L., Chang, J., Chevallier, F., Chini, L. P., Ciais, P., Fader, M., Feely, R. A., Gkritzalis, T., Harris, I., Hauck, J., Ilyina, T., Jain, A. K., Kato, E., Kitidis, V., Klein Goldewijk, K., Koven, C., Landschützer, P., Lauvset, S. K., Lefèvre, N., Lenton, A., Lima, I. D., Metzl, N., Millero, F., Munro, D. R., Murata, A., Nabel, J. E. M. S., Nakaoka, S., Nojiri, Y., O'Brien, K., Olsen, A., Ono, T., Pérez, F. F., Pfeil, B., Pierrot, D., Poulter, B., Rehder, G., Rödenbeck, C., Saito, S., Schuster, U., Schwinger, J., Séférian, R., Steinhoff, T., Stocker, B. D., Sutton, A. J., Takahashi, T., Tilbrook, B., van der Laan-Luijkx, I. T., van der Werf, G. R., van Heuven, S., Vandemark, D., Viovy, N., Wiltshire, A., Zaehle, S., and Zeng, N.: Global Carbon Budget 2015, Earth Syst. Sci. Data, 7, 349–396,, 2015. a

Le Quéré, C., Andrew, R. M., Friedlingstein, P., Sitch, S., Hauck, J., Pongratz, J., Pickers, P. A., Korsbakken, J. I., Peters, G. P., Canadell, J. G., Arneth, A., Arora, V. K., Barbero, L., Bastos, A., Bopp, L., Chevallier, F., Chini, L. P., Ciais, P., Doney, S. C., Gkritzalis, T., Goll, D. S., Harris, I., Haverd, V., Hoffman, F. M., Hoppema, M., Houghton, R. A., Hurtt, G., Ilyina, T., Jain, A. K., Johannessen, T., Jones, C. D., Kato, E., Keeling, R. F., Goldewijk, K. K., Landschützer, P., Lefèvre, N., Lienert, S., Liu, Z., Lombardozzi, D., Metzl, N., Munro, D. R., Nabel, J. E. M. S., Nakaoka, S., Neill, C., Olsen, A., Ono, T., Patra, P., Peregon, A., Peters, W., Peylin, P., Pfeil, B., Pierrot, D., Poulter, B., Rehder, G., Resplandy, L., Robertson, E., Rocher, M., Rödenbeck, C., Schuster, U., Schwinger, J., Séférian, R., Skjelvan, I., Steinhoff, T., Sutton, A., Tans, P. P., Tian, H., Tilbrook, B., Tubiello, F. N., van der Laan-Luijkx, I. T., van der Werf, G. R., Viovy, N., Walker, A. P., Wiltshire, A. J., Wright, R., Zaehle, S., and Zheng, B.: Global Carbon Budget 2018, Earth Syst. Sci. Data, 10, 2141–2194,, 2018. a, b

Leslie, A., Mencuccini, M., and Perks, M.: Preliminary growth functions for Eucalyptus gunnii in the UK, Biomass and Bioenergy, 108, 464–469,, 2018. a

Luyssaert, S., Jammet, M., Stoy, P. C., Estel, S., Pongratz, J., Ceschia, E., Churkina, G., Don, A., Erb, K., Ferlicoq, M., Gielen, B., Grünwald, T., Houghton, R. A., Klumpp, K., Knohl, A., Kolb, T., Kuemmerle, T., Laurila, T., Lohila, A., Loustau, D., McGrath, M. J., Meyfroidt, P., Moors, E. J., Naudts, K., Novick, K., Otto, J., Pilegaard, K., Pio, C. A., Rambal, S., Rebmann, C., Ryder, J., Suyker, A. E., Varlagin, A., Wattenbach, M., and Dolman, A. J.: Land management and land-cover change have impacts of similar magnitude on surface temperature, Nat. Clim. Change, 4, 389–393, 2014. a

Martínez-Vilalta, J., Vanderklein, D., and Mencuccini, M.: Tree height and age-related decline in growth in Scots pine (Pinus sylvestris L.), Oecologia, 150, 529–544,, 2007. a

Mauritsen, T., Bader, J., Becker, T., Behrens, J., Bittner, M., Brokopf, R., Brovkin, V., Claussen, M., Crueger, T., Esch, M., Fast, I., Fiedler, S., Fläschner, D., Gayler, V., Giorgetta, M., Goll, D. S., Haak, H., Hagemann, S., Hedemann, C., Hohenegger, C., Ilyina, T., Jahns, T., Jimenez de la Cuesta Otero, D., Jungclaus, J., Kleinen, T., Kloster, S., Kracher, D., Kinne, S., Kleberg, D., Lasslop, G., Kornblueh, L., Marotzke, J., Matei, D., Meraner, K., Mikolajewicz, U., Modali, K., Möbis, B., Müller, W. A., Nabel, J. E. M. S., Nam, C. C. W., Notz, D., Nyawira, S.-S., Paulsen, H., Peters, K., Pincus, R., Pohlmann, H., Pongratz, J., Popp, M., Raddatz, T., Rast, S., Redler, R., Reick, C. H., Rohrschneider, T., Schemann, V., Schmidt, H., Schnur, R., Schulzweida, U., Six, K. D., Stein, L., Stemmler, I., Stevens, B., von Storch, J.-S., Tian, F., Voigt, A., de Vrese, P., Wieners, K.-H., Wilkenskjeld, S., Winkler, A., and Roeckner, E.: Developments in the MPI-M Earth System Model version 1.2 (MPI-ESM 1.2) and its response to increasing CO2, J. Adv. Model. Earth Sy., 11, 998–1038,, 2019. a, b, c, d

Max Planck Institute for Meteorology: MPI-ESM 1.2.01p1, available at:, last access: 20 December 2019. a

Max Planck Society: Publication repository, available at:, last access: 20 December 2019. a

Meinshausen, M., Vogel, E., Nauels, A., Lorbacher, K., Meinshausen, N., Etheridge, D. M., Fraser, P. J., Montzka, S. A., Rayner, P. J., Trudinger, C. M., Krummel, P. B., Beyerle, U., Canadell, J. G., Daniel, J. S., Enting, I. G., Law, R. M., Lunder, C. R., O'Doherty, S., Prinn, R. G., Reimann, S., Rubino, M., Velders, G. J. M., Vollmer, M. K., Wang, R. H. J., and Weiss, R.: Historical greenhouse gas concentrations for climate modelling (CMIP6), Geosci. Model Dev., 10, 2057–2116,, 2017. a

Melton, J. R. and Arora, V. K.: Competition between plant functional types in the Canadian Terrestrial Ecosystem Model (CTEM) v. 2.0, Geosci. Model Dev., 9, 323–361,, 2016. a

Myneni, R., Hoffman, S., Knyazikhin, Y., Privette, J., Glassy, J., Tian, Y., Wang, Y., Song, X., Zhang, Y., Smith, G., Lotsch, A., Friedl, M., Morisette, J., Votava, P., Nemani, R., and Running, S.: Global products of vegetation leaf area and fraction absorbed PAR from year one of MODIS data, Remote Sens. Environ., 83, 214–231,, 2002. a, b

Nabel, J. E. M. S.: Upscaling with the dynamic two-layer classification concept (D2C): TreeMig-2L, an efficient implementation of the forest-landscape model TreeMig, Geosci. Model Dev., 8, 3563–3577,, 2015. a

Naudts, K., Ryder, J., McGrath, M. J., Otto, J., Chen, Y., Valade, A., Bellasen, V., Berhongaray, G., Bönisch, G., Campioli, M., Ghattas, J., De Groote, T., Haverd, V., Kattge, J., MacBean, N., Maignan, F., Merilä, P., Penuelas, J., Peylin, P., Pinty, B., Pretzsch, H., Schulze, E. D., Solyga, D., Vuichard, N., Yan, Y., and Luyssaert, S.: A vertically discretised canopy description for ORCHIDEE (SVN r2290) and the modifications to the energy, water and carbon fluxes, Geosci. Model Dev., 8, 2035–2065,, 2015. a, b, c

Naudts, K., Nabel, J. E. M. S., Sabot, M., and Pongratz, J.: Impact of age-dependent wood harvest on land surface properties, in preparation, 2020. a, b

Nyawira, S. S., Nabel, J. E. M. S., Don, A., Brovkin, V., and Pongratz, J.: Soil carbon response to land-use change: evaluation of a global vegetation model using observational meta-analyses, Biogeosciences, 13, 5661–5675,, 2016. a

Otto, J., Berveiller, D., Bréon, F.-M., Delpierre, N., Geppert, G., Granier, A., Jans, W., Knohl, A., Kuusk, A., Longdoz, B., Moors, E., Mund, M., Pinty, B., Schelhaas, M.-J., and Luyssaert, S.: Forest summer albedo is sensitive to species and thinning: how should we account for this in Earth system models?, Biogeosciences, 11, 2411–2427,, 2014. a

Pan, Y., Chen, J. M., Birdsey, R., McCullough, K., He, L., and Deng, F.: Age structure and disturbance legacy of North American forests, Biogeosciences, 8, 715–732,, 2011. a, b, c, d

Pan, Y., Birdsey, R. A., Phillips, O. L., and Jackson, R. B.: The Structure, Distribution, and Biomass of the World's Forests, Annu. Rev. Ecol. Evol. S., 44, 593–622,, 2013. a, b

Pongratz, J., Reick, C., Raddatz, T., and Claussen, M.: A reconstruction of global agricultural areas and land cover for the last millennium, Global Biogeochem. Cy., 22, GB3018,, 2008. a

Pongratz, J., Dolman, H., Don, A., Erb, K.-H., Fuchs, R., Herold, M., Jones, C., Kuemmerle, T., Luyssaert, S., Meyfroidt, P., and Naudts, K.: Models meet data: Challenges and opportunities in implementing land management in Earth system models, Glob. Change Biol., 24, 1470–1487,, 2017. a

Poorter, L., Bongers, F., Aide, T. M., Almeyda Zambrano, A. M., Balvanera, P., Becknell, J. M., Boukili, V., Brancalion, P. H. S., Broadbent, E. N., Chazdon, R. L., Craven, D., de Almeida-Cortez, J. S., Cabral, G. A. L., de Jong, B. H. J., Denslow, J. S., Dent, D. H., DeWalt, S. J., Dupuy, J. M., Durán, S. M., Espírito-Santo, M. M., Fandino, M. C., César, R. G., Hall, J. S., Hernandez-Stefanoni, J. L., Jakovac, C. C., Junqueira, A. B., Kennard, D., Letcher, S. G., Licona, J.-C., Lohbeck, M., Marín-Spiotta, E., Martínez-Ramos, M., Massoca, P., Meave, J. A., Mesquita, R., Mora, F., Muñoz, R., Muscarella, R., Nunes, Y. R. F., Ochoa-Gaona, S., de Oliveira, A. A., Orihuela-Belmonte, E., Peña-Claros, M., Pérez-García, E. A., Piotto, D., Powers, J. S., Rodríguez-Velázquez, J., Romero-Pérez, I. E., Ruíz, J., Saldarriaga, J. G., Sanchez-Azofeifa, A., Schwartz, N. B., Steininger, M. K., Swenson, N. G., Toledo, M., Uriarte, M., van Breugel, M., van der Wal, H., Veloso, M. D. M., Vester, H. F. M., Vicentini, A., Vieira, I. C. G., Bentos, T. V., Williamson, G. B., and Rozendaal, D. M. A.: Biomass resilience of Neotropical secondary forests, Nature, 530, 211–214,, 2016. a

Poulter, B., Aragão, L., Andela, N., Bellassen, V., Ciais, P., Kato, T., Lin, X., Nachin, B., Luyssaert, S., Pederson, N., Peylin, P., Piao, S., Saatchi, S., Schepaschenko, D., Schelhaas, M., and Shivdenko, A.: The global forest age dataset (GFADv1.0), link to NetCDF file,, 2018. a, b, c, d, e, f, g, h, i, j

Puettmann, K. J., Wilson, S. M., Baker, S. C., Donoso, P. J., Drössler, L., Amente, G., Harvey, B. D., Knoke, T., Lu, Y., Nocentini, S., Putz, F. E., Yoshida, T., and Bauhus, J.: Silvicultural alternatives to conventional even-aged forest management – what limits global adoption?, Forest Ecosystems, 2, 8,, 2015. a, b

Reick, C. H., Raddatz, T., Brovkin, V., and Gayler, V.: Representation of natural and anthropogenic land cover change in MPI-ESM, Journal of Adv. Model. Earth Sy., 5, 459–482,, 2013. a, b, c

Sato, H., Itoh, A., and Kohyama, T.: SEIB–DGVM: A new Dynamic Global Vegetation Model using a spatially explicit individual-based approach, Ecol. Model., 200, 279–307,, 2007. a, b

Shevliakova, E., Pacala, S. W., Malyshev, S., Hurtt, G. C., Milly, P. C. D., Caspersen, J. P., Sentman, L. T., Fisk, J. P., Wirth, C., and Crevoisier, C.: Carbon cycling under 300 years of land use change: Importance of the secondary vegetation sink, Global Biogeochem. Cy., 23, GB2022,, 2009. a, b, c, d, e, f, g

Smith, B., Prentice, I. C., and Sykes, M. T.: Representation of vegetation dynamics in the modelling of terrestrial ecosystems: comparing two contrasting approaches within European climate space, Global Ecol. Biogeogr., 10, 621–637,, 2001. a, b

Soja, A. J., Shugart, H. H., Sukhinin, A., Conard, S., and Stackhouse, P. W.: Satellite-Derived Mean Fire Return Intervals As Indicators Of Change In Siberia (1995–2002), Mitig. Adapt. Strat. Gl., 11, 75–96,, 2006. a

Stocker, B. D., Feissli, F., Strassmann, K. M., Spahni, R., and Joos, F.: Past and future carbon fluxes from land use change, shifting cultivation and wood harvest, Tellus B, 66, 23188,, 2014a. a, b

Stocker, B. D., Spahni, R., and Joos, F.: DYPTOP: a cost-efficient TOPMODEL implementation to simulate sub-grid spatio-temporal dynamics of global wetlands and peatlands, Geosci. Model Dev., 7, 3089–3110,, 2014b. a, b, c

Sun, G., Noormets, A., Gavazzi, M., McNulty, S., Chen, J., Domec, J.-C., King, J., Amatya, D., and Skaggs, R.: Energy and water balance of two contrasting loblolly pine plantations on the lower coastal plain of North Carolina, USA, Forest Ecolo. Manage., 259, 1299–1310,, 2010. a

Tramontana, G., Jung, M., Schwalm, C. R., Ichii, K., Camps-Valls, G., Ráduly, B., Reichstein, M., Arain, M. A., Cescatti, A., Kiely, G., Merbold, L., Serrano-Ortiz, P., Sickert, S., Wolf, S., and Papale, D.: Predicting carbon dioxide and energy fluxes across global FLUXNET sites with regression algorithms, Biogeosciences, 13, 4291–4313,, 2016. a, b

van Mantgem, P. J., Stephenson, N. L., Byrne, J. C., Daniels, L. D., Franklin, J. F., Fulé, P. Z., Harmon, M. E., Larson, A. J., Smith, J. M., Taylor, A. H., and Veblen, T. T.: Widespread Increase of Tree Mortality Rates in the Western United States, Science, 323, 521–524,, 2009. a

Yang, X., Richardson, T. K., and Jain, A. K.: Contributions of secondary forest and nitrogen dynamics to terrestrial carbon uptake, Biogeosciences, 7, 3041–3050,, 2010. a, b, c

Yao, Y., Piao, S., and Wang, T.: Future biomass carbon sequestration capacity of Chinese forests, Sci. Bull., 63, 1108–1117,, 2018. a

Yue, C., Ciais, P., and Li, W.: Smaller global and regional carbon emissions from gross land use change when considering sub-grid secondary land cohorts in a global dynamic vegetation model, Biogeosciences, 15, 1185–1201,, 2018a. a, b

Yue, C., Ciais, P., Luyssaert, S., Li, W., McGrath, M. J., Chang, J., and Peng, S.: Representing anthropogenic gross land use change, wood harvest, and forest age dynamics in a global vegetation model ORCHIDEE-MICT v8.4.2, Geosci. Model Dev., 11, 409–428,, 2018b. a, b, c, d

Zaehle, S., Sitch, S., Prentice, I. C., Liski, J., Cramer, W., Erhard, M., Hickler, T., and Smith, B.: The importance of age-related decline in forest NPP for modeling regional carbon balances, Ecol. Appl., 16, 1555–1574, 2006. a

Short summary
Models need to account for forest age structures when investigating land use influences on land–atmosphere feedbacks. We present a consolidated scheme to introduce forest age classes, combining age-dependent simulations of important processes with the possibility to trace forest age, and describe its implementation in JSBACH4, the land surface model of the ICON Earth system model. We evaluate simulations with and without age classes demonstrating the benefit of forest age classes in JSBACH4.