Mirco Migliavacca
email: m.migliavacca1@campus.unimib.it

Markus ; Reichstein

Andrew D Richardson

Roberto Colombo

Mark A ; Sutton

Gitta Lasslop

Georg Wohlfahrt

Enrico Tomelleri

Nuno Carvalhais

Alessandro Cescatti

Semi-empirical modeling of abiotic and biotic factors controlling ecosystem respiration across eddy covariance sites

Keywords: Ecosystem Respiration, Productivity, FLUXNET, Eddy Covariance, Leaf Area Index, Inverse Modeling Ecosystem Respiration, Productivity, FLUXNET, Eddy Covariance, Leaf Area 92 Index, Inverse Modeling 93 94

come

In this study we examined ecosystem respiration (RECO) data from 104 sites belonging to FLUXNET, the global network of eddy covariance flux measurements. The main goal was to identify the main factors involved in the variability of RECO: temporally and between sites as affected by climate, vegetation structure and plant functional type (PFT) (evergreen needleleaf, grasslands, etc.).

We demonstrated that a model using only climate drivers as predictors of RECO failed to describe part of the temporal variability in the data and that the dependency on gross primary production (GPP) needed to be included as an additional driver of RECO. The maximum seasonal leaf area index (LAIMAX) had an additional effect that explained the spatial variability of reference respiration (the respiration at reference temperature Tref=15°C, without stimulation introduced by photosynthetic activity and without water limitations), with a statistically significant linear relationship (r2=0.52 p<0.001, n=104) even within each PFT. Besides LAIMAX, we found that the reference respiration may be explained partially by total soil carbon content. For undisturbed temperate and boreal forest a negative control of the total nitrogen deposition on the reference respiration was also identified.

We developed a new semi-empirical model incorporating abiotic factors (climate), recent productivity (daily GPP), general site productivity and canopy structure (LAIMAX) which performed well in predicting the spatio-temporal variability of RECO, explaining >70% of the variance for most vegetation types. Exceptions include tropical and Mediterranean broadleaf forests and deciduous interactions (e.g. [START_REF] Cox | Acceleration of global warming 851 due to carbon-cycle feedbacks in a coupled climate model[END_REF][START_REF] Houghton | Missing sinks, feedbacks, and understanding 946 the role of terrestrial ecosystems in the global carbon balance[END_REF]. It has been hypothesized that relatively 99 F o r R e v i e w O n l y 4 small climatic changes may impact respiration with the effect of rivalling the annual fossil fuel 100 loading of atmospheric CO 2 [START_REF] Jenkinson | Model Estimates of Co2 Emissions from Soil in 967 Response to Global Warming[END_REF][START_REF] Raich | The global carbon dioxide flux in soil respiration and its 1082 relationship to vegetation and climate[END_REF]. 101

Recently, efforts have been made to mechanistically understand how temperature and other 102 environmental factors affect ecosystem and soil respiration, and various modeling approaches have 103 been proposed (e.g. Davidson et al., 2006a[START_REF] Lloyd | On the temperature dependence of soil respiration[END_REF][START_REF] Reichstein | Soil respiration across scales: The importance of a model-data 1090 integration framework for data interpretation[END_REF], 104 Reichstein et al., 2003a). Nevertheless, the description of the conceptual processes and the complex 105 interactions controlling R ECO are still under intense research and this uncertainty is still hampering 106 bottom-up scaling to larger spatial scales (e.g. regional and continental) which is one of the major 107 challenges for biogeochemists and climatologists. 108

Heterotrophic and autotrophic respiration in both data-oriented and process-based 109 biogeochemical models are usually described as a function of air or soil temperature and 110 occasionally soil water content (e.g. [START_REF] Lloyd | On the temperature dependence of soil respiration[END_REF], Reichstein et al., 2005, Thornton et al., 111 2002), although the functional form of these relationships varies from model to model. These 112 functions represent the dominant role of reaction kinetics, possibly modulated or confounded by 113 other environmental factors such as soil water content or precipitation, which some model 114 formulations include as a secondary effect (e.g. [START_REF] Carlyle | Abiotic controls of soil respiration beneath an eighteen-year-old 839 Pinus radiata stand in south-eastern Australia[END_REF], Reichstein et al., 2003a, 115 Richardson et al., 2006). 116

A large number of statistical, climate-driven models of ecosystem and soil respiration have been 117 tested and compared using data from individual sites (Del [START_REF] Grosso | 874 Modeling soil CO2 emissions from ecosystems[END_REF], Janssens & 118 Pilegaard, 2003[START_REF] Richardson | Statistical modeling of ecosystem respiration using eddy 1115 covariance data: Maximum likelihood parameter estimation, and Monte Carlo simulation of 1116 model and parameter uncertainty, applied to three simple models[END_REF][START_REF] Savage | Three scales of temporal 1129 resolution from automated soil respiration measurements[END_REF], multiple sites (Falge et al., 119 2001[START_REF] Rodeghiero | Main determinants of forest soil respiration along an 1119 elevation/temperature gradient in the Italian Alps[END_REF], and from a wide range of models compared across different 120 ecosystem types and measurement techniques [START_REF] Richardson | Comparing simple respiration models 1112 for eddy flux and dynamic chamber data[END_REF]. 121

Over the course of the last decades, the scientific community has debated the role of productivity 122 in determining ecosystem and soil respiration. Several authors [START_REF] Bahn | Soil Respiration in European Grasslands 808 in Relation to Climate and Assimilate Supply[END_REF], Curiel Yuste et 123 al., 2004, Davidson et al., 2006a[START_REF] Janssens | Productivity overshadows temperature in 958 determining soil and ecosystem respiration across European forests[END_REF], Reichstein et al., 2003a, Valentini et al., 124 2000) have discussed and clarified the role of photosynthetic activity, vegetation productivity and 125 their relationship with respiration. 126 Linking photosynthesis and respiration might be of particular relevance when modelling R ECO 127 across biomes or at the global scale. Empirical evidence for the link between GPP and R ECO is 128 reported for most, if not all, ecosystems: grassland (e.g. [START_REF] Bahn | Soil Respiration in European Grasslands 808 in Relation to Climate and Assimilate Supply[END_REF][START_REF] Bahn | Does photosynthesis affect 810 grassland soil-respired CO2 and its carbon isotope composition on a diurnal timescale?[END_REF], Craine 129 et al., 1999[START_REF] Hungate | Evapotranspiration and soil water content in a 952 scrub-oak woodland under carbon dioxide enrichment[END_REF], crops (e.g. [START_REF] Kuzyakov | Photosynthesis controls of rhizosphere respiration and organic 983 matter decomposition[END_REF][START_REF] Moyano | Response of mycorrhizal, rhizosphere and soil basal 1044 respiration to temperature and photosynthesis in a barley field[END_REF], 130 boreal forests [START_REF] Gaumont-Guay | Biophysical controls on 899 rhizospheric and heterotrophic components of soil respiration in a boreal black spruce stand[END_REF][START_REF] Hogberg | Large-scale forest girdling shows that current 939 photosynthesis drives soil respiration[END_REF] and temperate forests, both 131 deciduous (e.g. [START_REF] Curiel-Yuste | Annual Q 10 of soil respiration reflects 855 plant phenological patterns as well as temperature sensitivity[END_REF], Liu et al., 2006) and evergreen (e.g. [START_REF] Irvine | Coupling of canopy gas exchange with root and rhizosphere 954 respiration in a semi-arid forest[END_REF]. 132 This time lag depends to the vegetation structure it is related to the translocation time of assimilates 134 from aboveground to belowground organs through the phloem. Although the existence of a time lag 135 is still under debate, it has been found to be a few hours in grasslands, and croplands and a few 136 days in forests [START_REF] Baldocchi | How switches and lags in biophysical regulators affect 821 spatial-temporal variation of soil respiration in an oak-grass savanna[END_REF][START_REF] Knohl | Partitioning the net CO2 flux of a deciduous forest into respiration 976 and assimilation using stable carbon isotopes[END_REF][START_REF] Moyano | Soil respiration fluxes in relation to photosynthetic 1041 activity in broad-leaf and needle-leaf forest stands[END_REF], Savage et 137 al., 2009). 138

While the link between productivity and respiration appears to be clear, to our knowledge, few 139 model formulations include the effect of productivity or photosynthesis as a biotic driver of 140 respiration and these models are mainly developed for the simulation of soil respiration using a 141 relatively small data set of soil respiration measurements (e.g. [START_REF] Hibbard | An analysis of soil respiration across 934 northern hemisphere temperate ecosystems[END_REF], Reichstein et 142 al., 2003a). 143 In this context, the increasing availability of ecosystem carbon, water and energy flux 144 measurements collected by means of the eddy covariance technique (e.g. [START_REF] Baldocchi | Breathing of the terrestrial biosphere: lessons learned from a global network of 816 carbon dioxide flux measurement systems[END_REF] over 145 different plant functional types (PFTs) at more than 400 research sites, represents an useful tool for 146 understanding processes and interactions behind carbon fluxes and ecosystem respiration. These 147 data serve as a backbone for bottom-up estimates of continental carbon balance components (e.g. 148 [START_REF] Ciais | Europe-wide reduction in primary productivity caused 844 by the heat and drought in 2003[END_REF][START_REF] Papale | A new assessment of European forests carbon exchanges by eddy 1061 fluxes and artificial neural network spatialization[END_REF], Reichstein et al., 2007) and for ecosystem model 149 development, calibration and validation (e.g. [START_REF] Baldocchi | Measuring and modelling carbon dioxide and water vapour exchange over a 813 temperate broad-leaved forest during the 1995 summer drought[END_REF][START_REF] Hanson | Oak forest carbon and water simulations: 929 Model intercomparisons and evaluations against independent data[END_REF], Law et al., 150 2000[START_REF] Owen | Linking flux network measurements to 1055 continental scale simulations: ecosystem carbon dioxide exchange capacity under non-1056 water-stressed conditions[END_REF], Reichstein et al., 2003b, Reichstein et al., 2002[START_REF] Verbeeck | Parameter sensitivity and uncertainty of the 1180 forest carbon flux model FORUG: a Monte Carlo analysis[END_REF]. 151

The database includes a number of added products such as gap-filled net ecosystem exchange 152 (NEE), gross primary productivity (GPP), ecosystem respiration (R ECO ) and meteorological drivers 153 (air temperature, radiation, precipitation etc) aggregated at different time-scale (e.g. half-hourly, 154 daily, annual) and consistent for data treatment [START_REF] Papale | Towards a standardized processing of Net 1058 Ecosystem Exchange measured with eddy covariance technique: algorithms and uncertainty 1059 estimation[END_REF], Reichstei et al., 2005) 155

In this paper we analyze with a semi-empirical modeling approach the R ECO at 104 different sites 156 belonging to the FLUXNET database with the primary objective of synthesizing and identifying the 157 main factors controlling i) the temporal variability of R ECO , ii) the between-site (spatial) variability 158 and iii) to provide a model which can be used for diagnostic up-scaling of R ECO from eddy 159 covariance flux sites to large spatial scales. 160 Specifically, the analysis and the model development followed these two steps: 161 1. we developed a semi-empirical R ECO model site by site (site-by-site analysis) with the aim of 162 clarifying if and how GPP should be included into a model for improving the description of 163 R ECO and which factors are best suited for describing the spatial variability of reference 164 respiration (i.e. the daily R ECO at the reference temperature without moisture limitations). 165

We follow these three steps: 166

F o r R e v i e w O n l y
o the analysis of R ECO data was conducted by using a purely climate driven model: 'TP 167

Model' [START_REF] Raich | Interannual variability in global soil respiration, 1980-1080 94[END_REF]. The accuracy of the model and the main bias were 168 The data used in this analysis is based on the dataset from the FLUXNET (www.fluxdata.org) 187 eddy covariance network [START_REF] Baldocchi | Breathing of the terrestrial biosphere: lessons learned from a global network of 816 carbon dioxide flux measurement systems[END_REF], Baldocchi et al., 2001). The analysis was restricted to 188 104 sites (cf. Table in Appendix I and II) on the basis of the ancillary data availability (i.e. only 189 sites containing at least both leaf area index (LAI) of understorey and overstorey were selected) and 190 of the time series length (all sites containing at least one year of carbon fluxes and meteorological 191 data of good quality data were used). Further, we only analyzed those sites for which the relative 192 standard error of the estimates of the model parameters E 0 (activation energy) and reference 193 respiration (R 0 ) (please see further sections for more details on the meaning of parameters) were 194 less than 50% and where E 0 estimates were within an acceptable range (0-450 K). 195

The latitude spans from 71.32° at the Alaska Barrow site (US-Brw) to -21.62° at the Sao Paulo 196 Cerrado (BR-Sp1). The climatic regions include tropical to arctic. 197 All the main PFTs as defined by the IGBP (International Geosphere-Biosphere Programme) 198 were included in this study: the selected sites included 28 evergreen needleleaf forests (ENF), 17 199 At each site data are storage corrected, spike filtered, u * -filtered according to Papale et al. (2006) 209 and subsequently gap-filled and partitioned as described by [START_REF] Reichstein | On the separation of net ecosystem exchange into 1096 assimilation and ecosystem respiration: review and improved algorithm[END_REF]. Only days 210 containing both meteorological and daily flux data with a percentage of gap-filled half hours below 211 15% were used for this analysis. The median of the u * threshold applied in the FLUXNET database 212 for the site-years used in the analysis are listed in the Appendix II. The average of the median u * 213 values are lower for short canopies (e.g. for grasslands 0.075±0.047 ms -1 ) and higher for tall 214 canopies (e.g. for evergreen needleleaf forests 0.221 ±0.115 ms -1 ). 215

Along with fluxes and meteorological data, main ancillary data such as maximum ecosystem 216 LAI (overstory and understory for forest sites) (LAI MAX ), LAI of overstory (LAI MAX,o ), stand age 217 for forests (StandAge), total soil carbon stock (SoilC) and the main information about disturbance 218 (date of cuts, harvesting) were also downloaded from the database. Total atmospheric nitrogen 219 deposition (N depo ) is based on the atmospheric chemistry transport model TM3 [START_REF] Rodhe | The global distribution of acidifying wet deposition[END_REF] proposed by [START_REF] Raich | Interannual variability in global soil respiration, 1980-1080 94[END_REF] and further modified by (Reichstein et al., 2003a). Here we used the 230 'TP Model' for the simulation of R ECO at the daily time-step using as abiotic drivers daily T A and P: 231

232 ) ( ) ( P f T f R R A ref ECO ⋅ ⋅ = (
        - - - = 0 0 0 1 1 ) ( T T T T E A A ref e T f
(2) 242 243

We refine the approach of Reichstein et al. (2003) and propose a reformulation of the response 244 of R ECO to precipitation (Eq. 3), where k (mm) is the half saturation constant of the hyperbolic 245 relationship and α is the response of R ECO to null P. 246 247

( ) ( ) α α α - + - + = 1 1 ) ( P k P k P f (3) 248 249
Although soil water content is widely recognized as the best descriptor of soil water availability, 250 we preferred to use precipitation since the model developed is oriented to up-scaling and soil water 251 maps are more affected by uncertainty than precipitation maps. 252

The model parameters -R REF , E 0 , α, k -were estimated for each site in order to evaluate the 253 accuracy of the climate-driven model. At each site the Pearson's correlation coefficient (r) between 254 'TP Model' residuals (R ECO observed minus R ECO modelled ) and GPP was also computed. 255 256

Site-by-site analysis -Effect of productivity on the temporal variability of R ECO 257 258

The role of GPP, as an additional biotic driver of R ECO that has been included into Eq. 1, was 259 analysed at each site using three different formulations of the dependency of ecosystem respiration 260 on productivity f(GPP): 261 Linear response: According to different authors (e.g. [START_REF] Hibbard | An analysis of soil respiration across 934 northern hemisphere temperate ecosystems[END_REF], Reichstein et al., 2007) we hypothesized 266 that respiration might saturate at high productivity rates in a similar way to the Michaelis-Menten 267 enzyme kinetics. This saturation can also occur by a transition of carbon limitation to other 268 limitations. The exponential curve was used as another formulation of a saturation effect. 269

GPP k GPP f ⋅ = 2 ) ( (4) 262 Exponential response: ( ) GPP k e R GPP f ⋅ - - ⋅ = 2 1 ) ( 2 (5) 263 Michaelis-Menten: GPP h GPP R GPP f R + ⋅ = max max ) ( (6) 264
We tested two different schemes for the inclusion of f(GPP) (Eqs. 4,5,6) in the 'TP Model' 270 (Eq.1): 271 272 1) f(GPP) was included by replacing the reference respiration at reference temperature 273 (R ref in Eq. 1) with the sum of a new reference respiration (R 0 ) and the f(GPP): 274 ( )

GPP f R R ref + = 0 (7) 275
2) f(GPP) was included as an additive effect into the 'TP Model'. In this case one part 276 of ecosystem respiration is purely driven by biotic factors (e.g. independent from 277 temperature) and the other one by abiotic ones. 278

279

In Table 1, R 0 is the new reference respiration term (i.e. ecosystem respiration at T ref , when the 280 GPP is null and the ecosystem is well watered). This quantity is considered to be an indicator of the 281 ecosystem respiration of the site, strictly related to site conditions, history and characteristics, while 282 k 2 , R 2 , R max and h Rmax describe the assumed functional response to GPP. 283

[TABLE1] 285 286

The model parameters -R 0 , E 0 , α, k and the parameters of ƒ(GPP) -were estimated for each site 287 in order to evaluate which model formulation best describes the temporal variability of R ECO . 288 With the aim of confirming the existence of a time lag between photosynthesis and the 289 respiration response we ran the model with different time lagged GPP time-series (GPP lag,i ), starting 290 from the GPP estimated on the same day (GPP lag,0 ), and considering daily increments back to GPP 291 estimated one week before the measured R ECO (GPP lag,7 ). 292 GPP and R ECO estimated with the partitioning method used in the FLUXNET database are 293 derived from the same data (i.e. GPP=R ECO -NEE) and this may to some extent introduce spurious 294 correlation between these two variables. In literature two different positions on that can be found: 295 [START_REF] Vickers | Self-correlation between assimilation and 1186 respiration resulting from flux partitioning of eddy-covariance CO2 fluxes[END_REF] argue that there is a spurious correlation between GPP and R ECO when these 296 component fluxes are jointly estimated from the measured NEE (i.e. as estimated in the FLUXNET 297 database). [START_REF] Lasslop | Comment on Vickers et 988 al.: Self-correlation between assimilation and respiration resulting from flux partitioning of 989 eddy-covariance CO2 fluxes[END_REF] and R ECO-LASS ). The method by [START_REF] Lasslop | Separation of net ecosystem exchange into 992 assimilation and respiration using a light response curve approach: critical issues and global 993 evaluation[END_REF] do not compute GPP as a difference, but 302 derive R ECO and GPP from quasi-disjoint NEE data subsets. Hence, if existing, spurious correlations 303 is minimized. 304

To understand whether our results are affected or not by the 'spurious' correlation between GPP 305 and R ECO estimated in FLUXNET, we also performed the analysis using the GPP and R ECO 306 estimated by the partitioning method of [START_REF] Lasslop | Separation of net ecosystem exchange into 992 assimilation and respiration using a light response curve approach: critical issues and global 993 evaluation[END_REF]. The details of the analysis are 307 described in the Appendix IV. The results obtained confirmed (Appendix IV) that the data 308 presented and discussed in follow are not influenced by the possible 'spurious' correlation between 309 R ECO and GPP reported in the FLUXNET data set. 310 311 Site-by-site analysis -Spatial variability of reference respiration (R 0 ) 312 313 Once the best model formulation was defined, we analyzed the site-by-site (i.e. spatial) 314 variability of R 0 : the relationships between the estimated R 0 at each site and site-specific ancillary 315 data were tested, including LAI MAX , LAI MAX,o , N depo , SoilC and Age. Leaf mass per unit area and 316 aboveground biomass were not considered because these are rarely reported in the database for the 317 sites studied and poorly correlated with spatial variability of soil respiration, as reported by 318 Reichstein et al. (2003a). In this analysis the sites with incomplete site characteristics were removed 319 (Age was considered only for the analysis of forest ecosystems). On the basis of this analysis the 320 model was reformulated by adding the explicit dependency of R 0 on the site characteristics that best 321 explained its variability. 322

323

PFT-Analysis 324 325

In this phase we tried to generalize the model parameters in order to obtain a parameterization 326 useful for diagnostic PFT-based up-scaling. For this reason model parameters were estimated 327 including all the sites for each PFT at the same time. The dependency of R 0 was prescribed as a 328 function of site characteristics that best explain the spatial R 0 variability within each PFT class. 329

The model was corroborated with two different cross-validation methods: 330 analysis that optimize model parameters finding the minimum of a defined cost function. The cost 352 function used here is the sum of squared residuals weighted for the uncertainty of the observation 353 (e.g. [START_REF] Richardson | Statistical modeling of ecosystem respiration using eddy 1115 covariance data: Maximum likelihood parameter estimation, and Monte Carlo simulation of 1116 model and parameter uncertainty, applied to three simple models[END_REF]. The uncertainty used here is an an estimate of the random error 354 associated with the night-time fluxes (from which R ECO is derived). 355

Model parameter standard errors were estimated using a bootstrapping algorithm with N=500 356 random re-sampling with replacement of the dataset. As described by [START_REF] Efron | An Introduction to the Bootstrap[END_REF], 357 the distribution of parameter estimates obtained provided an estimate of the distribution of the true 358 model parameters. 359 360

Best model formulation selection 361 362

For the selection of the 'best' model from among the six different formulations listed in Table 1 363 and the 'TP Model' we used the approach of the information criterion developed by Akaike (1973) 364 which is considered a useful metric for model selection (Anderson et al., 2000, Richardson et al., 365 2006). In this study the Consistent Akaike Information Criterion (cAIC, eq. 8) was preferred to the 366 AIC because the latter is biased with large datasets [START_REF] Shono | Is model selection using Akaike's information criterion appropriate for catch per 1135 unit effort standardization in large samples[END_REF] tending to select more 367 complicated models (e.g. many explanatory variables exist in regression analysis): 368 369

F o r R e v i e w O n l y 12
( ) ( ) [ ] 1 log log 2 + + Θ - = n p L cAIC (8) 370 371
where L(Θ) is the within samples residual sum of squares, p is the number of unknown parameters 372 and n is the number of data (i.e. sample size). Essentially, when the dimension of the data set is 373 fixed, cAIC is a measure of the trade-off between the goodness of fit (model explanatory power) 374 and model complexity (number of parameters), thus cAIC selects against models with an excessive 375 number of parameters. Given a data set, several competing models (e.g different model 376 formulations proposed in Table 1) can be ranked according to their cAIC, with the formulation 377 having the lowest cAIC being considered the best according to this approach. 378

For the selection of the best set of predictive variables of R 0 we used the stepwise AIC, a 379 multiple regression method for variable selection based on the AIC criterion (Venables & Ripley, 380 2002[START_REF] Yamashita | A stepwise AIC method for variable selection in 1205 linear regression[END_REF]. The stepwise AIC was preferred to other stepwise methods for 381 variable selection since can be applied to non normally distributed data [START_REF] Yamashita | A stepwise AIC method for variable selection in 1205 linear regression[END_REF]. 382 383

Evaluation of model accuracy 384 385

Model accuracy was evaluated by means of different statistics according to Janssen and 386 Heuberger (1995): RMSE (Root Mean Square Error), EF (modelling efficiency), determination 387 coefficient (r 2 ) and MAE (Mean Absolute Error). In particular EF is a measure of the coincidence 388 between observed and modelled data and it is sensitive to systematic deviation between model and 389 observations. EF can range from -∞ to 1. An EF of 1 corresponds to a perfect agreement between 390 model and observation. An EF of 0 (EF = 0) indicates that the model is as accurate as the mean of 391 the observed data, whereas a negative EF means that observed mean is a better predictor than the 392 model. In the PFT-analysis for each validation set the cross-validated statistics were calculated. The 393 average of cross-validated statistics were calculated for each PFT both for the training/evaluation 394 splitting (EF cv , RMSE cv , r 2 cv ) and for the k-fold cross-validation (EF kfold-cv , RMSE kfold-cv , r 2 kfold-cv ). 395 The RMSE and EF obtained with 'TP Model' fitting (Table 2) showed a within-PFT-average EF 403 ranging from 0.38 for SAV to 0.71 for ENF and an RMSE ranging from 0.67 for SHB to 1.55 gC 404 m -2 d -1 for CRO. 405 406

[TABLE 2] 407 408
The importance of productivity is highlighted by residual analysis. A significant positive 409 correlation between the 'TP Model' residuals (z) and the GPP was observed with a systematic 410 underestimation of respiration when the photosynthesis (i.e. GPP) was intense. 411

In Fig. 1a, the mean r between the residuals and GPP for each PFT as a function of the time lag 412 is summarised. 413

The lowest correlation was observed for wetlands (r=0.29±0.14). The mean r is higher for 414 herbaceous ecosystems such as grasslands and croplands (0.55±0.11 and 0.63±0.18, respectively) 415 than for forest ecosystems (ENF, DBF, MF, EBF) which behaved in the same way (Fig. 1a), with a 416 r ranging from 0.35±0.13 for ENF to 0.45±0.13 for EBF. No time lag was observed with the 417 residuals analysis. 418 419

Gross Primary Production as driver of R ECO 420 421

The effect of GPP as an additional driver of R ECO was analyzed at each site by testing 6 different 422 models with the three different functional responses (Eqs. 4, 5 and 6) of respiration to GPP (Tab. 1). 423

The model ranking based on the cAIC calculated for each different model formulation at each site 424 showed agreement in considering the models using the linear dependency of R ECO on GPP 425 ('LinGPP') as the best model formulation (Tab. 2), since the cAICs obtained with 'LinGPP' were 426 lower than those obtained with all the other formulations. This model ranking was also maintained 427 when analysing each PFT separately, except for croplands in which the 'addLinGPP' formulation 428 provided the minimum cAIC although the difference between the average cAIC estimated for the 429 'addLinGPP' and 'LinGPP', respectively) and the standard errors of parameter estimates were 431 lower for the 'LinGPP' formulation. In general, the cAIC obtained at all sites with the 'LinGPP ' 432 model formulation (39.50 [37.50 -42.22], in squared parentheses the first and third quartile are 433 reported) were lower than the ones obtained with the 'TP Model ' (41.08 [39.02 -44.40]), although the 434 complexity of the latter is lower (one parameter less). On this basis we considered the 'LinGPP' as 435 the best one model formulation. 436

F o r R e v i
The statistics of model fitting obtained with the 'LinGPP' model formulation are reported in 437

Table 2. The model optimized site by site showed a within-PFT-average of EF between 0.58 for 438 EBF to 0.85 for WET with an RMSE ranging from 0.53 for SAV to 1.01 gC m -2 day -1 for CRO. On 439 average EF was higher than 0.65 for all the PFTs except for EBF. In terms of improvement of 440 statistics, the use of 'LinGPP' in the 'TP Model' led to a reduction of the RMSE from 13.4 % for 441 shrublands to almost one third for croplands (34.8%), grasslands (32.5%) and savanna (32.0%) with 442 respect to the statistics corresponding to the purely climate driven 'TP Model'. 443 444

[FIGURE 1] 445 446
No time lag between photosynthesis and respiration response was detected. In fact using GPP lag,-i 447 as a model driver we observed a general decrease in mean model performances for each PFT (i.e. 448 decrease of EF and increase of RMSE) for increasing i values (i.e. number of days in which the 449 GPP was observed before the observed R ECO ). The only exception were DBFs in which we found a 450 time lag between the GPP and R ECO response of 3 days as shown by the peak in average EF and by 451 the minimum in RMSE in Fig. 1b, although the differences were not statistically significant. 452 453

Spatial variability of reference respiration rates 454 455

The reference respiration rates (R 0 ) estimated site by site with the 'LinGPP' model formulation 456 represent the daily ecosystem respiration at each the site at a given temperature (i.e. 15°C), without 457 water limitation and carbon assimilation. Hence, R 0 can be consider the respiratory potential of a 458 particular site. R 0 assumed highest values for the ENF (3.01±1.35 gC m -2 day -1 ) while the lowest 459 values were found for SHB (1.49±0.82 gC m -2 day -1 ) and WET (1.11 ±0.17 gC m -2 day -1 ), possibly 460 reflecting lower carbon pools for shrublands or lower decomposition rates due to anoxic conditions 461 or carbon stabilization for wetlands. 462

F o r R e v i e w O n l y
By testing the pairwise relationship between R 0 and different site characteristics we found that 463 the ecosystem LAI MAX showed the closest correlation with R 0 (R 0 =0.44(0.04)LAI MAX +0.78(0.18), 464 r 2 =0.52, p<0.001, n=104, in parentheses standard errors of model parameters estimates were 465 reported), thus LAI MAX was the best explanatory variable of the retrieved R 0 variability (Fig 2a). 466

Conversely, LAI MAX,o correlated weakly (r 2 =0.40, p<0.001, n=104) with R 0 (Fig. 2b) indicating 467 that, for forest sites, understorey LAI must be also taken into account. A very weak correlation was 468 found with SoilC (r 2 =0.09; p<0.001, n=67) and no significant correlation with Age, N depo and 469 T MEAN were found for forest sites (Fig. 2 c-f). 470 471

[FIGURE 2] 472 473
The multiple regression analysis conducted with the stepwise AIC method including 474 simultaneously all sites, showed that the two best predictors of R 0 were LAI MAX and SoilC 475 (Multiple r 2 =0.57; p<0.001; n=68) which were both positively correlated with R 0 (Tab. 3). LAI MAX 476 was the best predictor of spatial variability of R 0 for all sites confirming the results of the pairwise 477 regression analysis above mentioned but the linear model which included the SoilC as additional 478 predictor led to a significant, though small, reduction in the AIC during the stepwise procedure. 479

Considering only the undisturbed temperate and boreal forest sites (ENF, DBF, MF), the 480 predictive variables of R 0 selected were LAI MAX and N depo . (Multiple r 2 =0.67; p<0.001; n=23). For 481 these sites both LAI MAX , which was still the main predictor of spatial variability of R 0 , and N depo 482 controlled the spatial variability of R 0 , with N depo negatively correlated to R 0 (Tab. 3). This means 483 that for these sites, once removed the effect of LAI MAX , N depo showed a negative control on R 0 with 484 a reduction of 0.025 gC m -2 day -1 in reference respiration for an increase of 1 kg N ha -1 year -1 . 485

Considering only the disturbed forest sites we found that SoilC and T MEAN were the best predictors 486 of spatial variability of R 0 (Multiple R 2 = 0.80, p<0.001, n=10). 487

In On the basis of the aforementioned results, the GPP as well as the linear dependency between R 0 499 and LAI MAX were included into the 'TP Model' leading to a new model formulation (Eq 9). The 500 final formulation is basically the 'TP Model' with the addition of biotic drivers (daily GPP and 501 LAI MAX ) and hereafter referred to as 'TPGPP-LAI Model', where the suffixes GPP and LAI reflect 502 the inclusion of the biotic drivers in the climate-driven model: 503 504

( ) ( ) α α α - + - + ⋅ ⋅           + ⋅ + =         - - - = 1 1 0 0 0 0 1 1 2 0 P k P k e GPP k LAI a R R T T T T E R MAX LAI LAI ECO A ref 4 4 4 3 4 4 4 2 1
(9) 505 506 where the term, R LAI=0 + a LAI LAI MAX , describes the dependency of the basal rate of respiration (R 0 507 in Table1) on site maximum seasonal ecosystem LAI. Although we found that SoilC and N depo may 508 help to explain the spatial variability of R 0 , in the final model formulation we included only the 509 LAI MAX . In fact the model is primarily oriented to the up-scaling and spatial distributed information 510 of SoilC, N depo and disturbance may be difficult to be gathered and usually are affected by high 511 uncertainty. 512

The parameters R LAI=0 and a LAI listed in Table 4 were introduced as fixed parameters in the 513 'TPGPP-LAI Model'. For wetlands and mixed forests the overall relationship between LAI MAX and 514 R 0 was used. For wetlands, available sites were insufficient to construct a statistically significant 515 relationship while for mixed forests the relationship was not significant (p=0.146). 516 PFT specific model parameters (k 2 , E 0 , k, α) of 'TPGPP-LAI Model' were then derived using all 517 data from each PFT contemporarily and listed with their relative standard errors in Table 4. No 518 significant differences in parameter values were found when estimating all the parameters 519 simultaneously (a LAI, R LAI=0, k 2 , E 0 , k, α). 520

The scatterplots of the observed vs modelled annual sums of R ECO are shown in Figure 3, while 521 results and statistics are summarized in Table 5. The model was well able to describe the 522 interannual and intersite variability of the annual sums over different PFTs, with the explained 523 variance varying between 40% for deciduous forests and 97% for shrublands and evergreen 524 broadleaved forests. Considering all sites, the explained variance is 81%, with a mean error of about 525 17% (132.99 gCm The r 2 cv ranges from 0.52 (for EBF) to 0.80 (for CRO) while the r 2 cv,kfold ranges from 0.58 (for 537 DBF) to 0.81 (for GRA). The cross-validated statistics averaged for each PFT are always higher for 538 the k-fold than for the training/evaluation splitting cross-validation. 539

The analysis of model residuals time series of the deciduous broadleaf forest (Fig. 4) showed a 540 systematic underestimation during the springtime development phase and, although less clear, on 541 the days immediately after leaf-fall. A similar behaviour was also found for croplands and 542 grasslands during the days after harvesting or cuts (Fig. 5). 543 The results obtained with the purely climate-driven model ('TP Model') and the best model 551 formulation selected in the site-by-site analysis (i.e. 'LinGPP', Tab. 1) confirm the strong 552 relationship between carbon assimilation and R ECO highlighting that this relationship must to be 553 included into models aimed to simulate temporal variability of R ECO . 554

Respiration appears to be strongly driven by the GPP in particular in grasslands, savannas and 555 croplands as already pointed out by several authors in site-level analysis [START_REF] Bahn | Soil Respiration in European Grasslands 808 in Relation to Climate and Assimilate Supply[END_REF], Moyano 556 et al., 2007, Wohlfahrt et al., 2008a[START_REF] Xu | Seasonal variation in carbon dioxide exchange over a Mediterranean 1203 annual grassland in California[END_REF]. For croplands and grasslands growth 557 respiration is controlled by the amount of photosynthates available and mycorrhizal respiration, 558 which generally constitutes a large component of soil respiration (e.g. [START_REF] Moyano | Response of mycorrhizal, rhizosphere and soil basal 1044 respiration to temperature and photosynthesis in a barley field[END_REF], 559 Kuzyakov & Cheng, 2001). 560

For wetlands instead the weak relationship between respiration and GPP can be explained by the 561 persistence of anaerobic conditions, decomposition proceeds more slowly with an accumulation of 562 water table depth rather than to other factors [START_REF] Lloyd | Annual carbon balance of a managed wetland meadow in the Somerset Levels, 1009 UK[END_REF]. 564

The lower correlation observed for forest ecosystems than for grasslands and croplands may be 565 due to the higher time for translocation, in trees, of substrates from canopy to roots, related to the 566 rates of phloem carbon transport [START_REF] Nobel | Physicochemical and Environmental Plant Physiology[END_REF], which affect the reactivity of the respiration and the 567 release of exudates or assimilates from roots as response to productivity (Mencuccini & Höltta, 568 2010). This is very often cause of time lags between photosynthesis and respiration response but 569 may justify the reduction of correlation between model residuals and GPP estimated at the same 570 day. 571

A clear time lag between GPP and R ECO response was not detected. In fact both the residual 572 analysis (Fig. 1a) and the analysis conducted with the 'LinGPP' model formulation (Fig. 1b) 573 confirmed the general absence of a time lag with the only exception of DBF where a time lag of 3 574 days was observed although the results were not statistically significant. However, in our opinion, 575 these results do not help to confirm or reject the existence of a time lag for several reasons: i) in 576 some studies (e.g. [START_REF] Baldocchi | How switches and lags in biophysical regulators affect 821 spatial-temporal variation of soil respiration in an oak-grass savanna[END_REF][START_REF] Tang | Spatial-temporal variation in soil respiration in an oak-grass 1155 savanna ecosystem in California and its partitioning into autotrophic and heterotrophic 1156 components[END_REF] a lag on the sub-daily time scale 577 was identified and the lags on the daily time scale were attributed to an autocorrelation in weather 578 patterns (i.e. cyclic passage of weather fronts with cycles in temperature or dry and humid air 579 masses) which modulates the photosynthetic activities, since our analysis focused on daily data we 580

were not able to identify the existence of sub-daily time lags; ii) lag effects may be more 581 pronounced under favorable growing conditions or during certain periods of the growing season, the 582 analysis of which analysis is out of scope of present study. 583 584

Spatial variability of reference respiration rates 585 586

The relationship between reference respiration rates (R 0 ) derived by using the 'LinGPP' model 587 formulation, and LAI MAX (Fig. 2a) is particularly interesting considering that the productivity was 588 already included into the model (i.e. daily GPP is driver of 'LinGPP'). While daily GPP describes 589 the portion of R ECO that originates from recently assimilated carbon (i.e. root/rhizosphere 590 respiration, mychorrizal and growth respiration), LAI MAX is a structural factor which has an 591 additional effect to the short-term productivity and allows to describe the overall ecosystem 592 respiration potential of the ecosystem. For instance, high LAI means increased autotrophic 593 maintenance respiration costs. Moreover LAI MAX can be considered both as an indicator of the 594 general carbon assimilation potential and as an indicator of how much carbon can be released to soil 595 yearly because of litterfall (in particular for forests) and leaf turnover which are directly related to 596

F o r R e v i e w O n l y
basal soil respiration [START_REF] Moyano | Response of mycorrhizal, rhizosphere and soil basal 1044 respiration to temperature and photosynthesis in a barley field[END_REF]. At recently disturbed sites, this equilibrium between 597 LAI MAX and soil carbon (through litter inputs) may be broken, for example thinning might lead to a 598 reduction of LAI MAX without any short-term effect on the amount soil carbon, while ploughing in 599 crops or plantations leads solely to a reduction in soil carbon content and not necessarily in LAI. 600 Also in cut or grazed grasslands maximum LAI does not correspond well with litter input because 601 most of this carbon is exported from the site and only partially imported back (as organic manure). 602

This explains why the multiple linear model including LAI MAX and SoilC was selected as the best 603 by the stepwise AIC regression using all the sites contemporarily and why considering only 604 disturbed forest ecosystems we SoilC was selected as best predictor of R 0 (Tab. 3). 605

Particularly interesting is also the negative control of N depo on R 0 with a reduction of 0.025 gC m - 606 2 day -1 in R 0 for an increase of 1 kg N ha -1 year -1 . The reduction of heterotrophic respiration in sites 607 with high total nitrogen deposition load was already described in literature and in some site-level 608 analysis and attributable to different processes. For instance soil acidification at high N depo loads 609 may inhibit litter decomposition suppressing the respiration rate [START_REF] Freeman | Export of dissolved organic carbon from peatlands 893 under elevated carbon dioxide levels[END_REF], Knorr et al., 610 2005) and increasing in N depo can increase N concentration in litter with a reduction of litter 611 decomposition rates [START_REF] Berg | Effect of N deposition on decomposition of plant litter and soil organic 827 matter in forest ecosystems[END_REF][START_REF] Persson | Carbon mineralization in 1072 European forest soils[END_REF] and the consequent reduction of 612 respiration. The latter process is more debated in literature because increased N supply may lead to 613 higher N release from plant litter, which results in faster rates of N cycling and in a stimulation of 614 litter decomposition (e.g. [START_REF] Tietema | 1165 Nitrogen cycling in acid forest soils subject to increased atmospheric nitrogen input[END_REF]. However this process is not always clear (e.g. Aerts 615 et al., 2006): in litter mixtures, N-rich and lignin-rich litter may chemically interact with the 616 formation of very decay-resistant complexes [START_REF] Berg | Litter mass-loss rates in pine forests of europe and 824 eastern united states -some relationship with climate and litter quality[END_REF]. In addition, litter with a high 617 concentration of condensed tannins may interact with N-rich litter reducing the N release from 618 decomposing litter as described in [START_REF] Hattenschwiler | The role of polyphenols in terrestrial ecosystem nutrient 932 cycling[END_REF]. Thus, the supposed 619 stimulating effects of N addition on N mineralization from decomposing litter may be counteracted 620 by several processes occurring in litter between N and secondary compounds, leading to chemical 621 immobilization of the added N (e.g. [START_REF] Pastor | Little bluestem litter dynamics in Minnesota old fields[END_REF][START_REF] Vitousek | Heterotrophic nitrogen fixation in decomposing litter: Patterns and 1189 regulation[END_REF] 622

Although the absolute values are a matter of recent debate (De [START_REF] Vries | Ecologically implausible carbon response? 869[END_REF], Magnani et al., 623 2007[START_REF] Sutton | 1146 Uncertainties in the relationship between atmospheric nitrogen deposition and forest carbon 1147 sequestration[END_REF], it is agreed that N depo stimulates net carbon uptake by temperate and 624 boreal forests. As net carbon uptake is closely related to respiration, once the effect of age is 625 removed, it can be seen that increased N depo has the potential to drive R ECO in either directions. The 626 stimulation of GPP as consequence of the increasing N depo is already include in the model since 627 GPP is a driver. Additionally our analysis suggests that overall an increased total N depo in forests 628 tends to reduce reference respiration. Without considering the effects introduced by N depo in our 629 models we may overestimate R ECO , with a consequent underestimation of the carbon sink strength 630 to other anthropogenic nitrogen inputs (fertilizers, animal excreta) (e.g. [START_REF] Galloway | Transformation of the nitrogen cycle: 895 Recent trends, questions, and potential solutions[END_REF]. 632 However, considering i) that LAI MAX is the most important predictor of R 0 , ii) that the uncertainty 633 in soil carbon and total nitrogen deposition maps is usually high, iii) that the spatial information on 634 disturbance is often lacking and finally iv) that our model formulation is oriented to up-scaling 635 issues, we introduced LAI MAX as the only robust predictor of the spatial variability of R 0 in the final 636 model formulation. 637

The use of LAI MAX is interesting for an up-scaling perspective (e.g. at regional or global scale) 638 since can be derived by remotely sensed vegetation indexes (e.g. normalized vegetation indexes or 639 enhanced vegetation indexes) opening interesting perspectives for the assimilation of remote 640 sensing products into the 'TPGPP-LAI Model'. 641

The intercepts of the PFT-based linear regression between R 0 and LAI MAX (Tab.4) suggest that, 642 when the LAI MAX is close to 0 ('ideally' bare soil), the lowest R 0 takes place in arid (EBF,SHB and 643 SAV) and agricultural ecosystems,. The frequent disturbances of agricultural soils (i.e. ploughing 644 and tillage), as well as management, reduce soil carbon content dramatically. In croplands, the 645 estimated R 0 is very low in sites with low LAI. However, with increasing LAI MAX , R 0 shows a rapid 646 increase, thus resulting in high respiration rates for crop sites with high LAI. For EBF, SHB and 647 SAV the retrieved slopes are typical of forest ecosystems, while the intercepts are close to zero 648 because of the lower soil carbon content usually found in these PFTs [START_REF] Raich | The global carbon dioxide flux in soil respiration and its 1082 relationship to vegetation and climate[END_REF]. 649

Because of the few available sites representing and on similarity in terms of climatic characteristics, 650 savannas, shrublands were grouped. 651

In grasslands, the steeper slope (a LAI ) value found (1.14 ± 0.33) suggests that R 0 increases 652 rapidly with increasing aboveground biomass as already pointed out in literature (Wohlfahrt et al., 653 2008a, Wohlfahrt et al., 2005a, Wohlfahrt et al., 2005b), i.e. an increase in LAI MAX leads to a 654 stronger increase in R 0 than in other PFTs. 655

In forest ecosystems, and in particular in evergreen needleleaf and deciduous broadleaf forests, 656 the physical meaning of the higher intercept may be found in less soil disturbance. In boreal forests, 657 the soil carbon stock is generally high even at sites with low LAI MAX , thus maintaining an overall 658 high R 0 which is less dependent on the LAI MAX . 659 training/evaluation splitting in fact, the excluded site for each PFT is modelled using a 666 parameterization derived from the other sites within the same PFT. However, the k-fold is more 667 optimistic than training/evaluation splitting cross-validation because the data set is less disturbed 668 and the calibration and validation datasets are statistically more similar. In the training/evaluation 669 splitting, instead, we exclude one site which is completely unseen by the training optimization 670 procedure. 671

The derived parameterization of the 'TPGPP-LAI Model' reported in Table 4 may be considered 672 as an optimized parameterization for the application of the model at large scale (e.g. continental or 673 global). For this application is necessary to link of the developed model with a productivity model 674 and remote sensing products necessary for the estimation of LAI. One of the main advances 675 introduced by this model formulation is the incorporation of GPP and LAI as driver of the 676 ecosystem respiration, which importance in modeling Reco is above discussed. These variables are 677 necessary to improve the description of both the temporal and spatial dynamics or R ECO . These 678 results imply that empirical models used with remote sensing (e.g. [START_REF] Reichstein | Reduction of ecosystem productivity and respiration 1093 during the European summer 2003 climate anomaly: a joint flux tower, remote sensing and 1094 modelling analysis[END_REF], 679 Reichstein et al., 2003a[START_REF] Veroustraete | Estimation of carbon mass fluxes over Europe using the 1184 C-Fix model and Euroflux data[END_REF] underestimate the amplitude of R ECO an might 680 lead to wrong conclusions regarding the interpretation of seasonal cycle of the global CO 2 growth 681 rate and annual carbon balance. 682

The values of the 'TPGPP-LAI Model' parameters (Tab. 4) related to the precipitation (k, α) 683 indicated a much stronger nonlinearity in the response of R ECO to precipitation for shrublands, 684 wetlands and croplands than for forest ecosystems (Fig. 6). Wetlands and croplands reached 685 saturation (no limitation of water on respiration) after a small rain event underlying their 686 insensitivity to precipitation owing to the presence of water in wetland soils and irrigation in 687 croplands. Grasslands are very sensitive to rain pulse as described in [START_REF] Xu | Seasonal variation in carbon dioxide exchange over a Mediterranean 1203 annual grassland in California[END_REF], 688 while savannas and evergreen broadleaved forests showed a strong limitation when rainfall was 689 scanty and f(P) saturation exceed 50 mm month -1 . The parameters related to GPP dependency (k 2 ) 690 estimated at PFT level confirm all the results obtained at site level indentifying a clear sensitivity of 691 grasslands and savannah to GPP. 692

[FIGURE 6] 693

However, when comparing these parameterizations, it is very likely that a background 694 correlation between precipitation, short-term productivity and soil respiration confused the apparent 695 response of respiration to water availability in the 'TPGPP-LAI Model '. 696 Despite the good accuracy, some criticisms and limitations of the 'TPGPP-LAI Model' were 697 identified, in particular for the deciduous broadleaf forests. The systematic underestimation during 698 due to the intense activity of vegetation during bud burst not described by the model. This 700 hypothesis is confirmed by different authors. For instance, [START_REF] Davidson | A distinct seasonal 860 pattern of the ratio of soil respiration to total ecosystem respiration in a spruce-dominated 861 forest[END_REF] pointed out that 701 during spring development, specific root respiration increases with increasing soil temperature and 702 the concomitant root growth increases the amount of respiring tissue. Moreover, during bud burst 703 also leaf growth, starch mobilisation and increased phloem transport may contribute to this pulse in 704 respiration as shown by [START_REF] Knohl | Large carbon uptake by an unmanaged 250-978 year-old deciduous forest in Central Germany[END_REF]. A systematic underestimation was also observed 705 immediately after the leaf-fall, in which the increase in heterotrophic respiration stimulated by the 706 decomposition of fresh litter was not completely described by the model. These results are in 707 accordance with [START_REF] Davidson | Soil water content and temperature as independent or 863 confounded factors controlling soil respiration in a temperate mixed hardwood forest[END_REF] whose showed that the sensitivity of respiration to 708 temperature derived using long-term data input is different from short-term sensitivity because it is 709 confused with other seasonally varying factors. At some DBF sites (US-HA1, Fig 4) the 710 observed fluxes are lower than the modelled ones during the foliaged period. Also the overall plot 711 for DBF in Fig 4 shows that model values are generally higher than observations. These 712 considerations suggest that the link between phenological models describing overall foliar 713 development [START_REF] Jolly | A generalized, bioclimatic index to predict foliar 969 phenology in response to climate[END_REF][START_REF] Migliavacca | European larch phenology in the Alps: can 1028 we grasp the role of ecological factors by combining field observations and inverse 1029 modelling[END_REF] and semi-empirical carbon flux models 714 may be useful for the correction of the long-term sensitivity in active spring or summer periods. 715

Another option is the assimilation of remotely-sensed time series from which the main phenological 716 phases may be derived (e.g. derivative methods) and used for instance for the correction of the 717 temporal variability of model parameters. 718

We also found a similar behaviour of croplands and grassland during the days after harvesting or 719 cuts, when respiration increased because of the decomposition of organic residues (e.g. grass or 720 crop residues) as depicted for example in Fig. 5. In this case, the model was unable to describe 721 increased respiration following the harvest. 722 In this study we proposed a model ('TPGPP-LAI Model') for the simulation of R ECO which 727 include the explicit dependency of the respiration to the productivity. We demonstrated that the 728 dependency of respiration on some measure of short-term productivity (e.g. GPP) needs to be 729 included in models simulating ecosystem respiration at regional and global scale in order to 730 improve the description of carbon fluxes and feedbacks between respiration and productivity. 731 In addition, the general site productivity (using maximum seasonal LAI as a proxy) is another 732 important additional variable which accounts for the spatial variability of reference respiration 733 within different plant-functional types. In other words, the LAI MAX can be used as an indicator of 734 the potential respiration for a specific site related to long-term respiration (i.e. low frequencies of 735 the modelled respiration) while GPP and climate drive the short-term respiration response (i.e. the 736 high frequencies of the modelled respiration). This opens interesting perspectives for assessing 737 properties related to respiration using remote sensing products. Soil carbon content and total 738 atmospheric nitrogen deposition may represent under certain circumstance additional parameters 739 enhancing and suppressing, respectively, reference respiration rates. 740

We demonstrated that variables related to productivity and site structure are necessary to 741 improve the description of both the temporal and spatial dynamics or R ECO . These results imply that 742 empirical models driven only by climate underestimate the amplitude of R ECO and might lead to 743 wrong conclusions regarding the interpretation of seasonal cycle of the global CO 2 growth rate and 744 annual carbon balance. 745

We provided a parameterization of the 'TPGPP-LAI Model' for a PFT-based application of the 746 model at large scale (e.g. continental or global). We have shown that the temporal, spatial and 747 interannual variability of ecosystem respiration can be captured quite well by the proposed model. 748

For this application is necessary a link of the developed model with a productivity model (for GPP 749 estimation) and remote sensing products (necessary for the estimation of LAI). One interesting 750 perspective is the integration of the proposed model formulation into the MODIS-GPP/NPP data 751 stream (e.g MOD17 Light Use Efficiency model) for regional and global estimates of R ECO . 752 Finally, we observed that a part of ecosystem respiration variance not explained by the model 753 may be related to phenology in forests and to management in grasslands and croplands. For these 754 reasons we consider the link between phenological models and/or remotely-sensed time series of 755 vegetation indexes and respiration models as well as the inclusion of total nitrogen deposition as an 756 additional driver for improving the description of ecosystem respiration in both space and time. 757

Model

Formula LinGPP ( ) ( ) ( )

α α α - + - + ⋅ ⋅ + =         - - - 1 1 0 0 0 1 1 2 0 P k P k e GPP k R R T T T T E ECO A ref ExpGPP ( ) [ ] ( ) ( ) α α α - + - + ⋅ ⋅ - + =         - - - 1 1 1 0 0 0 2 1 1 2 0 P k P k e e R R R T T T T E GPP k ECO A ref MicMenGPP ( ) ( ) α α α - + - + ⋅ ⋅       + + =         - - - 1 1 0 0 0 1 1 0 P k P k e hR GPP GPP R R R T T T T E MAX MAX ECO A ref addLinGPP ( ) ( ) GPP k P k P k e R R T T T T E ECO A ref 2 1 1 0 1 1 0 0 0 + - + - + ⋅ ⋅ =         - - - α α α addExpGPP ( ) ( ) α α α - + - + ⋅ ⋅ =         - - - 1 1 0 0 0 1 1 0 P k P k e R R T T T T E ECO A ref + ( ) GPP k e R 2 1 2 - addMicMenGPP ( ) ( ) MAX MAX T T T T E ECO hR GPP GPP R P k P k e R R A ref + + - + - + ⋅ ⋅ =         - - - α α α 1 1 0 0 0 1 1 0 F o r R e v i e w O n l y

  modeling of abiotic and biotic factors controlling ecosystem respiration across eddy covariance sites

  Part of the variability in respiration that could not be described by our model could be attributed to a range of factors, including phenology in deciduous broadleaf forests and management practices in grasslands and croplands. ecosystems (R ECO ) is one of the major fluxes in the global carbon cycle 97 and its responses to environmental change is important for understanding climate-carbon cycle 98

  authors have found a time lag between productivity and respiration response. 133

  analyzed and discussed; 169 o we evaluated the inclusion of biotic factors (i.e. GPP) as drivers of R ECO . A range of 170 different model formulations, which differ mainly in regard to the functional 171 responses of R ECO to photosynthesis, were tested in order to identify the best model 172 formulation for the daily description of R ECO at each site; 173 o we analyzed variability of the reference respiration estimated at each site with the 174 aim of identifying, among the different site characteristics, one or more predictors of 175 the spatial variability of this crucial parameter. This can be extremely useful for the 176 application of the model at large spatial scale; 177 2. we optimized the developed model for each PFT (PFT analysis) with the aim of generalizing 178 the model parameters in a way that can be useful for diagnostic, PFT-based, up-scaling of 179 R ECO . The accuracy of the model was assessed by a cross-validation technique and the main 180 weak points of model were critically evaluated and discussed.

  (DBF), 16 grasslands (GRA), 11 croplands (CRO), 8 mixed forests 200 (MF), 5 savannas (SAV), 9 shrublands (SHB), 7 evergreen broadleaved forests (EBF) and 3 201 wetlands (WET). Due to limited number of sites and their similarity, the class SAV included both 202 the sites classified as savanna (SAV) and woody savannas (WSA), while the class SHB included 203 both the open (OSH) and closed (CSH) shrubland sites. For abbreviations and symbols refer to 204 Appendix III. 205 Daily R ECO , GPP and the associated uncertainties of NEE data, together with daily 206 meteorological data such as mean air temperature (T A ) and 30-day precipitation running average 207 (P), were downloaded from the FLUXNET database. 208

  220and calculated at 1°x1° resolution. These data are grid-average downward deposition velocities and 221 do not account for vegetation effects. The data used for the selected sites are shown in the Appendix of R ECO we started from a widely used climate-driven model: 'TP Model' (Eq. 1) 229

  (gC m -2 day -1 ) is the ecosystem respiration at the reference temperature (T ref , K) 235 without water limitations. f(T A ) and f(P) are functional responses of R ECO to air temperature and 236 precipitation, respectively. 237Here temperature dependency f(T A ) is changed from the Q 10 model to an Arrhenius type equation 238 (Eq. 2). E 0 (K) is the activation energy parameter and represents the ecosystem respiration 239 sensitivity to temperature, T ref is fixed at 288.15 K (15°C) and T 0 is fixed at 227.13 K (-46.02°C): 240 241

  Beside the linear dependency the exponential and Michaelis-Menten responses were tested. 265

  /evaluation splitting cross-validation: one site at a time was excluded using the 332 remaining subset as the training set and the excluded one as the validation set. The model 333 was fitted against each training set and the resulting parameterization was used to predict the 334 R ECO of the excluded site. 3352) k-fold cross-validation: the whole data set for each PFT was divided into k randomly 336 selected subsets (k=15) called a fold. The model is fitted against k-1 remaining folds 337 (training set) while the excluded fold (validation set) was used for model evaluation. The 338 cross-validation process was then repeated k times, with each of the k folds used exactly 339 once as the validation set. 340 341 For each validation set of the cross-validated model statistics were calculated (see 'Statistical 342 Analysis' section). Finally, for each PFT we averaged the cross-validated statistics to produce a 343 single estimation of model accuracy in prediction. estimated using the Levenberg-Marquardt method, implemented in the 350 data analysis package "PV-WAVE 8.5 advantage"(Visual Numerics, 2005), a non-linear regression 351

  top of the mineral soil layer and respiration is closely related to temperature and 563

  ecosystems. It is also clear that, in managed sites, such interactions apply equally 631

660

  Final formulation of the model and weak points 661 662These results obtained with the 'TPGPP-LAI Model' cross-validation indicate that the developed 663 model describes the R ECO quite well. In particular results indicate a better description of the 664 R ECO rather than the spatial variability (or across-site variability). In the 665

  phase (Fig 4) is very likely related to the peak in autotrophic respiration 699

  like to thank all the PIs of eddy-covariance sites, technicians, postdoctoral 761 fellows, research associates, and site collaborators involved in FLUXNET who are not included as 762 co-authors of this paper, without whose work this analysis would not have been possible. This work 763 is the outcome of the La Thuile FLUXNET Workshop 2007, which would not have been possible 764 without the financial support provided by CarboEuropeIP, FAO-GTOS-TCO, iLEAPS, Max Planck 765 Institute for Biogeochemistry, National Science Foundation, University of Tuscia and the US 766

Figure 1

 1 Figure 1 -a) Pearson's correlation coefficients (r) for the residual of observed minus modelled R ECO versus measured GPP and a function of time lag; b) average model performances (EF and RMSE) for deciduous broadleaf forests as a function of the time lag between GPP and R ECO response.Results obtained running the 'LinGPP' formulation with different GPP time series, from the GPP measured at the same day up to the GPP measured one week before the R ECO . Error bars represent the standard deviation of model statistics calculated at each site. The definitions of different PFTs are: evergreen needleleaf forest (ENF), deciduous broadleaf forest (DBF), grasslands (GRA), croplands (CRO), savannah (SAV), shrublands (SHB), evergreen broadleaf forest (EBF), mixed forest (MF), wetland (WET).

Figure 2 -

 2 Figure 2 -Correlation between reference respiration (R 0 ) and a) seasonal maximum leaf area index (LAI MAX ) of understorey and overstorey, b) overstorey peak leaf area index (LAI MAX,o ), c) total soil carbon content (SoilC), d) stand age for forest ecosystems (Age), e) total atmospheric nitrogen deposition for forest sites (N depo ) and f) mean annual temperature. In panels a), b), c), d) and f) different symbols represent different PFT. In panel e) full circles represent disturbed sites while open circles the undisturbed ones. The r 2 , p and number of sites (n) were reported. The regression line and the 95% confidence interval are given if the relationship is significant. The definitions of different PFTs are: evergreen needleleaf forest (ENF), deciduous broadleaf forest (DBF), grasslands (GRA), croplands (CRO), savannah (SAV), shrublands (SHB), evergreen broadleaf forest (EBF), mixed forest (MF), wetland (WET).

Figure 3 -

 3 Figure 3 -Scatterplots of annual observed vs modelled R ECO obtained using the 'TPGPP-LAI Model'.Each panel represent a different plant functional type (PFT). The definitions of different PFTs are: evergreen needleleaf forest (ENF), deciduous broadleaf forest (DBF), grasslands (GRA), croplands (CRO), savannah (SAV), shrublands (SHB), evergreen broadleaf forest (EBF), mixed forest (MF), wetland (WET).

Figure 4 -

 4 Figure 4 -Time series of average monthly model residuals for different deciduous broadleaf forest (DBF) sites. The vertical grey dashed lines represent the phenological dates. Average phenological dates were derived for US-Ha1 from literature (Jolly et al. 2005) while for other sites they were retrieved from the FLUXNET database. Average phenological dates, bud-burst and end-of-growing season are respectively: US-Ha1 ( 115-296),DE-Hai (126-288), FR-Hes (120-290), FR-Fon (125-292), IT-Ro1 (104-298) and CA-Oas (146-258)..

Figure 5 -Figure 6 -

 56 Figure 5 -Time series of observed (open circles) and modeled (black circles) for the IT-MBo site (a,b) and for the ES-ES2 site (c, d), grey dashed lines represent the dates of cuts indicated in the database (the date may be indicative), the model underestimation of fluxes in the days after each cut is clear.

Figure 1 Figure 2 -Figure 3 -Figure 4 -Figure 5 -Figure 6 -

 123456 Figure 1 -a) Pearson's correlation coefficients (r) for the residual of observed minus modelled RECO versus measured GPP and a function of time lag; b) average model performances (EF and RMSE) for deciduous broadleaf forests as a function of the time lag between GPP and RECO response. Results obtained running the 'LinGPP' formulation with different GPP time series, from the GPP measured at the same day up to the GPP measured one week before the RECO. Error bars represent the standard deviation of model statistics calculated at each site. The definitions of different PFTs are: evergreen needleleaf forest (ENF), deciduous broadleaf forest (DBF), grasslands (GRA), croplands (CRO), savannah (SAV), shrublands (SHB), evergreen broadleaf forest (EBF), mixed forest (MF), wetland (WET). 191x108mm (300 x 300 DPI)

  was almost negligible (cAIC was 38.22 ± 2.52 and 38.26 ± 2.45 for 430

	two model formulations
	e
	w
	O n l
	y
	14

Table 5 (

 5 left column) the statistics of the pairwise regression analysis between R 0 and LAI MAX 488 for each PFT are reported. The best fitting was obtained with the linear relationship for all PFTs 489

	490	except for deciduous forests for which the best fitting was obtained with the exponential
	491	relationship R 0 =R LAI=0 (1-e -aLAI ).
	492
	493	[TABLE 3 AND TABLE 4]
	494
	495	PFT-Analysis
	496

  -2 yr -1 ) of the annual observed R ECO .

	529
	530	Evaluation of model predictions accuracy and weak points
	531
	532	The results obtained with the k-fold and training/evaluation split cross-validation are listed in
	533	Table 6.
	534
	535	[TABLE 6]
	536
		F o r
		R
		e v
		i e
		w
		O n l
		y
	526
	527
	528	[TABLE 5, FIGURE 3]
		17

Table 1 -Different model formulation of the dependency of ecosystem respiration (R ECO ) on Gross Primary Productivity (GPP) used in this analysis.

 1

Table 2 -Statistics of fit for the climate-driven model ('TP Model') and the best model selected among the models listed in Tab. 1 according to the consistent Akaike Information Criterion (cAIC). Statistics are averaged per Plant Functional Type (PFT). Except for croplands (CRO), 'LinGPP' is selected as the best model formulation. EF is the modelling efficiency while RMSE is the root mean square error (Jannsens and Heuberger, 1995). The definitions of different PFTs are: evergreen needleleaf forest (ENF), deciduous broadleaf forest (DBF), grasslands (GRA), croplands (CRO), savannah (SAV), shrublands (SHB), evergreen broadleaf forest (EBF), mixed forest (MF), wetland (WET). The list of acronyms is also provided in Appendix II. Values in brackets are the standard deviations.

 2

		'TP Model'	'LinGPP Model'	Best Model Selected
	PFT	EF	RMSE	EF	RMSE
	ENF	0.71(0.14)	1.02 (0.35) 0.78 (0.14) 0.83 (0.21)	LinGPP
	DBF	0.63 (0.17) 1.15 (0.51) 0.72 (0.13) 0.98 (0.41)	LinGPP
	GRA	0.62 (0.18) 1.35 (0.43) 0.83 (0.07) 0.91 (0.33)	LinGPP
	CRO	0.55 (0.18) 1.55 (0.53) 0.82 (0.08) 1.01 (0.33)	addLinGPP
	SAV	0.38 (0.16) 0.78 (0.24) 0.72 (0.06) 0.53 (0.15)	LinGPP
	SHB	0.59 (0.29) 0.67 (0.50) 0.66 (0.29) 0.58 (0.51)	LinGPP
	EBF	0.42 (0.27) 1.11 (0.55) 0.58 (0.23) 0.91 (0.49)	LinGPP
	MF	0.67 (0.18) 0.96 (0.72) 0.82 (0.13) 0.78 (0.50)	LinGPP
	WET	0.67 (0.18) 0.96 (0.51) 0.85 (0.48) 0.79 (0.07)	LinGPP

Table 3 -Results of the model selection conducted with the Stepwise AIC method for the sites belonging to all the PFT (All PFTs) and for undisturbed temperate and boreal forests identified in the Appendix II (Undisturbed Forests). Coefficients (a 1 ,a 2 , const), their significance and the statistics of the best model selected are reported. In parenthesis the standard error of the coefficients are reported. The significance of coefficients is also reported

 3

	(*** p<0.001, ** p<0.01, * p<0.05, . p<0.1).
	Model	Best Model Selected	a 1		a 2		const		r 2	r 2 adj.	p	n
	All PFTs	R 0 =a 1 LAI MAX +a 2 SoilC + const	0.412	***	0.045	**	0.582	*	0.58	0.57	<0.001	68
			(0.048)		(0.015)		(0.251)
	Undisturbed Forest	R 0 = a 1 LAI MAX + a 2 N depo + const	0.469	***	-0.025	.	0.948	*	0.70	0.67	<0.001	23
	(MF+DBF+ENF)		(0.069)		(0.017)		(0.377)
	Disturbed Forests	R 0 = a 1 SoilC + a 2 T MEAN + const	0.211	**	-0.188	**	3.487	*	0.85	0.80	<0.001	10
			(0.051)		(0.059)		(0.982)

Table 4 -Parameters of the relationships between reference respiration (R 0 ) defined at 15°C and seasonal maximum LAI for each Plant Functional Type (PFT). The standard errors of model parameters are reported in parenthesis. Determination coefficients and statistical significance are also shown.-'TPGPP-LAI Model' parameters estimated for each Plant Functional Type (see Appendix II). Standard errors estimated with the bootstrap algorithm are reported in parentheses. Model statistics are also given. ''TPGPP-LAI Model' is defined in Eq. 9. The definitions of different PFTs are: evergreen needleleaf forest (ENF), deciduous broadleaf forest (DBF), grasslands (GRA), croplands (CRO), savannah (SAV), shrublands (SHB), evergreen broadleaf forest (EBF), mixed forest (MF), wetland (WET).

 4

		Parameters and statistics	Final Model Parameters				Fitting statistics
	PFT	R LAI=0	(R0 vs LAI MAX ) a LAI r 2	p	k 2	E 0 [K]	Α	K [mm]	r 2	EF	RMSE [gCm -2 day -1 ]	MAE [gCm -2 day -1 ]
	ENF	1.02	0.42	0.50 <0.001	0.478	124.833	0.604	0.222	0.79 0.70	1.072	0.788
		(0.42)	(0.08)			(0.013)	(4.656)	(0.065)	(0.070)
	DBF	1.27	0.34	0.46	<0.01	0.247	87.655	0.796	0.184	0.65 0.52	1.322	0.899
		(0.50)	(0.10)			(0.009)	(4.405)	(0.031)	(0.064)
	GRA	0.41	1.14	0.60 <0.001	0.578	101.181	0.670	0.765	0.82 0.80	1.083	0.838
		(0.71)	(0.33)			(0.062)	(6.362)	(0.052)	(1.589)
	CRO	0.25	0.40	0.52 <0.001	0.244	129.498	0.934	0.035	0.80 0.79	0.933	0.659
		(0.66)	(0.11)			(0.016)	(5.618)	(0.065)	(3.018)
	SAV	0.42	0.57	0.54 <0.005	0.654	81.537	0.474	0.567	0.65 0.60	0.757	0.535
		(0.39)	(0.17)			(0.024)	(7.030)	(0.018)	(0.119)
	SHB	0.42	0.57	0.54 <0.005	0.354	156.746	0.850	0.097	0.73 0.60	0.618	0.464
		(0.39)	(0.17)			(0.021)	(8.222)	(0.070)	(1.304)
	EBF	-0.47	0.82	0.87 <0.001	0.602	52.753	0.593	2.019	0.55 0.41	1.002	0.792
		(0.50)	(0.13)			(0.044)	(4.351)	(0.032)	(1.052)
	MF	0.78	0.44	0.52 <0.001	0.391	176.542	0.703	2.831	0.86 0.79	0.988	0.723
		(0.18)	(0.04)			(0.068)	(8.222)	(0.083)	(4.847)
	WET	0.78	0.44	0.52 <0.001	0.398	144.705	0.582	0.054	0.87 0.86	0.403	0.292
		(0.18)	(0.04)			(0.013)	(8.762)	(0.163)	(0.593)

Table 5 -Statistics of the modelled (x-axis) vs measured (y-axis) annual R ECO with the 'TPGPP-LAI Model'. Number of site-years for each PFT are also reported. The definitions of different PFTs are: evergreen needleleaf forest (ENF), deciduous broadleaf forest (DBF), grasslands (GRA), croplands (CRO), savannah (SAV), shrublands (SHB), evergreen broadleaf forest (EBF), mixed forest (MF), wetland (WET).

 5

	Statistics

Table 6 -Results of Training/Evaluation splitting and k-fold cross-validation of the 'TPGPP-LAI Model' averaged per plant functional type as defined in the Appendix II. The definitions of different PFTs are: evergreen needleleaf forest (ENF), deciduous broadleaf forest (DBF), grasslands (GRA), croplands (CRO), savannah (SAV), shrublands (SHB), evergreen broadleaf forest (EBF), mixed forest (MF), wetland (WET).

 6

			Training/Evaluation Splitting			k-fold Cross-Validation
	PFT	r 2	EF	RMSE	MAE	r 2	EF	RMSE	MAE
				[gCm -2 day -1 ]	[gCm -2 day -1 ]			[gCm -2 day -1 ]	[gCm -2 day -1 ]
	ENF	0.74	0.74	1.170	0.854	0.76	0.76	1.145	0.827
	DBF	0.54	0.48	1.443	1.017	0.58	0.50	1.374	0.967
	GRA	0.79	0.79	1.227	0.881	0.81	0.80	1.174	0.819
	CRO	0.80	0.80	1.208	0.889	0.80	0.79	1.254	0.926
	SAV	0.57	0.54	0.831	0.623	0.60	0.59	0.717	0.515
	SHB	0.71	0.58	0.954	0.720	0.68	0.67	1.180	0.790
	EBF	0.52	0.28	1.350	0.985	0.70	0.69	0.957	0.928
	MF	0.71	0.71	1.326	0.927	0.75	0.74	1.254	0.871
	WET	0.79	0.75	0.566	0.320	0.83	0.82	0.490	0.312

  Table3were used to draw the curves. The definitions of different PFTs are: evergreen needleleaf forest (ENF), deciduous broadleaf forest (DBF), grasslands (GRA), croplands (CRO), savannah (SAV), shrublands (SHB), evergreen broadleaf forest (EBF), mixed forest (MF), wetland (WET). Blue boxes) partitioning. The median of the differences of parameters governing the response to GPP (k 2 ) estimated at each site with the two different datasets are not statistically different from 0 except for ENF and DBF (for both p<0.05). No statistical differences were found for model statistics. Data were grouped in box-plots for each PFT. The definitions of different PFTs are: evergreen needleleaf forest (ENF), deciduous broadleaf forest (DBF), grasslands (GRA), croplands (CRO), savannah (SAV), shrublands (SHB), evergreen broadleaf forest (EBF), mixed forest (MF), wetland (WET).

	Figure AI -Box-plot of the differences at each site between the Pearson's correlation coefficient
	between 'TP Model' residuals and GPP computed using FLUXNET partitioning (r TPModel-GPPFLUX )
	and Lasslop's partitioning (r TPModel-GPPLasslop ). Data were grouped in box-plots for each PFT. The
	definitions of different PFTs are: evergreen needleleaf forest (ENF), deciduous broadleaf forest
	(DBF), grasslands (GRA), croplands (CRO), savannah (SAV), shrublands (SHB), evergreen
	broadleaf forest (EBF), mixed forest (MF), wetland (WET)
	Figure AII -Box-plot of the parameters a) R 0 , b) k 2 , c) EF and d) RMSE estimated using
	FLUXNET (red boxes) and Lasslop's (

Table A II -Site characteristics derived from the FLUXNET database. R 0 is the reference respiration estimated with the LinGPP model formulation, LAI is the maximum seasonal leaf area index of the ecosystems (understorey and overstorey), LAI MAX,o is the maximum leaf area index of the solely overstorey, SoilC is the total soil carbon content, Age is the stand age, T mean is the annual average mean temperature, Ndepo is the total atmospheric nitrogen deposition derived as described in the method section, u* is the median of the yearly friction velocity threshold identified at each site by using the method described in Papale et al., (2006). Sites with (*) in the column dist (disturbance) represent sites with recent disturbance according to what reported in the FLUXNET database.

 A

	FI-Hyy US-IB1	Hyytiala Fermi National Accelerator Laboratory-Batavia (Agricultural site) 1.90 3.63	7.00 5.25	6.7 5.3	5.60 6.30	2.87 14.95	47	* *	4.47 13.83 0.010 0.296
	FI-Kaa FI-Sod US-KS2 US-Los APPENDIX III -Kaamanen wetland Sodankyla Kennedy Space Center (scrub oak) Lost Creek	1.27 2.09 1.92 1.94	0.70 1.20 2.50 4.24	0.7 1.2 2.5 4.2	3.14 3.60 4.50	1.30 1.07 7.00 3.02	11	*	0.20 1.10 22.26 0.053 0.089 0.211 4.72 0.140
	FR-Fon US-LPH	Fontainebleau Little Prospect Hill	2.20 3.19	5.05 5.00	5.1 5.0	10.20 3.70	23.38 12.27			11.50 0.163 8.82 0.221
	FR-Gri US-Me2	Grignon (after 6/5/2005) Metolius-intermediate aged ponderosa pine	2.16 2.15	3.34 2.80	3.3 2.7	7.90	21.09 3.45	95	*	11.25 0.100 6.82 0.601
	FR-Hes US-Me3	Hesse Forest-Sarrebourg Metolius-second young aged pine	3.17 0.90	6.70 0.52	7.3 0.5	7.17 10.00	26.30 3.45	43 21		10.37 0.152 8.47 0.064
	FR-LBr US-Me4	Le Bray (after 6/28/1998) Metolius-old aged ponderosa pine	3.51 1.28	4.00 2.10	2.5 2.1	10.90 5.56	14.30 3.45	39	*	13.66 0.206 8.32 0.034
	FR-Lq2 US-MMS	Laqueuille extensive Morgan Monroe State Forest	3.26 2.83	3.00 4.62	3.0 4.6	6.60	18.23 18.27		*	7.66 12.28 0.342 0.146
	FR-Pue US-MOz	Puechabon Missouri Ozark Site	2.66 2.09	3.90 4.20	1.9 4.2	6.10	14.46 17.17	66		13.64 0.229 14.87 0.224
	SITE ID AT-Neu AU-How BE-Lon BE-Vie BR-Sp1 BW-Ma1 CA-Ca1 CA-Ca3 CA-Gro CA-Let CA-Mer CA-NS1 CA-NS3 CA-NS6 CA-Oas CA-Ojp CA-Qfo CA-TP4 CA-WP1 CH-Oe1 CN-HaM CN-Ku1 CN-Ku2 CN-Xi2 DE-Bay DE-Hai DE-Kli DE-Tha DK-Ris ID-Pag IL-Yat IT-Amp IT-BCi IT-Cpz IT-MBo IT-Noe IT-Non IT-PT1 IT-Ren IT-Ro1 IT-Ro2 IT-SRo JP-Tef NL-Loo PT-Esp PT-Mi1 RU-Cok UK-EBu UK-Gri US-ARb US-ARM US-Aud US-Bar US-Bkg US-Bn1 US-NC2 US-Ne1 US-Ne2 US-NR1 US-Oho US-PFa US-SO2 US-SO3 US-SP1 US-SP2 US-Syv US-Ton US-UMB US-Var US-WCr US-Wi4 VU-Coc US-Bn2 US-Bn3 US-Bo1 US-Brw US-Dk3 US-FPe	Tower Name Neustift/Stubai Valley Howard Springs Lonzee Vielsalm Sao Paulo Cerrado Maun-Mopane Woodland British Columbia-Campbell River -Mature Forest Site British Columbia-Campbell River -Young Plantation Site Ontario-Groundhog River-Mat. Boreal Mixed Wood Lethbridge Eastern Peatland-Mer Bleue UCI-1850 burn site UCI-1964 burn site UCI-1989 burn site Sask.-SSA Old Aspen Sask.-SSA Old Jack Pine Quebec Mature Boreal Forest Site Ontario-Turkey Point Mature White Pine Western Peatland-LaBiche-Black Spruce/Larch Fen Oensingen1 grass Haibei Alpine Tibet site Kubuqi_populus forest Kubuqi_shrubland Xilinhot grassland site (X03) Bayreuth-Waldstein/WeidenBrunnen Hainich Klingenberg Tharandt-Anchor Station Risbyholm F o r R e v i e w R 0 LAI MAX LAI MAX,o gCm -2 day -1 m 2 m -2 m 2 m -2 4.83 6.50 6.5 1.84 2.40 0.9 2.23 5.62 5.6 2.47 5.10 5.1 3.54 3.50 3.5 0.67 1.10 1.1 2.77 8.40 7.1 3.84 6.70 3.0 4.88 4.30 4.3 1.05 0.80 0.8 0.94 1.30 1.3 3.43 5.68 5.2 6.10 9.81 5.3 2.40 2.97 3.0 3.70 5.10 2.1 1.76 2.60 2.6 2.14 3.70 3.7 3.56 8.00 8.0 0.76 2.61 1.3 3.83 4.85 4.9 2.97 2.78 2.8 0.23 0.23 0.2 0.61 0.20 0.2 0.88 0.25 0.3 5.04 5.60 5.3 2.93 6.08 6.1 4.42 9.73 5.5 5.64 7.60 5.2 2.77 6.00 6.0 O n l y SoilC N depo kgCm -2 kgN year -1 ha -1 years Age dist T mean °C 4.25 18.97 * 6.79 15.10 1.09 * 25.86 0.136 u * m s -1 0.035 3.70 23.12 * 10.88 0.134 3.82 25.22 96 8.31 0.459 8.00 8.32 22.70 0.263 0.50 3.54 * 22.83 0.159 1.51 60 8.67 0.295 1.65 21 * 9.97 0.102 1.82 78 * 3.84 0.408 3.01 * 6.66 5.79 6.69 0.039 16.53 0.69 159 -1.32 0.270 3.64 0.69 45 -1.04 0.192 4.40 0.69 20 -0.25 0.261 1.63 1.28 85 2.10 0.346 1.58 1.18 93 1.75 0.243 3.50 1.45 102 2.66 0.273 3.70 12.17 70 8.95 0.316 1.15 136 3.63 0.017 18.30 23.67 * 9.21 0.043 8.60 2.26 * -5.18 0.065 3.14 8 * 11.09 0.080 3.14 * 11.57 5.88 * 5.96 17.02 13.65 45 7.00 0.353 12.20 17.80 140 8.23 0.519 9.70 14.79 * 8.34 0.099 16.00 14.79 118 * 8.52 0.279 8.51 * 7.47 0.082 Palangkaraya 4.53 5.60 5.6 2.19 * 26.55 Yatir 0.68 2.50 2.5 3.70 7.18 42 * 18.68 0.338 Amplero (after 6/28/2004) 2.49 2.00 2.0 19.30 10.41 * 10.21 0.029 Borgo Cioffi 2.28 5.80 5.8 8.98 * 16.29 0.091 Castelporziano 1.31 3.50 3.5 4.31 11.25 14.82 0.096 Monte Bondone 4.82 2.82 2.8 35.00 18.78 * 5.09 0.075 Sardinia/Arca di Noè 2.84 2.10 2.1 10.00 10.22 45 16.87 0.091 Nonantola 1.27 1.70 1.7 4.80 16.96 17 * 13.91 0.080 Zerbolò-Parco Ticino-Canarazzo 2.65 4.45 2.2 4.59 18.91 14 * 14.53 0.185 Renon/Ritten (Bolzano) 1.79 5.11 4.6 15.20 18.78 * 4.71 0.119 Roccarespampani 1 2.97 4.30 3.0 11.30 13.72 7 * 15.64 0.218 Roccarespampani 2 2.46 4.08 3.9 11.84 13.72 17 14.79 0.095 San Rossore 2.89 4.20 4.2 2.15 16.10 57 15.44 0.201 Teshio Experimental Forest 4.76 7.50 4.5 1.83 * 6.30 0.130 Loobos 4.23 3.50 2.0 2.40 12.24 10.42 0.224 Espirra 2.06 2.80 2.8 5.62 16 16.03 0.231 Mitra (Evora) 1.10 2.30 0.7 5.62 15.86 0.228 Chokurdakh 1.20 1.50 1.5 4.35 0.20 2.62 Easter Bush-Scotland 2.00 3.90 3.9 22.95 6.27 * 9.00 Griffin-Aberfeldy-Scotland 3.72 7.00 7.0 15.00 4.54 25 7.61 0.175 ARM Southern Great Plains burn site-Lamont 2.66 3.25 3.3 13.51 10.71 * 16.97 0.195 ARM Southern Great Plains site-Lamont 0.84 2.10 2.1 11.52 * 15.57 0.075 Audubon Research Ranch 1.28 1.00 1.0 2.55 * 17.28 0.038 Bartlett Experimental Forest 3.91 4.70 5.1 15.50 6.98 70 7.15 0.050 Brookings 1.63 3.00 3.0 8.57 * 8.05 0.098 Bonanza Creek, 1920 Burn site near Delta Junction 1.73 3.50 3.5 0.62 89 NC_Loblolly Plantation 3.66 3.00 3.0 14.33 15 * 15.86 0.147 Mead -irrigated continuous maize site 3.82 6.30 6.3 18.40 13.20 * 11.36 0.098 Mead -irrigated maize-soybean rotation site 2.40 3.75 3.8 21.10 13.20 * 11.43 0.107 Niwot Ridge Forest (LTER NWT1) 3.04 5.60 5.1 16.00 3.77 102 2.46 0.308 Oak Openings 1.57 4.70 4.0 13.49 46 * 11.16 0.136 Park Falls/WLEF 3.31 4.10 4.1 20.20 4.32 4.59 0.211 Sky Oaks-Old Stand 1.15 3.00 3.0 0.87 3.56 78 * 13.77 0.038 Sky Oaks-Young Stand 0.66 1.10 1.1 3.56 4 * 15.87 0.104 Slashpine-Austin Cary-65yrs nat regen 3.04 4.50 4.5 8.00 9.15 65 21.04 0.186 Slashpine-Mize-clearcut-3yr,regen 3.60 3.88 2.9 9.15 9 * 20.56 0.050 Sylvania Wilderness Area 2.80 3.80 3.8 10.47 2.55 350 5.20 0.406 F o r R e v i e w Tonzi Ranch 1.88 2.00 0.6 4.85 1.87 17.36 0.143 Univ. of Mich. Biological Station 3.17 3.95 3.6 3.60 3.83 90 7.35 Vaira Ranch-Ione 2.15 2.50 2.5 1.87 * 15.93 0.047 Willow Creek 2.60 5.40 4.5 9.47 4.32 74 5.77 0.419 Mature red pine (MRP) 1.17 2.80 1.8 4.18 69 10.19 0.162 CocoFlux 4.44 5.65 3.0 4.25 0.39 24 * 24.76 0.188 F o r R O n l y e v -0.82 0.075 Bonanza Creek, 1987 Burn site near Delta Junction 0.88 2.50 2.5 0.62 22 * -0.29 0.071 Bonanza Creek, 1999 Burn site near Delta Junction 0.69 1.10 1.1 0.62 10 * -0.29 0.075 Bondville 2.57 5.25 5.3 16.50 * i e 11.14 0.108 Alaska -Barrow 1.12 1.50 1.5 16.50 0.15 -1.38 0.071 Duke Forest -loblolly pine 1.39 5.20 4.7 9.00 15.07 26 * 14.68 Fort Peck 1.25 2.50 2.5 3.74 * 5.74 0.060 w
	ES-ES1 ES-ES2 ES-LMa ES-VDA US-Fwf US-Ha1 US-Ho1 US-Ho2	El Saler El Saler-Sueca Las Majadas del Tietar Vall d'Alinya Flagstaff -Wildfire Harvard Forest EMS Tower (HFR1) Howland Forest (main tower) Howland Forest (west tower)	3.28 1.04 1.57 1.66 0.80 3.26 3.71 3.59 O n l 49 3.63 5.80 2.00 1.35 0.60 5.20 6.50 5.60 50	2.6 5.8 0.5 1.4 0.6 5.2 6.5 5.6	3.32 3.30 8.80 11.00 12.00	7.68 7.68 6.85 12.02 2.47 12.27 4.19 4.19	75 120 140 140	* * * *	17.41 0.255 18.01 0.070 16.16 0.153 6.51 0.069 12.26 0.082 8.16 0.392 6.60 0.224 6.51
			y
		51

List of acronyms and abbreviations Table AIII -Acronyms and abbreviations

  Precipitation model, from Raich et al. (2000) and modified by Reichstein et al. (2003) TP b Model TP biotic model, containing both the dependency on GPP and ecosystem LAI (Final model formulation)

	T Ref	Reference temperature (15 °C)
	y mod	Modeled data as a function of parameter vector
	y obs	Observed data
	Α	Response of R ECO to null precipitation
	Acronyms CRO Θ DBF Σ	Description Parameter vector Croplands Deciduous Broadleaf Forest Weight of cost function
	E 0 Ω LS	Activation Energy [K] Cost function
	EBF	Evergreen Broadleaf Forest
	EF	Modeling Efficiency from Jannssen and Heuberger (1995)
	ENF	Evergreen Needleleaf Forest
	GPP	Gross Primary Production
	GPP lag,i	GPP measured i days before the observation day of ecosystem respiration
	GRA hRmax IGBP K k 2 LAI MAX LAI MAX,o MAE MDS MF N depo NEE NEE mid NEE night P PFT R 0 R 2 R ECO Rmax RMSE R ref SAV SHB SoilC SWC T 0 T A TP Model	Grasslands GPP value at half saturation International Geosphere Biosphere Programme Half saturation constant of the hyperbolic relationship between R ECO and precipitation Parameter governing the linear and exponential response of R ECO to GPP Maximum Leaf Area Index (Understorey + Overstorey) Maximum Leaf Area Index (Overstorey) Mean Absolute Error from Jannssen and Heuberger (1995) Marginal Distribution Sampling Mixed Forest Total Nitrogen Depositions Net Ecosystem Exchange NEE mid-day NEE night-time 30-day Precipitation running average Plant Functional Type Respiration at reference temperature for TP Model with GPP dependency added Parameter of exponential dependency between GPP and R ECO Ecosystem Respiration Plateau of the R ECO response to GPP Root Mean Square Error from Jannssen and Heuberger (1995) Respiration at reference temperature for TP Model Savanna Shrublands Total soil stock (0-50 cm) Soil Water Content Constant temperature from Lloyd and Taylor (1994) at 46.02°C Air temperature o r R e v i e w O n l F o r R e v i e w O n l F o r R e v i e w O n l Temperature and F y y y
		52 54

Table A IV-Statistics of the sign test between the Pearson's correlation coefficient calculated between residuals of TP Model and GPP computed using FLUXNET partitioning (Reichstein et al., 2005) and Lasslop's partitioning (Lasslop et al., 2010). In the third colums NS means that the median is not significantly different to 0 while * means a significance level of p<0.05. Median of diff. represent the median of differences of two populations, p the level of significance, df the degree of freedom (i.e. number of sites (n) -1). The definitions of different PFTs are: evergreen needleleaf forest (ENF), deciduous broadleaf forest (DBF), grasslands (GRA), croplands (CRO), savannah (SAV), shrublands (SHB), evergreen broadleaf forest (EBF), mixed forest (MF), wetland (WET).

 A

	PFT	p	Median of		df
			Diff
	ENF	0.678	0.007	NS	25
	DBF	0.774	0.001	NS	14
	GRA	0.424	-0.015	NS	14
	CRO	<0.05	-0.050	*	8
	SAV	0.063	-0.064	NS	4
	SHB	0.999	0.015	NS	4
	EBF	0.688	0.046	NS	6
	MF	0.999	-0.022	NS	7
	WET	0.999	0	NS	2
	All	0.1875	-0.009	NS	92

Global Change Biology

Appendix List

APPENDIX I -Site Table . ID, Name, country, belonging network, coordinates PFT, climate and LAI MAX of the sites used in the analysis. Climate abbreviations follow the Koeppen classification [START_REF] Peel | Updated world map of the Köppen-Geiger climate 1068 classification[END_REF]. Networks are described in www.fluxdata.org APPENDIX II -Site characteristics derived from the FLUXNET database. R 0 is the reference respiration estimated with the LinGPP model formulation, LAI is the maximum seasonal leaf area index of the ecosystems (understorey and overstorey), LAI MAX,o is the maximum leaf area index of the solely overstorey, SoilC is the total soil carbon content, Age is the stand age, T mean is the annual average mean temperature, Ndepo is the total atmospheric nitrogen deposition derived as described in the method section. Sites with (*) in the column dist (disturbance) represent sites with recent disturbance according to what reported in the FLUXNET database. [START_REF] Roupsard | Partitioning energy and evapo-transpiration 1123 above and below a tropical palm canopy[END_REF] APPENDIX IV -Discussion of the 'spurious' correlation between R ECO and GPP.

APPENDIX I -Site table

To understand whether our results were affected by the 'spurious' correlation between GPP and R ECO as reported in FLUXNET (GPP FLUX ) we also perform the analysis using a 'quasi'-independent

Reco and GPP estimates as described by [START_REF] Lasslop | Separation of net ecosystem exchange into 992 assimilation and respiration using a light response curve approach: critical issues and global 993 evaluation[END_REF] (R ECO-LASS and GPP LASS, ). The method by [START_REF] Lasslop | Comment on Vickers et 988 al.: Self-correlation between assimilation and respiration resulting from flux partitioning of 989 eddy-covariance CO2 fluxes[END_REF] do not compute GPP as a difference, but derive R ECO and GPP from quasi-disjoint NEE data subsets. Hence, if existing, spurious correlations was minimized. The 'TP Model' was optimized against R ECO-LASS and GPP LASS and the Pearson's correlation coefficient between 'TP Model' residuals and GPP LASS was calculated (r TPModel-GPPLASS ) at each site and for each PFT.

At each site we compared the correlation between 'TP Model' residuals and GPP derived exploiting the FLUXNET database (r TPModel-GPPFLUX ) with the r TPModel-GPPLASS . The comparison was conducted by using the two sample paired sign test [START_REF] Gibbons | Nonparametric Statistical Inference, 4 th Edition, Marcel 902[END_REF]. We test the null hypothesis that the median of the difference between two samples is zero, for a 5% significance level. The sign test was selected instead the t-test because avoids: (i) the normal distribution assumption; and (ii) distribution symmetry.

The paired sign test between r TPModel-GPPFLUX and r TPModel-GPPLASS indicates that the median for the differences of the populations is not statistically different from 0 (p = 0.187) confirming that the bias observed in the purely climate driven model it is not imputable to a 'spurious' correlation between Reco and GPP introduced by the partitioning method used in the FLUXNET database. The differences are negligible also if we consider each PFT separately as depicted by the box-plot in Fig.

A-I and in Tab.

A-IV.

Once the best model formulation including GPP as driver is selected we also compared the parameters of the 'LinGPP' model formulation (i.e. best model selected by the consistent Akaike Information Criterion, cAIC in Table 1) estimated using the GPP and R ECO from FLUXNET and R ECO- LASS and GPP LASS . The statistics in fitting were also compared. The results are summarized in the boxplot in Fig. AII in which k 2 , R 0 and the main statistics in fitting (EF and RMSE) were schematically reported. These results showed that using the two different datasets the results are similar and the overall picture drawn using the Lasslop's method and the FLUXNET database is the same.