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

. Natural and anthropogenic disturbances, in particular forest management, affect forest age-structures all around the globe. Forest age-structures in turn inﬂuence 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 5 introduce forest age-classes in hierarchical tile-based DGVMs combining beneﬁts of recently applied approaches. The ﬁrst being a computationally efﬁcient age-dependent simulation of all relevant processes, such as photosynthesis and respiration, using a restricted number of age-classes. 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- 10 classes. We describe how we implemented this scheme in JSBACH4, the LSM of the ICON-ESM. Subsequently, we compare simulation output against 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 beneﬁt of the introduction of age-classes, with the optimal 15 number of age-classes being a compromise between computation costs and error reduction.


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 M km 2 of forest area, i.e. about 15% of global ice-free land, are under some kind of management , with 65% being under regular harvest schemes and another 7% being intensive plantations (Erb et al.,20 2017). Often, management practices make use of rotation cycles, as common in shifting-cultivation (Boserup, 1966) or evenaged forest management strategies that historically were common in temperate forests and are still the dominant management ferent 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 5 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 has been 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/disturbances. Comparably, ORCHIDEE-MICT (Yue et al., 2018b) introduced a fixed number of six tiles per woody PFT, with tile merging upon exceeding pre-defined woody biomass 10 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 of 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 15 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 (FFM). FFM takes the tile wood increment calculated 20 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/size classes of trees competing for soil resources and light. For each forest tile, POP gets the stem 25 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 30 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.

JSBACH4
The DGVM JSBACH4 is used as 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.,5 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, movement of area fractions from one 10 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., in prep.): 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. (in prep.) introduced a dependency of the maximum LAI on the available leaf biomass. Such a dependency is a 15 prerequisite for simulating forest re-growth and thus for the introduction of age-classes.

JSBACH4-FF
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 to introduce different processes and associated state variables on 5 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 pre-set in the configuration file for all forest PFTs (Fig. 1). In addition, a to be pre-defined upper age bound per age-class AC K (maxA K ) as well as a total maximum age (maxAge) were introduced. maxAge determines the oldest age up to which the age of an area is tracked, i.e. the length of the f ractP erAge 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 length 150 to track the exact age of the entire forest area up to 150 year old forest. Each AC K covers a certain interval of years [maxA K−1 , maxA K ) (Fig. 3), with the youngest AC (AC 1 ) always covering the range of year 0 to 1, and the oldest AC N covering all forest older than maxA N −1 , i.e. [maxA N −1 , IN F ), with maxA N −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 f ractP erAge vector.
Ageing The newly implemented process of forest "ageing" happens annually: upon ageing each tracked forest fraction gets one 5 year older. Yet, a shift from one age-class to the next age-class only happens for the area of the oldest age (maxA K−1 -1) of an age-class AC K-1 , i.e. only the forest area which upon getting one year older exceeds the upper age bound maxA K−1 of the AC K-1 needs to be shifted into AC K . Thanks to the tracking of the age in the f ractP erAge vector, the exact area fraction with age maxA K−1 -1 is known.
Harvest In the current version, we implemented harvest as a clear-cut of a certain fraction of an AC, which can happen 10 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 f ractP erAge vector, age-based harvest rules can be specified.
Disturbances 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 15 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 (AC 1 ).
Since JSBACH4 so far does not provide the infrastructure for the horizontal exchange of properties among tiles (Section 2.1), we had to implement such an infrastructure in JSBACH-FF, redistributing area fractions and associated state variables such as 20 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 Eq. 1 and Eq. 2. Here, f a is the area fraction moved from one AC to another, e.g. upon ageing; V S is the value of the state variable of the source AC S and delta_V S and delta_V T are the scheduled changes of the state variable on the source AC S and the target AC T , respectively.
In the second step, the redistributions are applied on the ACs (apply_state_changes in Fig. 2) according to Eq. 3. Here, V AC represents the affected state variable of the age-class AC, delta_V AC the scheduled change, f a is the incoming and f c 30 the current area of the age-class AC.

Simulation set-up and measures of model performance
The main purpose of JSBACH are 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 5 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 (Section 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 against GPP, LAI and AGB data based on observations using a normalised root mean squared error (NRMSE; Section 2.3.5). In addition, we created Taylor diagrams for 10 each variable, season and region (see Figures S5. 5-S5.11 in the supplementary).

Observation-based data
We used 2010 MODIS LAI (Myneni et al., 2002) and GPP data obtained from machine learning methods trained on fluxtower measurements (Tramontana et al., 2016). These LAI and GPP datasets had already been mapped to JSBACH's forests PFTs in a previous study (see supplementary of 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 (Sec-5 tion 2.3.2). The AGB per forest area (Avitabile et al., 2016) was downloaded from the GEOCARBON data portal (www.bgcjena.mpg.de/geodb/projects/Data.phd) 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 supplementary show maps of the pre-processed observation-based data.
2.3.2 General simulation set-up 10 We conducted simulations with JSBACH4 (4.20p7) feature/forest in standalone mode hosted within the MPI-ESM environment (see supplementary material 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 15 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 x 96 latitudes; 1.9 • × 1.9 • ), a climate forcing based on GSWP3 (Kim et al., 2012) and CO 2 from the collection of greenhouse gas concentrations for CMIP6 (Meinshausen et al., 2017). To obtain forest age distributions 20 comparable to those observed for 2010 (Poulter et al., 2018), harvest was conducted following prescribed maps (Section 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 25 distributions of age-classes (Section 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. 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 30 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.

Simulated number of age-classes and age distribution schemes
The two applied age distribution schemes are defined as follows: EAS The equal age-spacing (EAS) distribution scheme spreads the age-classes evenly over the age space.  Poulter et al. (2018). 5 IAS 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 (uLim K ) is defined following Eq. 5. 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. Table 1. Conducted 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. 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, independent 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 0.5 • resolution of the global forest age distribution of four forest PFTs: needleleaf evergreen (NE) and needleleaf deciduous (DE), as well as broadleaf evergreen (BE) and broadleaf deciduous (BD).

5
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 pre-processing 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: 10 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 that PFT k is evergreen) or both are deciduous PFTs (in case that PFT k is deciduous).
From these age maps, we derived harvest maps for each simulation year such that the simulated age distribution in simula- 15 tions 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 ten years covered by each age-class (see supplementary material 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, a clear-cut 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 5 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.

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 (Section 2.3.1), and 'SIM', the results of each simulation (see Table 1). The RMSE was calculated  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).
To more easily assess changes in performance when increasing the number of age-classes the different NRMSE Max-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.

Computational costs
In addition to determining the NRMSE for different variables, we also determined the computation costs of the different setups. 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 (Section 3.1). 10 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 (Section 3.2). Finally, we discuss assets and drawbacks of alternative schemes introducing age-classes in tile-based DGVMs (Section 3.3).

15
In this section we use the NRMSE Max-Min for different regions/seasons as 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 Section 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 20 the observation-based data than the simulations with age-classes (Fig. 5c). For most comparisons, the NRMSE Max-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.
Averaging the NRMSE Max-Min , giving each region and each season the same weight, results in an NRMSE Max-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 25 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 NRMSE Max-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 NRMSE Max-Min for the boreal AGB is due to an increased underestimation when accounting for more young forest, as is also discussed below  -TMP TROP  BOR NH-TMP TROP  BOR NH-TMP TROP   PFT  EAS03  IAS03  EAS06  IAS06  EAS08  IAS08  EAS11  IAS11  EAS13  IAS13 EAS16 IAS16 Figure 5. Evaluation of the conducted simulations (Table 1) (Table 2); and for LAI and GPP also for each of the four seasons.
This decrease in NRMSE Max-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, albeit the benefit of increasing the number of age-classes is slowly saturating towards larger numbers of age-classes (Fig. 6). 5 We performed the averaging of the NRMSE Max-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 lead to an increasing weight of the tropical region, and thus to an earlier saturation of the NRMSE Max-Min with increasing age-classes. 10 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 nearlinear 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 15 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, is not dominating 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. Nabel, 2015;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 manuscript, one set-up has been selected as an illustrative example: IAS11 (see Table 1), i.e. the simulation with 10+1 ageclasses and the age distribution scheme with increasing age space. This set-up is a compromise between the error reduction for 5 GPP and LAI comparisons on 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.

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  Fig. 7 we compare results of a JSBACH4-FF simulation with age-classes, a simulation only using PFTs, as representative for 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 15 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-cuts 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 one hand and young managed forest on the other hand in India (Fig. 7c); a grid-point with intensive harvest/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 pre-dominantly old-grown forests in central Africa (Fig. 7g). In general, the simulation with age-classes results in smaller GPP, LAI and AGB values ( Fig. 7 and Figs. S4.2-S4.4, column 3), which is 5 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).  15 replaced by a slight underestimation (Fig. 7 panel d and f). Globally, reduced overestimations get 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. S4.d)  S4.4). These are particularly areas with pre-dominantly old-grown forests, i.e. without a distinct age-structure, such as central Africa, 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 NRMSE Max-Min in areas where the PFT simulation was biased high and an increase in the NRMSE Max-Min in areas where the PFT simulation was biased low. Thus, if such a forest 25 age-structure would be implemented in a DGVM being predominately biased low, the difference to the observation-based data could increase. In this context, also caveats regarding the observation-based data themselves 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, 30 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   5 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 10 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 a clear-cut.

Limitations and alternative schemes
In the previous sections we compared a JSBACH4-FF simulation when only using PFTs with simulations including forest ageclasses and discussed associated trade-offs and benefits. In this section we discuss limitations and advantages of the applied 15 and of alternative schemes.
Since JSBACH is a tile-based DGVM, the introduction of an individual/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, 20 is one major aspect and purpose of JSBACH, which historically has been part of the MPI-ESM (Mauritsen et al., 2019) and now also is 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 one hand, this approach shares the advantage with individual/cohort-based DGVMs that it provides a 25 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 (as e.g. depicted for spring LAI in Fig. 7). This limitation can only be avoided by increasing 30 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  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. Center 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 timespan 2001-2010 (this is not the case for the age-classes 9 and 10 in panel 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. (Kuusela, 1994;Pan et al., 2011;Puettmann et al., 2015;Kuuluvainen and Gauthier, 2018). 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 to shift the actually affected parts of an age-class, and not entire age-classes/"cohorts" as 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. 5 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 tilehierarchy 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. 10 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 ageclasses/"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 15 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 20 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 indicated in Fig. 7b or described in Zaehle et al. (2006) andBellassen 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 25 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/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 30 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.
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.

5
JSBACH4 itself is still highly under development regarding infrastructure and processes integrated from JSBACH3. In the version used for this paper (4.20p7), particularly the representation of natural and anthropogenic land cover change 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 cuts, following the implementation of other disturbances in JSBACH3 (see Brovkin et al., 2009). For a representation of different 10 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 15 provides a valuable tool to study forest management effects, particularly due to its integration with the ICON-ESM.