Biogeosciences Constraining a global ecosystem model with multi-site eddy-covariance data

Assimilation of in situ and satellite data in mechanistic terrestrial ecosystem models helps to constrain critical model parameters and reduce uncertainties in the simulated energy, water and carbon fluxes. So far the assimilation of eddy covariance measurements from flux-tower sites has been conducted mostly for individual sites (“single-site” optimization). Here we develop a variational data assimilation system to optimize 21 parameters of the ORCHIDEE biogeochemical model, using net CO 2 flux (NEE) and latent heat flux (LE) measurements from 12 temperate deciduous broadleaf forest sites. We assess the potential of the model to simulate, with a single set of inverted parameters, the carbon and water fluxes at these 12 sites. We compare the fluxes obtained from this “multi-site” (MS) optimization to those of the prior model, and of the “single-site” (SS) optimizations. The model-data fit analysis shows that the MS approach decreases the daily root-mean-square difference (RMS) to observed data by 22 %, which is close to the SS optimizations (25 % on average). We also show that the MS approach distinctively improves the simulation of the ecosystem respiration ( Reco), and to a lesser extent the gross primary productivity (GPP), although we only assimilated net CO2 flux. A process-oriented parameter analysis indicates that the MS inversion system finds a unique combination of parameters which is not the simple average of the different SS sets of parameters. Finally, in an attempt to validate the optimized model against independent data, we observe that global-scale simulations with MS optimized parameters show an enhanced phase agreement between modeled leaf area index (LAI) and satellite-based observations of normalized difference vegetation index (NDVI).


Introduction
Terrestrial ecosystem models have been a tool of growing importance in order to understand and simulate the behavior of land ecosystems and their response to various disturbances, be it natural or anthropogenic.Mechanistic terrestrial ecosystem models are widely used to assess the current land carbon balance (Sitch et al., 2008) and to predict its future evolution under climate change (Friedlingstein et al., 2006;Cox et al., 2000), as a major driver of the future climate itself.In this context, there has been a growing effort to evaluate and validate the simulated carbon fluxes and stocks against in situ or remote sensing observations.In this study, we investigate the potential of eddy-covariance flux measurements from a dozen FLUXNET sites (http://fluxnet.ornl.gov/;Baldocchi et al., 2001Baldocchi et al., , 2008) ) to improve a global process-based ecosystem model, ORCHIDEE (Krinner et al., 2005).These data provide near-continuous in situ measurements of carbon dioxide, water and energy fluxes; measurements are currently conducted at more than 500 sites, spanning a wide range of climate regimes and vegetation types worldwide.Here, we focus on modeling deciduous broadleaf forests, which are particularly well represented in the FLUXNET database.

S. Kuppel et al.: Constraining a global ecosystem model
The discrepancies between the fluxes simulated by terrestrial biosphere models and the observations have five main origins: errors in flux measurements, errors in meteorological forcings, error in structural equations of the model (including missing processes), inadequate calibration of the model parameters, and erroneous initial state of the model.While the first two types of error are foreign to the biosphere model itself, the last three items are crucial to improve model simulations.Most global biosphere models describe the terrestrial ecosystems with a small number of categories, referred to as plant functional types (PFT, usually on the order of 10-15).With this classification, each PFT is considered sufficiently homogeneous to be described by a single set of equations and parameters at the global scale.In this context, the choice of a representative value for each parameter becomes a critical step that might add significant error to the simulated fluxes for a given PFT.
Numerous studies in various types of ecosystem have illustrated the capacity of data assimilation to provide optimized sets of parameters that significantly improve the model-data fit (see review by Williams et al., 2009).Many of these efforts have used eddy covariance measurements from flux towers to improve their model.Initial studies focused on individual sites to investigate the potential for carbon, latent heat and sensible heat flux measurements to serve as model constraints (Braswell et al., 2005;Knorr and Kattge, 2005;Santaren et al., 2007;Wang et al., 2001Wang et al., , 2007)).However, the small spatial footprint of each flux tower (a few hectares) often resulted in model parameters overly tuned to the specificities of a particular site.A few recent studies have sought independent evaluation using parameters optimized at one location when simulations are run at other sites (Medvigy et al., 2009;Verbeeck et al., 2011).This last approach is useful to evaluate the potential of a given model structure to simulate the spatial heterogeneity in the fluxes (its genericity).Complementarily, Groenendijk et al. (2011) optimized the parameters of a photosynthesis model using groups of sites within and across PFTs.Their results show a large intra-PFT variability of model parameters, which suggests that perhaps the PFT concept needs to be re-defined.
The present study further investigates the potential of the simultaneous assimilation of carbon net ecosystem exchange (NEE) and latent heat (LE) flux measurements from an ensemble of flux towers situated in temperate deciduous broadleaf forests (DBF).We focus on the following questions: -Can we find a single set of parameters that improves the model-data fit at all sites (multi-site optimization)?
-How does the multi-site optimization perform in comparison to independent optimizations at each site (single-site optimization)?
-Which parameters are well constrained by the NEE and LE flux measurements?
-What are the typical time scales of model-data mismatches that are improved by the optimization, and which processes remain poorly captured by the model after the optimization?
-What is the impact of an optimized set of parameters on a global simulation?Can we evaluate changes in the simulated vegetation activity (using the new parameters) with remote sensing data?
In order to investigate these different questions, we use the following methodology.
We selected a set of representative DBF sites for which we would assimilate daily NEE and LE (see Sect. 2.2).We present in Sect.2.3 the optimization approach and the set of parameters that are estimated.Section 2.4 describes the different optimizations/sensitivity tests that are performed.Finally, we introduce in Sect.2.5 the additional data used for the posterior evaluation of the optimization, at the local and the global scale.
Throughout the analysis, the results of the model after the multi-site optimization (hereafter referred as MS) are compared with those of the model optimized independently at each measurement site (single-site optimization, hereafter referred as SS).The results section is divided into five parts.The overall performance of the optimization (i.e., modeldata mismatch) is first evaluated at different time scales (Sect.3.1).Then, we analyze the optimized parameters for each biophysical process separately (Sect.3.2).To this end, the Bayesian inversion scheme allows us to use and interpret statistical information regarding both prior and posterior parameter error correlations.The relevancy of using MS parameters rather than extrapolating a combination of SS parameters is discussed in Sect.3.3, while the choice of the assimilated data is evaluated in Sect.3.4.We further compare the gross primary productivity (GPP) and ecosystem respiration (R eco ) of the optimized model against GPP and R eco derived from in situ observations at the flux towers (Sect.3.5.1).Finally, we investigate the impact of the optimized parameters at the global scale, using remote sensing observations of the vegetation activity (Sect.3.5.2).

ORCHIDEE model
The biogeochemical vegetation model used in this study is ORCHIDEE (Krinner et al., 2005).It calculates the water, energy and carbon fluxes between the land surface and the atmosphere at a half-hourly time step.The exchange of carbon and water during photosynthesis and the energy balance are treated at the smallest time scale (30 min), while carbon allocation, autotrophic respiration, foliar onset and senescence, mortality and soil organic matter decomposition are computed on a daily time step.The equations involving the parameters optimized in the present study can be found in the results section.More extensive descriptions of ORCHIDEE are given elsewhere (Ducoudre et al., 1993;Krinner et al., 2005;Santaren et al., 2007;Verbeeck et al., 2011).
As in most biogeochemical models, the vegetation is grouped into several PFTs, 13 in the case of ORCHIDEE.Except for the modeled phenology (Botta et al., 2000), the equations governing the different processes are generic, but with specific parameters for each PFT.In the present study, ORCHIDEE is mainly used in a "grid-point mode" at one given location, forced with the corresponding local halfhourly gap-filled meteorological measurements obtained at the flux towers.Only in Sect.3.5.2,ORCHIDEE is run at the global scale, forced by the global ERA-Interim meteorology (http://www.ecmwf.int/research/era/do/get/era-interim).Importantly, the modeled carbon pools are initially brought to equilibrium by cycling the available meteorological forcing over a long period (1300 years), with the prior parameterization of the model.It ensures a net carbon flux close to zero over annual-to-decadal time scales.

Eddy covariance flux data
The present study focuses on one single type of ecosystemtemperate deciduous broadleaved forests (DBF).We selected 12 measurement sites (Fig. 1 and Table 1), with a "footprint" that corresponds to a vegetation cover represented by at least 70 % of DBF, the rest being mostly C3 grasslands or evergreen needleleaf trees.The eddy-covariance flux data used are part of the FluxNet network, with standard flux data processing methodologies (correction, gap-filling and partitioning) of the La Thuile dataset (Papale et al., 2006;Reichstein et al., 2005).The gap-filled measured fluxes of net ecosystem exchange (NEE) and latent heat fluxes (LE) at half hourly time step are used to compute daily means.We choose to assimilate daily mean observations and not half hourly measurements in order to focus the optimization on time scales ranging from synoptic to seasonal and inter-annual varia-tions, and to avoid the complicated treatment of the error correlation between half hourly data (Lasslop et al., 2008).The remaining gaps in observations are distributed somewhat evenly along the course of the day.Note that individual days with more than 20 % of missing half-hourly observations were not included in the assimilation; over a total of 43 site-years used in this study, the missing data represent 196 and 2 days in the NEE and LE data streams, respectively.

Data assimilation system
The model parameters are optimized using a variational data assimilation method.Within this Bayesian inversion framework, we account for uncertainties regarding the model, the observations, and the prior parameters.The approach is based on Santaren et al. (2007).Assuming Gaussian distribution for errors on both the parameters and the observations, the optimized set of parameters corresponds to the minimization of the following cost function J (x) (Tarantola, 1987): which contains both the misfit between modeled and observed fluxes, and the misfit between a priori and optimized parameters.
x is the vector of unknown parameters, x b the background parameter values, H (x) the model output, and y the vector of observed fluxes.P b describes the prior parameter error variances/covariances, while R contains the prior data error variances/covariances, as described in Sect.2.2.Both matrices are diagonal since uncertainties are supposed to be uncorrelated between parameters/observations, and all parameters/observations are identically weighted in the cost function.The latter is minimized iteratively using a gradientbased algorithm called L-BFGS, which provides the possibility to prescribe boundaries for each parameter (Byrd et al., 1995).The standard deviation for each parameter used for P b variances is equal to 40 % of the range between lower and higher boundaries, which have been carefully specified following the physical and empirical expertise of ORCHIDEE modelers, based on literature reviews or databases (such as TRY, Kattge et al., 2011) providing estimated parameter values from in situ or laboratory measurements.
Regarding the observational error statistics, the error covariance matrix R should include both the error on the measurements and the error on the model process representation.On the one hand, the random measurement error on the observed fluxes is known not to be constant and can be estimated as the residual of the gap-filling algorithm (Richardson et al., 2008).On the other hand, model errors are rather difficult to assess and may be much larger than the measurement error itself.Therefore, we chose to focus on the model error whose correlations cannot be neglected (Chevallier et al., 2006).The difficulty of evaluating the structure of the model error leads us to keep R diagonal and, as compensation, www.biogeosciences.net/9/3757/2012/Biogeosciences, 9, 3757-3776, 2012  1999-2004Cook et al. (2004) ) artificially inflate the variances (Chevallier, 2007).Firstly, the variances in R are defined as constant in time for each type of observation (NEE and LE).Secondly, their values are chosen based on the mean squared difference between the prior model and the observations, multiplied by the inflation factor k σ , which represents our estimation of the characteristic autocorrelation time length (in days) of the model error.
The value of k σ is fixed to 30, which for the error propagation is equivalent to assimilating one observation every 30 days.At each iteration, the gradient of the cost function J (x) is computed, with respect to all the parameters.For most of the parameters, we use the tangent linear (TL) model of OR-CHIDEE to compute the gradient, generated with the automatic differentiator tool TAF (Transformation of Algorithms in Fortran; see Giering et al., 2005).Some processes in the model are described by functions that are not smooth.Examples include the threshold functions controlling the temperature dependence of foliage onset and senescence (parameters K pheno,crit and T senes ; see Table 2).In both cases we use a finite difference approach with prescribed perturbation steps respectively equal to 4 % and 2 % of the allowed variation range.These values were found in prior sensitivity tests as the smallest possible without provoking numerical errors.Once the cost function reaches the minimum, the posterior parameter error variance/covariance matrix P a is explicitly calculated from the prior error variance/covariance matrices (P b and R) and the Jacobian of the model at the minimum of J (H ∞ ), using the linearity assumption (Tarantola, 1987): Large absolute values of error correlation indicate that the observations do not provide enough information to distinguish between the effects of each corresponding parameter throughout the optimization.

Performed optimizations: single-site and multi-site
One goal is to optimize a mean set of parameters using information from several measurements sites simultaneously (MS approach) and to compare the results with optimizations conducted separately for each site (SS approach).We thus have extended the site-scaled inversion algorithm used in Santaren et al. (2007) and Verbeeck et al. (2011) to include the observations from all sites and to share a common set of parameters in the optimization procedure.Concerning the parameters, we distinguish two cases: the site-specific parameters that are only relevant for a particular site (i.e., that cannot be applied to other sites) and the generic parameters that apply to all sites.The initial state of the model is optimized with the only parameter chosen as site-specific: a multiplier of the different soil carbon pool contents, which are closely related to the local land-use history (Carvalhais et al., 2008).The rest of the parameters are considered as generic across sites in this study.We acknowledge that this assumption brings some limitations given the potential inter-site variability of some parameters (e.g., soil water availability), which will be kept in mind in the results analysis.The list of optimized ORCHIDEE parameters is given in Table 2. Besides, two other data assimilation experiments are conducted for the purpose of this study.Firstly, a series of multisite optimizations are performed in order to evaluate the individual impact of the parameters related to heterotrophic respiration, notably regarding the initial soil carbon content (see Sect. 3.2.1.).Second, we use the multi-site procedure to separately assimilate each one of the data stream (NEE and LE), so as to evaluate their respective information content.In the case of LE, non-sensitive parameters were left out the optimization.
In the end, we performed the following optimizations, including sensitivity tests:  -1 MS optimization with only LE data (14 generic parameters)

Photosynthesis and respiration
The model is evaluated at the sites using the two components of the NEE flux: the gross ecosystem productivity (GEP) and the ecosystem respiration R eco , estimated via the fluxpartitioning method described in Reichstein et al. (2005).This method extrapolates nighttime measurement of NEE, representing R eco , into daytime R eco using a short-termcalibrated temperature response function.GEP is then derived as the difference between R eco and NEE.Due to the small quantitative difference between GEP and the gross primary productivity (GPP) (Stoy et al., 2006), from now on we use the term "GPP" instead.Similarly to NEE and LE, we also use daily-averaged data.We acknowledge that GPP and R eco are not fully independent data (with respect to the assimilated NEE) and are essentially model-derived estimates that are to some degree conditional on our underlying assumptions, but note that these concerns have been largely addressed in previous analyses (e.g., Desai et al., 2008).

MODIS remote sensing NDVI
In order to evaluate if the optimized set of parameters improves the phenology of the model at the global scale, we use the normalized difference vegetation index (NDVI) estimated from measurements made by the MODIS instrument, aboard the Terra satellite (see Sect. 3.5.2).MODIS measures irradiances (converted to reflectances) which are first corrected for atmospheric absorption and scattering (Vermote et al., 2002), then from directional effects (Vermote et al., 2009), before NDVI is calculated (Maignan et al., 2011).The result is a daily product with 5 km spatial resolution.Spatial averaging is used to match the ERA-Interim resolution.Note that NDVI was preferred to other satellite products (such as LAI or FAPAR) as it is directly calculated from the surface reflectances, contrary to LAI and FAPAR which require intermediate models to be generated, possibly adding significant uncertainty to the retrieved data (Garrigues et al., 2008).

Model-data fit: multi-site versus single-site optimization
Throughout this section, the results of the MS optimization are systematically compared to those given by the SS optimization.

Overall model-data fit
To evaluate the performance of the different optimizations, we first look at the misfit between the model output and the assimilated observations for three cases: prior model, MS optimization and SS optimization.Figure 2 shows the seasonal cycles of NEE and LE at two of the 12 sites used in this study, where only two years of data are shown.Plots with all years at all the "sites" can be found in the Supplement.Note that, for the sake of clarity, data have been smoothed with a 15-day moving average window in all the figures showing seasonal cycles, but not in any of the optimizations.The brackets give the average annual carbon budget (gC m −2 yr −1 ) and the averaged LE flux (W m −2 ).Lastly, the error bars on the left part represent for both fluxes the total uncertainty averaged over all the periods: where σ param is the parameter error contribution in the observation space, while σ model represents the structural model error.The former is calculated using the parameter error covariance matrix and Jacobian matrix of the model, similarly to Eq. ( 25) in Rayner et al. (2005) (assuming linearity at the minimum of the cost function).The model error is reported at each site as a standard deviation from the statistical analysis of the prior and the posterior residuals (model minus observations), following Eq.( 3) in Desroziers et al. (2005).The all-site average values are estimated to be 1.8 gC m −2 d −1 (NEE) and 24.1 W m −2 (LE).NEE shows a significant seasonal cycle with large negative values (uptake) in summer and positive values (release) during winter, which is captured by the prior model relatively well.There are however three major types of mismatch: -recurrent overestimation of the winter carbon release; -underestimation of the summer carbon uptake at most sites; and -significant phase shift at some sites (onset and senescence).
The first, and to a certain extent the second, item partly results from the model being spun up before each simulation, thus forcing the annual carbon budget to be close to zero (see Sect. 2.1).The optimization generally manages to correct the aforementioned winter bias, and this is clearly visible in the shift of annual carbon budget after both MS and SS optimizations.The summer uptake is also increased after the optimizations, but does not always reach the amplitude of the measured NEE.The phase shift is corrected at some sites, but not consistently at all sites.
The large LE seasonal cycle, with a peak of evapotranspiration in summer and values close to zero in winter, is relatively well reproduced by the prior model.The model tends to overestimate LE at the end of winter and in spring (see also the Supplement).This discrepancy is more pronounced at sites with snow cover (such as JP-Tak and most of the American sites), most likely because the current model overestimates snow sublimation (Slater et al., 2001).The optimizations manage to correct the summer misfit, although less significantly than for NEE.The sublimation-related misfit of LE is not corrected, as we did not include specific parameters of snow buildup and sublimation.Note, however, that sensitivity tests have shown no significant impact of this discrepancy of LE upon the optimization of the other energy balance parameters.
In the end, we found that (1) the MS parameter set significantly improves the model-data fit at most sites, and (2) the results of the MS optimization are often comparable to those of the SS optimizations.

Model-data fit as a function of time scale
To further evaluate and quantify the optimization performances, we analyze the model-data misfit for different time scales: daily, monthly anomaly, monthly average, seasonal average, and yearly average (Fig. 3).The model-data misfit is quantified using the root-mean-square difference (RMS), calculated for each site from the following quantities: x daily = x x monthly anomaly = x month − x month all years x monthly average = x month − x year x seasonal average = x season x yearly average = x year where x is either the observation or the model output vector, x month the average over a month, x month all years the average of x month over all the years available at a given site, x season the average over one season (DJF, MAM, JJA or SON), and x year is the average over the year.On the daily time scale, the RMS is calculated for each site (first row in Fig. 3), while at other time scales all the site-years are used together in one single RMS per case (second and third rows) to increase statistical power.
A large RMS reduction is found on yearly average, where the RMS of NEE is reduced by half from the prior (49 % decrease in MS case, 59 % in SS case).More specifically, the best relative improvement of the fit happens in winter (Fig. 3, second row), where the RMS is reduced by an even larger amount (MS: -55 %, SS: -69 %).The all-site averaged net annual NEE shifts from -39 gC m −2 yr −1 (prior model) to -260 gC m −2 yr −1 (MS) and -251 gC m −2 yr −1 (SS), which is much closer to the observed value (-344 gC m −2 yr −1   parameter sets result in significant improvements in annual carbon and water budgets. On the monthly time scales, we distinguish the average and the anomaly.The RMS reduction on monthly average is significant mostly for NEE (MS: -16 %, SS: -24 %), becoming smaller in the case of LE (MS: -8 %, SS: -18 %).This indicates that only a small improvement of the NEE mean seasonal cycle is possible with the current model structure, and even less for LE.Regarding the monthly anomaly, we see that neither the SS nor the MS optimizations bring any significant improvement.However, the inner monthly interannual variations (RMS between the quantities x month and x month all years introduced in Eq. ( 4)) of NEE and LE are already similar between observations and prior model: 0.63 gC m −2 d −1 and 6.88 W m −2 against 0.41 gC m −2 d −1 and 7.21 W m −2 , respectively.We deduce that the current model structure captures a significant part of the inter-annual variations but no further improvement is given by optimizations.Therefore, we suggest that the remaining gap could be bridged by taking into account the impact of biotic factors and climate anomalies on photosynthesis/respiration in the model structure and/or the parameterization.
On the daily time scale (raw data used in the optimization), the RMS reduction is 22 % on average for NEE after the MS optimization, while it is 25 % in the SS case.Concerning LE, the difference between the two optimizations is slightly more significant, with an average of 16 % RMS reduction versus 23 % in the SS case.The larger SS/MS discrepancy observed for LE reflects the fact that the RMS is significantly higher for the MS case at 3 sites (DK-Sor, FR-Fon and US-UMB).
Overall, the NEE fit is more improved than the LE fit, partly reflecting the larger set of optimized parameters that are exclusively related to carbon assimilation and respiration processes (the focus of this study).

Level of constraint on the different processes
In the following subsections, we analyze the values and the associated errors of the relevant parameters for the main simulated processes and we investigate the strengths and weaknesses of the current model structure.Figure 4 displays the optimized values for each parameter, and the error correlations between all parameters in the MS case can be found in Fig. S13.For each process, the relevant ORCHIDEE equations are shown, with the optimized parameters highlighted by bold fonts.

Parameters of the initial soil carbon pools and heterotrophic respiration
As described in Sect.3.1.1,the optimization strongly reduces the winter carbon release to match the observations, at most sites.In winter, the carbon exchange for DBF is mainly driven by heterotrophic respiration R h , here modeled as the sum of the respirations from four litter pools (metabolic and structural, both above and below ground) and three soil organic matter pools (active, passive and slow): where C p , α p , τ p , c H and c T are the size of the carbon pool, a pool-specific partitioning coefficient based on Parton et al. (1988), the pool-specific carbon residence time also based on Parton et al. (1988), and the temperature and humidity control factors, respectively.The initial sizes of the carbon pools integrate the past input of carbon to the soil and thus all land use history at the site.These impacts are difficult to account for, and we choose to scale all initial soil carbon pool sizes through the optimization procedure with one parameter, K soilC (Table 2), as in Santaren et al. (2007) and Carvalhais et al. (2010): c H and c T represent the slowing down of respiration in too dry soils or at low temperatures, respectively: where H is the relative humidity of the above-ground litter or the soil, and T is the surface/soil temperature for the above/below-ground pools.HR H,b , HR H,c , and Q 10 are critical parameters that are optimized (Table 2).At most sites, both the MS and SS optimizations lead to similar results, with a smaller initial pool size (K soilC lower than 1).This is linked to the spin-up procedure (see Sect. 2.1) that brings the ecosystem to near-equilibrium (no net release or uptake of carbon).However, most of the selected sites are young growing forests (less than 90 year old).Soil carbon content is probably lower in these forests than it would be for a mature forest.As proposed by Carvalhais et al. (2010), the optimization of soil carbon pools can alone explain most of the RMS reduction at yearly time scale for NEE.We have further investigated such a hypothesis with 6 MS experiments where, as compared to the reference MS optimization, one or more of the four parameters related to the heterotrophic respiration (K soilC , Q 10 , HR H,b and HR H,c ) are in turn left out the assimilation procedure.Table 3 shows that without K soilC there is a 38.5 % reduction of the model-data RMS at yearly time scale, a value significantly lower than in the standard case (49.3 %).A similar degraded performance is found when the three other parameters are simultaneously left out of the optimization, whereas we observe a slight improvement (as compared to the standard MS case) when each one of them is individually excluded.When the four parameters are not considered, the optimization procedure becomes significantly less efficient with a 16.6 % RMS reduction on yearly average.With these tests, we conclude that the initial sizes of the carbon pools have the highest individual leverage upon the simulation of the annual carbon mean flux, but the parameters controlling the temperature and humidity dependences of the respiratory processes also play an important role when combined together.
The MS optimization results in an increased value of Q 10 , which represents a middle ground between the spread of SS values.Below 30 • C, the temperature control factor c T is a decreasing function of Q 10 (Eq.8), and so is the heterotrophic respiration.Regarding the humidity control factor c H , both parameters, HR H,b and the offset HR H,c , are decreased after the MS optimization (Fig. 4), reflecting the trend of most SS optimizations.Overall, we observe a significant reduction of R h at 10 out of 12 sites after the MS optimization (not shown).Note that errors of HR H,b and HR H,c are strongly anti-correlated, and also correlated with Q 10 errors (Fig. S13).This indicates that using only NEE and LE measurements does not allow full separation of temperature and humidity impacts on R h .

Parameters of the autotrophic respiration
The autotrophic respiration R a is computed as the sum of the growth respiration R g and the maintenance respirations R m,i from the various biomass pools i: The maintenance respiration follows from , (10) where T i , B i , and LAI are respectively the soil or surface temperature, the biomass content and the leaf area index, while c 0,i (g g −1 day −1 ) is the maintenance respiration coefficient at 0 • C, which is prescribed depending on the PFT and the biomass pool i: 0 (heartwood), 1.19×10 −4 (sapwood, fruits, and carbohydrate reserve), 1.67×10 −3 (roots) and 2.62×10 −3 (leaf).MR a and MR b are two critical parameters that are optimized (Table 2).The growth respiration is calculated as a fraction of the remaining total biomass: where B a is the total biomass, t the time step (one day), and GR frac a fraction to be optimized (Table 2).For all three parameters, the SS optimizations tend to reduce the prior values (Fig. 4).Surprisingly, the MS values correspond approximately to the lowest values in the spread of SS optimizations (and not a median or mean value).This result has to be related to the cross-dependence between parameters and the nonlinearity of the model.Figure S13 shows that MR a and MR b errors are anti-correlated (-0.28), and similarly between MR a and GR frac (-0.21), but also that GR frac /V cmax , MR b /Q 10 , and MR a /Q 10 error pairs are correlated.This indicates that we might face an equifinality problem: a range of different parameter sets may yield similar model performance, and similar model predictions.
Overall, the inversion system always reduces R g and R m , and consequently R a .The main reason for this is the improved fit to the summer carbon fluxes: the prior model always underestimates the magnitude of summer uptake (see Sect. 3.1.1).Carbon assimilation is unchanged or even reduced by the optimization (see next section), and R h is already drastically reduced to fit the winter NEE (Fig. 3).Hence, the optimization further decreases the autotrophic respiration.

Parameters of the photosynthesis
The photosynthesis model developed by Farquhar et al. (1980) and Collatz et al. (1991) is used for C3 plants in ORCHIDEE.GPP is an increasing function of the maximum carboxylation capacity V cmax , which is one of the parameters optimized here.The effective maximum carboxylation capacity V cmax,effective is modulated by several coefficients: where ε leaf , ε temp , and ε water are the leaf efficiency, the dependence on temperature and the dependence on soil water availability, respectively.The leaf efficiency is determined by the relative leaf age, which is a fraction of the critical leaf age parameter L age,crit , following the law shown in Fig. A1 in Krinner et al. (2005).The temperature efficiency varies between 1 at the optimal temperature T opt , and 0 at the minimum and maximum temperatures T min and T max , respectively: In this study, only the parameters T opt and T min are optimized, as we found little sensitivity to the value of T max (fixed to 38 • C) in preliminary tests.The dependence factor on soil water availability ε water is an increasing function of the water fraction f w available for the plant in the root zone, represented by a double bucket scheme (Ducoudre et al., 1993).f w is calculated based on the exponential root profile Hum cste along the soil depth Dpu cste , both of which are optimized: + 1 − x top exp −Hum cste Dpu cste a deep , where x top , a top and a deep are a normalized coefficient related to the wetness of the top soil water reservoir, and the dry fraction of the top/deep soil water reservoirs, respectively.The maximum leaf area index LAI max is also optimized, and the stomatal conductance g s is expressed as following Ball et al. (1987): where G s,slope is optimized.h r , C a and G s,c are the air relative humidity (%), the atmospheric CO 2 concentration, and an offset fixed to 0.01 for the current PFT, respectively.The optimized values of the photosynthesis parameters can be classified in three categories: -those reducing carbon assimilation: decrease of G s,slope , LAI max , Dpu cste and the increase of values for restriction-related parameters such as T opt and T min (Fig. 4); -those slightly increasing GPP: increased SLA (with a large associated error); and -those with no significant trend, such as V cmax .
Overall, this combination of optimized parameters results in a decrease of the carbon assimilation, as discussed in Sect.3.5.1.Note that V cmax has strong error correlations with other parameters: G s,slope (-0.48),T opt (0.73), and L age,crit (-0.97) (Fig. S13).The same applies to the two parameters T opt and T min involved in the temperature dependence of the photosynthesis efficiency, whose errors are significantly anticorrelated (-0.45).These error correlations indicate that various combinations of these parameter values may result in decreases in the carbon assimilation.We checked whether leaving R a parameters out of the optimization would lead to an increase of GPP (as the NEE constraint remains the same in the inversion system, but with an higher modeled R a ), yet the effect is rather weak: the average annual assimilated carbon changes by less than 2 % (not shown), with a substantially worse fit to summertime NEE.We deduce that climate dependencies of photosynthesis and autotrophic respiration are different enough so that the two processes can be distinguished in the ORCHIDEE model with only daily mean NEE data.

Parameters of the evapotranspiration
In this study, we optimize three parameters directly related to the energy balance: G s,slope , Z0 overheight , and K albedo,veg .For a given GPP, G s,slope modulates the stomatal conductance and thus LE (see Eq. 15).Z0 overheight is a characteristic rugosity length used to calculate the aerodynamic resistance to mass transfer.K albedo,veg is a multiplying factor of the vegetation albedo at each time step.G s,slope is generally decreased by the optimizations with a significant spread among SS values (see Sect. 3.2.3),while K albedo,veg is more homogeneously increased.An enhanced albedo decreases the amount of radiation absorbed by the vegetation, hence reducing the evapotranspiration (LE).Besides, a decreased value of G s,slope induces a lower stomatal conductance (for a constant GPP) and thus reduces the transpiration from the leaves.These two factors combine to result in less LE after optimization (Fig. 3 and the Supplement): the prior annually averaged LE flux, equal to 39 W m −2 , is reduced by 23 % and 22 % after the MS and the SS optimizations, respectively.This is consistent with the fact that the prior model overestimates LE throughout the year at most sites (Sect.3.1.1)except for DK-Sor, US-LPH, and US-UMB.At these three sites, LE is underestimated by the prior model; the SS optimization consistently improves the fit by increasing G s,slope and decreasing K albedo,veg (Fig. 4), while the MS trend toward a reduction of LE further enhances the underestimation.Note finally that changes in the values of Z0 overheight are more difficult to interpret, due to the large spread among SS values; the associated errors remain large and are correlated with those from G s,slope (Fig. S13).

Parameters of the phenology
We optimize three critical phenology parameters: K pheno,crit , T senes and L age,crit .K pheno,crit is a multiplicative factor of the growing degree-days (GDD) threshold initiating spring leaf onset.In temperate deciduous broadleaf forests, T senes directly gives the threshold temperature T crit triggering leaf senescence.The higher the value of K pheno,crit (resp.T senes ) is, the later (resp.the earlier) in the year leaf onset (resp.leaf senescence) occurs.For L age,crit the smaller it is, the sooner the leaf loses its photosynthetic efficiency.K pheno,crit is generally increased by the optimization, and so is T senes (Fig. 4).It means that the growing season is delayed and ends earlier, thus being shortened, as compared to the prior model.L agecrit is almost unchanged in MS case, while the shift after optimization oscillates between -40 days and +35 days depending on the SS optimization considered.To further quantify these changes, for each site we computed the smooth seasonal cycle of NEE, based on the signal decomposition proposed in Thoning et al. (1989), and we used the days when the smooth curve crosses the zero line to define the boundaries of the growing season (Fig. 5).Regarding the beginning of the growing season, the prior simulation is too early by 2 to 19 days, except at DK-Sor (Fig. 5, first row).The optimizations improve the agreement with the data with average positive shifts of 8 days in the SS case and 5 days after MS optimization, although overcorrection is apparent at 3 sites (DK-Sor, FR-Fon and FR-Hes).Regarding leaf senescence (Fig. 5, second row), the prior model is generally too late in ending the growing season, by 7 (± 12) days, and both SS and MS optimizations tend to advance the date by 3.5 days on average.Overall, the length of the growing season is significantly improved in both MS and SS cases: the mean prior overestimation by 16 days is reduced to respectively 8 and 5 days.
The seasonality of evapotranspiration is tightly linked to that of photosynthesis.The reduction of the latent heat flux after optimization (see Sect. 3.2.4)generally delays the spring increase of LE and advances the subsequent autumnal decrease, in general agreement with the observations (Fig. 2 and S1 to S12).

Is the multi-site set of parameters the most generic one?
We have shown that the MS optimization is able to provide significant improvement in the model results.However, one can wonder if deriving across-site information directly from SS optimizations would not prove to be as efficient.Two methods have been considered: directly using SS sets of parameters at the other sites, or deriving a mean set of parameters calculated as the average of the 12 SS sets of parameters.In Fig. 6, we show the model-data RMSs for NEE and LE at each site, resulting from applying 15 sets of parameters: prior model (grey), the 12 SS sets of parameters (yellow, and red for the parameters optimized at the given site), the MS set of parameters (blue), and the mean set of SS parameters (black).Note that, in the mean set of parameters, we keep using the local SS-optimized carbon pool scaling factor K soilC , instead of an averaged value.Figure 6a first shows that transposing SS set of parameters most often results in NEE RMSs significantly higher than in the corresponding SS and MS cases (yellow bars higher than red and blue).This is less significant for LE, where one or more foreign SS sets of parameters give RMS reductions similar to the local SS set of parameters (Fig. 6b).We can deduce that in general the SS sets of parameters are not generic enough to be trans-posed at other sites, and we can link the small difference for LE to the small number of LE-related parameters.Besides, the transposition of SS parameters to other sites also gives hints regarding which sites do not "fit" in the group here studied: for example, the Sorø forest site (DK-Sor) shows a high RMS for NEE whenever optimized parameters from a different site or MS parameters are used.It shows that this site is in some way atypical relative to the other deciduous broadleaf forest sites in our analysis, and Fig. S2 suggests that this is because only at this site both GPP and R eco fluxes are underestimated by the prior model, calling for corrections opposed to the general trend observed throughout this study.Second, the MS optimization generally results in a better RMS reduction than when using the mean set of SS parameters, with averaged RMS differences of 0.097 gC m −2 d −1 (NEE) and 1.83 W m −2 (LE).These discrepancies are larger than the averaged RMS difference between MS and SS optimization (0.056 gC m −2 d −1 and 1.26 W m −2 ).It suggests that the non-linearity of the model implies that a potentially generic set of parameters cannot simply be the average of site-specific values and that the MS optimization brings additional information.In addition, we checked that using an averaged value also for K soilC in the mean parameter set (instead of the local SS value) results in much poorer results (not shown), most likely because of the strong dependence of this parameter on the land-use history at each site.Finally, there are a few sites where the mean set of parameters does better than the MS optimization for one, but not both, fluxes: UK-Ham for NEE, FR-Fon and US-UMB for LE.We can notably notice that UK-Ham is one of the very few sites where the prior GPP is underestimated (Fig. S6), while at FR-Fon and US-UMB the prior LE is either commensurate with observations or even underestimated -as compared to a general overestimation.MS corrections might thus be inconsistent with the local peculiarities, and Fig. 3 shows indeed that, in these three cases, the MS RMS for the corresponding flux is significantly larger than the SS one.We conclude that these sites might not fit in the multi-site group with respect to the mentioned fluxes, so that a site-specific parameterization would be needed here.

Information content of the different observations
In order to estimate the respective contribution NEE and LE to the assimilation, we conducted the MS optimizations using either NEE or LE, but not both, fluxes.The impact on the RMS is shown in Fig. 7, which compares the results of the different cases, with and without each type of data.Regarding the fit to NEE, the performance of the optimization is very similar at most sites whether LE data are used or not, whereas using only LE data results in a significant degradation of the NEE model-data fit from the prior in most cases (averaged RMS increased by 22 %).This degradation stems from LE having no leverage on the modeled R eco , while the assimilation of LE still decreases GPP from the prior model, via the reduced stomatal conductance.This difference leads to higher NEE values than for the already-overestimated prior NEE, thus degrading the fit.Regarding the modeled LE, the situation is almost symmetrical: using or not NEE barely affects the performance of the model after optimization.One can however notice that the assimilation of just NEE does not degrade the fit, and even improves it: only

Evaluation of the optimized model
A crucial step following any optimization procedure is to assess whether the new set of parameters improves the overall model performance, using additional data (i.e., not used in the assimilation).

Optimized GPP and R eco
We first compare the model GPP and R eco fluxes after the optimization with estimates derived directly from the observations (see Sect. 2.5.1).Although not independent, these "data-oriented" estimates provide valuable insights into the model performance.Figure 8 shows the seasonal cycle of GPP and R eco at two of the sites for a 2-year time period (the full period at each site can be found in the Supplement).In general, the optimizations decrease the seasonal amplitude of GPP and R eco at these sites, with a reduction in the growing season length, confirming the patterns shown in Fig. 5 and discussed in Sect.3.2.Besides, we observe that the model-data fit is generally improved both for R eco and GPP, although more significantly for R eco .This is quantified in Fig. S14 where the RMSs for the different cases are shown for both fluxes at different time scales, in the same way as in Fig. 3. Unsurprisingly, the most important RMS reduction happens for R eco on the yearly time scale (an average decrease of 54 % after MS optimization, and 59 % in the SS case), consistent with the changes in initial carbon pools discussed in Sect.3.2.1.More specifically, the amplitudes of the seasonal cycles (see Fig. 8 and S1 to S12) show a strong decrease of R eco at most sites, compared with the overestimation of this flux by the prior model.In some cases, however, the correction after optimization is either too small (DE-Hai, US-Bar, and US-LPH) or too large (DK-Sor).
Regarding GPP, there is also an overestimation by the prior model at most sites (except DK-Sor and UK-Ham): the average annual GPP is equal to 1754 gC m −2 yr −1 , as compared to 1463 gC m −2 yr −1 given by the flux-partitioning estimates.
After MS and SS optimizations, the reduction in carbon assimilation respectively leads to annual GPPs of 1352 and 1479 gC m −2 yr −1 .In-depth comparisons with site-specific gross flux estimates at each site (e.g., Granier et al., 2008) are beyond the scope of this paper but deserve further attention for a more precise evaluation of the optimization procedure at all sites.Overall, our optimization scheme is able to provide a set of parameters that fairly improves the simulation of assimilation and respiration processes in the ORCHIDEE model, although we have chosen to assimilate daily NEE and not to separate between nighttime and daytime values.

Global-scale evaluation: use of MODIS data
One objective of the data assimilation system described in this study is to local information provided by flux tower measurements order to improve continental-to globalscale simulations of the carbon and water balance, with an optimized set of parameters.Such improvement is not guaranteed, and the evaluation of the simulated fluxes at large scales for the DBF ecosystem considered in this study after the MS optimization is thus a necessary step.We use below the information on temporal variations of the vegetation activity retrieved globally by MODIS and not on the absolute values of these measurements (see Sect. 2.5.2).We thus correlate the satellite-derived NDVI time series against those of the Fraction of Absorbed Photosynthetically Active Radiation (FAPAR) modeled by ORCHIDEE.The method has been extensively described in Maignan et al. (2011); only a brief summary and its adaptation to our study are given here.FAPAR is estimated from modeled LAI with a simple Beer's law:

Biogeosciences
Correlations are computed using weekly time series over the period 2000-2008.Note that the NDVI data, originally at the 5-km spatial resolution, were aggregated at the resolution of the vegetation model (determined by the meteorological forcing, here the ERA-Interim, i.e., 0.7 degree).The mesh cells where no clear NDVI annual cycle is visible are ignored in the calculation (i.e., when observed time series have a standard deviation lower than 0.04).The analysis is made only for the pixels with at least 50 % of DBF (the PFT considered in the optimization).We use a high-resolution vegetation map such as CORINE over Europe (Heymann et al., 1993) and UMD (Hansen et al., 2000) to assign a PFT to each satellite pixel, while model boxes are assigned to the PFT whose fraction exceeds 50 % (box unused otherwise).Finally, the spatial averaging of satellite NDVI data is made over pixels where the assigned PFT matches the model box PFT.
Changes in NDVI/FPAR correlations obtained either with the prior model or the MS-optimized model are shown in Fig. 9 for several regions (median value).Note that the boxes used for the calculation are shown in orange on the background map.First one should notice that the correlations with the prior model are relatively high for the three Northern Hemisphere regions (North America, Europe, Asia), with values over 0.88.On the other hand, in the Southern Hemisphere, the prior model performances are much lower, with correlations below 0.5.If we now consider the changes brought by the MS optimization, we find that the correlations are improved everywhere except in Oceania.In the Northern Hemisphere, the MS optimization improves the phase correlation factor between NDVI and FPAR by nearly 4 %, with posterior correlations above 0.91.In the southern areas, the MS optimization leads to contrasting effects.While it does not affect much South America (+7 %), there is a significant improvement in Africa (+36 %) and a strong degradation in Oceania (-28 %).Globally, we observe an improvement brought by the MS optimization with a median correlation factor going from 0.83 to 0.88 (not shown).As mentioned in Maignan et al. (2011), there is no expectation of a rigorous correspondence between NDVI and FPAR temporal variations because NDVI is impacted by other variables than FPAR (such as vegetation geometry, soil reflectance, fractional cover, mixture of grass understory with trees, measurement noise, or leaf spectral signature), so that a perfect correlation should not be taken as a target.The much poorer results in the Southern Hemisphere could either reveal a problem of the phenology model (purely temperature-dependent for the DBF PFT in ORCHIDEE), or simply point out the limit of using a single generic temperate DBF classification at the global scale.

Conclusions
In the framework of eddy-covariance flux data assimilation in biosphere models, we have sought to expand the geographically limited approach of site-scaled model optimizations.To this end, we have built a data assimilation system able to simultaneously integrate the information given by several measurements sites in order to derive a unique set of optimized model parameters.This so-called multi-site (MS) optimization procedure has been here focused on one type of ecosystem -the temperate deciduous broadleaved forests (DBF).
The MS optimization is able to provide daily model-data RMS reductions (with respect to the prior model) that are often as good as the single-site (SS) optimizations.This consistency is also true at yearly time scale where the NEE misfit is reduced by half for both the single-site and multi-site cases.The major contribution to the yearly improvement is the scaling of the initial carbon pools, governed by the only parameter purposely kept site-specific in the MS case.This scaling is crucial because of the discrepancy between the state of the modeled ecosystem after the initial spin-up procedure, corresponding to a mature ecosystem, and the young forests mostly used in this study (with lower soil carbon content).Note however that this first-order correction remains linked to the other respiration parameters, and overall the www.biogeosciences.net/9/3757/2012/Biogeosciences, 9, 3757-3776, 2012 assimilation of net carbon fluxes does not allow to fully separate between pool size and turnover rate effects in the calculation of the respiratory fluxes (see error correlations in Fig. S13).Additional measurements of soil carbon pool content could be used in the future to obtain cross constraints on all factors controlling the heterotrophic respiration.The autotrophic respiration (R a ) is also generally reduced after optimization.Using total ecosystem respiration derived from NEE partitioning techniques, we validated the reduction in model respiration after optimization.However, we could not assess whether the autotrophic reduction may compensate for an insufficiently reduced heterotrophic component.
In parallel, the carbon assimilation is slightly reduced at most sites following optimization.Comparisons with estimations of GPP derived from NEE indicate that this correction is somewhat relevant, but the summer carbon uptake remains underestimated at half of the sites after optimization, suggesting model structural errors.The latter could be related to the fact that for example the temperature and the soil-water control on the photosynthesis are still simply parameterized, and that no biotic effect is taken into account.Besides, we have observed that the growing season is consistently shortened as compared to the prior overestimation, but remains too long at most sites.Similarly, in an evaluation of more than a dozen different terrestrial biosphere models, Richardson et al. (2012) found that most models tended to overestimate the growing season length at five North American deciduous broadleaf forests, resulting in misrepresentation of the seasonality of leaf area index as well as photosynthetic uptake.Furthermore, most models could not successfully predict interannual variability in spring or autumn phenology either.
Regarding the water cycle, the prior model generally overestimates the latent heat flux, and both MS and SS optimizations generally improve the model-data fit with a reduced stomatal conductance and a larger vegetation albedo.This result should however be tempered given the lack of energy balance closure at the sites, a potential source of bias in the measurement values of LE.Additionally, we observed highly fluctuating modeled values of LE in winter and spring at some sites in contradiction with observations, most likely caused by an inconsistency in the simulation/parameterization of the snow sublimation in OR-CHIDEE, as the parameters associated with sublimation were not optimized here.
In general, we have observed comparable parameter changes between SS and MS optimizations, even if the former provides (not surprisingly) somewhat better model data fits.Performing MS optimization is more complex in terms of optimization code, and we thus investigated its actual benefit as compared to a series of SS optimizations.It turns out that, for more than half of the parameters, the MS values are close to the means of the SS values.This was not the case for all the parameters we optimized, presumably because of nonlinearities in the model structure.Overall the MS set of parameters provides a significantly better fit than the mean of the SS sets, while using of a single SS set at others sites is most likely to degrade the performances of the model.This consistency between individual and grouped optimizations contrasts with the results of Groenendijk et al. (2011) which showed a significant degradation of the model-data agreement when changing from site-specific to PFT-generic parameters, in an analysis carried out with a different model and optimization scheme.Besides, we acknowledge that the optimized parameter sets might still correspond to local minima.Ensemble methods could possibly provide a better modeldata fit, but the relatively large number of optimized parameters makes our variational method much more affordable in terms of computational time.
Regarding the assimilated data, we observed that the amount of information brought by both fluxes (NEE and LE) is complementary: it greatly improves the fit to each of the two outputs almost independently.
Finally, we have evaluated the performance of the phenology in the ORCHIDEE model at the global scale via the correlation coefficient between modeled FAPAR and measured NDVI, restricted to the DBF ecosystem.Our analysis shows that the prior model does far better in the Northern than in the Southern Hemisphere.From this starting point, the MS optimization brings a slight improvement in the Northern Hemisphere, and contrasting results in the Southern Hemisphere, where none of the sites used in the optimization is located: significant improvement in South Africa but degradation in Australia.At the global scale, the correlation median shifts from 0.83 to 0.88.The degradation in Australia might reflect the limits of the phenological scheme of deciduous forests in the model, solely based on a temperature criterion.The tree species in the arid Australian forests, although classified as DBF, are likely to have a phenology strongly controlled by the available soil moisture, a feature much less prevalent at the sites used in this study.At present, we can only suggest the need for further investigations regarding the formulation of the DBF phenology in the model towards a refinement of the PFT classification.
This work should be considered as a guideline for the assimilation of eddy-covariance data at several sites, within carbon cycle data assimilation system (CCDAS) (Rayner et al., 2005) together with complementary data streams.It emphasizes the strengths but also the limitations of using sitespecific daily NEE.Although the temporal variations of NEE bring specific information on the gross fluxes, assimilating additional data would help to better distinguish between processes, and potentially better constraint sensitivities to environmental drivers.The separation between daytime and nighttime NEE, as well as data on within-tree carbon allocation, leaf area index and/or phenology, soil respiration fluxes, biomass and litter/soil C pools, could also be used as model constraints, where such data are available.A multipleconstraint approach can be used to improve the internal dynamics of the modeled system, and reduce bias and model error when the model is used for extrapolation or projection in time or space (Richardson et al., 2010).We also plan to extend our method to other types of ecosystems and to expand the evaluation process using additional independent data at the global scale.For example, the seasonal cycle contained in the atmospheric CO 2 concentration reflects primarily that of the terrestrial biosphere and could be used through an atmospheric transport model.Simulated LE flux can be spatially integrated and evaluated against measured stream flows for large basins.Finally, inter-annual variations of the model phenology could be compared for instance to more than 20 years of NDVI data from the AVHRR instrument.

Fig. 2 .
Fig. 2. Seasonal cycle of NEE and LE at (a) Hainich and (b) Harvard Forest sites, smoothed with a 15-day moving average window.The observations (black) are compared with the prior model (grey), the MS optimization (blue) and the SS optimization (orange).The error bars give the total flux uncertainties.The annual carbon budget (gC m −2 yr −1 ) and the annually averaged LE flux (W m −2 ) are given between brackets.

Fig. 3 .
Fig. 3. Model-data RMSs at different time scales for (a) NEE and (b) LE.Prior model is in grey, MS optimization in blue and SS optimization in orange.

Fig. 5 .
Fig. 5. Starting day, ending day, and length of the growing season (average over all the years available at each site).Measurements are in black, prior model in grey, MS-optimized model in blue, and SS-optimized model in orange.

Fig. 6 .Fig. 7 .
Fig. 6.Daily model-data RMSs for (a) NEE and (b) LE.For each site, one bar corresponds to the results from one set of the 21 parameters: prior model (grey), parameters from SS optimizations at the other sites (yellow), local SS optimization in orange, MS optimization (blue), mean of the SS sets of parameters (black).

Fig. 8 .
Fig. 8. Seasonal cycle of GPP and R eco at (a) Hainich and (b) Harvard Forest sites, smoothed with a 15-day moving average window.The estimations derived from flux-partitioning of NEE (black) are compared with the prior model (grey), the MS optimization (blue) and the SS optimization (orange).The averaged annual fluxes in gC m −2 yr −1 are given between brackets.

Fig. 9 .
Fig. 9.Continental medians of NDVI/FPAR correlation of prior model (grey) and after MS optimization (blue), using weekly time series for the 2000-2008 period and the ERA-I simulation.Correlations are only calculated for boxes with dominant DBF ecosystem and where cycles in NDVI and FPAR are detected (orange boxes on the map).

Table 1 .
List of the sites.

Table 2 .
ORCHIDEE parameters optimized in this study.

Table 3 .
All-site yearly reduction of the model-data NEE RMS from the prior model, for various combinations of R h parameters left out the optimization.
R h parameter(s) left out None K soilC Q 10 HR b HR c Q 10 +HR H,b +HR H,c All