A vertically discretised canopy description for ORCHIDEE (SVN r2290) and the modiﬁcations to the energy, water and carbon ﬂuxes

Since 70 % of global forests are managed and forests impact the global carbon cycle and the energy exchange with the overlying atmosphere, forest management has the potential to mitigate climate change. Yet, none of the land surface models used in Earth system models, and therefore none of today’s predictions of future climate, 5 account for the interactions between climate and forest management. We addressed this gap in modelling capability by developing and parametrizing a version of the land surface model ORCHIDEE to simulate the biogeochemical and biophysical e ﬀ ects of forest management. The most signiﬁcant changes between the new branch called ORCHIDEE-CAN (SVN r2290) and the trunk version of ORCHIDEE (SVN r2243) are 10 the allometric-based allocation of carbon to leaf, root, wood, fruit and reserve pools; the transmittance, absorbance and reﬂectance of radiation within the canopy; and the vertical discretisation of the energy budget calculations. In addition, conceptual changes towards a better process representation occurred for the interaction of radiation with snow, the hydraulic architecture of plants, the representation of forest management and 15 a numerical solution for the photosynthesis formalism of Farquhar, von Caemmerer and Berry. For consistency reasons, these changes were extensively linked throughout the code. Parametrization was revisited after introducing twelve new parameter sets that represent speciﬁc tree species or genera rather than a group of unrelated species, as is the case in widely used plant functional types. Performance of the new model was 20 compared

the allometric-based allocation of carbon to leaf, root, wood, fruit and reserve pools; the transmittance, absorbance and reflectance of radiation within the canopy; and the vertical discretisation of the energy budget calculations. In addition, conceptual changes towards a better process representation occurred for the interaction of radiation with snow, the hydraulic architecture of plants, the representation of forest management and 15 a numerical solution for the photosynthesis formalism of Farquhar, von Caemmerer and Berry. For consistency reasons, these changes were extensively linked throughout the code. Parametrization was revisited after introducing twelve new parameter sets that represent specific tree species or genera rather than a group of unrelated species, as is the case in widely used plant functional types. Performance of the new model was 20 compared against the trunk and validated against independent spatially explicit data for basal area, tree height, canopy strucure, GPP, albedo and evapotranspiration over Europe. For all tested variables ORCHIDEE-CAN outperformed the trunk regarding its ability to reproduce large-scale spatial patterns as well as their inter-annual variability over Europe. Depending on the data stream, ORCHIDEE-CAN had a 67 to 92 % librium due to their slow fill rates, an approach known as model spin-up (Thornton and Rosenbloom, 2005;Xia et al., 2012). For a long time, spin-ups have been performed by brute force, i.e., running the model iteratively over a sufficiently long period which allows even the slowest carbon pool to reach equilibrium. This native approach is reliable but slow (in the case of ORCHIDEE it takes 3000 years) and thus comes with Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | horizontal location following a Poisson distribution across the stand. Each PFT contains a user-defined number of model trees, each one corresponding to a circumference class. Model trees are replicated to give realistic stand densities. Following tree growth, canopy dimensions and stand density are updated (Fig. 1, arrow 3). This formulation results in a dynamic canopy structure that is exploited in other parts of the model, 5 i.e., precipitation interception, transpiration, energy budget calculations, albedo (Fig. 1, arrow 4) and absorbed light for photosynthesis ( Fig. 1, arrow 5). In the trunk version these processes are driven by the big-leaf canopy assumption. The introduction of an explicit canopy structure is thought to be a key development with respect to the objectives of the ORCHIDEE-CAN branch, i.e., quantifying the biogeochemical and 10 biophysical effects of forest management on atmospheric climate. The radiation transfer scheme at the land surface benefits from the introduction of canopy structure. The trunk version of ORCHIDEE prescribes the vegetation albedo solely as a function of LAI. In the ORCHIDEE-CAN branch each tree canopy is assumed to be composed of uniformly distributed single scatterers. Following the as- 15 sumption of a Poisson distribution of the trees on the land surface, the model of Haverd et al. (2012) calculates the transmission probability of light to any given vertical point in the forest. This transmission probability is then used to calculate an effective LAI, which is a statistical description of the vertical distribution of leaf mass that accounts for stand density and horizontal tree distribution. The complexity and computational Introduction Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | flected by the canopy at vertically discretized levels, which are then used for the energy budget ( Fig. 1, arrow 6) and photosynthesis calculations (Fig. 1, arrow 5). The canopy radiative transfer scheme of (Pinty et al., 2006) separates the calculation of the fluxes resulting from downwelling direct and diffuse light, with different scattering parameters available for near-infrared (NIR) and visible (VIS) light sources. The snow 5 albedo scheme in the trunk does not distinguish between these two shortwave bands. Therefore, the snow albedo scheme of the Biosphere-Atmosphere Transfer Scheme (BATS) for the Community Climate Model (Dickinson et al., 1986) was incorporated into the ORCHIDEE-CAN branch, since it distinguishes between the NIR and VIS radiation. The radiation scheme of Pinty et al. (2006) requires snow to be put on the soil below 10 the tree canopy instead of on the canopy itself. The calculation of the snow coverage of a PFT therefore had to be revised according to the scheme of Yang et al. (1997), which allows for snow to completely cover the ground at depths greater than 0.2 m. The parameter values of Yang et al. (1997) were used in the ORCHIDEE-CAN branch.
The ORCHIDEE-CAN branch differs from any other land surface model by the in- 15 clusion of a newly developed multi-layer energy budget. There are now subcanopy wind, temperature, humidity, longwave radiation and aerodynamic resistance profiles, in addition to a check of energy closure at all levels. The energy budget represents an implementation of some of the characteristics of detailed single site, iterative canopy models (e.g., Baldocchi, 1988;Ogee et al., 2003) within a system that is coupled im- 20 plicitly to the atmosphere. Contrary to the trunk version of ORCHIDEE (see Table 1), the new approach generates a leaf temperature, using the same vegetation and radiation profile generated in the radiation scheme above, which will be fully available when parametrisation of the scheme has been completed across test sites corresponding to the species within the model. As with the trunk version, the new energy budget is 25 calculated implicitly (Polcher et al., 1998;Best et al., 2004), so as to allow for, given the 15 min time-step, a computationally efficient and stable coupling to the atmospheric model LMDz. Parameters were derived by optimizing the model against observations from short-term field campaigns. The new scheme may also be reduced to the existing single layer case so as to provide a means of comparison and compatibility with the trunk version of ORCHIDEE. The combined use of the new energy budget and the hydraulic architecture of plants required changes to the calculation of the stomatal conductance and photosynthesis ( Fig. 1, arrow 7). When water supply limits transpiration, stomatal conductance 5 is reduced and photosynthesis needs to be recalculated. Given that photosynthesis is among the computational bottlenecks of the model, the semi-analytical procedure as available in previous trunk versions is replaced by an adjusted implementation of the analytical photosynthesis scheme of Yin and Struik (2009), which is also implemented in the latest ORCHIDEE-trunk version. In addition to an analytical solution for 10 photosynthesis the scheme includes a modified Arrhenius function for the temperature dependence that accounts for a decrease of k Vcmax and k Jmax at high temperatures and a temperature dependent k Jmax/Vcmax ratio (Kattge and Knorr, 2007). The temperature response of k Vcmax and k Jmax was parametrized with values from reanalysed data in literature (Kattge and Knorr, 2007), whereas k Vcmax and k Jmax at a reference 15 temperature of 25 • C were derived from observed species-specific values in the TRY database (Kattge et al., 2011). As the amount of absorbed light varies with height (or canopy depth), the absorbed light computed from the albedo routines is now directly used in the photosynthesis scheme resulting in full consistency between the top of the canopy albedo and absorption. This new approach replaces the old scheme which Introduction Species parameters were extracted from a wide range of sources including original observations, large databases, primary research and remote sensing products (see Sect. 4). The use of age classes is introduced through externalization of the PFT parameters as well. Age classes are used during land cover change and forest management to simulate the regrowth of a forest. Following a land cover change, biomass and 5 soil carbon pools (but not soil water columns) are mixed. The number of age classes is user defined. Contrary to typical age classes, the boundaries are determined by the tree diameter rather than the age of the trees. Finally, the forest management strategies in the ORCHIDEE-CAN branch were refined from the original forest management branch (Bellassen et al., 2010). Self-thinning 10 was activated for all forests regardless of human management, contrary to the original FM branch. The new default management strategy thus has no human intervention but includes self-thinning, which replaces the fixed 40 year turnover time for woody biomass. Three management strategies with human intervention have been implemented: (1) "High stands", in which human intervention is restricted to thinning op- 15 erations based on stand density and diameter, with occasional clearcuts. Aboveground stems are harvested during operations, while branches and belowground biomass are left to litter. (2) "Coppices" involve two kinds of cuts. The first coppice cut is based on stem diameter and the aboveground woody biomass is harvested whereas the belowground biomass is left living. From this belowground biomass new shoots sprout, which 20 increases the number of aboveground stems. In subsequent cuts the amount of shoots is not increased, although all aboveground wood biomass is still harvested. (3) "Short rotation coppices", where rotation periods are based on age and are generally very short (3-6 years). The different management strategies can occur with or without litter raking, which reduces the litter pools and has a longterm effect on soil carbon (Gimmi Introduction analytical spin-up method (see Sect. 2.1). The semi-analytical spin-up is now run for nine C pools.

Allocation
Following bud burst, photosynthesis produces carbon that is added to the labile carbon 5 pool. Labile carbon is used to sustain the maintenance respiration flux (F rm ), which is the carbon cost to keep existing tissue alive (Amthor, 1984). Maintenance respiration for the whole plant is calculated by summing maintenance respiration of the different plant compartments, which is a function of the nitrogen concentration of the tissue (Zaehle and Friend, 2010, their Eqs. 6 and 7) and subtracted from the whole-plant labile pool (up to a maximum of 80 % of the labile pool). The remaining labile carbon pool is split into an active and none-active pool. The size of the active pool is calculated as a function of plant phenology and temperature and was formalized following Ryan (1991);Sitch et al. (2003); Zaehle and Friend (2010). The remaining non-active pool is used to restore the labile and carbohydrate reserves 15 pools according to the rules proposed in Zaehle and Friend (2010). The labile pool is limited to 1 % of the plant biomass or 10 times the actual daily photosynthesis. Any excess carbon is transferred to the non-respiring carbohydrate reserve pool. The carbohydrate reserve pool is capped to reflect limited starch accumulation in plants, but carbon can move freely between the two reserve pools. After accounting for growth 20 respiration (F rg ), i.e., the cost for producing new tissue excluding the carbon required to build the tissue itself (Amthor, 1984), the total allocatable C used for plant growth is obtained (M totinc ).
New biomass is allocated to leaves, roots, sapwood, heartwood, and fruits. Allocation to leaves, roots and wood respects the pipe model theory (Shinozaki et al., 1964)  of sapwood to transport water from the roots to the leaves as well as a proportional fraction of roots to take up the water from the soil. The scaling parameter between leaf and sapwood mass is derived from: where d l is the one-sided leaf area of an individual plant, d s is the sapwood area of an 5 individual plant, k ls a parameter linking leaf area to sapwood area and, m w is the water stress as defined in Sect. 3.2. Alternatively, leaf area can be written as a function of leaf mass (M l ) and the specific leaf area (k sla ): Sapwood mass M s can be calculated from the sapwood area d s as follows: where d h is the tree height and k ρs is the sapwood density. Following substitution of Eqs.
(2) and (3) into Eq. (1), leaf mass can be written as a function of sapwood mass: where, where, k ls is calculated as a function of the gap fraction as supported by site-level observations (Simonin et al., 2006): using the gap fraction as a control of k ls more carbon will be allocated to the leaves until canopy closure is reached. Following Magnani et al. (2000), sapwood mass and root mass (M r ) are related as follows: where the parameter k sar is calculated according to Magnani et al. (their Eq. 17): where k rcon is the hydraulic conductivity of roots, k scon is the hydraulic conductivity of sapwood, k τs is the longevity of sapwood and k τr is the root longevity. Following substitution of Eq. (4) into Eq. (7) and some rearrangement, leaf mass can be written 10 as a function of root mass: where, Parameter values used in Eqs.
(1) to (9), i.e., k lsmax , k lsmin , k sar , k sla , k ρs , k rcon , k scon , 15 k τs and k τr , are based on literature review (Tables 3 and 4). The allometric relationships between the plant components and the hydraulic architecture of the plant are both based on the pipe model theory, hence, the same parameter values for the hydraulic conductivities of the plant components were used in their calculations. In this version of ORCHIDEE, forests are modelled to have k ncirc circumference 20 classes with d ind identical trees in each one. Hence, the allocatable biomass (M totinc ) needs to be distributed across l diameter classes: where M inc(l ) is the biomass that can be allocated to diameter class l . Mass conservation thus requires: where M linc(l ) , M rinc(l ) and, M sinc(l ) are the increase in leaf, root and wood biomass for a tree in diameter class l , respectively. Eqs. (4) and (8) can be rewritten as An allometric relationship is used to describe the relationship between tree height and basal area: The change in height is then calculated as where d ba(l) and d bainc(l ) are the basal area and its increment, respectively. The distribution of C across the l diameter classes depends on the basal area of the model tree within each diameter class. Trees with a large basal area are assigned more carbon for 15 wood allocation than trees with a small basal area, according to the method of Deleuze et al. (2004).
where k m is a parameter, f γ and g σ are calculated from parameters and d circ(l ) is the circumference of the model tree in diameter class l . g σ is a function of the diameter 20 distribution of the stand at a given time step.  (10) to (16) need to be simultaneously solved. An iterative scheme was avoided by linearising Eq. (15), which was found to be an acceptable numerical approximation as allocation is calculated at a daily time step, and hence the changes in height are small and the relationship is locally linear: where f s is the slope of the locally linearised Eq. (15) and is calculated as: Equations (10)- (14) and (16)-(18) are then solved for f γ . f γ distributes photosynthates across the different diameter classes and as such controls the intra-species competition within a stand. f γ thus depends on the total allocatable carbon and needs 10 to be optimised at every time step. Once f γ has been calculated, M linc(l ) , M rinc(l ) and M sinc(l ) can be calculated.
The different biomass pools have different turnover times, and therefore at the end of the daily time step, the actual biomass components may no longer respect the allometric relationships. Consequently, at the start of the time step carbon is first allocated 15 to restore the allometric relationships before the remaining carbon is allocated in the above manner.

Hydraulic architecture
The representation of the impact of soil moisture stress on water, carbon and energy fluxes has been identified as one of the major uncertainties in land surface models Introduction C assimilation through a constrain on k Vcmax in the ORCHIDEE-trunk, by a constrain based on the amount of water plants can transport from the soil to their leaves. The model calculates plant water supply according to the implementation of hydraulic architecture by Hickler et al. (2006). Plant water supply is the amount of water the plant can transport from the soil to its stomata, accounting for the resistances to water 5 transport in the roots, sapwood and leaves. If transpiration rate exceeds plant water supply, the stomatal conductance is reduced until equilibrium is reached.
The water flow from the soil to the leaves is driven by a gradient of decreasing water potential. Using Darcy's law (Slatyer, 1967;Whitehead, 1998), the supply of water for transpiration through stomata can be described as: where p δ is the pressure difference between the soil and the leaves; and R r , R sap and R l are the hydraulic resistances of fine roots, sapwood and leaves, respectively. p delta is calculated following Whitehead (1998): 15 where k ψl is a PFT-specific minimal leaf water potential, which means that plants are assumed to maximise water uptake by lowering their k ψl to the minimum, if transpiration exceeds F Trs (Tyree and Sperry, 1989). The product of d h , k ρw and k g accounts for the loss in water potential by lifting a mass of water from the soil to the place of transpiration at height d h , k ρw is the density of water, and k g is the gravitational constant. The soil 20 water potential in the rooting zone (p ψsr ) was calculated by adding a modulator (m ψ ) to the bulk soil water potential, which was calculated as the sum of the soil water potential in each soil layer weighted by the relative share of roots (d rd ) in the individual soil layer: The soil water potential for each layer p ψsl is calculated from soil water content according to Van Genuchten (1980).
where M swc is the volumetric soil water content, k swcr and k swcs are respectively the residual and saturated soil water content and k av , k mv and k nv are parameters.

5
Root resistance is related to the root mass and thus can be expressed as Weatherly (1982): where k rcon is the fine root hydraulic conductivity per unit biomass. Sapwood resistance is calculated according to Magnani et al. (2000): where k scon is the sapwood specific conductivity, which is decreased when cavitation occurs. The loss of conductance as a result of cavitation is a function of p ψsr and was implemented by using an s-shaped vulnerability curve: 15 where k ψ50 is the p ψsr that causes 50 % loss of conductance; and k c is a shape parameter. R l is related to the specific leaf conductivity per unit leaf area (k l ) and the leaf area index: The response of water viscosity to low temperatures increases the resistance (Cochard et al., 2000). The relationship is described as: where k α1v and k α2v are empirical parameters (Cochard et al., 2000), R temp is the temperature adjusted R l , R sap or R r , T is air temperature for R l and R sap and T is soil 5 temperature for R r . If, for any time step, the transpiration calculated by the energy budget exceeds the amount of water the plant can transport from the soil to its stomata, transpiration is limited to the plant water supply. As the transpiration is now reduced, the initial calculations of the energy budget and photosynthesis, solely based on atmospheric infor-10 mation, are no longer valid. As a result the energy budget and photosynthesis must be recalculated for the time step in question. For this recalculation, stomatal conductance at the canopy level is calculated such that transpiration equals the amount of water the plant can transport. Owing to the feedback between stomatal conductance, leaf surface temperature and transpiration this calculation may require up to 10 iterations to 15 converge. Canopy level stomatal conductance is then decomposed to obtain the stomatal conductance at each canopy layer assuming that each layer is equally restricted by drought stress. Finally, the restricted stomatal conductance is used to calculate CO 2 assimilation rate according to the photosynthesis model by Farquhar, von Caemmerer and Berry (see Sect. 3.6).

Canopy structure
Stand structure controls the amount of light that penetrates to a given depth in the canopy. For example, the amount of light reaching the forest floor will be higher for a stand with few mature trees compared to many young trees even if both stands have the same leaf area index. Where a big leaf approach assumes a homogeneous 25 block shaped canopy (as in the trunk version of ORCHIDEE) and can therefore rely 8584 Introduction  Sect. 3.5) and the amount of light reaching the forest floor (see for example Sect. 3.1). The gap fraction, which is the basic information in calculating light penetra-5 tion at different depths in the canopy, is calculated following the approach presented by Haverd et al. (2012) and formalized in their semi-analytical model. Rather than a spatially explicit approach, Haverd et al. (2012) follow a statistical approach which reduces the memory requirements for the simulations and limits the space requirements for storing the model output files. 10 The model of Haverd et al. (2012) represents the canopy by a statistical height distribution with varying crown sizes and stem diameters for each height class. The crown canopies are treated as spheroids containing homogeneously distributed single scatterers. Although this f Pgap model can explicitly include trunks, we made the decision to exclude them, as the spectral parameters for our radiation model (see Sect. 3.4) 15 are extracted from remote sensing data (see Sect. 4.8) without distinguishing between leafy and woody masses. This gives the gap probability for trees as a function of height (z) and solar zenith angle (θ z ): where d λ is the inverse of the tree density, d c is the projected crown area (for an opaque 20 canopy), and f Pwc is the mean crown porosity. The overbar depicts the mean over the tree distribution as a function of tree height or, in our case, the mean over the l circumference classes. Following minor adaptations, the implementation of Haverd and Lovell  was incorporated in ORCHIDEE-CAN. As there also exist crops, grasses, and bare soil in the model, f Pgap was adjusted for these situations 25 as well. For grasses and crops, the same formulation is used: where d LAIabove is the total amount of LAI above height z, and m LAIcorr is a correction factor to account for the fact that grasses and crops are treated as homogeneous blocks of vegetation with no internal structure and is often referred to as a clumping factor. Here it is treated as a tunable parameter and therefore the term correction factor was used. For bare soil, there is no vegetation to intercept radiation, and therefore 5 f bs Pgap (θ z , z) is always unity.

Multi-layer two-way albedo for tall canopies
Species-specific radiation absorbance, reflectance and transmittance by the forest canopy were calculated from a radiation transfer model (Pinty et al., 2006) which was parametrized by satellite-derived species-specific scattering values (see Sect. 4.8).
Given the complexity of radiation transfer, it remains challenging to accurately simulate radiation transfer through structurally and optically complex vegetation canopies without using explicit 3-D models. The applied 1-D model belongs to the family of twostream models (Meador and Weaver, 1980) and thus calculates transmittance, absorbance and reflectance of both the incoming and outgoing radiation. The calculation 15 of the reflectance at the top of the canopy due to a collimated source (i.e., the sun) is divided in three components: 1. scattering of radiation between the vegetated elements with a black background 2. scattering of radiation by the background with a black canopy 3. multiple scattering of radiation between the canopy and the background Term (1) is widely used in cloud reflectance calculations, and depends on the cosine of the solar zenith angle (θ mu ), the reflectance and transmittance of the single leaves (f rl and f tl , respectively), the leaf orientation function (g G ), and the effective LAI (d LAIeff ). The exact definition of this term is given in Eq. (B2) in Pinty et al. (2006). In term (2), f Rbgd is the reflectance of the ground beneath the canopy and f T UnColl, veg is the transmit-5 ted fraction of light to the ground which has not collided with any canopy elements. In term (3), f f R Coll, bgd,1 is the fraction of light which has struck vegetation and collided with the background a single time, while f f R Coll, bgd,n is the fraction which has collided multiple times (n) with the background. The sum of the three components results in the canopy albedo (Pinty et al., 2006). Similar equations can be derived for light originating from 10 diffuse sources (e.g., clouds and other atmospheric scattering). Implementations of the calculations of the canopy fluxes for a single level are available from the JRC, and these implementations were used as the basis of the routines put into ORCHIDEE-CAN for both the single-and multi-level cases . This implementation relies on the use of the effective LAI, which is the LAI that needs to be used in a 1-D process 15 representation to obtain the same reflectance, absorbance and transmittance as would be obtained by a 3-D-canopy representation (Pinty, 2004). In this study, the effective LAI was calculated by first computing the canopy gap probability, i.e. the probability that light is transmitted to a specified height in the canopy at a given solar angle. The gap probability is then converted into the effective LAI by passing it as an input to the 20 inverted Beer-Lambert's law (with an assumed extinction coefficient of 0.5). cally applies the 1-D-two stream canopy radiation transfer model by Pinty et al. (2006) to each canopy level where the light transmitted by the overlaying level becomes the input for the lower level.
As the multi-level approach is built around the solution of the one-level scheme for each canopy level, no new equations are introduced. The method can be summarized 5 by the following algorithm for which the details are given in . First, three fluxes are calculated for each level independently: the fraction of light transmitted through the layer without striking vegetation, the fraction of light reflected after striking vegetation, and the fraction of light transmitted through the layer after striking vegetation. These three fluxes represent the only possible fate of light (any light not taking one of these paths must be absorbed for energy conservation). Next, an iterative approach is invoked which follows the path of a single photon entering the top level. Based on the solutions for each single level, probabilities can be calculated that the photon is transmitted to a lower level or reflected to a higher level. Any fraction which is reflected upwards from the top level is added to the total canopy albedo and not considered 15 further. The fraction which is transmitted through the top level enters the next highest level, and again the single level solutions determine where this light goes. Any fraction reflected upwards is considered in the next iteration as part of the light entering the upper level. The steps continue until the bottom canopy level is reached. Here, any fraction which is transmitted into the soil is removed from consideration and added to 20 the total transmittance through the canopy. The algorithm then proceeds to the above canopy level. Now the "transmitted" fluxes are moving in the upwards direction towards to the sky, while "reflected fluxes" are moving towards the ground. The code continues towards the top level, taking as input from below both the flux reflected by downwelling light from the level below the current level and the flux transmitted from the lower level 25 by upwelling light. After each iteration (moving from the top of the canopy to the bottom and back to the top), the total amount of light considered "active" has been reduced by light escaping to the sky or being absorbed by the canopy or ground. Eventually, this "active" light falls below a pre-defined threshold and the calculation is considered as converged.
Due to the iterative procedure, energy is not strictly conserved, although we have attempted to choose a threshold which minimizes this loss. The multilevel albedo calculation is currently the most expensive part of the model, due to the iterations and the fact that it must be performed over all canopy levels (currently set to 10), grid points, and PFTs at every physical time-step. Levels with no LAI are no less expensive to compute, either, although we have arranged our canopy levels to make sure no levels are empty in most cases. 10 The present generation of land surface models have difficulties in reproducing consistently the energy balances that are observed in field studies (Pitman et al., 2009;Jiménez et al., 2011;de Noblet-Ducoudré et al., 2012). The ORCHIDEE-CAN branch implemented an energy budget scheme that represents more than one canopy layer to simulate the effects of scalar gradients within the canopy for determining more ac- 15 curately the net sensible and latent heat fluxes that are passed to the atmosphere. As outlined in Polcher et al. (1998), the use of an implicit solution for coupling between the atmospheric model and the surface layer model is the only way to keep profiles of temperature and humidity synchronised across the two models when the coupled-model is run over large time steps (e.g., of 30 min). The difference between explicit and im-20 plicit schemes is that an explicit scheme will calculate each value of the variable (e.g., temperature and humidity) at the current time step entirely in terms of values from the previous time step. An implicit scheme requires the solution of equations written only in terms of those at the current time step.

Multi-layer energy budget
The modelling approach formalises three constraints that ensure energy conserva- 1. the energy balance at each layer is the sum of incoming and outgoing fluxes of latent and sensible heat and of shortwave and longwave radiation: where F LW,i is the sum total of long wave radiation, that is, the net LW radiation 5 absorbed into layer i and F SW,i is the net absorbed short wave radiation as calculated by the albedo scheme in Sect. 3.4. k shc is the specific heat capacity of air. The source sensible heat flux from the leaf at level i is the difference between the leaf temperature (T L,i ) and the atmospheric temperature at the same level (T a,i ), divided by R a,i , which is the leaf resistance to sensible heat flux (a combination 10 of stomatal and boundary layer resistance). Similarly, the source latent heat flux from the leaf at level i is the difference between the saturated humidity in the leaf (q L,i ) and that in the atmosphere at level i (q a,i ), divided by R s,i which is the leaf resistance to latent heat flux. R a,i is calculated based upon the leaf boundary layer resistance, and is described in the present model according to Baldocchi (1988). 15 R s,i is the stomatal resistance of the leaf that may be calculated using the model of Ball et al. (1987).
2. Transport of sensible heat flux between the vegetation ("the leaf") and the surrounding atmosphere at each level, and between adjacent atmospheric levels above and below is provided by the following expression: where z denotes the height above the soil suface. We have re-written the scalar conservation equation, as applied to canopies, in terms of the sensible heat flux, temperature and source sensible heat from the vegetation at each layer. 3. The transport of latent heat flux between the vegetation and surrounding atmosphere at each level, and between adjacent atmospheric levels above and below is described in a form that is analogous to Eq. (36), above: In addition to these three basic equations various terms had to be parameterised.

5
The 1-D second-order closure model of Massman and Weil (1999) was used to simulate the vertical transport coefficients k k,i within the canopy while accounting for the vertical and horizontal distribution of LAI (see Sect. 3.3). This set of equations were then written in an implicit form and solved by induction. More details on the implicit multi-layer energy budget and a complete mathematical documentation are given in Ryder et al. (2014).
To complete the energy budget calculations, the multi-layer 1-D canopy radiation transfer model (see Sect. 3.4) was used to calculate the net shortwave radiation at each canopy layer. Further, the canopy radiation scheme makes use of the Longwave Radiation Transfer Matrix (LRTM) (Gu, 1988;Gu et al., 1999). This approach separates 15 the calculation of the radiation distribution completely from the implicit expression. Instead, a single source term for the long wave radiation is added at each level. This means that the distribution of LW radiation is now explicit (i.e., makes use of information only from the "previous" and not the "current" time step) but the changes within the timestep were small enough to not affect the overall stability of the model). However, 20 an advantage of the approach is that it accounts for a higher order of reflections from adjacent levels than the single order assumed in the process above.

Analytical solution for photosynthesis
The photosynthesis model by Farquhar, von Caemmerer and Berry (Farquhar et al., 1980)  Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | rate of CO 2 assimilation and the electron transport-limited rate of CO 2 assimilation (Farquhar et al., 1980). The ORCHIDEE-CAN branch calculates net photosynthesis following an analytical algorithm as described by Yin and Struik (2009). In addition, the C4 photosynthesis is calculated by an equivalent version of the Farquhar, von Caemmerer and Berry model that was extended to account for noncyclic electron transport 5 (Yin and Struik, 2009). A detailed derivation of the analytical solution of the Farquhar, von Caemmerer and Berry model is given in Yin and Struik (2009).
Owing to the canopy structure simulated in this model version and the layering of the canopy, the amount of absorbed light now varies with canopy depth. This new approach replaces the old scheme which uses multiple levels based on the leaf area index, not the physical height within the canopy. Photosynthesis is now calculated at each vertically resolved canopy level independently, using the total amount of absorbed light calculated by the radiation transfer scheme, which means that radiation transfer inside the canopy and photosynthesis are now fully consistent. In the new photosynthesis scheme, photosynthesis thus indirectly depends on canopy structure.

Mortality
A variety of changes have been made to processes involving vegetation death. A whole PFT is now killed if, at the end of the day, there is no carbon available in leaf, carbohydrate reserve, or labile pools. In this situation, it will be impossible for the plant to assimilate new carbon from the atmosphere as it will not be able to grow new leaves 20 and thus initiate plant recovery. In addition, a forest can die if the density falls below a certain prescribed value. In the next time step a new young forest will be prescribed. Different age classes are distinguished to better account for the structural diversity and its possible effects on the element, energy and water fluxes. A clear hierarchy was established for the mortality processes when it comes time to actually kill the trees (i.e., 25 move their biomass to the litter or harvest pools). All of the processes determine first how much biomass they would remove in the absence of all the other processes. Afterwards, the killing is arranged in the most realistic way possible. A clear-cut event 8592 Introduction has the highest priority, followed by human thinning and finally natural mortality including self-thinning. If, for example, a forest is scheduled to be clear-cut, the entire forest biomass is subjected to the rules of the clear-cut and no other mortality occurs in that time-step. If a forest is thinned, it is assumed that the weakest trees will be thinned, and there-5 fore human thinning reduces or even eliminates the natural mortality for that time-step. Natural mortality still happens on a daily time-step, while human-induced mortality happens only at the end of the year. Self-thinning, as described below, takes priority over environmental mortality. Environmental mortality is calculated as a fraction of the total site biomass; if self-thinning is greater than or equal to this percentage, no environ-10 mental mortality occurs in this time-step. Otherwise trees are selected so that the total amount of biomass killed due to self-thinning and environmental mortality is equal to the total amount predicted by multiplying the stand biomass by the mortality percentage. The use of circumference classes adds a good deal of realism and flexibility to 15 the ORCHIDEE-CAN simulations, but it also raises additional questions. For example, which trees should be targeted by which mortality? Given that self-thinning reflects the outcome of continuous resource competition, the largest trees are expected to be most successful when competing for resources, and therefore we mainly kill the smallest trees to reduce the stand density. Conversely, larger trees are more likely to die be-20 cause of environmental stress factors, being more prone to cavitation, wind damage, lightening, and, heart rot. Therefore, we select more older trees to die from environmental mortality. While doing this also trees in the other diameter classes were killed based on the following recursive definition (cf. Bellassen et al., 2010): where k ddf is the death distribution factor, which is the factor by which the smallest and largest circumference classes differ (e.g., k ddf = 10 means that the largest circum-8593 Introduction Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | ference class will lose ten times as much biomass as the smallest as a result of the mortality), m Ndeath is a normalization factor so that sum of f icir death is unity, and f 1 death is set equal to unity before normalization. As the stands are very close to even-aged, we set the factor k ddf to be equal to 1. This means the same number of trees is killed in each circumference class If, for some reason, there is not enough biomass in a given class 5 to satisfy this distribution, the extra biomass is taken from the next smallest class (in the case the smallest class does not have enough, it is taken from the largest class).
Related to mortality is the question of the circumference class distribution. As mentioned above, trees in different circumference classes are preferentially killed by different processes. If the simulation is long enough (or if the morality is aggressive enough), eventually the number of trees in some circumference classes may become zero. This would reduce the numerical resolution of the allocation scheme. When only one circumference remains populated, the scheme effectively loses its meaning as all the newly produced biomass is now be allocated to the only remaining circumference class. In order to maintain the same level of detail through the simulation, the distribution of all 15 the circumference classes is recalculated at the end of each day. A normalized target distribution is specified as an input parameter (an exponential distribution is currently used), and this distribution is scaled to produce a target distribution for the current number of individuals. All of the current individuals are placed in these new classes until the target distribution is satisfied. The target distribution now contains, however, trees Introduction Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | a temperate and boreal needleleaf evergreen and a temperate and boreal broadleaved residual group. For use outside Europe, the original MTC classification of ORCHIDEE was kept. The parameters of the residual groups and MTCs are the mean of the parameters of the species-level PFTs that are in the MTC, with the exception of albedo parameters that could be extracted from remote sensing products. 5 Finally, separate PFTs were introduced for boreal grasses and croplands, which allowed for a boreal parametrization of phenology, senescence and growth. This approach, which distinguishes a total of 28 PFTs, allows a higher taxonomic resolution over Europe, better defines forest types compared to the more general MTC approach and facilitates the use of observations to derive parameters.

Sources of parameter values
Species parameters were extracted from a wide range of sources including original observations (i.e. POPFULL), large databases (national forest inventories, TRY), primary research reports and remote sensing products (JRC-TIP Pinty et al., 2011a, b). The mean values and standard deviations were calculated without weighting. For most 15 parameters, the mean values were used as the parameter value in ORCHIDEE without further processing. Using the mean parameter estimates at the species level avoids hidden model-tuning and largely reduces the likelihood that simulation results are biased by hidden calibration owing to a poor taxonomic definition of PFTs (Scheiter et al., 2013). The phenology-related parameters of the deciduous MTCs were optimised by Introduction Dependencies of k ls on tree height Novick et al., 2009), tree diameter following stand thinning (Simonin et al., 2006) and CO 2 (Pataki et al., 2006) have been reported. Most observations come from experiments where time was substituted by space which hampers teasing apart the sources of variability. Given the variation and uncertainty in the observations and the model sensitivity to this parameter, we 5 tuned its value within the observed range to jointly match European-wide observations of GPP, evapotranspiration and effective leaf area index. Both the trunk and ORCHIDEE-CAN branch reduce the definition of net primary production to biomass production; hence, carbon leaching from the roots, volatile organic emissions from the leaves, dissolved and particulate carbon losses through water 10 fluxes and carbon subsidies to mycorryhzae are not accounted for in the model. These fluxes are (incorrectly) accounted for in the modelled autotrophic respiration. Modelled autotrophic respiration should therefore be considered an effective rather than a true value. For this reason, the basal rate of autotrophic respiration was optimized against site observations of biomass production efficiency (i.e., the ratio between an-15 nual biomass production and annual photosynthesis), using an optimization scheme that minimizes the mismatch between the model and the observations in a rigorous statistical framework (see for example Tarantola, 2005).

Allocation
In addition to k ls the leaf to sapwood area, the allocation scheme required several new 20 parameter values that were estimated from fitting regression models to the national forest inventory data of Spain, France, Germany and Sweden or through literature search (Table 4).
Strictly speaking, the basal rate of autotrophic respiration is not an allocation parameter but here it was optimized against 126 site observations of the biomass production 25 efficiency (k cmaint , Table 3) calculated as the ratio between annual biomass production and annual photosynthesis (Vicca et al., 2012;Campioli et al., 2014). Where the model now simulates an acceptable GPP, it will also simulate an acceptable biomass 8597 Introduction production which justifies the inclusion of this respiration parameter in the section on allocation. Following this approach it remains untested how well the simulated effective autotrophic respiration represents the (rarely) observed autotrophic respiration. Note that in the cases of both the trunk and the ORCHIDEE-CAN branch of ORCHIDEE, a match between effective and observed autotrophic respiration should not be inter-5 preted as evidence of desired model behaviour because several components of net primary production are not modelled yet (see Sect. 4.3).

Hydraulic architecture
Initial choices of parameters for this scheme were based on the values and parameter sources listed by Hickler et al. (2006). All data sources were revisited and the search was extended to obtain values at the PFT rather than MTC level (Table 4). Given that plant hydraulogy is rather well studied, observed parameters were available for most of the species. Our implementation of hydraulic architecture required the introduction of a tuning parameter to account for processes that are currently absent in, e.g. plant water storage and soil-root resistance. A more process-based description of these pro-15 cesses (i.e., Sperry et al., 1998;Steppe et al., 2005) is being tested and should allow removal of this parameter.

Canopy structure
The relationship between diameter and projected crown surface area follows the model proposed by Pretzsch (2009) between d csa and d dbh was fitted to all available data for Pinus sylvestris. Next, all observations with a d csa that falls below the predicted d csa were selected as considered to represent a boreal subset. Given the importance of snow pressure on crown structure, selecting observations with sub average d csa is justifiable as a first approximation. Subsequently, the parameters were fitted to this subset of data. For Quercus ilex no data 10 were available and parameters were tuned such that the crown diameter was 0.85 m less than the tree height.

Multi-layer two-way albedo for tall canopies
The radiation transfer scheme makes use of parameters describing leaf and background properties, i.e., leaf single scattering and prefered scattering direction (for  Only parameter values for which the posterior standard deviation of the probability density functions were significantly smaller than the prior standard deviation were selected from the JRC-TIP since this condition ensures statistically significant values. Species and MTC specific values were derived from JRC-TIP by performing a multiple regression (Table 6). This methods determines, in an objective way, how the fractions of each MTC or species explain the JRC-TIP parameter. The multiple regression was performed separately for the six parameters: the single scattering of leaves (for both VIS and NIR), the scattering direction of leaves (VIS and NIR) and the background 5 albedo (VIS and NIR). Each JRC-TIP parameter was used as the dependent variable and the independent variables consisted of the fractions of each MTC (Poulter et al., 2014) or species (Brus et al., 2012). These fractions were used to find a linear function that best predicted each JRC-TIP parameter. The corresponding slope of a regression of each MTC or species fraction gives the MTC or species dependent JRC-TIP 10 value. The multiple regression was performed without an intercept. To avoid pollution by the seasonal cycle, the multiple regression was applied only for the pixels of the Northern Hemisphere. Only pixels that were less than 10 % covered by non-vegetative fractions where selected for the analysis and only significant results following an F-test and positive r 2 values were selected.

Analytical solution for photosynthesis
Three originally MTC-specific photosynthetic parameters (k Vcmax , k Jmax and k sla ) were derived at the species level by obtaining weighted site means for each species from the global leaf trait database TRY (Kattge et al., 2011) and additionally from Medlyn et al. (2002). Only k Vcmax and k Jmax standardized to a common formulation and 20 parametrization of the photosynthesis model by (Farquhar et al., 1980)  Depending on the species this resulted in 5 to 183 observations for k sla and 11 to 173 observations for k Vcmax,opt and k Jmax,opt . From these observations species-specific means were calculated, weighted for differences in the number of observations per site (Table 5). 5 Tree mortality is controlled by (see Sect. 3.7): (1) maximum tree diameter, (2) minimum stand density, (3) environmental mortality, (4) self-thinning and, (5) anthropogenic thinning. Maximum tree diameter was extracted from the French, Swedish, German and Spanish forest inventories as the observed 50 % quantile for diameter at breast height. The 50 % quantile rather than the observed maximum was used to account for the fact 10 that large scale land surface models are expected to reproduce large scale patterns rather than local extremes. Minimum stand density was estimated as the expected stand density for the maximum tree diameter for a stand under self-thinning. Although both criteria are related to each other through the observed self-thinning relationship, the minimum number of trees is used to decide when unmanaged forests should be 15 replaced, whereas both the maximum diameter and the minimum number are used for managed sites as criteria to initiate a clear cut. Self-thinning was parametrized based on the observed relationship between stand density and quadratic mean stand diameter. These observations include mortality due to intra-stand competition which strictly speaking is the process described by the self-20 thinning relationship as well as the mortality of individuals by insects, lightening, wind, drought, frost and heart rot; these are referred to here as environmental mortality. Because the national forest inventory data rarely distinguish between these causes of mortality we introduced a two step approach. First, self-thinning mortality is calculated. Where self-thinning is less than an assumed constant environmental mortality of 25 1/k t resid , self-thinning is complemented by additional mortality to reach the set environmental mortality. Where self-thinning mortality exceeds the set environmental mortality, simulated self-thinning is assumed to include environmental mortality. The fire module 8601 Introduction

Tables Figures
Back Close

Full Screen / Esc
Printer-friendly Version

Interactive Discussion
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | that is available for the trunk but not for the ORCHIDEE-CAN branch simulates stand replacing fires rather than individual-tree based mortality due to lightening. The approach implemented in the ORCHIDEE-CAN branch could therefore be extended with models that simulate stand replacing mortality from fire, insects and storms. Resource competition between trees in the same stand has been reported to result 5 in the so-called self-thinning relationship that relates the number of individuals within a stand to the stand biomass (Reineke, 1933;Kira et al., 1953;Yoda et al., 1963): where k α and k β are the constants of the self-thinning relationship. Furthermore, stem volume can be written as a function of tree diameter (d dbh ), tree height and stem form 10 factor (k α ) to account for the fact that the stem shape is not a perfect cylinder: Following the allometric relationship given in Eq. (14), tree height can be written as a function of tree diameter. Hence, the self-thinning relationship can be re-written to relate stand diameter to stand density: where, k β2 relates to k β1 (as in Eq. 14) as follows: k α1 and k β1 were estimated by fitting Eq. (14) to observed diameter and height of individual trees from NFI of Sweden, Germany, France and Spain. k β2 was calculated 20 from Eq. (43) and k α2 was estimated by fitting Eq. (42) to observations of the quadratic mean stand diameter and stand density from NFI data.

Validation
ORCHIDEE-CAN is designed as the land surface model to be coupled to the atmospheric model LMDz. As such future applications of ORCHIDEE-CAN are expected to be regional to global in the spatial domain and to span several years in the temporal domain. Given its anticipated uses, the ability of the model to reproduce large-scale 5 spatial patterns as well as their inter-annual variability is essential. The first applications of the model, both offline and coupled to the atmosphere, will focus on Europe. indices are calculated for winter, spring, summer and autumn and thus allow to evaluate the capacity of the model to reproduce observed annual cycles.
In addition to the root mean square error, a land performance index (LPI) based on the principles laid out for the Climate Performance Index (Murphy et al., 2004, their SI) was also calculated. LPI normalizes the root of the squared differences between the 15 simulations and observations by the observed spatial and temporal variance. The LPI was used to estimate the likelihood that the simulated variable belongs to the same population as the observed variable, defined as exp(−0.5LPI 2 ). An LPI equal to 1 indicates that the model correctly reproduces the mean observed value and implies a likelihood of 61 % (Murphy et al., 2004) that the simulations and observations come from 20 the same population. Similarly, an LPI of 2 reduces this likelihood to 13 %. An LPI of less than 0.32 has a likelihood of more than 95 % and therefore indicates a statistically significant result.
While developing ORCHIDEE-CAN, the numerical approaches that added functionality to the code were selected on the basis of their performance at the site-level (see 25 below). Rather than running the same site-level tests for our implementation, we performed a complementary large-scale validation. The strength of our approach lies not in the details, as is the case for site-level validation, but in its width by simultaneously  (Simard et al., 2011), biogeochemical fluxes such as GPP (Jung et al., 2011), biophysical fluxes such as albedo (Schaaf et al., 2002) and fluxes at the interface of biogeochemistry and biophysics such as evapotranspiration (Jung et al., 2011). The selection of variables was 5 limited by the availability of spatially explicit data-derived products for Europe. For the validation, both the trunk and ORCHIDEE-CAN branch were run from 1850 to 1900 using CRU-NCEP climate forcing from 1901-1950 at 0.5 • resolution. From 1901 until 2012, the corresponding CRU-NCEP forcing data for each year were used. Both versions used the 11-layer soil hydrology, the single-layer energy budget and the same 10 land cover map (Poulter et al., 2014). Given that no European-wide, spatially explicit and data-derived products were found for the validation of the net carbon flux, there was no need for a carbon spin-up. For the ORCHIDEE-CAN branch, the observed tree height and basal area were compared against the simulation values at the end of 2010 (the trunk does not simulate these variables). For both the trunk and the ORCHIDEE- 15 CAN branch, the observed GPP, evapotranspiration, effective LAI and VIS and NIR albedos were compared against monthly means between 2001 and 2010.

Allocation
In ORCHIDEE-CAN, functional relationships which vary by species and light stress are used to allocate carbon among the fine roots, foliage and sapwood. The allocation 20 scheme largely follows Zaehle and Friend (2010), who in turn was inspired by Sitch et al. (2003). Approaches simulating allocation based on functional relationships were found to outcompete allocation schemes based on constant fractions or resource limitation (De Kauwe et al., 2014). The ability of these schemes to reproduce foliage, fine root and sapwood reported in large observational data sets (for example, Luyssaert et al., 2007) demonstrates that these schemes capture the main observed features (Zaehle and Friend, 2010). In addition, allocation schemes making use of functional relationships were also capable of simulating the observed effect of elevated CO 2 on 8604 Introduction

Tables Figures
Back Close

Full Screen / Esc
Printer-friendly Version

Interactive Discussion
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | two mature forest ecosystems (De Kauwe et al., 2014). Despite these successes, the schemes were reported to be sensitive to their parametrization. Differences in parameters were reported to result in substantial differences in the simulated allocation. The parameters for the functional relationships used in ORCHIDEE-CAN are given in Table 4. The main conceptual difference between the allocation scheme by Zaehle and Friend (2010) and ORCHIDEE-CAN is that the latter was designed to simulate one or more diameter classes. Given that photosynthesis is still calculated at the stand level (and thus not at the tree level) the allocation rule of Deleuze et al. (2004) was integrated in the functional allocation scheme to account for light and resource competition within a stand (see 10 Sect. 5.6). Where the functional relationships are used to simulate carbon allocation within an individual tree of a given diameter, the rule of Deleuze and Dhote allocates carbon across the different diameter classes. The allocation rule which models the radial increment for individual trees in pure even-aged stands was successfully tested for Norway spruce and Douglas fir stands in France (Deleuze et al., 2004). A similar ap-15 proach for modelling radial increment has already been implemented in a version close to the trunk of ORCHIDEE (Bellassen et al., 2010) and was able to successfully simulate stand characteristics such as height, basal area and stand diameter (Bellassen et al., 2011). This previous implementation differs from the current implementation in its time resolution (which is now daily instead of yearly), its analytical solution and the 20 underlying allocation scheme (which is now based on functional relationships instead of resource limitation).
The aforementioned studies performed a detailed validation of the two approaches dealing with carbon allocation, which were combined in ORCHIDEE-CAN. Complementary to these studies we performed a European-wide validation of our implementa- Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | its capacity to simulate the foliage biomass, a correct simulation of height reflects the model's capacity to simulate aboveground woody biomass and its capacity to reproduce observed basal areas suggest that the interaction of stand density and individual tree diameter are well-captured.
The new implementation and parametrization of the within-tree and within-stand allo-5 cation schemes were found to have an 91, 68 and 72 % chance that the simulations reproduced the observations for GPP, tree height and basal area for Europe, respectively ( Table 7). Given that basal area and height are not available from the trunk version of ORCHIDEE, we could not compare the performance of model versions in this respect. With respect to GPP, the ORCHIDEE-CAN branch was found to outperform the trunk by 12 % and thus increased the likelihood that ORCHIDEE-CAN is an unbiased simulator of the spatial and temporal variability of GPP from 79 to 91 %. Improved performance of the ORCHIDEE-CAN branch compared to the trunk is observed for all regions in summer where the RMSE of GPP was halved from 2.5-5 to 1-2 gC m −2 day −1 (Figs. 2 and 3). 15 Although part of the high likelihood could be due to the fact that the observed GPP was upscaled making use of similar climatologies being used as the forcings of the models, this circularity could neither have contributed to the improved performance between the trunk and the ORCHIDEE-CAN branch nor to the decrease in RMSE. The improvements are thought to be due to structural changes to the model such as 20 allocation, hydraulic architecture and canopy structure as well as to the use of more consistent parametrization as the ORCHIDEE-CAN branch makes use of tree species rather than plant functional types.

Plant water supply
Our implementation of plant hydraulic architecture was largely based on the scheme 25 of Hickler et al. (2006), which was tested globally and at site level. Global simulation results for actual evapotranspiration were found to reproduce available data (Baumgartner and Reichel, 1975;Henning, 1989 with the magnitude and seasonality of eddy-covariance measurements of actual evapotranspiration for 15 European forests sites (EUROFLUX), with a tendency to slightly overestimate actual evapotranspiration for 6 sites (Hickler et al., 2006). The maximum amount of water that can be transported by a tree relies on the hydraulic architecture of the tree and therefore on the capacity of the model to simulate 5 tree and stand dimensions as well as on the model's capacity to simulate soil water content. As an additional test, our implementation of the model was compared against the upscaled eddy-covariance measurements for GPP and actual evapotranspiration (Jung et al., 2011). The capacity to jointly reproduce GPP and actual evapotranspiration is an indicator that the model successfully reproduces the coupling between CO 2 10 and water exchange. Model validation showed 91 and 87 % chance (compared to 79 and 45 % for the trunk) that ORCHIDEE-CAN reproduces the upscaled GPP and actual evapotranspiration data ( Table 7). The RMSE for actual evapotranspiration during summer dropped well below 1 mm day −1 for most regions (Fig. 2), whereas it never dropped below 1 mm day −1 for the trunk (Fig. 3). 15

Canopy structure
The canopy structure model by Haverd et al. (2012) was previously validated against ground-based LIDAR data for several test sites with varying density, structural complexity, layering and clumping . Model-derived canopy gap probabilities compared with observations using a one-sample t-test were significant for 11 out of 12 20 test sites. We considered this result as a sufficient proof to use this canopy structure model in the ORCHIDEE-CAN branch and added to its validation by comparing the simulated canopy structure model over Europe against a remote-sensing based map of tree height (see Sect. 5.1, Simard et al., 2011) and the JRC-TIP effective LAI product (Pinty et al., 2011a). The effective LAI value expresses the capability of the canopy to 25 intercept direct radiation, and is thus associated with the probability distribution function of the canopy gaps  branch, canopy structure is used to calculate the albedo, roughness length, absorbed light for photosynthesis and leaf area that is coupled to the atmosphere for, e.g., transpiration and interception of precipitation. The ORCHIDEE-CAN branch is the first branch of ORCHIDEE that makes use of an effective LAI to calculate the interaction between the canopy and the atmosphere. The 5 LPI and RMSE of the branch, therefore, cannot be compared against the trunk. Overall, the combined implementation of the allocation scheme and the canopy structure model shows a 67 % chance to reproduce the satellite-based estimates for effective LAI. Surprisingly, effective LAI is better simulated in spring and autumn when dynamics within the canopy are substantial due to leaf on-set and senescence. For the periods when the effective LAI is expected to be most stable, i.e., summer and winter, LPI approached and frequently exceeded 1 (data not shown). Part of this shortcoming may be due to the lack of shrubs in the land cover classification. In the model, shrublands are replaced by forest and/or grasslands, likely resulting in differences between the observed and simulated canopy structure. This lapse also appears in the RMSE of effective LAI 15 (RMSE higher than 0.8, Fig. 2).

Top of the canopy albedo
The radiation transfer model (Pinty et al., 2006) has been validated extensively against realistic complex three-dimensional canopy scenarios (Pinty et al., 2006) and as part of the RAdiation transfer Model Intercomparison (RAMI) project. The 1-D canopy radiation 20 transfer model by Pinty et al. (2006) was demonstrated to accurately simulate both the amplitude and the angular variations of all radiant fluxes with respect to the solar zenith angle (Widlowski et al., 2011). In addition, the radiation transfer model and its effective values extracted from the JRC-TIP data set were successfully applied to a single forest site (Pinty et al., 2011c). 25 Previously we reported on the capacity of the radiation transfer model to simulate the effects of forest management on albedo . For the latter, forest properties were prescribed and the radiation transfer model was validated against top-of- the-canopy albedo data from five observational sites. Differences in the spatial scales between the observed and simulated albedo values were accounted for by presenting the mean June albedo during 2001-2010. The simulated summertime canopy albedo falls within the range of observation. However, there occurs a slight overestimation in the near-infrared wavelength band compared to the single site mea-5 surement. Too high near-infrared single scattering albedo values for pine, as obtained from the JRC-TIP product, are the most likely cause. The observed deviation is not due to a shortcoming in the model itself but reflects the difficulties the JRC-TIP has with optimizing parameter values in the absence of field observations in the specific case of sparse canopies .
For the spatial validation we use the white-sky albedo (VIS and NIR) from Moderate Resolution Imaging Spectroradiometer (MODIS, Schaaf et al., 2002) at 0.5 • resolution (distributed in netCDF format by the Integrated Climate Data Center (ICDC, http://icdc.zmaw.de, University of Hamburg, Hamburg, Germany). Over large spatial and temporal domains the ORCHIDEE-CAN branch reproduces the observed VIS and 15 NIR albedo and its variability; LPI for the albedo in the visible light is especially satisfying with a likelihood of 92 % for the simulations to come from the same population as the observations (Table 7). This high overall performance index, however, hides performance issues over Scandinavia and the Alps during the snow season. The RMSE for VIS and NIR albedo without snow lies around 0.05, whereas during the snow season 20 the RMSE increases to 0.20 (VIS) and 0.18 (NIR) over these regions (Fig. 2). When the ORCHIDEE-CAN branch is coupled to an atmospheric model, however, these deviations will only have a minor effect on the climate, owing to low incoming radiation during most of the snow season, especially in Scandinavia.
Previous validation of the radiation transfer model showed that the largest discrepan-25 cies were occurring in the near-infrared domain with a snow covered background (Pinty et al., 2006). With the exception of the snow-covered season, the new albedo scheme, that relies on the simulated canopy structure, resulted in a substantial improvement of 0.05-0. 15 (Figs. 2 and 3). The European LPI-based likelihood that our model simulations come from the same populations as the MODIS albedo increased by a remarkable 11 and 23 % for, respectively, NIR and VIS albedo (from 61 and 69 % for the trunk to 72 and 92 % for the ORCHIDEE-CAN, Table 7). Given that the parametrization of the canopy radiation transfer model used in 5 ORCHIDEE-CAN relies on MODIS, the high likelihood may not come as a surprise. However, our implementation of the radiation transfer model also relies on the simulated absorbed light, simulated GPP, simulated allocation and simulated canopy structure (which depends on mortality and forest management). In the absence of all these processes our canopy radiation transfer model is expected to reproduce the MODIS data with a probability of 100 %. Hence, the likelihood of 72 and 92 % (for NIR and VIS, respectively) could also be interpreted as a verification of the aforementioned calculations; all calculations that determine the canopy structure reduce the reproducibility of the data by only 8-28 % (100 to 72 or 92 %). 15 The multilayer scheme is in the process of a detailed evaluation across a range of tests conditions , and further validation across a range of sites is on-going. The scheme is able to produce within-canopy temperature and humidity profiles, and successfully simulates the in-canopy radiation distribution, as well as the separation of the canopy from the soil surface. However, in order to preserve a measure 20 of continuity with previous evaluations of the model, the multilayer solution is here set to single layer operation mode, which includes the effects of hydraulic limitation (see Sect. 3.2) and canopy structure (see Sect. 3.3) on the energy budget. The single-layer set-up of the multi-layer solution makes use of an improved albedo estimation and is therefore expected to better simulate the net radiation that needs to 25 be redistributed in the canopy. This has been confirmed at a single site with a sparse canopy . Furthermore, the improvements in actual evapotranspira-

Tables Figures
Back Close

Full Screen / Esc
Printer-friendly Version

Interactive Discussion
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | tion in addition to the low RMSE (Fig. 2) are expected to be propagated in the performance of the energy budget.

Forest management strategies
Model comparison has previously demonstrated that explicitly treating thinning processes is essential to reproduce local and large scale biomass observations (Wolf 5 et al., 2011). This finding justifies the implementation of generic approaches to forest management despite the difficulties associated with defining and quantifying forest management and its intensity (Schall and Ammer, 2013). Although the use of so-called naturalness indices, in which the current state of the forest in referenced against the potential state of the forest, has been criticised because of difficulties in defining the 10 potential state of the forest (Schall and Ammer, 2013), such approaches were demonstrated to correctly rank different management strategies according to their intensity (Luyssaert et al., 2011). Naturalness indices making use of only diameter and stand density or the so-called Relative Density Index (RDI) have been previously implemented at the stand-level 15 (Fortin et al., 2012) and as well as in large scale models (Bellassen et al., 2010). This approach was shown to successfully reproduce the biomass changes during the life cycle of a forest (Bellassen et al., 2011;Fortin et al., 2012). The implementation of a forestry model based on the relative density index was reported to perform better than simple statistical models for stand-level variables such as stand density, basal 20 area, standing volume and height (Bellassen et al., 2011). Although the performance of the model was reported as less satisfying for tree-level variables, the approach is nevertheless considered reliable to model the effects of forest management on biomass stocks of forests across a range of scales from plot to country (Bellassen et al., 2011).
The forestry model implemented in ORCHIDEE-CAN is based on the RDI approach Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | served basal area from national forest inventories (de Rigo et al., 2014) and height from remote-sensing (Simard et al., 2011). With an RMSE of 3-7 for height and 7-15 for BA, and a chance of, respectively, 68 and 72 % to reproduce the data at the European scale (Table 7), our model is capable of correctly simulating the mean height and basal area but fails to capture much of the spatial variability (temporal variability 5 was not considered because the data products were only available for one time period). This finding could be due to the simulation protocol that started in 1850 with 2 to 3 m tall trees all over Europe. A longer simulation accounting for the major historical changes in forest management such as the reforestation in the 1700 s following an all time low in the European forest cover, the start of high stand management at the expense of coppicing in the early 1800s, and the reforestation programs following World War II (Farrell et al., 2000) is expected to improve the spatial variability in tree height and basal area. Regional deviations such as those observed in the Iberian Peninsula or over the entire Mediterranean (thus including part of the Iberian Peninsula) may be due to the lack of shrubs in the land cover map and parametrization of the ORCHIDEE- 15 CAN branch. Therefore the models simulates a higher stand density and higher basal area for regions where in reality shrubs occur. The parametrization of the forestry module strongly depends on the national forest inventories from Spain, France, Germany and Sweden. Therefore verification against the same data contains little information about the model quality. Nevertheless, no time-20 dependent relationships were used in the ORCHIDEE-CAN branch thus the model's capacity to reproduce the relationship between basal area and stand age, diameter and stand age or wood volume and stand age could be considered as largely independent test of the model quality. These tests were performed over 8 bioclimatic regions of France and the ORCHIDEE-CAN branch was found to largely capture the time de-Introduction

Conclusions
ORCHIDEE-CAN (SVN r2290) differs from the trunk version of ORCHIDEE (SVN r2243) by the allometric-based allocation of carbon to leaf, root, wood, fruit and reserve pools; the transmittance, absorbance and reflectance of radiation within the canopy; and the vertical discretisation of the energy budget calculations. Conceptual changes 5 towards a better process representation were made for the interaction of radiation with snow, the hydraulic architecture of plants, the representation of forest management and a numerical solution for the photosynthesis formalism of Farquhar, von Caemmerer and Berry. Furthermore, these changes were extensively linked throughout the code to improve the consistency of the model. By making use of observation-based parameters 10 the physiological realism of the model was improved and significant reparametrization was done by introducing twelve new parameter sets that represent specific tree species or genera rather than a group of phylogenetically often unrelated species, as is the case in widely used plant functional types (PFT). As PFTs have no meaning outside the modelling community, the species level parametrization of the ORCHIDEE-
The scheme overlooks the effect of vegetation shading bare soil for sparse canopies and gives the ground in all PFTs the same reflectance properties as bare soil.
Soil hydrology Vertical water flow in the soil is based on the Fokker-Planck equation that resolves water diffusion in non-saturated conditions from the Richards equation (Richards, 1931). The 4 m soil column consists of eleven moisture layers with an exponentially increasing depth (D'Orgeval et al., 2008).

No change
Soil temperature The soil temperature is computed according to the Fourier equation using a finite difference implicit scheme with seven numerical nodes unevenly distributed between 0 and 5.5 m (Hourdin, 1992).

No change
Energy budget The coupled energy balance scheme, and its exchange with the atmosphere, is based on that of Dufresne and Ghattas (2001). The surface is described as a single layer that includes both the soil surface and any vegetation.
A big leaf approach does not account for within canopy transport of carbon, water and energy. Further, it is inconsistent with the current multilayer photosynthesis approach and the new multilayer albedo approach.
Photosynthesis C3 and C4 photosynthesis is calculated following Farquhar et al. (1980) and Collatz et al. (1992), respectively. Photosynthesis assigns artificial LAI levels to calculate the carbon assimilation of the canopy. These levels allow for a saturation of photosynthesis with LAI, but have no physical meaning.
The scheme uses a simple Beer's law transmission of light to each level, which is inconsistent with the albedo scheme.

Autotrophic respiration
Autotrophic respiration distinguishes maintenance and growth respiration. Maintenance respiration occurs in living plant compartments and is a function of temperature, biomass and, the prescribed carbon/nitrogen ratio of each tissue (Ruimy et al., 1996). A prescribed fraction of 28 % of the photosynthates allocated to growth is used in growth respiration (McCree, 1974). The remaining assimilates are distributed among the various plant organs using an allocation scheme based on resource limitations (see allocation).   Table 3 coeff_maint_init Fraction of allocatable photosynthates that is consumed for maintenance and growth respiration Ryan (1991) k ls ,k lsmin , k lsmax m