Accounting for forest management in the estimation of for est carbon balance using the dynamic vegetation model LPJ-GUESS (v4.0, r9333): Implementation and evaluation of simulations for Europe

Global forests are the main component of the land carbon sink, which acts as a partial buffer to CO 2 emissions into the atmosphere. Dynamic vegetation models offer an appr oach to making projections of the development of fo rest carbon sink capacity 10 in a future climate. Forest management capabilities n dynamic vegetation models are important to incl ude the effects of age and species structure and wood harvest on carbon stocks and carbon storage potential. This article describ es the introduction of a forest management module in the dynamic vegetation model L PJ-GUESS. Different ageand species-structure setu p strategies and harvest alternatives are introduced. The model is u sed to represent current European forests and an au tom ted harvest strategy is applied. Modelled carbon stocks and fluxes are eval uated against observed data at the continent and co untry levels. Including wood 15 harvest in simulations increases the total European carbon sink by 32 % in 1991-2015 and improves the fi to the reported European carbon sink, growing stock and net annual increment (NAI). Growing stock (156 m3 ha-1) and NAI (5.4 m3 ha-1 y-1) densities in 2010 are close to reported values, while the carbon sink density in 2000-2007 (0.085 kgC m -2 y-1) is 63 % of reported values. The fit of modelled values and observations for individ ual European countries vary, but NAI is generally c loser to observations when including wood harvest in simulations. 20


Introduction
Forests globally provide ecosystem services including provision of timber, fuel and water, regulation of local climate and hydrology, carbon sequestration, support of biodiversity and recreation (Bonan 2008;Mori et al., 2017). The effects of climate change on forest productivity and biodiversity may be predicted to be negative due to increased evapotranspiration and reduced rainfall in many forested areas, an increase in extreme events like drought, wild-fires, storms and insect attacks and local or regional 25 extinctions of plant and animal species (Easterling et al., 2000;Seidl et al., 2011;Anderegg et al., 2015;Urban 2015). On the other hand, productivity may increase due to the fertilising effect of increased nitrogen deposition and higher atmospheric CO2 levels (Zaehle and Dalmonech, 2011;Luyssaert et al., 2008) as well as shifts in tree species composition and longer growing seasons at high latitudes caused by higher temperatures (Sitch et al., 2015;Morin et al., 2017).
Forests make up the largest portion of the current land carbon sink, and are estimated to have absorbed 20-50 % of CO2 emitted by 30 fossil fuel combustion and industry during the first decade of this century (Pan et al., 2011;Le Quéré et al., 2018;Pugh et al., 2019). The suggested basis for this carbon uptake is the recent history of the drivers increasing productivity mentioned above, especially increased CO2, and the recovery of carbon pools in regrowth forests (Pugh et al., 2019). The size of the forest carbon sink has been estimated by either book-keeping methods (Pan et al., 2011;Houghton et al., 2012) or by global vegetation models Shevliakova et al., 2009;Pugh et al., 2019) but this sink is associated with relatively large uncertainties, 35 resulting in differing estimates using different approaches and models. Key uncertainties include the magnitude of CO2 fertilisation https://doi.org/10.5194/gmd-2020-440 Preprint. Discussion started: 27 January 2021 c Author(s) 2021. CC BY 4.0 License.
-which may be limited by soil availability of nutrients such as N and P (Zaehle and Dalmonech, 2011;Jiang et al., 2020) -and the extent of shifting cultivation in the tropics (Heinimann et al., 2017). While the net atmosphere-to-land flux (FL) is relatively well constrained by measurements, large uncertainties in the net land-use and land-cover flux (FLULCC) make the size of the residual (land) sink (FRL) itself uncertain (FL= FRL-FLULCC) (Arneth et al., 2017).

40
Forests cover 33 % of the Europe's land area (Forest Europe, 2015) and store approximately 13 PgC in vegetation and 28 PgC in soils (Pan et al., 2011). The carbon sink of European forests in [2000][2001][2002][2003][2004][2005][2006][2007] has been estimated at 0.27 PgC y -1 or about 12 % of the global carbon sink of established forests (Pan et al., 2011). Europe has been identified as a region where regrowth forests dominate carbon sequestration (Pugh et al., 2019) and has a history of thousands of years of human impact on forest structure and species composition (Perlin, 2005). Forest management practices of the past few hundred years are relatively well documented (McGrath 45 et al., 2015). Depending on the region, different management strategies are applied (Cardellini et al., 2018). The relatively young tree age and the removal of wood in managed forests influence carbon stocks and fluxes e.g. by increasing productivity and reducing self-thinning and age-related mortality and litter production compared to pristine forests (Zaehle et al., 2006). In addition to the effects on atmospheric CO2 , forest management influences local climate by changing albedo, evapotranspiration and surface roughness (Luyssaert et al., 2014).

50
Dynamic global vegetation models (DGVMs) provide a potential framework for predicting the combined effects of climate and forest management scenarios on forest ecosystem structure and carbon balance. Based on such information, the potential of forest landscapes to contribute to climate change mitigation by maintaining or enhancing carbon sinks, and to climate adaptation through sustained production of forest products and other ecosystem services in the face of climate change, can be assessed. Applications of DGVMs to represent climate responses of potential natural vegetation (PNV) have been shown in the past, for example as a 55 basis for nature conservation planning (Hickler et al., 2012). Human management of land, including cropland, pasture and managed forest, has been introduced in a number of DGVMs (Bondeau et al., 2007;Bellassen et al., 2010;Lindeskog et al., 2013;Arneth et al., 2017). Key elements required to represent managed forests in a DGVM framework include the ability to initialise a simulation with historical land use, to represent age/size structure of forests stands and their change over time, to account for tree species composition and to apply silvicultural treatments that modify stand composition and structure, like planting, thinning and 60 harvesting.
LPJ-GUESS (Smith et al., 2001;Smith et al., 2014) is a second-generation DGVM incorporating a detailed representation of forest ecosystem composition and stand dynamics, suitable for the implementation of a forest management scheme. It captures the distribution of European PNV at species level and can make projections of vegetation shifts under future climate scenarios (Hickler et al., 2012). The model has been shown to represent vegetation growth and succession successfully (Smith et al., 2014). It has 65 been used to estimate forest vulnerability to climate change (Seiler et al., 2015) and carbon mitigation potential of regrowth forests and forests under alternative management scenarios (Pugh et al., 2019;Krause et al., 2020). Earlier versions of LPJ-GUESS have been modified to enable analysis of clearcut forest management and the effects of wind damage and insect outbreaks Jönsson et al., 2012). In this study, we describe the implementation of forest management capabilities in LPJ-GUESS v.4.0, which considers, in addition to detailed carbon-and water-cycle processes, nitrogen-cycling and nitrogen-limitation (Smith 70 et al., 2014). Model alternatives for forest stand initialisation (land-use history and species-and age-distribution) and silvicultural management (detailed and automated harvest strategies) are presented in detail. Using an automated thinning and clearcut approach for European forests, we compare modelled carbon stocks and fluxes with observational data and explore the performance under a changing climate.

75
2.1 General description of LPJ-GUESS and overview over simulated processes LPJ-GUESS (Smith et al., 2001;Smith et al., 2014) simulates the dynamics of terrestrial vegetation and soils across a regional or global grid, forced by meteorological inputs and soil physical properties. In the original implementation, each grid cell encompasses a landscape of natural, climatically determined vegetation (PNV). A number (5-100) of replicate patches, nominally 0.1 ha in size, represent disturbance-related variation in stand age across the wider landscape of a grid cell. In each patch, age cohorts of tree and 80 shrub plant functional types (PFTs) or species and grass PFTs compete for light, water, nitrogen and space (Fig. 1). In its original version, the model only simulates PFTs that capture the major vegetation zones globally. The parameter set of these PFTs has been extended to simulate the most important tree species in northeastern USA (Hickler et al., 2004) and Europe (Koca et al., 2006;Hickler et al., 2012) as distinct PFTs. The new functionality defined in this paper can operate equally on individual tree species or more generalised PFTs. Hereinafter 'species' is thus used synonymously with 'PFT'. The canopy is represented as a multi-layered 85 structure. Leaves, fine roots and stem heartwood and sapwood are represented as dynamic pools for each age cohort of each PFT.
Branches and course roots are not explicitly discriminated but are implicit in the wood biomass pool. The patches are subject to stochastic vegetation-destroying disturbance events (representing e.g. wind-storms or landslides) with a prescribed return interval (e.g. 100-400 years). Disturbance results in the loss of vegetation in a patch, after which a secondary succession of grass and tree PFTs follows. Establishment is affected by forest floor light conditions and is subject to PFT-specific environmental envelopes 90 defined by bioclimatic limits. A slightly different set of bioclimatic limits govern survival (Table A1). Growth-efficiency-, selfthinning-, background-(age-related) and fire mortality are applied to individual cohorts. Establishment and mortality have a stochastic component. Soil carbon and nitrogen cycling are based on the CENTURY model (Parton et al., 1993) and soil hydrology on a two-layered "leaking bucket" model. A soil mineral nitrogen pool is provided by atmospheric deposition, biological nitrogen fixation and gross nitrogen mineralisation of soil organic matter. Plant nitrogen uptake is driven by the demand from photosynthesis 95 and biomass growth, and is limited by the supply from the soil mineral nitrogen pool. The nitrogen dynamics in the model are described in detail by Smith et al. (2014). Different land-use/land-cover types in addition to PNV are represented in the model by stand types with different management, e.g. cropland, pasture and managed forest (Lindeskog et al., 2013, Fig. 1). Transitions between different stand types may occur at 105 any point in time, according to historic land-use data, to recreate land-use history. When a potentially forested stand type area expands, new stands are created, keeping the soil history from the previous stand type intact and allowing vegetation succession to proceed from bare ground (in most cases, but cf. 2.2). In modelled wood harvest events 66 % of wood biomass and 30 % of leaf biomass is removed from the stand and the rest remains as litter. Removed leaf biomass and part of wood biomass (by default 67 %) is oxidised the same year. The remaining wood biomass is put into a product pool with a turnover rate of 4 % per year.

Forest structure initialisation routines
Forest stand age-and species distributions can be achieved in the model by several alternatives, utilising the structure of a previous PNV stand or defining a new age-and species structure at various levels of detail.

Stand creation
A managed forest stand may be created in the model by two different options (Fig. 2, B1). By cloning the stand of origin, the 115 complete state with all patches intact is inherited by the secondary stand. If the origin is previous woodland (PNV or secondary forest), a cutting scheme may start with the existing tree structure, optionally cutting unwanted species. In the other alternative, tree growth starts from bare ground after an initial clearcut or when expanding on former cropland or pasture. In this case (with an even-age stand and if disturbance and fire are turned off), the secondary stand can in many cases be modelled by a smaller number of replicate patches since there is usually no random variation in the timing of management events. 120

Secondary forest age structure
Managed forest stands with an uneven age structure can be represented in the model by different options. An age structure may be created in individual patches by thinning (enabling regeneration by increased light at the forest floor) at defined intervals during an initialisation period, allowing for competition between both cohorts and species (Fig. 3a). To achieve an age structure among patches within a stand, the age structure of PNV, achieved during the spinup by patch-distroying disturbance events, may be 125 conserved after the conversion to managed forest if the cloning functionality is used, copying the PNV stand with the semirandomised age structure intact (Fig. B1). Alternatively, multiple patches in a secondary stand may be clearcut successively at regular intervals during an initialisation period (Fig. 3b). In the final approach, a prescribed age structure, either representing a specific moment in time, or a historical development, may be created among stands representing a stand type, using land-cover change input data (Fig. 3c).

Secondary forest species composition
Species mixtures may be defined either at the management type level (Fig. 1), using predefined planting densities for individual species and/or later cuttings to achieve prescribed relative abundances of the different species within a patch (Fig. 4a, see below), or at the landscape level, using land-cover input data to achieve predefined mixtures of monocultures (Fig. 4b), or a combination of both of these options.

155
Depending on the level of detail in historic forest management input data or, in simulations of future scenarios, whether the management should be able to adapt to a changing climate or other factors, various model alternatives are available.

Simplified clearcut forestry
A simplified method to represent forestry using global wood harvest input data (e.g. harvested area) is achieved by creating secondary forest stands after clearcutting either a PNV stand or other secondary forest stands, representing cutting of primary or 160 secondary forest, respectively. In cutting events, looping through the stands, these are cut according to stand age rules (cut oldest or youngest stands first, avoid cutting stands younger than 15 years old), allowing the allocation of wood harvest to primary forest and mature or young secondary forest. This method was used with LUH2 input data by Pugh et al., (2019).

Detailed forestry
A number of forest management options can be selected at the stand type or management type level in the LPJ-GUESS instruction text file required to run a simulation (Table 1).

Species selection
A forest stand may contain a full selection of tree species (as in PNV) or a selection of species defined in the management type.

175
After a clearcut event, or after creating a new forest stand from bare ground/grassland, selected species may be planted at defined sapling densities with or without the additional need to fall within the envelope of the bioclimatic limits that govern PFT establishment in PNV mode (Table A1). Alternatively, the standard establishment method can be used. After the initial planting/establishment, re-establishment can be optionally enabled or disabled for selected and unselected species. If several tree species are selected, it is possible to define a target relative abundance for each species (relative biomass) and apply selective 180 cutting (Fig. 4a). Start and end calendar years for this treatment may also be defined.

Clearcut forestry
A fixed rotation period is defined, at the end of which a clearcut takes place (Fig. 5a). Alternatively, clearcut may be triggered by attainment of a prescribed stand density limit (Fig. 5b). The timing of a number of thinning events (default 5) may be defined as https://doi.org/10.5194/gmd-2020-440 Preprint. Discussion started: 27 January 2021 c Author(s) 2021. CC BY 4.0 License.
fractions of the rotation period in the case of a fixed rotation period. The harvest amount (strength) for such thinning events is 185 defined as a fraction of current biomass, with the option of different settings for selected and unselected species. At each thinning event, trees may be cut using alternative strategies. Size/age criteria: (1) old or big trees first; (2) young or small trees first; (3) a specified harvest amount pertaining to trees above a specified diameter limit only; (4) all sizes/ages cut equally, may be combined with species criteria: (1) selected species first; (2) unselected species first; (3) separatly defined harvest amounts for selected and unselected species; (4) shrubs and shade-intolerant species first; (5) all species cut equally (Fig. 5a). young/old first small/big first selected/unselected species first shrubs and shade-intolerant first diameter limit (cut only trees above a diameter limit) 2. Automated: Self-thinning rule-based method Clearcut 1. Fixed rotation time 2. Automated: Tree density limit Cut to species distribution target Relative biomass of selected species monitored at 5-year intervals; if value deviates more than 10 %, cut dominant species to reach target. Start and end of cutting period may be defined. N fertilisation kgN ha -1 y -1 , evenly distributed the whole year Irrigation Bypass water stress in photosynthesis Fire/disturbance suppression Switch off fire and disturbance

Management change
Change management type a specific calender year (optionally wait for clearcut) 1 1 All management options in this table except re-establishment can be defined in separate management types, which may be selected in a stand type rotation scheme at pre-defined calender years.

Continuous cutting
When modelling continuous cutting, it is possible to define the same harvest parameters and cutting priority settings as described above for the clearcut case, for two different periods: the first for a specified "regeneration" time following a clearcut, and the second for a "continuous" phase, in which the cutting cycle is repeated indefinitely (Fig. 5c). During the continuous phase, the minimum diameter limit in tree size selection option (3) above can be adapted to low productivity by automatically lowering the  https://doi.org/10.5194/gmd-2020-440 Preprint. Discussion started: 27 January 2021 c Author(s) 2021. CC BY 4.0 License.

Automated wood harvest
As an alternative to specifying the thinning rules in clearcut forestry in detail, a thinning scheme based on Reineke's self-thinning rule may be chosen (Fig. 5b)

, 5 225
where 5 is the initial tree density and *+, is the density limit for clearcut (se below). As an alternative to imposing a specified rotation length in clearcut forestry, clearcut may be triggered by stand density when it is below *+, as in Bellassen et al. (2010). ( *+, and *+, were selected and further adjusted to give rotation times around 100 years in the early 2000s in LPJ-GUESS simulations (Table A3).

Nitrogen fertilisation and irrigation
A specified amount of nitrogen may be applied to the soil every year (Fig. B2). With irrigation enabled, the amount of water required to avoid water stress is calculated and applied every year.

Management change
To capture management changes, a new silvicultural treatment of a stand type can be prescribed any specified calender year, 235 changing from one specified management type to another with the next harvest event as an optional trigger (Fig. 6).

Demonstration simulation protocol
To demonstrate the implemented forest management functionality and its effects on simulated stand structure, composition and 245 productivity, we performed demonstration simulations for a representative location (grid cell) in Europe, and across Europe as a whole. PNV stands were modelled using 25 replicate patches and a disturbance return time of 400 years. Managed forest stands contained only one patch except where explicitly stated, disturbance and fire were turned off and mortality was deterministic. In managed forest stands created after clearing the previous vegetation, this setup saves computational time and produces almost identical results compared to using multiple patches and adding the stochastic component to establishment and mortality. The 250 European species as described by Hickler et al. (2012) were used with updated parameters (Table A1- In future climate scenario simulations, monthly temperature, radiation and precipitation data for 1850-2100 were adopted from the In future forest projections, either the historic environmental drivers were recycled after 2015 or future climate, CO2 and nitrogen projections were used to demonstrate model behaviour under a time-span of several forest rotations.

285
The EFI Tree species map describes the spatial distribution (fraction of land area) of 20 tree species groups at 1 x 1 km resolution (Brus et al., 2011). The map is based on ICP-Forest Level-I plot data combined with National Forest Inventory (NFI) data of 18 countries. In areas with NFI data, spatial interpolation of the plot data was used, whereas in areas without NFI data, statistical relationships between tree species and covariates (soil, biogegraphy and bioindicators) were used (Brus et al., 2011). The EFI Tree species map was aggregated to 0.5 o x 0.5 o resolution in this study and was used to further refine the species distribution derived 290 from GFAD.
The structure of European forests in 2010 was reconstructed by using a combination of the the GFAD age database and the EFI Tree species map. For each gridcell, the most common species or species group within the GFAD NE and BD forest types were obtained from the EFI Tree species map and these species groups were then mapped to LPJ-GUESS tree species/species groups (Table C1, Fig. C2). In the multi-species LPJ-GUESS groups, species compete with each other for resources (cf. above, 2.1). BE 295 was mapped to Quercus ilex and ND to Larix decidua, the only available PFTs in the model to represent these two functional tree classes.

Modelling current and future European managed forests
Secondary forest stands were created in the model from 1871 to 2010 to obtain the GFAD age (1-140 years) distribution in 2010 and species selections were planted (without climate restrictions for NE and ND stands to bypass establishment temperature limits 300 used in PNV). The oldest forest class in GFAD (>140 years) was modelled as PNV and was not subject to any management. In secondary stands, automated thinning and clearcut (cf. 2.3.2) were implemented using the parameters in Table A3  is defined as the difference in the sum of vegetation and soil carbon pools between two consecutive years plus the removed 310 harvested carbon, not taking into account any reductions in wood products and residues following removal from the site. Similarly, NAI is defined as the difference in growing stock volume between two consecutive years plus the removed harvested volume.
Harvested carbon is not included in the total carbon pool and includes both wood products and removed wood residues. The forested area in 2010 as defined by GFAD and Forest Europe (2015)

Implications of secondary forest initialisation and land-use history
Secondary forest stand initialisation and land-use history have long-term effects on the development of tree species distribution, 320 productivity and carbon fluxes in the model (Fig. 2). When the age distribution and species composition from spinup is retained in each patch (i.e. cloning PNV), both the warming climate in the 20th century and the prevention of fires and other disturbances result in an increase in tree biomass and a tree species shift from a Q.robur-P.sylvestris-dominated forest landscape to a forest increasingly dominated by the shade-tolerant species P.abies and F.sylvatica in an example forest simulated at a Southern Swedish site (Fig. 2a). Older patches contribute to an early-onset stagnation of the carbon sink. A forest stand created after a clearcut of 325 PNV displays a mixed broad-leaf forest with a late establishment and dominance by P.abies (Fig. 2b&c). Leaving harvested biomass on site results in an extended litter-induced carbon source (Fig. 2b). When the previous land-use history is grassland, the initial dominance by shade-intolerant species is more pronounced and the slow accumulation of the litter pool results in a stronger and more protracted carbon sink (Fig. 2d&e). Soil carbon and nitrogen depletion due to intensive harvest of the previous grassland influences productivity, succession of tree species and carbon sink capacity of the secondary forest. Initial tree growth is delayed 330 by several decades, the dominant shade-intolerant species is P.sylvestris rather than B.pubescens and Q.robur competes more successfully than under normal soil nitrogen (Fig. 2e). Also, the long-term carbon sink is larger than in any other option. The big differences in tree species succession and the timing and magnitude of the carbon sink using the different stand creation options illustrate the importance of land-use history for modelling secondary forest stands.

335
The choice between the different age-and species structure setup options depends on whether competition between species and cohorts within patches is required or not (Fig. 3-4). Also, the desired level of detail of the age structure might decide whether to use a simplified setup or a detailed structure with many separate stands, increasing computation time. Setups using separate stands https://doi.org/10.5194/gmd-2020-440 Preprint. Discussion started: 27 January 2021 c Author(s) 2021. CC BY 4.0 License.
for each species-age combination offer the possibility of reflecting regional distributions based on inventory data, but might not represent competition correctly e.g. in mixed forests.

340
Although management changes during the course of a simulation may be prescribed, using detailed, but static, harvest methods would not reflect foresters' choice of gradual adaptation of harvest parameters under changing CO2/climate conditions in future scenarios. In these cases, the simplified dynamic harvest methods might be a better option (Fig. 5b).

European-wide simulations of managed forest
Dominant tree species in managed forests based on the EFI species map differ from PNV simulations in large parts of Europe. In 345 central and eastern Europe, broadleaved species are to a large degree replaced by needleleaved species in managed forests, especially by P.sylvestris, but since old-growth forest is modelled as PNV in this study, the dominance by needleleaves in this region seen in the original EFI data is moderated in the total forest landscape (Fig. C3, C4).        carbon sink seen in a simulation with thinning (0.12 kgC m -2 y -1 ) (Fig. 10) is correlated with a strong reduction of natural mortality (-80 % in 1991-2015) in regrowth forest stands, induced by thinning and, after 2010, rejuvenation of regrowth forest stands resulting from clearcuts (Fig. D6). The reduced natural mortality following thinning results in a lower soil respiration (Fig. D7).

Robustness of automated harvest methods under future climates
To demonstrate the automated harvest methods, in which thinning intensity and rotation times are adjusted to maintain standing 470 stock when stand productivity changes in response to forcing conditions, we used CO2/climate projections in extended simulations with an otherwise identical setup as in the European-wide historic simulations. A significant modelled increase in NAI is accompanied by shorter rotation periods (Fig. D8), while a stable vegetation pool in managed forest is maintained (Fig. D9). The mean thinning fraction of the total harvest over the rotation for NE and BD stands increased over the the 21st and 22nd centuries from 0.50 to 0.53 and 0.40 to 0.46, respectively, for both the RCP4.5 and RCP8.5 simulations (not shown).

4 Discussion
LPJ-GUESS representations of unmanaged forest have previously been compared favourably with observed forest vegetation succession, growth, stand structure, biomass and regrowth timescales (Smith et al., 2001;Smith et al., 2014, Pugh et al., 2019 and land use and land-cover change (LULCC) functionality has been included in the model since version 4.0 (Lindeskog et al., 2013).
In a recent global study that used the model to analyse the carbon stocks of old-growth and regrowth forests (modelled as primary 480 and secondary forest stands, respectively), without wood harvest in regrowth stands (Pugh et al., 2019), the total forest carbon sink https://doi.org/10.5194/gmd-2020-440 Preprint. Discussion started: 27 January 2021 c Author(s) 2021. CC BY 4.0 License.
was found to be about 50 % of values reported by Pan et al. (2011). The absence of wood harvest has been identified as an important factor for under-estimating carbon sinks in vegetation models (Zaehle et al., 2006;Ciais et al., 2008). In an effort to improve the ability to simulate carbon pools and fluxes on managed land, we here introduce new forest management options into LPJ-GUESS v4.0 and provide a comprehensive description of forest initialisation and wood harvest alternatives. The initialisation and harvest 485 alternatives used are constrained by the available forest inventory data and harvest information. Ideally, both age and species structure as well as land-use history and currrent wood harvest strategy should be taken into account, but this is not always possible for simulations with a large spatial extent because of limited data availability. To demonstrate a possible workaround, we used an automated thinning and clearcut alternative to represent European forests, initialised on the basis of inventory-based age-and species data, but without wood harvest-or LULCC data input.

490
The modelled mean vegetation carbon density in European forests in 2000-2010 is close to observations from several published sources (Pan et al., 2011;Liu et al., 2015;FOREST EUROPE, 2015). Including thinning in the simulation has a rather small impact on vegetation carbon (<5 %), but after clearcutting starts in regrowth stands after 2010, simulations with and without harvest in regrowth stands diverge strongly (Fig. D7). Also the modelled mean European growing stock is close to observations. Modelled carbon sink density (= -NEE) for European forests in simulations without thinning in the present study is about 46 % of the 2000-  value, respectively. Thinning reduces natural mortality due to relaxed competition between trees, and since a large part of harvested biomass is removed from forest stands, litter input to the soil, and the resulting heterotrophic respiration, is also reduced (Fig. D6-D7), increasing the carbon sink.
Details in the simplified European setup might explain the remainder of the "missing" carbon sink, relative to reported values. One cause is that old-growth forests are represented by unmanaged PNV (with a much lower carbon sink, cf.  (2019)), which is most likely inappropriate for Europe. Including wood harvest in old-growth forests would be expected to increase the carbon sink. Wildfires also contribute to a lower carbon sink in modelled PNV. A further likely cause of the discrepancy between the modelled and reported carbon balance is that secondary forests are created from PNV stands, without taking land-use history into account. Reforestation of cropland, which generally has a much lower soil carbon content than forests in Europe (Guo et al., 2002), has a higher carbon storage potential than clearing of existing forests. Also, soils of existing European 510 forests have probably been depleted of carbon historically because of higher harvest rates, fuel-wood collection and litter raking McGrath et al., 2015). Higher initial soil carbon pools will increase the release of CO2 in regrowth forests, especially under rising temperatures. Alternative methods to initialise secondary forests (fate of cleared wood, land-use history) have large implications for simulated carbon pools and fluxes as seen in the example Swedish site in this study, e.g. a mean carbon sink over 150 years spanning from 0.078 to 0.188 kgC m -2 y -1 (Fig. 2). This has also been shown at the global scale (Pugh et al.,  Our simulation results using LPJ-GUESS are consistent with results from the ORCHIDEE DVM, which was applied with the same automated thinning method (Bellassen et al., 2010). The ORCHIDEE simulation with automated thinning, compared to a simulation without thinning, gave a similar modest vegetation reduction (7 %), thinning fraction (0.55), reduced heterotrophic respiration (ca. 20 %) and carbon sink increase (67 %). The forest NPP reduction over time in ORCHIDEE simulations (ca. 10 %) is also seen in the average value for unharvested regrowth forests in European simulations with LPJ-GUESS, while in simulations 525 with clearcuts in regrowth forests, a balance between stands with different age is seen after clearcut starts in 2011 (Fig. D7b). The decline of NPP directly after thinnings in ORCHIDEE is not included in this version of LPJ-GUESS, but both models display a short-lived increase in heterotrophic respiration after thinnings (not shown). The recovery time after a clearcut (when the stand turns into a carbon sink) is 6 years in the example south Swedish site with a standard harvest removal, but 18 years if the harvested biomass is left on site (Fig. 2). This is similar to the ORCHIDEE results with a stand recovery time of 10-20 years after a clearcut.
The automated thinning/clearcut modelling strategy applied in the model in the present study is intended as an example for demonstrating the new forest management capabilities and an improvement on the age-structure setup of Pugh et al. (2019) and does not include all available possibilities in the model. In addition to the shortcomings in the setup already noted concerning land-535 use history, many central European forests are managed by continuous wood harvest and not by clearcutting and also consist of species mixes. Estimating the effect of such different wood harvest strategies and monoculture/mixed-species alternatives on carbon stocks and fluxes is now possible and will be done in a further study. The self-thinning and tree density-based harvest method is less successful in the northernmost and southernmost parts of Europe, where productivity is strongly limited by temperature and precipitation, respectively, and the self-thinning relationship between biomass and tree density is weaker. The

540
poor productivity of forests in the Mediterranean probably reflects a requirement for a revision of the parameterisation of tree species to better reflect Mediterranean managed forests or possibly the introduction of tree species that are not currently represented in the model (Fig. D8). While the model shows a good fit of mean values for Europe's vegetation carbon and productivity, the correlation between modelled results and observations for the individual countries show a large spread with no simple pattern for the deviations (Fig. D1-D5). However, it is obvious that countries in the Balkans, except Albania and Greece, have modelled 545 thinning fractions higher than the reported total harvest fractions. These countries also show a worse fit to observed NAI values in a simulation with thinning compared to a simulation without thinning. In any case, including thinning in simulations improves the fit to observed national NAI values in most other countries.
New forest management functionallity in LPJ-GUESS includes the most important requirements for the improvment of modelling carbon pools and fluxes as well as the development of forest stands under future climates, but a few important additions will be 550 desirable to include in the future. These include e.g. automated continuous wood harvesting and coppice management. For a good representation of coppicing, the model should also be improved to include plant carbohydrate storage. For better representations of European forests, land-use history, including litter raking, should be included in the setup to generate more realistic soil carbon pools, using functionality already available in the model.    (Smith et al., 2001); kla:sa = leaf area to sapwood crosssectional area ratio; gmin = minimum canopy conductance; fAWC = minimum growing-season (daily temperature > 5°C) 580 fraction of available soil water holding capacity in the first soil layer; CAmax = maximum woody crown area; z1 = fraction of roots in first soil layer; rfire = fraction of individuals surviving fire; aleaf = leaf longevity; aind = maximum, non-stressed longevity; fnstorage: fraction of sapwood (root for herbaceous pfts) that can be used as a nitrogen longterm storage scalar   Competing interests. The authors declare that they have no conflict of interest.