E Dufrêne

H Davi

C Franc ¸ois

G Le Maire

V Le Dantec

A Granier

Ecological Modelling xxx (2005) xxx-xxx Modelling carbon and water cycles in a beech forest Part I: Model description and uncertainty analysis on modelled NEE

Keywords: Model, Carbon, Water, Growth, Eddy fluxes, NEE, NPP, Soil, Beech forest, Uncertainty analysis, Sensitivity analysis

A forest ecosystem model (CASTANEA) is developed with the aim to bridge the gap between soil-vegetation-atmosphere (SVAT) and growth models. A physiologically multi-layer process-based model is built, completed with a carbon allocation model and coupled with a soil model. CASTANEA describes canopy photosynthesis and transpiration, maintenance and growth respiration, seasonal development, partitioning of assimilates to leaves, stems, branches, coarse and fine roots, evapotranspiration, soil heterotrophic respiration, water and carbon balances of the soil. Net primary productivity (NPP) is calculated as the difference between gross photosynthesis and plant respiration. The net ecosystem exchange (NEE) between soil-plant system and atmosphere is calculated as the difference between gross photosynthesis and total respiration (soil + plants). The meteorological driving variables are global radiation, rainfall, wind speed, air humidity and temperature (either half-hourly or daily values). A complete description of the model parameterization is given for an eddy flux station in a beech stand (Hesse, France). A parametric sensitivity analysis is carried out to get a classification of the model parameters according to their effect on the NEE. To determine the key input parameters, a +10% or -10% bias is applied on each of the 150 parameters in order to estimate the effect on simulated NEE. Finally 17 parameters, linked to photosynthesis, vegetative respiration and soil water balance, appear to have a significant effect (more than 2.5%) on the NEE. An uncertainty analysis is then presented to evaluate the error on the annual and daily NEE outputs caused by uncertainties in these input parameters. Uncertainties on these parameters are estimated using data collected in situ. These uncertainties are used to create a set of 17,000 simulations, where the values of the 17 key parameters are randomly selected using gaussian random distributions. A mean uncertainty of 30% on the annual NEE is obtained. This uncertainty on the simulated daily NEE does not totally explain the discrepancies with the daily NEE measured by the eddy covariance technique (EC). Errors on daily measurements by EC technique and uncertainty on the modelling of several processes may partly explain the discrepancy between simulations and measurements.

Introduction

Long-term accurate estimates of (i) carbon and water fluxes between forest ecosystems and the atmosphere and (ii) carbon storage in trees and soil, are crucial to assess the role of forested areas in the global carbon cycle and in the continental water balance.

During the last century, several types of forest models have been developed with different objectives. Among the oldest, the empirical "forestry" models [START_REF] Schober | Ertragstafeln wichtiger Baumarten bei vershiedener Durchforstung[END_REF][START_REF] Dhôte | Modèles de la dynamique des peuplements forestiers: articulation entre les niveaux de l'arbre et du peuplement. Application à la sylviculture des hêtraies[END_REF][START_REF] Dhôte | Modélisation de la croissance des peuplements réguliers de hêtre: dynamique des hiérarchies sociales et facteurs de production[END_REF], with a time step of one year or more, predict stem growth without considering seasonal and inter-annual climatic changes. These models use empirical rules based on large data sets from field plots. They are able to reproduce tree growth over a century, according to forest management, age of stocking assuming no change in climate trend. However, climate change probably explains the increasing trend in radial growth [START_REF] Becker | Short-term variations and long term changes in oak productivity in northeastern France. The role of climate and atmospheric CO 2[END_REF] and in tree height [START_REF] Dhôte | Changements de productivité dans quatre forêts de chêne sessiles depuis 1930: une approche au niveau du peuplement[END_REF].

During the last two decades, several biophysical carbon/water fluxes models have been developed, considering the canopy as a single layer (i.e. "big-leaf"), a multi-layer or a three dimensional volume [START_REF] Wang | Influence of crown structural properties on PAR absorption, photosynthesis and transpiration in Sitka spruce: application of a model (MAESTRO)[END_REF][START_REF] Aber | A generalized, lumped-parameter model of photosynthesis, evapotranspiration and net primary production in temperate and boreal forest ecosystems[END_REF][START_REF] Amthor | Scaling CO 2 -photosynthesis relationships from the leaf to the canopy[END_REF][START_REF] Baldocchi | Scaling carbon dioxide and water vapour exchange from leaf to canopy in a deciduous forest. II. Model testing and application[END_REF][START_REF] De Pury | Simple scaling of photosynthesis from leaves to canopies without the errors of big-leaf models[END_REF]. They assess mass and energy exchange between the canopy and the atmosphere by coupling the fluxes of carbon dioxide and water vapour. The leaf physiologically based photosynthesis model by [START_REF] Farquhar | A biochemical model of photosynthetic CO 2 assimilation in leaves of C 3 species[END_REF], linked with a stomatal conductance model [START_REF] Jarvis | The interpretation of the variations in leaf water potential and stomatal conductance found in canopies in the field[END_REF][START_REF] Ball | A model predicting stomatal conductance and its contribution to the control of photosynthetis under different environmental conditions[END_REF][START_REF] Collatz | Physiological and environmental regulation of stomatal conductance, photosynthesis and transpiration: a model that includes a laminar boundary layer[END_REF] provides a theoretical framework for spatially integrating fluxes from leaf to canopy level. These process-based models integrate accurately the canopy functioning over time from seconds to days [START_REF] Caldwell | Canopy structure, light microclimate and leaf gas exchange of Quercus coccifera L. in a Portuguese macchia: measurements in different canopy layers and simulations with a canopy model[END_REF][START_REF] Baldocchi | A Lagrangian random walk model for simulating water vapour, CO 2 , and sensible heat flux densities and scalar profiles over and within a soybean canopy[END_REF][START_REF] Leuning | Leaf nitrogen, photosynthesis, conductance and transpiration: scaling from leaves to canopy[END_REF][START_REF] Williams | Modelling the soil-plantatmosphere continuum in a Quercus-Acer stand at Harvard Forest: the regulation of stomatal conductance by light, nitrogen and soil/plant hydraulic properties[END_REF][START_REF] Wang | A two-leaf model for canopy conductance, photosynthesis and partitioning of available energy I: model description and comparison with a multi-layered model[END_REF]. They are not designed, however, to predict the seasonal and interannual variations of tree growth and stand biomass increment. Forest growth models, based on ecological and biogeochemical principles, focus on how carbon and water fluxes vary from daily to annual and decadal time scales. They are able to predict changes in plant carbon pools (i.e. organs) by understanding there respective size and turnover time [START_REF] Mohren | Simulation of forest growth, applied to Douglas fir stands in the Netherlands[END_REF][START_REF] Mcmurtrie | Modelling the yield of Pinus radiata on a site limited by water and nitrogen[END_REF]. Some of these models consider litter and soil mineralization processes with the aim to predict soil organic matter dynamics [START_REF] Running | A general model of forest ecosystem processes for regional applications. I. Hydrologic balance, canopy gas exchange and primary production processes[END_REF][START_REF] Korol | Testing a mechanistic carbon balance model against observed tree growth[END_REF][START_REF] Running | Generalization of a forest ecosystem process model for other biomes, BIOME-BGC, and an application for global-scale models[END_REF][START_REF] Hoffmann | FAGUS, a model for growth and development of beech[END_REF][START_REF] Bossel | TREEDYN3 forest simulation model[END_REF]. However, they tend to ignore the effects of microclimate spatial variability within the plant canopy and, as they generally use a daily time step, they need empirical (i.e. not physiologically based) leaf and canopy photosynthesis sub-models. As hydrologic processes control drought effect on photosynthesis [START_REF] Schulze | Carbon dioxide and water vapour exchange in response to drought in the atmosphere and in the soil[END_REF] soil carbon and nitrogen dynamics [START_REF] Parton | Analysis of factors controlling soil organic matters levels in great plains grasslands[END_REF], some of these models also couple the carbon budget with a model simulating the water cycle. The rainfall reaching the ground is shared out into soil evaporation, transpiration, interception, infiltration or runoff. According to the application domain, the hydrology models are more or less sophisticated. For example, the evapotranspiration can be calculated as function of a potential evaporation (Kim et al., 1996) or estimated following [START_REF] Monteith | The state and movement of water in living organisms[END_REF] or [START_REF] Shuttleworth | Evaporation from sparse crops-an energy combination theory[END_REF]. The soil can be divided into numerous layers [START_REF] Braud | A simple soil-plant atmosphere transfer model (SiS-PAT) development and field verification[END_REF] or parameterized into one or several buckets [START_REF] Eagleson | Climate, soil and vegetation 3. A simplified model of soil moisture movment in the liquid phase[END_REF]. Nevertheless, detailed SVAT models are difficult to use for the investigation of the spatial and temporal variability of land surface fluxes. The large number of parameters they involve requires detailed field studies and experimentation to derive parameter estimates [START_REF] Boulet | A simple water and enrgy balance model designed for regionalization and remote sensing data utilization[END_REF]. Simple water balance models using simple soil and stand parameters and basic climatic data are often sufficient to predict temporal variation in soil water content [START_REF] Granier | A lumped water balance model to evaluate duration and intensity of drought constraints in forest stands[END_REF].

In this paper we focus on the complete description of a model simulating carbon, water and energy fluxes and we precisely describe the parameterization in the case of a beech forest stand (Hesse Euroflux site, France). A precise description of the CASTANEA structure and parameterization is given to account for the complexity and the large number of variables and parameters of this physiologically multi-layer process-based model. An effort concerning the parameterization has been done to have confidence in the model validation at eddy flux stations. Given the large number of input parameters, a sensitivity study, to determine the key parameters and their effect on the simulated output variables (uncertainty analysis), is presented. The uncertainty on the daily and annual simulated net ecosystem exchange (NEE) is then estimated using Monte-Carlo simulations.

Model description

General description

CASTANEA is a physiologically process-based model aiming at predicting the carbon and the water balances of an even-aged monospecific forest stand. Flux densities variables concerning CO 2 , and H 2 O are simulated and allow the calculation of the state variables, i.e. water and carbon contents.

Light interception and photosynthesis submodels are implemented in CASTANEA by using a multilayer canopy description. On the other hand, growth and allocation submodels, are treated using six compartments. Five compartments correspond to the main organs (see Fig. 1), while the last one, without explicit location, corresponds to the carbohydrate storage in the whole tree. One averaged tree is considered as representative of the whole stand (i.e the variability between trees is not taken into account). The soil water balance is simulated with a bucket model including three layers. The soil organic carbon (SOC) model is similar to [START_REF] Parton | Analysis of factors controlling soil organic matters levels in great plains grasslands[END_REF] approach except that two soil layers are considered (instead of one). Two time steps, half-hourly and daily, are used in the model. Most variables involving fluxes are simulated half-hourly; all state variables plus growth and phenology are daily simulated. Input meteorological driving variables can be either half-hourly or daily. If only daily values are available, the diurnal course (half-hourly values) of weather variables is generated according to the approach proposed by [START_REF] Mcmurtrie | Modelling the yield of Pinus radiata on a site limited by water and nitrogen[END_REF].

Radiation interception

The canopy is assumed to be horizontally homogeneous and vertically subdivided into a variable number of layers (i.e. multi-layer canopy model). Each layer contains the same amount of leaf area (typically less than 0.1 m 2 m -2 ). For each layer, the leaf inclination distribution function (LIDF) is represented through an ellipsoidal function [START_REF] Campbell | Extinction coefficients for radiation in plant canopies calculated using an ellipsoidal inclination angle distribution[END_REF], discretized in eighteen steps between 0 • and 90 • and the same average leaf angle is assumed for all canopy layers.

Three different radiative balances are performed, in the PAR (400-700 nm), in the NIR (700-2500 nm), and in the thermal infrared. In the PAR and global regions, incident light is split into direct and sky diffuse radiation using equations given by [START_REF] Spitters | Separating the diffuse and direct component of global radiation and its implications for modelling canopy photosynthesis. Part II. Calculation of canopy photosynthesis[END_REF] and [START_REF] Spitters | Separating the diffuse and direct component of global radiation and its implications for modelling canopy photosynthesis[END_REF]. In the thermal infrared, the diffuse atmospheric radiation is computed from air temperature according to [START_REF] Idso | A set of equations for full spectrum and 8 to 14 m and 10.5 to 12.5 m thermal radiation from cloudless skies[END_REF].

Radiative balance in the thermal infrared

In the thermal infrared the radiative balance coefficients are based on the general formalism given in Franc ¸ois (2002) (hereafter referred to as FR02). The soil and leaf emissivity are assumed to be constant, equal to 0.97. The canopy and the surface soil temperature are calculated iteratively by closing the energy budget (see later). The interception coefficient is computed according to [START_REF] Goudriaan | Crop micrometeorology: a simulation study[END_REF]. The soil-vegetation radiative balance coefficients are written according to the third model (MOD3) presented in FR02, including multiple scattering parameterization (Eqs. ( 11), ( 14), ( 16) and ( 17) in FR02) and adapted for diffuse radiation. These coefficients are τ v , ρ v (resp. canopy simple transmittance and reflectance), τ t , ρ t , (resp. multiple scattering canopy transmittance and reflectance), and ω t and ε t (resp. multiple scattering vegetation emissivity and canopy emissivity). The details of the expression for the coefficients and the expression for the radiative transfer balance are given in Appendix A.

Radiative transfer in the PAR and NIR waveband

The leaf reflectance and transmittance in the PAR (resp. the NIR) are assumed to be constant inside the canopy (see Table 1). The diffuse and direct absorbed PAR (aPAR d and aPAR s ) and the absorbed NIR are computed for each layer. The radiation extinction and diffusion are based on the SAIL model [START_REF] Verhoef | Light scattering by leaf layers with application to canopy reflectance modeling: the SAIL model[END_REF][START_REF] Verhoef | Earth observation modelling based on layer scattering matrices[END_REF]. The SAIL model is used to compute the five elementary reflectances (ρ) and transmittances (τ) for an elementary layer: ρ dd , ρ sd , τ dd , τ sd , τ ss (where the subscript s stands for direct radiation, and d for the diffuse radiation; first position is for incident radiation and second for transmitted or scaterred radiation). The calculation of the multi-layer radiative balance is described in Appendix B. The soil is treated as a particular layer, with no transmittance and having a reflectance based on the surface wetness Finally, a classical radiative balance equation is written combining the radiative balance in the PAR, NIR and thermal infrared region, to calculate the net radiation Rn veg for the vegetation and Rn soil for the soil. The leaf (T leaf ) and soil temperatures (T soil ), are then computed by solving iteratively the equations of the energy balance.

Canopy clumping is taken into account in the model by using a clumping factor (agreg) in the radiative transfer model. For each leaf layer this factor reduces the leaf area used in the SAIL sub-model for radiation intercepting. The way to estimate the clumping factor is specified in §3-in situ measurements.

Leaf photosynthesis and stomatal conductance

The leaf photosynthesis (A) is represented using a C 3 plants biochemical process-based model, according to [START_REF] Farquhar | A biochemical model of photosynthetic CO 2 assimilation in leaves of C 3 species[END_REF]. The CO 2 demand is determined as the minimum between Rubisco carboxylation and RuBP regeneration, while CO 2 supply depends on the difference in CO 2 concentration between the air outside the leaf and the carboxylation sites. Base temperature for forcing budburst (Granier et al., 2000).

Three main Eqs. (1), ( 2), (4) determine the assimilation rate (A):

First, the carbon dioxide demand Eq. (1) writes:

A = Vc -Rd (1)
where Vc is the carboxylation rate and Rd the dark respiration during the day. Rd is calculated assuming a light inhibition phenomenon (see respiration module). Second, the carbon dioxide supply Eq. ( 2) is similar to the Ohm's law in electricity where stomata conductance for carbon dioxide (g s CO 2 ) acts like an inverse of a resistance, while the gradient of CO 2 between the evaporative site (C i ) and the leaf surface (C s ) acts like a potential difference.

A = g s CO 2 × (C s -C i ) (2)
C s depends on atmospheric CO 2 (Ca), on A and on the leaf boundary layer conductance (g b CO 2 )

A = g b CO 2 (C a -C s ) ⇔ C s = C a -A g b CO 2 (3)
Third, the carbon dioxide control Eq. ( 4), which allows the stomatal conductance calculation according to [START_REF] Ball | A model predicting stomatal conductance and its contribution to the control of photosynthetis under different environmental conditions[END_REF]:

g s H 2 O = g 0 + g 1 × A × RH C s (4)
where g 1 depends on soil water stress (see Eq. ( 55)), g 0 corresponds to cuticular conductance of the leaf and RH is the relative humidity in the surrounding air.

Stomatal conductance for water g s H 2 O is proportionnal to g s CO 2 according to the molecular properties of water vapour and carbon dioxide.

The equations system (1), ( 2) and ( 4) is solved using an analytical solution proposed by [START_REF] Baldocchi | An analytical solution for coupled leaf photosynthesis and stomatal conductance models[END_REF] to calculate A, g s and C i .

The carboxylation rate (V c ) is limited either by RU-BISCO activity (W c ) at low C i or by electron transport rate (W j ) at higher C i . Note that we have neglected the gradual transition from electron transport to Rubisco limitation [START_REF] Kirschbaum | Temperature dependence of whole-leaf photosynthesis in Eucalyptus pauciflora Sieb. ex Spreng[END_REF][START_REF] Collatz | Physiological and environmental regulation of stomatal conductance, photosynthesis and transpiration: a model that includes a laminar boundary layer[END_REF] which have very small effect, as only a small fraction of leaves are near the transition to light saturation (De Pury and Farquhar, 1997).

Effect of temperature on photosynthesis

The biochemical basis of the temperature dependence of photosynthesis has been included in the model following [START_REF] Nolan | Multi-temperature effects on Hill reaction, activity of barley chloroplasts[END_REF] for electron transport rate (V jmax ), [START_REF] Long | Modification of the response of photosynthetic productivity to rising temperature by atmospheric CO 2 concentration: has its importance been underestimated?[END_REF] for the carbon dioxide (C i ) and oxygen (O i ) concentrations at the evaporative sites, and [START_REF] Bernacchi | Improved temperature response functions for models of Rubisco-limited photosynthesis[END_REF] for all others.

Leaf nitrogen effect on photosynthesis

Leaf nitrogen effect on photosynthesis is taken into account assuming a linear relationship between the maximal carboxylation rate (V cmax ) and leaf nitrogen content per unit area N a (constant along the leafy season, see Table 1), and a fixed ratio β between V cmax and the potential rate of electron flow (V jmax ) (Wullschleger, 1993).

V cmax = α N a N a (5) V jmax = V cmax × β (6)

Canopy photosynthesis

The leaf mass per area profile (LMA L ) is simulated using the following equation:

LMA L = LMA sun exp(-k LMA L) (7)
where LMA sun is the leaf mass per area for sun leaves at the top of the canopy, k LMA is an exponential decrease coefficient, and L the leaf area index above the considered level inside the canopy. LMA sun varies during leaf growth to reach a maximal value LMA Sunmax (see Eq. ( 21)).

Leaf nitrogen concentration (per mass) N mleaves (g N g dm -1 ) is assumed to be constant inside the canopy despite the fact that a slight increase is generally measured from the top to the bottom. Leaf nitrogen per unit leaf area N a (g N m -2 leaf ) is then calculated for each layer:

N a L = LMA L × N mleaves (8)
Photosynthesis is calculated layer per layer (typically less than 0.1 m 2 m -2 ), on sunlit (leaf area intercepting diffuse and direct PAR) and shaded (leaf area intercepting diffuse PAR only) foliage separately. Both contributions are summed and photosynthesis is finally summed for all layers over the canopy.

Phenology

All phenological events and growth are calculated on a daily time step. The phenological events (budburst, full leaf area development, full leaf maturity, start of leaf yellowing, complete canopy yellowing) depend on day-degrees and day duration. Leaf fall is assumed to occur simultaneously with senescence (i.e. yellowing). The phenology model has been calibrated in an additional site (Fontainebleau forest, near Paris), independent of the validation data set (see §3-in situ measurements).

Budburst

Budburst is simulated following the Eq. ( 9)-( 11):

R frcBB = T if T > T 2 and N > N start1 0 ifT ≤ T 2 or N < N start1 (9)
where R frcBB is the rate of forcing for bud break, T the mean daily temperature, T 2 the base temperature, N the day of year and N START1 the date of onset of rest.

S frcBB = N N START1 R frcBB if S frcBB < F critBB (10) N BB = N if S frcBB = F critBB (11)
with S frcBB the state of forcing, F critBB the critical value of state of forcing for the transition from quiescence to the active period and N BB the day when bud break occurred.

Leaf area development

Leaf area development is then computed following the Eqs. ( 12)-( 14):

S frcLA = N LAmax N BB T if T > T 3 (12) dLA dt = (LA max × T )/S frcLA if LA < LA max 0 i f L A = LA max (13) LN = L max LA max (14) L = LA × LN (15)
where S frcLA is the state of forcing for leaf growth, T 3 the base temperature for leaf growth, LA the mean area per leaf for the day considered, LA max the maximum area per leaf, L the leaf area index, L max the maximum leaf area index, N LAmax the day when L reaches L max and LN the number of leaves in the canopy. Note that LA max and LN are input parameters, consequently the maximum LAI (L max ) is forced according to the year and site.

Leaf fall (i.e. yellowing)

Finally, in the same way the yellowing and the leaf fall are simulated (Eq. ( 16)-( 19)):

R frcfall = T 4 -T if T > T 4 and N > N start2 0 i f T ≥ T 4 or N < N start2 ( 16
)
with R frcLfall the rate of forcing for leaf fall, T 4 the base temperature and N Start2 the date of onset of leaf ageing.

S frcLfall = N N start2 R frcLfall if S frcLfall < F critLfall (17) N Lfall = N if S frcLfall = F critLfall (18)
with S frcLfall the state of forcing, F critLfall the critical value of state of forcing for the transition from full leaf development to the leaf fall period and N Lfall the day when leaf fall started.

dL dt =    L (N-1) {[(Dl N -Dl min )/ (Dl NLfall -Dl min )] Pfall -1} if N > N Lfall 0 i f L = 0    (19)
with L (N -1) the leaf area index of the previous day, Dl N the day duration (min), Dl min the shortest day duration of the year and Dl NLfall the duration of the day when leaf started to fall. P fall is an empirical coefficient.

Allocation

Allocation of assimilates (daily time step) is based on constant coefficients coupled with a system of priorities varying along the year according to the phenological stage. The use of assimilates for maintenance respiration has priority over growth and storage allocation all over the year.

For temperate deciduous forest, the year can be divided into three periods:

(1) from budburst to maximum leaf maturity: sinks are maintenance respiration and leaves formation.

The carbon previously stored is used as a source in addition to gross photosynthesis. The other organs can only grow if gross photosynthesis exceeds the sum of maintenance respiration and leaf growing sink.

(2) from leaf maturity to complete yellowing (end of senescence): If net photosynthesis (i.e. gross photosynthesis minus autotrophic respiration) is positive, the allocation to all parts of the tree is allowed according to constant coefficients for each organ. Otherwise, no growth is allowed and storage is used for maintenance respiration. (3) from the end of leaf senescence to budburst: source is the storage compartment, sink is the maintenance respiration and no growth occurs.

Vegetative growth of forests may be temperature limited in temperate environment [START_REF] Cannell | Light use efficiency and woody biomass production of poplar and willow[END_REF]. However, except for leaves (see above), no direct effect of temperature on growth has been incorporated in the model. The water stress effect operates only indirectly by reducing gross photosynthesis. Moreover, there are no age-related effects on carbon allocation besides the indirect effect of increasing living biomass.

The partitioning model is somehow simplistic, since we assume constant partionning coefficients for aerial wood, fine roots, coarse roots and reserves (following [START_REF] Mcmurtrie | Above and below ground growth of forest stands: a carbon budget model[END_REF]. It corresponds to the empirical approach described in [START_REF] Thornley | Parttionning during vegetative growth[END_REF]. But it requires very few parameters compared with more mechanistic approaches as the transportresistance mechanistic models. Moreover, it seems quite efficient, for a given site and one year (see Part II). However, if one want to take into account the age [START_REF] Magnani | Age-related decline in stand productivity: the role of structural acclimation under hydraulic constraints[END_REF] and fertility [START_REF] Mcmurtrie | Forest productivity in relation to carbon partitioning and nutrient cycling: a mathematical model attributes of trees as crop plants[END_REF] effects on the allocation processes, the allocation coefficients may be calculated by adding empirical or functional constraints on the sinks (e.g. fine roots to leaf ratios) as in [START_REF] Barbaroux | Analyse et modélisation des flux de carbone de peuplements forestiers pour la compréhension de la croissance des espèces feuillues Quercus petraea et Fagus sylvatica[END_REF] and [START_REF] Davi | Développement d'un modèle générique simulant les flux et les stocks de carbone et d'eau dans le cadre des changements climatiques[END_REF].

Sunlit leaf mass per area (LMA) evolution

The evolution of sunlit leaf mass per area (LMA Sun ) is simulated in the same way as the leaf area index (L), with a dependency on air temperature. A higher critical value of state of forcing F crit LA is used to calculate the end of LMA growth.

S frcLMA = N LMAmax N BB T if T > T 3 (20) LMA Sun = S frcLMA F critLMA × LMA Sunmax ( 21
)
S frcLMA is the state of forcing for leaf mass growth, LMA Sunmax the maximum value of sunlit leaf mass per area reached during the year and N LMAmax the day when LMA Sun reaches LMA Sunmax . Note that LMA Sunmax is an input parameter and then it varies according to the site and the species.

The foliar biomass is then calculated as the integral over the canopy of leaf mass per area times the leaf area index.

B leaf = L max 0 LMA L dL (22)
Then the foliar growth (GB leaf ) (which, if negative, turns out to be the mortality MB leaf ), is calculated as the difference between foliar biomass (B leaf ) of the day and foliar biomass of the previous day.

Carbon allocation to the different compartments

After accounting for leaf growth and need of carbohydrate for maintenance respiration, the available carbon is allocated to the four other compartments (aerial woody organs, storage, coarse roots and fine roots). The carbon growth of the different organs but leaves (GB organ ), is a proportion (allocation coefficients) of the available carbon during the same day:

GB organ = AG organ CR organ ×(A canopy -RM-RG leaf -GB leaf ) (23)
AG organ is the allocation coefficient per organ type (the sum of the four coefficients is one), A canopy the gross canopy photosynthesis, CR organ the organs construction cost (see growth respiration), RM the total maintenance respiration (all organs), RG leaf the leaf growth respiration and GB leaf the leaf growth.

During the leafy period, if maintenance respiration (RM) is greater than gross photosynthesis, the required carbon is taken from storage (MB storage , "storage emptying"). During the leafless period, all the carbon needed for maintenance respiration contribute to storage emptying.

dB storage dt = GB storage -MB storage (24)
The fine roots mortality depends on a constant turn over (TO fine roots ) and on the fine roots biomass (B fine roots ). The fine roots biomass is thus given by: dB fine roots dt = GB fine roots -TO fine roots × B fine roots (25)

The mortality of woody organs (branches, trunk and coarse roots) is not simulated. Hearthwood formation (i.e. duramenisation) is not taken in consideration seasonally but only annually. The equations of daily biomass evolution are therefore:

dB organ dt = GB organ (26)
The distribution of aerial woody biomass between stems and branches depends on an allometric coeffi-cient P branch which is function of stand age and species.

B branch = P branch × B aerial wood (27)
The living biomass of woody organs (B alive organ , i.e. trunk, branches and coarse roots) is calculated each day after growth calculation using the proportion of living wood (P alive organ ). This living biomass is used to estimate maintenance respiration and changes with species and age.

B alive organ = P alive organ × B organ (28)

Maintenance respiration

Maintenance respiration is calculated from the nitrogen content of living biomass and assuming an exponential relationship to account for temperature dependence [START_REF] Ryan | Effects of climate change on plant respiration[END_REF]. For each compartment (i.e. type of organ) maintenance respiration is simulated every half hour.

RM organ = B alive organ × MRN × N morgan ×Q (Tsuf-T MR /10) 10 organ (29)
where T surf is the surface temperature, B alive organ the biomass of living tissue actually respiring, N morgan the amount of nitrogen per alive mass unit and MRN the respiration rate per nitrogen unit. Q 10 of aerial woody parts (i.e. trunk and branch) is different from other compartments.

During the night, the foliar respiration is calculated using Eq. ( 29). During the day, light is assumed to inhibit this respiration according to a constant ratio P leaf inhib . The degree of inhibition ranges between 17 and 66% depending on species [START_REF] Sharp | Kok effect and the quantum yield of photosynthesis[END_REF][START_REF] Brooks | Effect of temperature on the CO 2 /O 2 specificity of ribulose-1,5-bisphosphate carboxylase\oxygénase and the rate of respiration in the light. Estimates from gas-exchange measurements on spinach[END_REF][START_REF] Kirschbaum | Investigation of the CO 2 dependance of quantum yield and respiration in Eucalyptus pauciflora[END_REF]. [START_REF] Villar | Dark leaf respiration in light and darkness of an evergreen and a deciduous plant species[END_REF] give a mean rate of 51% for evergreen tree species and 62% for deciduous tree species.

Growth respiration

Organ growth is calculated on a daily time step and, consequently, growth respiration is first calculated on a daily basis (RG d organ ), assuming that daily growth respiration depends on both growth rate (GB organ ) and construction cost (CR organ ) of the considered organ. The construction cost can be calculated from gas exchange measurements or following the approach of Penning de Vries et al. (1974) and Penning de Vries (1975aVries ( , 1975b)), based on the knowledge of the organ biochemical composition (i.e. nitrogen, lignin, lipids, organic acids, carbohydrates, total mineral content).

RG d organ = GB organ (CR organ -1) (30)

During a second step, growth respiration is calculated half-hourly assuming that surface temperature is the driving variable (Q 10 -based relationship). We can then calculate the net ecosystem CO 2 flux (NEE) on an half-hourly time step.

Water fluxes 2.9.1. Rainfall interception by the canopy

The canopy is regarded as having a surface storage capacity, recharged by rainfall and discharged by evaporation and drainage. The storage capacity of leaves and woody parts (i.e. stems and branches) are considered separately. During the leafy period, leaves intercept a fraction of rainfall, depending on both leaf area index and gap fraction. If the amount of water intercepted exceeds the leaf capacity storage, the excess is lost by throughfall and can either be intercepted by woody parts or reach the soil surface. If the amount of rainfall intercepted by woody parts exceeds the storage capacity of the bark then water can both run out along branches and stems or throughfall on to the litter layer. Evaporation rate occurs according to the Penman-Montheith equation [START_REF] Monteith | The state and movement of water in living organisms[END_REF] assuming a zero stomatal resistance.

R leaf max = L × RA leaf (31) R bark max = 2 × WAI × RA bark (32) R can max = R bark max + R leaf max (33)
R leaf max , R bark max , and R can max are respectively the bark, leaf and canopy water storage capacity. L and WAI are respectively leaf and wood area indices. RA leaf and RA bark are respectively the water storage capacity per unit of leaf and bark area.

The gross water interception of leaves (I leaf ) and barks (I bark ) depends on precipitation (Pi) and interception coefficients for leaves (C leaf ) and bark (C bark ) respectively: 34)

I leaf = C leaf × Pi (
I bark = C bark × Eg leaf ( 35
)
during leafy period (if L > WAI). Eg leaf is the drainage occurring, when the water on the leaf exceeds the leaf storage capacity.

I bark = C bark × Pi (36)
during leafless period (if L < WAI)

C leaf = L L × C ia + C ib ( 37
)
C bark = WAI WAI × C ia + C ib (38)
C ia and C ib are empirical coefficients. R leaf , R bark and R can are respectively the current bark, leaf and canopy (bark + leaves) water reserve, which are calculated half hourly by summing gross water interception (respectively I leaf and barks I bark ).

If the amount of water on bark (R bark ) exceeds the storage capacity (R bark max ) then drainage occurs and the excess of water is separated into two parts: the first one flowing along the bark (Ec) and the other one dripping over the soil (Eg bark ). Allocation between flows is determined with an empirical coefficient (PROP EC ).

Then the net rainfall interception of the canopy (IN), the rain falling directly on the soil through the canopy (PSd soil ) and the total amount of precipitation reaching the soil (PS soil ) are calculated:

IN = (I leaf -Eg leaf ) + (I bark -Eg bark -Ec) (39) PS soil = Pi -IN (40) PSd soil = PS soil -Eg leaf × (1 -C bark ) -Eg bark -Ec (41)
2.9.2. Canopy evapotranspiration Following [START_REF] Rutter | A predictive model of rainfall interception in forests. I. Derivation of the model from observations in a plantation of Corsican pine[END_REF] and Chassagneux andChoisnel (1986, 1987), the canopy evapotranspiration (ETR canopy ) is calculated by adding water evaporated (EP) from the wet parts of the canopy and transpiration (Tr) from the dry parts:

ETR canopy = R can R canmax × EP + 1 - R can R canmax × Tr (42)
R can represents the current quantity of water retained by the canopy (i.e. bark + leaves) and R canmax the storage capacity. The Penman-Monteith equation [START_REF] Monteith | The state and movement of water in living organisms[END_REF]) is applied at the canopy level to calculate both transpiration (Tr) and evaporation (EP).

Tr = h × Rn veg + (ρ cp × C ph × VPD)/R ac λh( h + γ × (1 + 1/(g c × R ac ))) (43)
Heat storage and soil heat flux are assumed to be negligible. The canopy conductance g c is calculated by averaging leaf stomatal conductance over the canopy. For water evaporation calculation, an infinite conductance g c is assumed (no resistance) and Eq. ( 43) is simplified to the Penman formulation [START_REF] Penman | Natural evaporation from open water, bare soil and grass[END_REF] modified by [START_REF] Van Bavel | Potential evaporation: the combination concept and its experimental verification[END_REF]. Net radiation (Rn veg ) is driven by energy balance (see radiative balance in the global domain), air vapour pressure deficit (VPD) is an input meteorological variable and the canopy aerodynamic resistance (R ac ) is calculated assuming a logarithmic wind profile above the canopy and an exponential decrease inside the canopy (adapted from [START_REF] Lafleur | Application of an energy combination model for evaporation from sparse canopies[END_REF]. λh is latent heat; h is the slope of the relation between saturate vapor pressure and surface temperature, γ is the psychometric constant; ρ cp is the air density and C ph the air specific heat. All of these parameters vary according to the temperature except C ph .

Soil evaporation

Soil evaporation (EP soil ) is calculated in the same way as canopy transpiration using Penman-Monteith equation. It differs from canopy evaporation only by the available energy for evaporation (Rn soil ), the aerodynamic resistance from soil to atmosphere (R as ) and the water vapour soil conductance (g soil )

Soil conductance to water vapour depends on water status both in the upper part of soil (RW top ) and in litter (RW lit ).

g soil = (g soil max -g soil min ) × RW top -RW top wilt RW top fc -RW top wilt + g soil min (44)
g soil max is the maximal soil conductance which is a site specific constant. Note that when leaf litter above top soil is wet, g soil is taken to g soil max value. RW top fc is the water content at field capacity of upper part of soil. RW top wilt , the water content at wilting point upper part of soil. When litter and top-soil layers are dry, evaporation is limited according to g soil min .

Soil water balance

The soil water balance model (see Fig. 2) is basically a bucket one with three layers (litter, a top-soil layer and a total-soil root zone including top-soil layer). Current water contents (RW lit , RW top , RW soil ) are the state variables corresponding to the three compartments. For each layer, water content is calculated daily as the difference between inputs (stemflow, throughfall and drainage from above layer) and outputs (evaporation, transpiration and drainage). The way to arrange in time the different fluxes is different for litter and soil compartments: -Output 1: evaporation of litter (E lit ).

E lit = min(EP soil ; RW lit int ) (46)
-Output 2: after evaporation and interception, if RW lit exceeds the water litter content at "field capacity" (RW lit fc ) then water drains down (D lit ) to the topsoil layer:

D lit = max(RW lit int -E lit -RW lit fc ; 0) (47) RW lit(t) = RW lit(t-1) + PSd soil + Eg leaf × (1 -C bark ) + Eg bark -D lit -E lit ( 48
)
• Soil: 2 steps are considered Note that two water balances are calculated independently; one for the top-soil layer and another one for the total soil root zone (including top soil layer).

Both are used to calculate the soil water stress (see below Eq. ( 56)), which controls both transpiration and photosynthesis.

-Input: stemflow (Ec) and litter drainage (D lit ).

For each layer (j), a temporary water content is calculated:

RW int (j) = RW (d-1) (j) + Ec + D lit ( 49
)
For each soil layer, the water content at day d RW (d) is calculated as the product of actual soil volumetric humidity (θ) and the layer thickness H(l).

RW fc (j) = θ fc (j) × H(j) (50)

RW fc (j) is the water content at field capacity of the layer (j). RW wilt (j), the water content at wilting point, is calculated in the same way.

-Output: soil water drainage (D(j)), soil evaporation (E(j)) and canopy transpiration TR(j), for the layer (j):

D(j) = max(RW int (j) -Ec -RW fc (j); 0) (51) RW (d) (j) = RW (d-l) (j) + Ec + D lit -D(j) -E(j) -TR(j) (52) E(j) = EP soil -E lit ( 53
)
TR top = prac × TR ( 54
)
where prac is the proportion of fine roots in top soil layer, depending on soil type and soil depth.

Effect of soil water status on canopy gas exchange

During water stress period the slope (g 1 ) of the relationship, proposed by [START_REF] Ball | A model predicting stomatal conductance and its contribution to the control of photosynthetis under different environmental conditions[END_REF] between leaf assimilation (A) and stomatal conductance (gs) (Eq. ( 4)), is assumed to decrease linearly when soil water storage decreases [START_REF] Sala | Simulations of canopy net photosynthesis and transpiration in Quercus ilex L. under the influence of seasonal drought[END_REF]. The effect of soil water stress on photosynthesis translates into g 1 through a reduction factor (reduc):

g 1 = (g 1max -g 1min ) × reduc + g 1min (55)
g 1min and g 1max are respectively the minimum and maximum values taken by g 1 . While the water content of one of the soil layer (j) is above a specific threshold based on proportion of soil extractable water content (S stress ×[RW fc (j) -RW wilt (j)] + RW wilt (j)) then no water stress occurs (i.e. reduc = 1). On the contrary when both soil layers (top and total) are below this threshold, the soil water stress increase linearly with decreasing soil water content until zero: reduc = RW soil -RW soil wilt S Stress (RW soil fc -RW soil wilt ) (56) RW(j), RW wilt (j) and RW fc (j) are respectively the water content, the water content at wilting point and the water content at field capacity for the considered soil layer (j). S Stress is a threshold parameter for soil water stress.

Heterotrophic respiration and soil organic matter

The soil organic carbon sub-model (SOC), based on the soil organic matter (SOM) sub-model of CEN-TURY [START_REF] Parton | Analysis of factors controlling soil organic matters levels in great plains grasslands[END_REF], simulates the dynamic of carbon in the soil system at the daily time step. Soil organic carbon is divided into three major components including active, slow and passive soil carbon. Active SOC includes live soil microbes plus microbial products, the slow pool includes resistant plant material (lignin-derived material), soil-stabilized plant and microbial material. The passive material is very resistant to decomposition and includes physically and chemically stabilized SOM. The model also includes a surface microbial pool, which is associated with decomposing surface litter (mainly leaf litter). Carbon flows between these pools are controlled by decomposition rate and microbial respiration loss parameters, both of which are a function of soil texture, soil temperature and soil water content.

The original version developed by [START_REF] Parton | Analysis of factors controlling soil organic matters levels in great plains grasslands[END_REF] has been adapted to be coupled with the forest plant sub-model (Le Dantec, 2000;[START_REF] Epron | Seasonal dynamics of soil carbon efflux and simulated rhizosphere respiration in a beech forest[END_REF]. In this soil model, all fractions are located in the soil vertical profile (see Fig. 2). The soil is divided into three layers according to the soil water balance sub-model. A deep soil layer, however, is considered explicitly rather than a total soil layer (see Fig. 1). The first layer is the surface litter one, containing both metabolic and structural fractions (calculated from leaf litter input) and the associated microbial pool. The top (0-30 cm) and the deepest layer (below 30 cm) include all carbon pools, i.e. the microbial population (active pool), both metabolic and structural fractions of root litter, and the slow and passive soil organic carbon.

The carbon inputs from plant residue (leaves and fine roots) are simulated by the plant phenology submodel. Surface layer receives the leaf litter and soil layers receive the root litter in proportion of root biomass (i.e. assuming a similar mortality rate for each layer). No carbon migration is assumed between the superficial and the deep soil layers. The "slow pool" of the superficial layer, however, is supplied with the organic carbon coming from both the structural fraction and the active SOC of the surface layer.

For the different pools (i), the decomposition rate (DR) is the output flow:

DR (i) (j) = K (i) (j) × Q (i) (j) × A w (j) × A t (j) × C (i) (j) ( 57
)
where j is the soil layer where the considered carbon pool i is located.

According to [START_REF] Parton | Analysis of factors controlling soil organic matters levels in great plains grasslands[END_REF], K is the maximum decomposition rate, C the carbon biomass, A w and A t , the impact function of soil water content and temperature respectively. Q depends on soil textural properties of the considered j layer for the soil active and slow carbon pools; it depends on lignin content of the considered litter (leaves or fine roots) for the litter structural material, and it equals 1 for the other carbon pools.

A w (j) = 1 (1 + 30e -8.5×RW(j)/[RW fc (j)-RW wilt (j)] ) (58)

where RW(j) is the soil water content of the layer (j), simulated by the soil water balance model. RW wilt(j) and RW fc(j) are the soil water content at wilting point and at field capacity respectively for the considered layer.

A T (j) = 45 -T (j) 45 -T b (j) 0.2 × e 0.076(1-[(45-T (j))/(45-T b (j))] h(j) ) (59)

where T(j) is the temperature of the considered layer, T b (j) a base temperature depending on soil layer and h(j) a parameter depending on soil layer. The model assumes that all C decompositions flows (DR) are associated with microbial activity and that microbial respiration occurs for each of these flows. Heterotrophic soil respiration (R heterotrophic ) is the sum of all these microbial respiration processes. Depending on the considered soil C pool (see Fig. 2), output flows can be separated into two parts according to lignin content (flows from structural fraction of litter) or with textural composition of the soil layer (flows from both active and slow soil pools).

Calculation of CO 2 fluxes: NPP (Net Primary Productivity) and NEE (Net Ecosystem Exchange)

NPP = A canopy -RM -RG (60) NEE = NPP -R heterotrophic ( 61
)
A canopy is the canopy gross phototosynthesis, RM the maintenance respiration of all organs, RG the growth respiration of all organs and R heterotrophic the soil heterotrophic respiration.

3. In situ measurements (Hesse, France, 1997)

Site characteristics

The experimental plot area is located in the State forest of Hesse (East of France; 48 • 40 N, 7 • 05 E), and is mainly composed of beech (Fagus sylvatica L.); understorey vegetation is very sparse. The elevation is 300 m and the slope is less than 2%. The studied plot covers 0.63 ha of a pole beech forest (30 years old in 1997) with a density of 3482 trees ha -1 and a dominant height close to 14 m. Three towers were installed: one (18 m high) is used for eddy covariance and microclimate measurements, the two other ones (15 m high) for physiological measurements. A hut containing the data acquisition systems is located near the first tower. Sixty sub-plots are delimited for geo-statistical studies of soil and vegetation. For more details, see Granier et al. (2000a).

Model parameterization data set

Meteorological driving variables

Micro-meteorological sensors are installed above the stand at a height of 17.5 m: a global radiometer (Mod 180, Cimel, France), a laboratory-made quantum sensor [START_REF] Pontailler | A simple red:far-red sensor using gallium arsenide phosphide detectors[END_REF], a ventilated psychrometer (INRA, with Pt 100), a raingauge (ARG 100, Campbell Scientific, Logan, Utah), and a switching anemometer (Vector Instruments, UK). Additionally, soil temperature sensors (copper/constantan thermocouples) are installed at -10 cm (five replicates) and in one profile (-5, -10, -20, -40, and -80 cm) in a central plot. Weather data were averaged/summed and recorded half-hourly using an acquisition system (Campbell Scientific) (see Table 3).

Species and site-specific parameters

Values of the main input parameters are given in Table 1 and initial values of state variables are given in Table 2.

Stand and canopy structure

Initial biomass of organs have been determined by inventories of aerial (see Granier et al., 2000a) and below-ground tree parts (Le [START_REF] Goff | Root biomass and biomass increment in a beech (Fagus sylvatica L.) stand in North-East France[END_REF] or by an estimation from literature for fine roots less than 2 mm diameter [START_REF] Bauhus | Fine-root growth in beech (Fagus sylvatica) forest gaps[END_REF]. The amount of carbohydrate stored was estimated from measurement on the different woody organs [START_REF] Barbaroux | Distribution of aboveground and belowground carbohydrate reserves in adult trees of two contrasting broad-leaved species (Quercus petraea and Fagus sylvatica)[END_REF]. Maximum leaf area index (L max ) was estimated by litter collection during leaf fall (42 square litter traps, 0.25 m 2 each) (Granier et al., 2000a). Foliage aggregation was calculated by dividing an apparent leaf area index measured using a plant canopy analyzer (LAI2000, Li-Cor, Lincoln, Nebraska), calculated with four rings, by the true leaf area index derived from litter collection. Wood area index was calculated using allometric relationships between both trunk and branch area with trunk girth measured at 1.30 m height [START_REF] Damesin | Stem and branch respiration of beech: from tree measurements to estimations at the stand level[END_REF].

A profile of leaf mass per area was measured at end of June in 1997 (Montpied, personaal communication) allowing us to calculate the coefficient k LMA , which is close to those observed in other canopies [START_REF] Aussenac | Etude bioclimatologique d'une futaie de feuillus (Fagus sylvatica L. et Quercus sessiliflora Salisb.) de l'Est de la France. I. Analyse des profils microclimatiques et des caractéristiques anatomiques et morphologiques de l'appareil foliaire[END_REF][START_REF] Ducrey | Etude bioclimatologique d'une futaie de feuillus (Fagus sylvatica L. et Quercus sessiliflora Salisb.) de l'Est de la France. III. Potentialités photosynthétiques à différentes hauteurs dans le peuplement[END_REF][START_REF] Rambal | Optimization of carbon gain in canopies of mediterranean evergreen oaks[END_REF].

Rainfall interception and canopy evaporation rate

Bark water storage capacity (R barkmax ) was estimated from literature data on smooth bark species, by dividing canopy water storage (0.3-0.65 mm) during leafless period [START_REF] Helvey | Canopy and litter interception of rainfall by hardwoods of eastern United States[END_REF][START_REF] Leyton | Rainfall interception in forest and moorland[END_REF] by wood area derived from allometric relationship [START_REF] Damesin | Stem and branch respiration of beech: from tree measurements to estimations at the stand level[END_REF]. We used leaf water storage capacity (R leafmax ) given by [START_REF] Nizinski | Mesures et modélisation de l'interception nette dans une futaie de chênes[END_REF] for Quercus petraea. We obtain water storage during leafless (R barkmax ) period of 0.45 mm and canopy water storage during leafy period of 1.6 mm (very similar to the threshold value of 1.7 mm given by Granier et al., 2000b).

Interception coefficients for leaves (C leaf ) and woody parts (C bark ) were estimated using gap canopy fractions measured over 0-20 • zenithal angle from hemispherical photography. Then parameters C ia and C ib were calculated by solving Eqs. ( 37) and ( 38).

Allocation between stemflow and drainage (PROP EC ) was calibrated using data measured during the leafy period of 1997 (Granier et al., 2000b).

Maintenance respiration

Temperature dependency (Q 10 ) for bole and branches was measured in Hesse forest by [START_REF] Damesin | Stem and branch respiration of beech: from tree measurements to estimations at the stand level[END_REF]. We used the same value of Q 10 for coarse roots as for bole wood. For other organs, fines roots and leaves, we used Q 10 values from literature, respectively from [START_REF] Epron | Fine root respiration in forest trees[END_REF] and from [START_REF] Vose | Challenges to modelling NPP in diverse eastern deciduous forests: species-level comparisons of foliar respiration responses to temperature and nitrogen[END_REF].

Bole and branches nitrogen concentrations (N m ) were measured in situ [START_REF] Ceschia | Spatial and seasonal variations in stem respiration of beech trees (Fagus sylvatica)[END_REF]. Profiles of N m in leaves were measured by P. Montpied (personal communication) showing little variation with leaf layer. N m in coarse roots was taken the same as for bole wood. For fine roots we used nitrogen concentration from [START_REF] Van Praag | Root turnover in a beech and spruce stand of the Belgian Ardennes[END_REF].

We used the nitrogen dependency given by [START_REF] Ryan | Effects of climate change on plant respiration[END_REF] to calculate the basal respiration (i.e. respiration at base temperature, here 15 • C).

Phenology

To derive phenological parameters, we used observations and measurements performed on beech plots in Fontainebleau forest (France). For budburst we used a 17-year data set, for leaf growth and leaf fall we used 4 years LAI values, measured at regular intervals along the leafy periods (Li-Cor LAI 2000) and 9 years of semi-quantitative leaf yellowing observations. Despite that budburst models generally do not work very well for beech species [START_REF] Kramer | Selecting a model to predict the onset of growth of Fagus sylvatica[END_REF], the accuracy of the obtained prediction is better than the null hypothesis which is not the case of different models tested by [START_REF] Kramer | Selecting a model to predict the onset of growth of Fagus sylvatica[END_REF]. This is specific to Fagus sylvatica because our model performs very well for both sessile oak and hornbeam phenology measured using the same protocol in the Fontainebleau forest (data not shown). The data set used to parameterise leaf growth and leaf fall is small. However, as the leaf growth is very fast (less than 3 weeks for beeches), the influence of the function and the parameters used is not crucial. Leaf yellowing is more variable from year to year according to late seasonal water stress influence. Consequently the parameterisation could possibly be improved by adding more observations. Few data are available for LMA dynamic during leaf maturation phase on mature beech trees: measurements in the Italian site of Collelongo in 1995 [START_REF] Matteucci | Bilancio del carbonio in una faggeta dell'Italia Centro-Meridionale: determinanti ecofisiolo-gici, integrazione a livello di copertura e simulazione dell'impatto dei cambiamenti ambientali[END_REF], a study from [START_REF] Schulte | Saisonale und interannuelle Variabilität des CO 2 -Gaswechsels von Buchen (Fagus sylvatica L.). Bestimmung von C-Bilanzen mit Hilfe eines empirischen Modells. Aachen: Shaker Schulze. E. -D. 1970. Der CO 2 -Gaswechsel der Buche (Fagus silvatica L.) in Abhängigkeit von den Klimafaktoren im Freiland[END_REF] in the German site of Solling (1986)(1987)(1988), an unpublished experiment on mature trees during spring 1998 in Orsay (France) and a monitoring in Hesse in 1998 (Montpied, personal communication). We determined the critical value of state of forcing (from budburst to leaf maturity) by averaging results from all experiments (424 • C/day ± 81).

Allocation of assimilates

Leaf allocation is predicted by the phenology submodel. We used the coefficient for storage allocation estimated by [START_REF] Barbaroux | Analyse et modélisation des flux de carbone de peuplements forestiers pour la compréhension de la croissance des espèces feuillues Quercus petraea et Fagus sylvatica[END_REF]. The coefficient for fine roots was estimated by inversing the model and assuming a constant biomass on a yearly basis. Finally the coefficient for coarse roots was deduced assuming a constant ratio of 0.2 between coarse roots and trunks.

Mortality of organs

Leaf mortality is simulated by the phenology submodel. Fine roots mortality is assumed to be proportional to the biomass with a turnover of 1 over the year [START_REF] Bauhus | Fine-root growth in beech (Fagus sylvatica) forest gaps[END_REF]. We assumed no bole or branches mortality over one year simulation.

Growth respiration

For aerial woody organs, the construction cost is calculated using both respiration and growth measurements [START_REF] Damesin | Stem and branch respiration of beech: from tree measurements to estimations at the stand level[END_REF]. The construction cost for coarse roots is considered to be the same as for bole. For leaves and fine roots, data from literature were used [START_REF] Niinemets | Energy requirement for foliage formation is not constant along canopy light gradients in temperate deciduous trees[END_REF]Ågren and Axelsson, 1980 respectively).

Soil water balance

During 1997, the soil volumetric water content (θ) was measured weekly with a neutron probe (NEA, Ballerup, Denmark) along eight aluminium access tubes of 160 cm long. The volumetric minimum water content (θ wilt ) and the water content at field capacity (θ fc ) were estimated for both top and deep soil layers. For each soil layer, θ fc was taken as the average observed θ during the winter when soil was refilled (after a delay of one day to take into account of rapid drainage) and θ wilt was estimated from retention curves.

The height of total soil root zone (H soil ) is assumed to be 160 cm in Hesse (Granier et al., 2000b). The height of the upper soil layer (H top ) is set to 30 cm, corresponding to the part of the soil which is more dynamic in terms of volumetric water content evolution and contributes more to the stomata control [START_REF] Dufrêne | Photosynthèse. In: consommation en eau et modélisation de la production chez le palmier à huile (Elaeis guineensis Jacq[END_REF]. The upper soil layer includes 63% of fine roots [START_REF] Quentin | Etude des sols de la forêt de Hesse (Lorraine) Contribution à l'étude du bilan hydrique[END_REF]. This value is similar to those given for temperate deciduous species by [START_REF] Jackson | A global buget for fine root biomass, surface area, and nutrient contents[END_REF].

Very few data regarding soil water vapour conductance of forest are available in the literature. No measurements were done on Hesse forest and consequently the minimum soil vapor conductance (g solmin ) was taken equal to the value given by [START_REF] Schaap | Forest floor water content dynamics in a Douglas fir stand[END_REF] and the maximal soil vapour conductance (g solmax ) was derived from [START_REF] Kelliher | Estimating the effects of understorey removal from Douglas fir forest using a twolayer canopy evapotranspiration model[END_REF]. These two values were measured in a Douglas fir stand.

The soil water stress threshold (S Stress ) is not very variable according to soil and species (see review from [START_REF] Granier | A lumped water balance model to evaluate duration and intensity of drought constraints in forest stands[END_REF] and the value was fixed to 0.4.

Soil heterotrophic respiration and organic matter

The soil mineral composition and the soil total carbon content were measured in situ [START_REF] Farque | Contribution des racines de hêtre (Fagus sylvatica L.) au bilan carboné d'un sol forestier[END_REF]. The initialization of the different carbon pools was set from the vertical distributions of carbon concentrations and of residence times in soil [START_REF] Elzein | Mechanistics simulation of vertical distribution of carbon concentrations and residence times in soils[END_REF], and from the vertical distribution of microbial biomass in soil [START_REF] Wolters | Microbial carbon turnover in beech forest soil at different stages of acidification[END_REF][START_REF] Ross | Microbial biomass, and C and N mineralization, in litter and mineral soil of adjacent montane ecosystems in a southern beech (Nothofagus) forest and a tussock grassland[END_REF]. Equilibrium values for carbon pools obtained after 20-year simulation periods were used to initialize the sub-model [START_REF] Epron | Seasonal dynamics of soil carbon efflux and simulated rhizosphere respiration in a beech forest[END_REF].

The eddy covariance data set

From January through December 1997, CO 2 fluxes were measured in a meteorological tower at 18 m (i.e 3 m above the canopy) using the eddy covariance method [START_REF] Leuning | Eddy-covariance CO 2 flux measurements using open-path and close path CO 2 analysers. Corrections for analyser water vapour sensitivity and damping of fluctuations in air sampling tubes[END_REF]. A complete description of the system used is given by Granier et al. (2000a).

Sensitivity analysis methodology

CASTANEA uses many input parameters; some are coming from literature, other ones are specifically fixed. Most of these parameters are estimated from field sampling, fitting curve technique, interpolation or theoretical evaluation. It is important to identify the input parameters having the greatest effect on model predictions, and to quantify the uncertainty in the model predictions due to the combined effects of uncertainties in the set of input parameters. In this study, the Net Ecosystem Exchange, which sums up most of the model processes, was chosen as the main output variable. The sensitivity study and the subsequent uncertainty study are based on the Monte Carlo technique. The principle was described in [START_REF] Spear | Eutrophication in Peel Inlet II. Identification of critical uncertainties via generalised sensitivity analysis[END_REF] and further developed in [START_REF] Franks | Bayesian estimation of uncertainty in land surface -atmosphere flux predictions[END_REF], and applied in [START_REF] Franks | Multi-objective conditioning of a simple SVAT model[END_REF]. In our study, each parameter is associated with a randomly set error and an ensemble of simulation is performed by a random selection of input parameter values following a Gaussian distribution. [START_REF] Haan | Effect of parameter distributions on uncertainty analysis of hydrologic models[END_REF] have evaluated the impact of parameter distribution assumptions on estimates of model output uncertainty and concluded that good estimates of the means and variances of the input parameters are of greater importance than the actual form of the distribution [START_REF] Haan | Effect of parameter distributions on uncertainty analysis of hydrologic models[END_REF].

Three successive steps are carried out: (i) identification of key parameters through a sensitivity analysis involving all parameters, (ii) determination of the real uncertainty on the key input parameters, and (iii) uncer-tainty analysis involving all the selected key parameters together.

The results of such a study could depend on the site, soil and climate that are considered. Our sensitivity and uncertainty analysis applies to the Hesse site for the year 1997, using the above described parameterisation.

Determination of the key input parameters

In the first step, our purpose was to identify key input parameters, i.e. the parameters whose variations have the largest effect on the model outputs. For this purpose, we applied a +10% or -10% bias on each of the 146 parameters and 23 initial conditions of state variables (Table 2).

For each parameter and initial state variables, we calculated the difference between the Net Ecosystem Exchange computed with +10% (respectively -10%) of a chosen parameter and the reference run no bias . Subsequently the percent of variation (VR +10% and VR -10% ) were calculated as:

VR +10% = 100 × |run no bias -run +10% | run no bias (62)
VR -10% = 100 × |run no biasrun -10% | run no bias (63)

The maximum percent of variation between VR +10% and VR -10% was used to represent the effect of the parameter (or of the initial condition) on the annual NEE.

Afterwards, the parameters and initial conditions were sorted in descending order according to the effect on NEE. All initial conditions and parameters with an effect of more than 2.5% on NEE were kept for the uncertainty analysis.

Estimation of the in situ uncertainty of key parameters

It is a difficult task to estimate the real uncertainty of model parameters. For each parameter or variable selected as described above, we have estimated an in situ error. For this purpose, measurements collected on the Hesse site (or on other beech forests near Paris) or data collected from literature have been used. We obtained estimates of input parameters uncertainties based on errors due to ground measurements sampling, instrumentation error, or curve fitting error. In each case a variation coefficient (in percent) was computed. If no measurement error estimation was possible, a variation coefficient of 10% was chosen.

Evaluation of errors on output variables caused by uncertainty on input parameters

A new set of Monte Carlo simulations was carried out with the uncertainties determined above. This time, the key parameters and variables were all randomly selected with the same procedure described in step (i), except that the standard deviations are not fixed to 10% but correspond to the uncertainty estimated in step (ii). All selected key parameters are processed together to obtain the whole uncertainty on the output variables. This process is repeated many times (e.g. 1000 × n, where n is the number of selected key parameters) to generate a wide range of input parameter datasets and to be able to get statistical information on the model outputs. The effect of uncertainties on the key parameters are tested on eight main annual output variables: the net ecosystem exchange (NEE), the gross primary production (GPP), the autotrophic respiration (RA), the heterotrophic respiration (RH), the aerial wood growth (DBwood) in gC m -2 year -1 and the transpiration (Tr) and evapotranspiration (ETR) in mm year -1 .

The mean and the standard deviation of each annual output variable (NEE, GPP, RA, RH, Tr, ETR) of the 1000 × n runs were calculated, as well as the variation coefficient in percent. The mean daily and standard deviation of NEE were also computed, with the aim to compare them with the daily CO 2 fluxes measured by the eddy covariance technique.

Results and discussion

Determination of key input parameters

Among the 169 parameters and initial conditions, 17 parameters with an effect greater than 2.5% on NEE have been selected as key input parameters (see Table 4). We can distinguish five types of key parameters:

-the dependency between V cmax and leaf nitrogen density, the quantum yield, the ratio between V Cmax and V Jmax , and the slope of the Ball relationship (g1 max ), which describe the response of photosynthesis and transpiration on meteorological variables and on biochemical parameters. -the leaf nitrogen content, the maximum sunlit leaf mass per area, the clumping factor and the leaf area index, which represent biochemical and structural input parameters. -the soil volumetric humidity at field capacity and the soil volumetric humidity at wilting point controlling the extractable water in the soil and therefore the drought effect. -the fine roots construction cost, the leaves construction cost, the nitrogen dependency for all organs and biomass of trunk and branches allowing the estimation of autotrophic respiration. -the critical values of state of forcing for budburst and the critical values of state of forcing for leaf fall, which control the phenology.

Among all parameters and initial conditions, only six key parameters have an effect greater than 10% on NEE when they are affected by an Gaussian noise of 10%: soil water content at field capacity, dependency between V cmax and leaf nitrogen density, soil water content at wilting point, leaf nitrogen content, leaf mass per area of sun leaves and the quantum yield.

For an improved sensitivity analysis, it would be interesting to realize this work on several years, to test the possible interannual effects and to have an estimation of the real in situ uncertainty of each parameter. Such a study, coupled with a process uncertainty study would allow us to reduce the number of input parameters and to simplify the model structure.

Evaluation of the in situ uncertainty of key parameters

The estimated in situ uncertainty (in percent) for the 17 key input parameters are given in Table 5. For four parameters this inference was not achievable. The estimated errors ranged between 3% and 22% with a mean value of 10.2%. The parameters which describe the response of photosynthesis on meteorological variables and on biochemical parameters had a lower uncertainty than the other ones. For 14 key input parameters we were unable to calculate the additional uncertainty due to measurement errors, which was neglected. For this reason the total error is probably underestimated.

The parameters uncertainties are not strictly comparable because they were estimated in different ways: instrumentation error, sampling error (spatial variability) and curve fitting error. Furthermore, different scales and sites are used to assess the spatial variability. 10.0 * * No estimation was available and a mean value of 10% was taken.

Effect of input parameters uncertainty on output variables

Annual results

The results of the standard simulation and of the ensemble of 17,000 simulations including the 17 key parameters uncertainties are given in Table 6 and Fig. 3. For all output variables except heterotrophic respiration (RH), the mean simulation over 17,000 simulations is smaller than the simulation with the mean input parameters. This result is produced when highly non linear responses f(x) are involved in the model (including threshold effect), and when the current parameter value x is close to the threshold. In our case, this effect is not negligible (the difference between f (x) and f (x) reaches 49 gC for the annual Fig. 3. Annual main output variables (see text for explanation) simulated with mean parameter (black chart). Mean over 15,000 simulations of the same output variables with a random Gaussian noise on key parameters (white chart). The standard deviation used for the random process corresponds to the values reported Table 4. Error bar represents the standard deviation of the output variables. The main output variables chosen are the opposite of net ecosystem exchange (NEE), the gross primary production (GPP), the autotrophic respiration (RA), the heterotrophic respiration (RH), the total wood growth (DB) in gC m -2 year -1 and the transpiration (Tr) and evapotranspiration (ETR) in mm year -1 . NEE). A detailed study shows that in our case, two parameters are the main cause of the observed difference between f (x) and f (x): θ fcsoil and θ wiltsoil , i.e. the soil volumetric humidities at field capacity and wilting point, both linked to the Soil Extractable Water (SEW = RW fc -RW wilt ). In our case an increase of the SEW above the prescribed value has no effect on the NEE, whereas a decrease in SEW strongly affects NEE. This indicates two things: first, the water stress was moderate in our site during the study (we are near the plateau of f(x): an increase of SEW has little impact on NEE). Second, there is a strong sensitivity to water stress (would there be one): when SEW decreases, NEE largely decreases (i.e. f(x) depends greatly on x). These particular conditions lead to a strong nonlinearity around the prescribed value of x, causing a bias in NEE: choosing a mean value for θ fcsoil and θ wiltsoil is not neutral, and the uncertainty in these parameters has to be taken into account in our site and year of study. Generally speaking, this result shows that using the mean value of input parameters instead of explicitly taking into account their spatial variation over the stand leads to significant errors. Therefore, Monte-Carlo approach, in addition to produce estimations of uncertainty in the output variables, corrects the biased results obtained when single values (mean values) of parameters are used, which is almost always the case in current simulation studies.

The output uncertainty is ∼10% on GPP, RV and DB and ∼30% on NEE (see Table 6 and Fig. 3). The large effect on NEE in percent is due to its small value compared to others fluxes and also because NEE was chosen as the key output variable. The standard deviation on NEE corresponds to only 120 gC and cannot fully explain the difference observed between the simulations (NEE = 414 gC m -2 year -1 ) and the measurements (NEE = 257 gC m -2 year -1 ).

On hydrological output variables the uncertainty is slightly smaller: respectively 6% and ∼12% for ETR and Tr (note, however, that we selected the parameters for their effect on the NEE only).

One must keep in mind that these model uncertainties correspond only to parameters (and initial conditions) errors and do not include the model uncertainty, i.e. uncertainties on processes and modelling errors.

Daily results

The seasonal course of NEE simulation is given in Fig. 4, together with the daily uncertainty estimated using the method described above, and compared to the measurements. The first, logical, observation is that the uncertainty on NEE is greater during the vegetative season. The simulated uncertainty is particularly important at the end of summer (August and September), when a drought occurred. The uncertainty due to errors in leaf fall simulation is also observed in October. One notes that the parameterization uncertainty could explain a part of the difference between eddycovariance measurements and simulations during the winter and the beginning of the leafy season. However, during summer, even if we take into account the uncertainty on simulations, the model seems not to be able to assess the extreme values of measured NEE. On the other hand, the decrease of NEE measured late in summer is fairly reproduced by the model, except in October. This NEE decline might be caused by a delayed action of the summer drought or by a decrease in photosynthetic capacities, which was not taken into account by the model [START_REF] Wilson | Spatial and seasonal variability of photosynthetic parameters and their relationship to leaf nitrogen in a deciduous forest[END_REF].

Conclusion

A model simulating the carbon and the water budget at stand scale has been developed. In this study, attention was paid to obtain a precise model parameterization. A comparison with independent fluxes data measured by eddy covariance allows an accurate validation of the model in the Hesse Euroflux site. An overestimation of the simulated NEE was observed and different leads could be investigated to better understand why. The first one is to evaluate the effect of the real uncertainties on key input parameters on the simulated NEE. Seventeen key input parameters lead to an uncertainty of 30% on simulated NEE. The mean of 17,000 simulations with input parameters randomly selected around their mean value is smaller than the simulation using the mean input parameters. The proposed procedure to estimate bias, which is due to non linearities in the model, is an important result and should be further studied and taken into account in future works.

Nevertheless, the difference between modelled and measured NEE is not completely explained by the uncertainty on the key input parameters and must be further understood. Two other hypotheses can be put forward to explain the fact that modelled fluxes, including uncertainties, still differ from measurements: either some processes are not correctly described in the model or some discrepancies exist between the scale at which the model is parameterized and the scale at which the model is validated. To study these issues, it is necessary to test the model for each process with independent measurements and to estimate the carbon and the water budget from the organ to the stand scales.

where ε v is the leaf emissivity.

τ t = τ v 1 -(1 -ε g )ρ v
where ε g is the soil emissivity

ρ t = ρ v + (1 -ε g )τ v 2 1 -(1 -ε g )ρ v ε t = (1 -ρ v ) [1 -(1 -ε g )ρ v ]
These formulations are the adaptation for diffuse radiation of Eqs. ( 11), ( 14), ( 16) and ( 17) published in FR02.

Appendix B. Multi-layer radiative balance in the PAR and the NIR

In CASTANEA the radiative balances in the PAR and NIR are performed through a multi-layer parameterization to account for the vertical distribution of absorbed radiation into the canopy. A forest with LAI = 6 is typically divided into 30 elementary layers.

The SAIL model is used to compute the five elementary reflectances (ρ) and transmittances (τ) for a single elementary layer k: ρ dd , ρ sd , τ dd , τ sd , τ ss (where the subscript s stands for direct radiation and d for diffuse radiation). The input data for the SAIL model are the leaf reflectance and transmittance in the PAR (ρ PAR and τ PAR ), the Leaf Inclination Distribution Function (LIDF), the sun elevation and the thickness of the elementary vegetation layer ∆ (in LAI units).

Then, for each layer k, the corresponding multiple scattering reflectances (R) and transmittances (T) for layer k,R k dd , R k sd , T k ss , T k dd , T k sd are computed:

T k ss = τ ss T k dd = τ dd (1 -R k-1 dd ρ dd )(1 -ρ dd R k+1 bot dd )(1 -R k+1 bot dd τ dd R k-1 dd τ dd ) T k sd = τ sd (1 -R k-1 dd ρ dd ) + τ sd R k+1 bot dd τ dd R k-1 dd τ dd (1 -R k-1 dd ρ dd )(1 -ρ dd R k+1 bot dd )(1 -R k+1 bot dd τ dd R k-1 dd τ dd ) + τ ss R k-1 sd ρ dd (1 -R k-1 dd ρ dd ) + ρ sd R k+1 bot dd τ dd (1 -R k-1 dd ρ dd )(1 -ρ dd R k+1 bot dd )(1 -R k+1 bot dd τ dd R k-1 dd τ dd ) + τ ss R k-1 sd ρ dd R k+1 bot dd τ dd R k-1 dd τ dd (1 -R k-1 dd ρ dd )(1 -ρ dd R k+1 bot dd )(1 -R k+1 bot dd τ dd R k-1 dd τ dd ) R k dd = ρ dd (1 -ρ dd R k+1 bot dd )(1 -R k+1 bot dd τ dd R k-1 dd τ dd ) R k sd = ρ sd
(1 -ρ dd R k+1 bot dd )

+ ρ sd R k+1 bot dd τ dd R k-1 dd τ dd (1 -R k-1 dd ρ dd )(1 -ρ dd R k+1 bot dd ) (1 -R k+1 bot dd τ dd R k-1 dd τ dd ) R k bot dd = ρ dd (1 -R k-1 dd ρ dd )(1 -R k+1 bot dd τ dd R k-1 dd τ dd ) T k
bot dd = T k dd R bot and T bot correspond to reflectances and transmittances involving the bottom of layer k (concerning upcoming radiation, instead of descending radiation).

R k dd , R k sd , T k ss , T k dd and T k sd include multiple scattering with the lower (k -1) and upper (k + 1) layer. The calculation of T dd , T sd , R dd , R sd , R botdd and T botdd are recursive (the calculations are done for the top layer first, then the second layer, etc., down to the soil layer, and then again). They generally converge after five iterations (five iterations are performed in the model). The subsequent iterations provide only negligible corrections. Note that the multiple reflections are taken into account to infinite order through the denominators. Iterations are required only because of the (k -1) reflectances and transmittances that appear in the equations (coefficients concerning the lower layer) and which are still unknown when coefficient for layer k are computed. They are initialized to their non multiple scattering value, ρ dd , ρ sd , τ dd , τ sd , and τ ss , respectively, for the first iteration.

The details for the four terms of multiple scattering R sd (direct incident radiation reflected to diffuse radiation) are given here for illustration:

ρ sd (1 -ρ dd R k+1
bot dd ) Single reflection. Subsequent multiple scattering between the top of layer k and the bottom of the upper layer k + 1 are taken into account through the denominator.

ρ sd R k+1 bot dd τ dd R k-1 dd τ dd (1-R k-1 dd ρ dd )(1 -ρ dd R k+1 bot dd )(1 -R k+1 bot dd τ dd R k-1 dd τ dd )
The direct incident radiation if first reflected and scattered on the top of layer k, then on the bottom of the upper layer (k + 1), then transmitted through layer k down to layer (k -1), reflected on the top of layer (k -1) and transmitted again upwards through layer k. Subsequent multiple scattering are taken into account through the denominator following three possibilities (more reflections between the top of layer (k -1) and the bottom of layer k, more reflections between the bottom of the layer (k + 1) and the top of layer k, more reflections between the bottom of the layer (k + 1) and the top of layer (k -1) through layer k).

In this example all possibilities are taken into account for a direct incident radiation to be reflected and scattered into diffuse radiation by layer k, including all multiple scattering with layers (k + 1) and (k -1): this is the definition of R sd .

The calculations of T dd , T sd , R dd , R sd , R botdd and T botdd follow the same pattern and may be easily reconstructed following the same approach.

Intop k dd = 1 (1 -ρ dd R k+1 bot dd ) + R k+1 bot dd τ dd R k-1 dd τ dd (1 -R k-1 dd ρ dd )(1 -ρ dd R k+1 bot dd )(1 -R k+1 bot dd τ dd R k-1 dd τ dd ) Intop k sd = ρ sd R k+1 bot dd
(1 -ρ dd R k+1 bot dd )

+ ρ sd R k+1 bot dd R k+1 bot dd τ dd R k-1 dd τ dd (1 -R k-1 dd ρ dd )(1 -ρ dd R k+1 bot dd )(1 -R k+1 bot dd τ dd R k-1 dd τ dd ) + τ sd R k-1 dd τ dd R k+1 bot dd (1 -R k-1 dd ρ dd )(1 -ρ dd R k+1 bot dd )(1 -R k+1 bot dd τ dd R k-1 dd τ dd )
We are then able to compute the incident direct I k s and diffuse I k d radiation above layer k: Then, using the layer absorption coefficients a s and a d , for diffuse and direct radiation respectively, we can compute the direct aPAR k s and diffuseaPAR k d absorbed radiation by layer k:

I k s =
aPAR k s = a s I s (K) aPAR k d = a d [(Intop k dd + Inbottom k dd )I d (k)Intop k sd + Inbottom k sd )I s (k)]
The absorption coefficients a d and a s for each layer are:

a d = 1 -ρ dd -τ dd a s = 1 -τ ss -ρ sd -τ sd
The incoming coefficients (Intop dd , Inbottom dd , Intop sd , Inbottom sd ) are computed as follows:

Inbottom k dd = τ dd R k-1 dd (1 -R k-1 dd ρ dd )(1 -ρ dd R k+1 bot dd )(1 -R k+1 bot dd τ dd R k-1 dd τ dd ) Inbottom k sd = τ sd R k-1 dd (1 -R k-1 dd ρ dd ) + τ sd R k+1 bot dd τ dd R k-1 dd τ dd (1 -R k-1 dd ρ dd )(1 -ρ dd R k+1 bot dd )(1 -R k+1 bot dd τ dd R k-1 dd τ dd ) + ρ sd R k+1 bot dd τ dd R k-1 dd (1 -R k-1 dd ρ dd )(1 -ρ dd R k+1 bot dd )(1 -R k+1 bot dd τ dd R k-1 dd τ dd )
The approach is very similar to the one that allowed the calculation of R k dd , R k sd , T k ss , T k dd , T k sd . The difference, since all R k dd , R k sd , T k ss , T k dd , T k sd are known, is that iterations are not required any more. For instance, Intop dd is the a factor (greater than 1) indicating that incident diffuse radiation on top of layer k is greater than I dd because of multiple scattering with the bottom of upper layer (k + 1) (first term) and multiple scattering, through layer k, with the top of layer (k -1) (second term). where I PAR s0 is the incident direct radiation in the PAR. Furthermore the factor (a s )/(∆(1 + τ ss )/2) is in fact independent of LAI and layer thickness ∆, and only depends on leaf reflectance ρ PAR , leaf transmittance τ PAR and LIDF.

Fig. 1 .

 1 Fig. 1. Flow diagram for both water (left) and carbon (right) sub-models. Right angle boxes are state variables (C and H 2 O), rounded angle boxes are flow variables (C and H 2 O) and elliptic boxes are forcing meteorological variables. Each arrow corresponds to a water or carbon flow from one compartment to the next one. Dotted arrows correspond to (i) the influence of soil water content on stomatal opening and (ii) the stomatal control of both transpiration and photosynthesis.

Fig. 2 .

 2 Fig. 2. Flow diagram for the soil carbon sub-model. Each boxes is a state variable for soil carbon (i.e. soil organic matter compartments). Except input from leaf and root litters, each arrow (simple) correspond to a carbon flow (decomposition rate) from one compartment to the next one with the double arrow output corresponding to the proportion of flow respired (heterotrophic respiration). γ Rh4 depends on soil textural properties.

•

  Litter: three steps are considered -Input: incoming water is the sum of throughfall and rain reaching directly the soil through the canopy. Note that the water flowing along the bark (Ec) goes directly to the top-soil. RW lit int = RW lit(t-1) + PSd soil + Eg leaf ×(1 -C bark ) + Eg bark (45) RW lit int is an intermediate variable (transitory water content before evaporation and drainage).

Fig. 4 .

 4 Fig.4. Daily NEE measured by eddy covariance technique (black circle) and simulated (white circle). Simulations data corresponds to the mean over 15,000 simulations with a yearly random Gaussian noise on key parameters and the error bar is the standard deviation. Here a net input of carbon to the atmosphere is positive.

  direct and diffuse radiation I s and I d above layer k, we have to compute the incoming coefficients (Intop dd , Inbottom dd , Intop sd , Inbottoms dd ) to obtain the incoming diffuse and direct radiation at each layer k (the incoming radiation takes into account the multiple scattering).

  aPAR k s and aPAR k d are then used to compute the radiation absorbed by the sunlit and shaded leaves, aPAR k sunlit and aPAR k shade layer thickness (in LAI units) and p k sun the mean proportion of sunlit leaves in layer k: ss is the total transmittance for the direct radiation at layer k:TT k ss = (T k ss ) n layer +1-kn layer is the total number of layers, about 30 in the model. aPAR k sunlit and aPAR k shade are the quantities that finally enter the photosynthesis module. Note: the direct radiation aPAR k s /∆ • p k sun absorbed by the sunlit leaves is in fact constant and equal to: aPAR k s ∆ • p k sun = a s (∆ • (1 + τ ss ))/2 I PAR s0

Table 1

 1 List of parameters

	Symbol	Description	Units				Value	Reference
	Radiation and canopy structure
	Long	Longitude	Degree			48.66
	Lat	Latitude	Degree			-7.08
	ρ PAR	Leaf reflectance for PAR radiation	Dimensionless	0.05	mean typical value (Franc ¸ois, pers. com.).
	τ PAR	Leaf transmittance for PAR radiation	Dimensionless	0.11	mean typical value (Franc ¸ois, pers. com.).
	ρ NIR	Leaf reflectance for NIR radiation	Dimensionless	0.32	mean typical value (Franc ¸ois, pers. com.).
	τ NIR	Leaf transmittance for NIR radiation	Dimensionless	0.26	mean typical value (Franc ¸ois, pers. com.).
	ρ soildNIR	Dry soil reflectance for NIR radiation	Dimensionless	0.10	Nagler et al. (2000)
	ρ soilwNIR	Wet soil reflectance for NIR radiation	Dimensionless	0.28	Nagler et al. (2000)
	α L	Average leaf angle inside the canopy	Radians			0.42	Planchais and Pontailler, (1999)
	LMA sunmax	Maximum sunlit leaf mass per area	g dm m -2			101	Montpied (pers. com.)
	k LMA	Leaf mass per area decrease coefficient	Dimensionless	0.187	Montpied (pers. com.)
	WAI	Wood area index	m wood	2 m soil	-2	0.75	Damesin et al. (2002)
	L max	Maximum leaf area index	m leaf	2 m soil	-2	5.6	Bréda (pers. com.)
	Agreg	Clumping factor						0.79	Soudani (pers. com., data 2000)
	Carbon
	Leaf and canopy photosynthesis
	Ca	Atmospheric carbon dioxide concentration	mol CO 2 mol air -1 370	Granier (pers. com.)
	α Na	Dependency between V Cmax and leaf	mol CO 2 g N	-1 s -1 20	Liozon et al. (2000)
		nitrogen density
	β	Ratio between V Cmax and V Jmax	Dimensionless	2.1	Liozon et al. (2000)
	α	Quantum yield	mol electrons (mol	0.292	Ehleringer and Björkman (1977)
			quanta) -1
	θ	Curvature of the quantum response of the	Dimensionless	0.1	fixed
		electron transport rate
	g bCO 2 g 0	Leaf boundary layer conductance Leaf cuticular conductance	mol CO 2 m -2 s -1 mol H 2 O m -2 s-1	1 0.001	calculated from leaf size fixed
	g 1max	Slope of the Ball relationship (maximum	Dimensionless	11.8	Medlyn et al. (2001)
		value)
	g 1min	Slope of the Ball relationship (minimum	Dimensionless	0.001	fixed
		value)
	Maintenance respiration
	Q 10 stem	Temperature effect for stem wood	Dimensionless	1.7	Damesin et al. (2002)
	Q 10 branches	Temperature effect for branches	Dimensionless	2.8	Damesin et al. (2002)
	Q 10 coarse roots Temperature effect for coarse roots	Dimensionless	1.7	fixed
	Q 10 fine roots	Temperature effect for fine roots	Dimensionless	2.2	Epron et al. (2001)
	Q 10 leaves	Temperature effect for leaves	Dimensionless	2.1	Vose and Bolstad (1999)
	T MR	Base temperature for maintenance	• C					1 5	Damesin et al. (2002)
		respiration
	N mleaves	Leaves nitrogen content	mg N g dm			24.2	Montpied (pers. com.)
	N mbranches	Branch nitrogen content	mg N g dm			5.50	Ceschia et al. (2002)
	N mtrunks	Trunk nitrogen content	mg N g dm			1.20	Ceschia et al. (2002)
	N mcoarseroots Coarse roots nitrogen content	mg N g dm			1.20	fixed
	N mfineroots	Fine roots nitrogen content	mg N g dm			9.90	Van Praag et al. (1988)
	MRN	Nitrogen dependency for all organs	mol CO 2 g N	-1 h -1	5.5 × 10 -4 Ryan (1991)
	P leaf inhib	Inhibition of leaf respiration at light	Dimensionless	0.62	Villar et al. (1995)
	Leaf phenology: budburst, growth and fall
	N Start1	The date of onset of rest	Julian day			25	Dufrêne (unpublished data)
	T 2

Table 1 (

 1 Continued )

	Symbol	Description	Units			Value	Reference
	T 3	Base temperature for leaf growth	• C			0	Fixed
	LA max	Maximum area reached by leaf (in	m 2			0.002	Le Dantec et al. (2000)
		average)
	N Start2	The date of onset of ageing	Julian day	213	Fixed
	T 4	Base temperature for forcing leaf fall	• C			20	Fixed
	F critLfall	Critical value of state of forcing (from	• C			225	calibrated
		N Start2 to leaf fall period)
	P fall	Leaf ageing parameter	Dimensionless	0.4	calibrated
	F critLMA	Critical value of state of forcing (from	• C			424	estimated (see text)
		leaf development to leaf maturity)
	Allocation
	AG aerialwood	Aerial wood allocation coefficient	Dimensionless	0.51	calculated as a difference
	AG coarseroots	Coarse roots allocation coefficient	Dimensionless	0.10	ratio to aerial (0.2)
	AG fineroots	Fine roots allocation coefficient	Dimensionless	0.22	estimated
	AG storage	Storage allocation coefficient	Dimensionless	0.18	Barbaroux (2002)
	TO fineroots	Fine roots turn over	jour -1			1/365	Bauhus and Bartsch (1996)
	P branch	Proportion of branches among aerial	Dimensionless	0.125	Damesin et al. (2002)
		wood
	P alive branch	Proportion of living cells in branches	Dimensionless	0.37	Ceschia et al. (2002)
	P alive trunk	Proportion of living cells in trunks	Dimensionless	0.21	Ceschia et al. (2002)
	Growth respiration
	CR leaves	Leaf construction cost	gC gC -1	organ	1.20	Niinemets (1999)
	CR wood	Wood construction cost	gC gC -1	organ	1.38	Damesin et al. (2002)
	CR fine roots	Fine roots construction cost	gC gC -1	organ	1.28	Ågren and Axelsson (1980)
	Heterotrophic (soil) respiration
	Rmob leaf	Fraction of N not remobilized before leaf	Dimensionless	0.6	Montpied (pers. com.)
		dead
	Rmob fineroots Fraction of N not remobilized before fine	Dimensionless	1.0	fixed
		root dead
	Lig leaf	Lignin concentration in leaves	mg g dm	-1	196	Le Dantec (2000)
	Lig fineroot	Lignin concentration in fine roots	mg g dm	-1	229	Chen et al. (2002)
	T b (top)	Base temperature for decomposition rate	• C			35.0	Parton et al. (1987)
		in the top soil layer
	h (top)	Parameter for temperature effect on DR	Dimensionless	2.63	Parton et al. (1987)
		in the top soil layer
	T b (deep)	Base temperature for decomposition rate	• C			25.0	Parton et al. (1987)
		in the deep soil layer
	h (deep)	Parameter for temperature effect on DR	Dimensionless	4.9	Parton et al. (1987)
		in the deep soil layer
	T B (surface)	Base temperature for decomposition rate	• C			25.0	Le Dantec (2000), Käärik (1974),
		at the surface					Rayner and Boddy (1988)
	h (surface)	Parameter for temperature effect on DR	Dimensionless	4.9	Le Dantec (2000), Dommergues and
		at the surface					Mangenot (1970), Käärik (1974), Rayner
							and Boddy (1988)
	T silt (top)	Fraction of silt in the top soil layer	Dimensionless	0.68	Farque (1997)
	T clay (top)	Fraction of clay in the top soil layer	Dimensionless	0.26	Farque (1997)
	T sand (top)	Fraction of sand in the top soil layer	Dimensionless	0.06	Farque (1997)
	T silt (deep)	Fraction of silt in the deep soil layer	Dimensionless	0.62	Farque (1997)
	T clay (deep)	Fraction of clay in the deep soil layer	Dimensionless	0.32	Farque (1997)
	T sand (deep)	Fraction of sand in the deep soil layer	Dimensionless	0.06	Farque (1997)
	K LLst	Maximum decomposition rate for the	Day -1			1.07 × 10 -2 Waelbroeck (1995)
		surface structural litter pool
	K LLm	Maximum decomposition rate for the	Day -1			4.05 × 10 -2 Waelbroeck (1995)
		surface metabolic litter pool

* Estimated from vertical profiles of soil water content

Table 2

 2 List of state variables

	Symbol	Description	Units		Initial value	Reference
	Carbon plant
	B leaf	Leaf biomass	g C m soil	-2	0
	B branch	Branch biomass	g C m soil	-2	519	Damesin et al. (2002)
	B trunk	Trunk biomass	g C m soil	-2	3691	Damesin et al. (2002)
	B coarseroot	Coarse roots biomass	g C m soil	-2	815	Le Goff and Ottorini (2001)
	B fineroot	Fine roots biomass	g C m soil	-2	171	Bauhus and Bartsch (1996)
	B storage	Carbohydrate storage biomass	g C m soil	-2	160	Barbaroux (2002)
	Carbon soil (three layers: surface litter, top and deep soil)
	C m (lit)	Metabolic fraction of litter (leaves)	g C m soil	-2	43	Epron et al. (2001)
	C St (lit)	Structural fraction of litter (leaves)	g C m soil	-2	335	Epron et al. (2001)
	C a (lit)	Active C pools (micro-organism) of litter	g C m soil	-2	35	Epron et al. (2001)
	C m (top)	Metabolic fraction of top soil (fine roots litter)	g C m soil	-2	22	Epron et al. (2001)
	C St (top)	Structural fraction of top soil (fine roots litter)	g C m soil	-2	320	Epron et al. (2001)
	C a (top)	Active C pools (micro-organism) of top soil	g C m soil	-2	230	Epron et al. (2001)
	C S (top)	Slow C pools of top soil	g C m soil	-2	2833	Epron et al. (2001)
	C P (top)	Passive C pools of top soil	g C m soil	-2	827	Epron et al. (2001)
	C m (deep)	Metabolic fraction of deep soil (fine roots litter)	g C m soil	-2	7	Epron et al. (2001)
	C St (deep)	Structural fraction of deep soil (fine roots litter)	g C m soil	-2	113	Epron et al. (2001)
	C a (deep)	Active C pools (micro-organism) of deep soil	g C m soil	-2	104	Epron et al. (2001)
	C S (deep)	Slow C pools of deep soil	g C m soil	-2	1442	Epron et al. (2001)
	C P (deep)	Passive C pools of deep soil	g C m soil	-2	2400	Epron et al. (2001)
	Water
	R leaf	Amount of water on leaves in the canopy	mm		0
	R bark	Amount of water on bark	mm		0
	RW lit	Amount of water in soil surface litter	mm		0
	RW top	Amount of water in top soil layer	mm		RW top fc
	RW soil	Amount of water in the deep soil layer	mm		RW soil fc

Table 4

 4 List of key parameters and their effect (max(VR +10% ,VR -10% )), Eq. (62) and (63) on output variables when a 10% bias is applied on these key parameters

	Description	Symbol	NEE	GPP	RH	RA	DBwood	Tr	ETR
	Soil volumetric humidity at field capacity	θ fcsoil	23.9	8.5	0.3	2.4	8.6	8.6	4.0
	Dependency between V cmax and leaf nitrogen density	α Na	13.5	5.2	0.1	2.4	6.0	4.0	1.8
	Soil volumetric humidity at wilting point	θ wiltsoil	12.9	4.6	0.4	1.6	5.9	6.5	3.1
	Leaf nitrogen content	N mleaves	10.9	5.4	0.1	4.3	4.8	3.6	1.6
	Leaf mass per area of sun leaves	LMA sunmax	10.3	5.2	0.2	4.2	5.7	3.4	1.5
	Quantum yield	α	10.2	4.0	0.1	1.9	4.7	2.7	1.2
	Wood construction cost	CR wood	8.6	0.0	0.0	5.4	4.1	0.0	0.0
	Nitrogen dependency for all organs	MRN	8.3	0.2	0.0	5.6	2.8	0.3	0.1
	Ratio between VC max and VJ max	β	7.9	3.1	0.1	1.5	3.6	2.1	0.9
	Critical value of state of forcing for budburst	F critBB	6.5	3.0	0.2	2.2	3.5	2.2	1.3
	Clumping factor	Agreg	5.2	2.0	0.3	0.9	2.3	1.4	0.8
	Leaf construction cost	CR leaves	4.1	0.0	0.0	2.6	0.6	0.0	0.0
	Critical value of state of forcing for leaf fall	F critLfall	3.7	1.5	0.5	0.9	1.4	0.6	0.2
	Maximum leaf area index	L max	3.5	2.3	0.8	2.3	2.5	1.9	0.6
	Fine roots construction cost	CR fine roots	3.0	0.0	0.2	2.0	1.6	0.0	0.0
	Trunk and branch biomass	Btrunk + Bbranch	3.0	0.0	0.0	1.9	0.8	0.0	0.0
	Slope of the ball relationship (maximum value)	g 1max	2.6	0.7	0.5	0.1	1.0	7.8	3.6

Table 5

 5 Estimation of in situ uncertainty (%) of the key parameters and type of error which are taken into account

		Stand variability	Measurement error	Fit error	Total uncertainty (%)	Sources
	θ fcsoil	X			1 2 .1	Granier et al. (2000b)
	α Na	X		X	6 .1	Liozon et al. (2000)
	θ wiltsoil				10.0 *
	N mleaves	X	X		11	Unpublished data
	LMA sunmax	X	X		8 .8	Unpublished data
	α	X		X	8	Unpublished data
	CR wood	X		X	5	Damesin et al. (2002)
	MRN			X	10	Ryan (1991)
	β	X		X	3	Liozon et al. (2000)
	F critBB	X			4 .1	Bouriaud (2003)
	Agreg	X			14	Soudani personal communication
	CR leaves				10.0 *
	F critLfall	X			2 2 .2	Bouriaud (2003)
	L max	X	X		1 1	Dufrêne and Bréda (1995)
	CR fine roots				10.0 *
	B trunk + B branch	X			1 7 .6	Bouriaud (2003)
	g 1max

Table 6

 6 Results of annual simulations: reference simulation compared to mean and standard deviation (S.D.) of the 17,000 ensemble simulations including key parameter uncertainties (The percentage uncertainty, or variation coefficient VC is given for the ensemble simulation (mean/S.D.))

	Output	Reference	Ensemble	Ensemble	Uncertainty
		simulation	mean	S.D.	(%)
	NEE (gC year -1 ) 463	414	120	29
	GPP (gC year -1 ) 1518	1456	151	10
	RA (gC year -1 )	735	722	75	10
	RH (gC year -1 )	321	321	5	2
	DBwood	363	348	48	14
	(gC year -1 )
	ETR (mm year -1 ) 596	590	33	6
	Tr (mm year -1 )	311	302	37	12

Acknowledgments

We are grateful to the Office National des Forêts (ONF, the French National Office for Forestry) for giving us facilities for in situ measurements. Financial support was provided by the European community through the Euroflux program (contract ENV4-CT95-0078 and CarboEurope-IP). The author acknowledge Bernard Saugier and Franz Badeck for valuable discussions and J.Y Pontailler for his comments. The authors thanks the anonymous reviewers for their constructive remarks that significantly helped strengthen the paper

Appendix A. Radiative balance in the thermal infrared

In the thermal infrared (TIR) the radiative balance coefficients are based on the general formalism given in Franc ¸ois (2002) (hereafter referred to as FR02). The soil-vegetation radiative balance coefficients are written according to MOD3 multiple scattering parameterization (Eqs. ( 11), ( 14), ( 16) and ( 17) in FR02) adapted for diffuse radiation. These coefficients are τ v , ρ v (resp. canopy simple transmittance, reflectance and emissivity), τ t , ρ t , (resp. multiple scattering canopy transmittance and reflectance), and ω t and ε t (resp. multiple scattering vegetation emissivity/absorbance and canopy emissivity/absorbance). Using these coefficients, the radiative transfer balance for the vegetation and the soil writes:

where Rn veg (resp. Rn soil ) are the net vegetation (resp. soil) absorbed TIR radiation, Rv (resp. Rg) are the blackbody TIR radiation emitted by the leaves (resp. the soil) and Ra the incident TIR atmospheric radiance. Ra is computed according to [START_REF] Iziomon | Downward atmospheric longwave irradiance under clear and cloudy skies: measurement and parameterization[END_REF], including a dependence on the cloud cover. Rg and Rv are computed according to Stefan equation using the soil and leaf temperatures. The soil and leaf emissivity are assumed to be constant, equal to 0.97.

The simple and mutiple scattering coefficients τ v , ρ v , τ t , ρ t , ω t and ε t are given by: τ v = 1 -σ f , where the hemispherical shielding factor (or interception coefficient) σ f is computed according to [START_REF] Goudriaan | Crop micrometeorology: a simulation study[END_REF].