Deciphering the components of regional net ecosystem fluxes following a bottom-up approach for the Iberian Peninsula

N. Carvalhais, M. Reichstein, G. J. Collatz, M. D. Mahecha, M. Migliavacca, C. S. Neigh, E. Tomelleri, A. A. Benali, D. Papale, and J. Seixas Departamento de Ciências e Engenharia do Ambiente, DCEA, Faculdade de Ciências e Tecnologia, FCT, Universidade Nova de Lisboa, 2829-516 Caparica, Portugal Max-Planck-Institut für Biogeochemie, P.O. Box 10 01 64, 07701 Jena, Germany NASA Goddard Space Flight Center, Greenbelt, Code 614.4, Greenbelt, MD 20771, USA Remote Sensing of Environmental Dynamics Laboratory, DISAT, University of Milano-Bicocca, Milano, Italy Science Systems Applications Inc., Lanham, Maryland, 20706, USA

nitude and sign of most of the trends.However, by removing the time series of pure model recovery from the time series of the overall fluxes, we are able to retrieve estimates of interannual variability and trends in net ecosystem fluxes that are quasi-independent from the initial conditions.This approach reduced the sensitivity of the net fluxes to initial conditions from 47% and 174% to −3% and 7%, for strong initial sink and source conditions, respectively.
With the aim to identify and improve understanding of the component fluxes that drive the observed trends, the net ecosystem production (NEP) trends are decomposed into net primary production (NPP) and heterotrophic respiration (R H ) trends.The majority (∼97%) of the positive trends in NEP is observed in regions where both NPP and R H fluxes show significant increases, although the magnitude of NPP trends is higher.Analogously, ∼83% of the negative trends in NEP are also associated with negative trends in NPP.The spatial patterns of NPP trends are mainly explained by the trends in f APAR (r = 0.79) and are only marginally explained by trends in temperature and water stress scalars (r = 0.10 and r = 0.25, respectively).Further, we observe the significant role of substrate availability (r = 0.25) and temperature (r = 0.23) in explaining the spatial patterns of trends in R H .These results highlight the role of primary production in driving ecosystem fluxes.

Introduction
The quantification of terrestrial net ecosystem fluxes is of significant importance to the understanding of the global carbon cycle (e.g., Heimann and Reichstein, 2008) and has been the subject of active research (e.g., Ciais et al., 2000;Piao et al., 2009b).In this regard, model-data synthesis approaches have focused on the improvement of ecosystem flux estimates through model structure development (e.g., Richardson et al., 2006), estimation of parameters (e.g., Knorr and Kattge, 2005) and initial conditions (e.g., Braswell et al., 2005;Carvalhais et al., 2008), adjustments in state variables (e.g., Jones et al., 2004) and sensitivity to forcing variables (e.g., Abramowitz et al., 2008).The wide assessment of uncertainties in the different modeling components can further be integrated in bottom-up approaches and should be a robust indicator of the overall methodological uncertainties (Raupach et al., 2005).The further integration of different models in ensemble approaches establishes the boundaries of prognostic scenarios for future climate conditions (e.g., Sitch et al., 2008;van den Berge, 2010).
In process-based biogeochemical modeling the estimation of carbon fluxes is dependent on the magnitude of prior ecosystem carbon pools.The initial estimates of ecosystem pools are usually prescribed by long initialization routines that drive models to equilibrium between C uptake and efflux from the ecosystem (e.g., Morales et al., 2005;Thornton et al., 2002).These routines, also called spin-up runs, rely on consecutive model iterations for periods varying from hundreds to thousands of years of simulations with average climate datasets (e.g., Thornton and Rosenbloom, 2005).More sophisticated approaches include additional transient runs in which model drivers embody the trajectories of land cover change, climate and atmospheric CO 2 concentrations since the beginning of the industrial revolution until present (Hurtt et al., 2002;McGuire et al., 2001;Zaehle et al., 2007) or incorporate post disturbance recovery dynamics (Masek and Collatz, 2006).But, ecosystem carbon pools are rarely initialized at non-steady state conditions for regional or larger scale simulations.In general, data availability is sparse and there is reluctance to constrain models with datasets that are not harmonized with the spatial resolution of models or to match conceptual carbon pools of models with in situ measurements (e.g.regarding soil carbon pools, Trumbore, 2006).For bottom-up approaches, site level evaluations help develop confidence in model structure and model parameters.However, the dependence of ecosystem pools on site history of past climate, management and disturbance regimes hampers our ability to regionalize the initial conditions of carbon pools.Consequently there is a certain degree of arbitrariness in the initial estimates of ecosystem carbon pools which contributes to additional uncertainties to net ecosystem fluxes estimates.
The response of net ecosystem production (NEP) to climate variability is a result of the separate responses of the component fluxes that yield NEP; photosynthesis and respiration.Environmental drivers can stress or boost simultaneously individual processes that remove or emit carbon into the atmosphere through photosynthesis or respiration, respectively.Hence, changes in net primary production (NPP) can be directly associated with changes in temperature or water availability conditions following mechanistic reasoning (e.g., Haxeltine and Prentice, 1996).Analogously, changes in temperature (e.g., Rey and Jarvis, 2006) and soil moisture (e.g., Orchard and Cook, 1983) have been shown to drive changes in heterotrophic respiration (R H ) patterns.Understanding how these components of faster and slower carbon fluxes respond to climate conditions provides the mechanistic understanding of NEP variability.The full accounting of the regional ecosystem carbon fluxes would imperatively include the effects of disturbance events (e.g.fire, insect outbreaks) and management regimes (e.g.grazing, logging); which summed up to NEP would yield the net biome production (NBP) (e.g., Chapin et al., 2006).Although our approach attempts to be comprehensive, we do not consider the dynamics of NBP, which would require further model evaluation and parameterization efforts.
This study builds directly on two previous studies that explored the impacts of initial steady state assumptions in inverse parameter optimization exercises conducted at site level (Carvalhais et al., 2008(Carvalhais et al., , 2010)).At the site level these studies showed that initial steady state assumptions tend to force neutral NEP estimates which tend to bias parameters and inflate parameter uncertainty.These issues can be circumvented through the relaxation of initial conditions in ecosystem pools using empirical approaches.Here, we follow a bottom-up approach to investigate the role of initial conditions on temporal trends and inter-annual variability in net ecosystem fluxes during a time span of 25 years for the Iberian Peninsula (IP).Ecosystem fluxes are estimated using the Carnegie Ames Stanford Approach (CASA) biogeochemical model (Potter et al., 1993;Randerson et al., 1996).NPP is estimated independently from vegetation pools (Monteith, 1972) and the effect of initial conditions on the ecosystem fluxes can be explored solely through the soil carbon pools (Carvalhais et al., 2008).The modeled time series of ecosystem fluxes are consequently a function of the assumptions about the initial conditions and the climate and phenology drivers.Here, we examine how modeled ecosystem responses exclusively induced by the model drivers can be independent from the initial conditions.We follow a mechanistic approach to remove the dynamics of recovery from non-equilibrium to explore the differences in driver induced trends for different initial conditions.Further, we decompose NEP into NPP and R H fluxes to evaluate the sensitivity of regional ecosystem fluxes to the model drivers.The optimized parameter vectors are then upscaled according to the plant functional type and phenology and climate regimes to the entire Iberian Peninsula.

Materials and methods
This section describes the bottom-up approach from site level optimization to regional flux estimations and analysis, which is summarized in Fig. 1.

Eddy flux sites and data
The European network of eddy-covariance measurements sites from the Carboeurope Integrated Project extends from Mediterranean to boreal ecosystems (http://www.carboeurope.org)and has been supporting broad research on processes of mass and energy transfer at the ecosystem level (Aubinet et al., 2000) throughout different vegetation types (Ciais et al., 2010;Luyssaert et al., 2009).Here, we rely on a selection of 33 sites that includes a significant diversity of ecosystems representing a significant range of climate regimes and net ecosystem flux magnitudes (Table 1).These sites were selected based on: (1) their representativeness of the ecosystems present in the Iberian Peninsula, though not necessarily located in this region; (2) a minimum data availability of bi-weekly flux integrals with more than 80% of the half-hourly original data or after gap-filling with high confidence (Category A in Reichstein et al., 2005).In situ measurements of climate variables and MODIS remotely sensed normalized difference vegetation index (NDVI, Huete et al., 2002) were used to drive site level simulations in the inverse model optimization.

The CASA model
CASA is a process-based biogeochemical model that estimates net ecosystem production (NEP) fluxes as the difference between NPP and R H (Field et al., 1995;Potter et al., 1993).The estimates of NPP follow the radiation use efficiency approach of Monteith (1972): where f APAR is the fraction of photosynthetically active radiation absorbed by vegetation; PAR is the amount of photosynthetically active radiation; and ε is the light use efficiency, which is calculated by down-regulating maximum light use Fig. 1.Workflow of the current analysis: from site level observations to regional dynamics of ecosystem C fluxes.The model-data integration followed at site level yields optimized parameter vectors used for regional simulations.The propagation of parameter vectors is dependent on the classification results of each IP gridcell according to site level characteristics, which support the sensitivity analysis to the initial conditions.Ultimately, the evaluation of trends in ecosystem fluxes quasi-independently from initial conditions allows the decomposition of modeled ecosystem fluxes into component fluxes and investigating the underlying dynamics of NEP trends.The numbers in gray boxes indicate the sections of the manuscript focusing on the different aspects of each step: materials and methods (2.x); results (3.x); and discussion (4.x).
efficiency (ε * ) via the effect of temperature and water availability stress factors (T ε and W ε , respectively): where T ε follows an asymmetric bell shaped curve that peaks at an optimum temperature (T opt ), dropping faster above than below T opt , and W ε responds linearly to the ratio between actual and potential evapotranspiration.Here, the actual net photosynthetic capacity (f APAR • ε) considers both canopy structural components -via f APAR -and instant climate effects -in ε.However, changes in the f APAR also implicitly include the effects of environmental conditions in the canopy's structure, which hampers the isolation of the pure climate effects on NPP.For site level and regional simulations we estimated f APAR as the mean of the two linear scaling methods with the simple ratio and the NDVI, individually per plant functional type (PFT), according to Los et al. (2000).We computed leaf area index (LAI) as varying exponentially with f APAR, or linearly for clustered PFTs (e.g.coniferous trees and shrubs) or as a combination of the previous two for clustered and evenly distributed vegetation types, following Sellers et al. (1996).The carbon assimilated by vegetation is partitioned between the different vegetation pools according to the dynamic allocation scheme of Friedlingstein et al. (1999) which adjusts the investments in leaves, stems and roots towards the most limiting resource (water, light and mineral nitrogen).Carbon is then transferred from living vegetation pools to soil pools through leaf litter fall, root and wood mortality (Potter et al., 1993;Randerson et al., 1996).The cycling of carbon between the different soil C pools follows a simplified version of the CENTURY model (Parton et al., 1987).R H is estimated as the integral of decomposition from the different soil C pools: where each pool i is characterized by a different turnover rate k i that is regulated by the effect of temperature (T s ) and water availability conditions (W s ) (Potter et al., 1993).The effect of temperature in decomposition (T s ) follows an exponential Q 10 -type response; while the effect of soil water conditions (W s ) is estimated as a non-linear function of the ratio between water availability (soil moisture plus precipitation) over the potential evapo-transpiration, following a linear positive response from 0 to 1, with a maximum set for values between 1 and 2, and then decreasing gradually as excess conditions of water supply increase (Potter et al., 1993).The carbon content of each pool (C i ) results from the integrated transfers between the different pools, which is further regulated by the carbon assimilation efficiency of microbes (M ε ).CASA assumes nine soil carbon pools (P = 9) that circulate carbon inputs from vegetation debris and have different turnover rates, ranging from monthly (metabolic litter pools) to centennial (old pools) scales.According to the model structure, temperature (T s ) and water availability (W s ) controls on the instantaneous turnover rates equally affect all soil C pools.Higher inertia is seen in the slower turnover soil carbon pools, which take longer to respond to changes in climate conditions and changes in vegetation states -and vice versa.Hence, these pools take longer to reflect a perturbation as well as to recover from it.In general, the robustness of the CASA model has been corroborated by its wide application in studies that range from ecosystem to global scales (e.g., Carvalhais et al., 2008;Randerson et al., 2002) and that focus on different ecological and biogeochemical issues (e.g., Potter et al., 2001;van der Werf et al., 2003).

Inverse model parameter optimization
The optimized parameters in the CASA model are responsible for governing the responses of NPP and R H to environmental conditions (temperature and water availability) and maximum energy mass conversion rates (maximum light use efficiency, Monteith, 1972;Potter et al., 1993).In addition, a parameter is introduced that adjusts the initial conditions of microbial and more recalcitrant C pools after the spin-up routinesη -and regulates the ecosystem's initial distance to equilibrium (Carvalhais et al., 2008).In this regard, an equilibrium situation would be easily achieved by the faster pools while non-equilibrium conditions would be sustained when lower turnover carbon pools are adjusted.Therefore the strategy to explore non-equilibrium conditions in the recalcitrant pools was chosen.The inclusion of microbial pools in the set of pools affected by η was motivated by the effective role on model performance and the structural dependencies from the slower pools.This relaxation of the steady state assumption gives a structural flexibility to the model that was shown to reduce parameter uncertainties and biases inherent from wrong structural model assumptions (Carvalhais et al., 2008).Additional non-equilibrium conditions can be observed in the vegetation pools which may influence directly gross primary production and autotrophic respiration (R A ) and indirectly heterotrophic respiration fluxes (through transfer of carbon from vegetation to soil reservoirs) (e.g., Nabuurs et al., 2004).CASA relies on the implicit calculations of autotrophic respiration as a fixed fraction of NPP, and estimates of NPP rely on inputs of f APAR (Potter et al., 1993), which minimizes the impacts of nonequilibrium conditions of vegetation on NPP.These implicit treatments of R A do not degrade the CASA model's ability to simulate net ecosystem fluxes (Carvalhais et al., 2010).The model performance is also not different between an exclusive consideration of nonequilibrium conditions in soil versus nonequilibrium conditions in soil and vegetation carbon pools, revealing that η could compensate for the indirect effects of nonequilibrium conditions in R H (Carvalhais et al., 2010).Given the amenability of both representations and the straightforward treatment of η in the current exercise we opted to solely consider nonequilibrium in soil level carbon pools.However, due to inter-annual variability in climate and ecosystem fluxes, the accurate estimation of non-equilibrium states may be hampered because of the short time series of fluxes, or periods that do not represent the development stage of the observed ecosystem, yielding a non representative η.In principle this problem in not exclusive to the determination of η since such representativeness issues can be in the basis of biases in the estimation of other parameters as well.Also the different historical backgrounds within the selection of sites imply different past dynamics behind the current ecosystem states and fluxes observations.However, the application of this empirical approach -through η -has shown robust model performances and reduced biases N. Carvalhais et al.: Upscaling ecosystem C fluxes for the Iberian Peninsula in parameter estimation for different sites with different historical backgrounds (Carvalhais et al., 2008).The parameter optimization method consisted on the minimization of the sum of residual squares between eddy-covariance measurements and model estimates of biweekly NEP fluxes using the Levenberg-Marquardt algorithm (Draper and Smith, 1981).

Upscaling of model parameters
The upscaling of model parameters aims at attributing parameter vectors optimized at site level, on a per-pixel basis, to the whole Iberian Peninsula.The conceptual idea here is that, within an ecosystem type or PFT, each pixel p in the map would be treated according to its similarity with eddy covariance site j , and parameterized accordingly.The assignment of a parameter vector to a given pixel p is supported by a nearest neighborhood classification of the climatic and phenological conditions: to a given pixel p, the parameter vector S corresponding to site j (S j ) is applied when the climate and phenological characteristics of p are closer to j 's than to any other site's from the same PFT.So, S p = S j when d * p,j = argmin j =1,...,N d p,j finds the minimum distances between climate and phenological characteristics between a pixel p and an eddy covariance site j calculated as d p,j = 1 − NS V p ,V j .Here, V p and V j are vectors containing the normalized biweekly time series of mean air temperature, precipitation and solar radiation for the period 1960to 1990;andmean NDVI between 1982 and2005.NS is the Nash-Sutcliffe coefficient: which quantifies the relative association between two vectors relative to the association between the reference vector and its mean (the nominal situation) (Janssen and Heuberger, 1995;Nash and Sutcliffe, 1970), where N is the length of the vectors V p and V p (N = 96) and i the column index of these vectors.In addition to providing a measure of association, the NS also measures the agreement between two vectors (proximity to the 1:1 line).Consequently, the climaticphenological distance measure (from here on identified as CP d ) reflects both the proximity in the magnitude and seasonality of climate and phenology drivers between site level observations and the regional datasets for the IP.The results associate one site of each PFT per pixel which means that if a given PFT is present in that pixel the respective parameter vector corresponds to the site where the climate and phenological characteristics are closer to the pixel's.

Data for spatial runs
CASA requires drivers related to climate (temperature, precipitation and solar radiation) as well as to biophysical properties of vegetation (NDVI or f APAR and LAI).Due to the diagnostic nature of CASA, the temporal extent of the model runs for the Iberian Peninsula was bounded by the longest remotely sensed NDVI record available: the Global Inventory Modelling and Mapping Studies (GIMMS) NDVI (Tucker et al., 2005).The option of using the biweekly MODIS NDVI products (2 by 2 window of 250 m × 250 m pixels) for site level evaluation of the model and the GIMMS NDVI datasets (8 km × 8 km pixels) for the regional runs was based on two assumptions: (1) the eddy-covariance footprint depends on the local conditions (Göckede et al., 2008) and height of the tower (Barcza et al., 2009) but generally is not larger than 1 km 2 , making the 8 km by 8 km areas too large for representation of the eddy-covariance data; and (2) the MODIS and GIMMS NDVI products are comparable (Tucker et al., 2005) (Rodell et al., 2004).Precipitation was extended using 0.5 degrees datasets from the University of Delaware (Matsuura and Willmott, 2007).Every dataset was spatially interpolated to the GIMMS NDVI 8 km by 8 km grid following a weighted average considering a non-linear distance to the four neighboring pixels Zhao et al. (2005).We used linear interpolation to downscale from monthly to biweekly periods.Research, 1985).The proportion of each PFT within every 8 km by 8 km pixel is defined by the CORINE land cover map (Bossard et al., 2000).The fraction of forest PFTs is defined by the tree cover value of the MODIS vegetation continuous fields at 1 km (Hansen et al., 2003) bounded by the class intervals defined in each CORINE class.

Regional model runs for a range of initial conditions
The CASA model estimates of ecosystem fluxes for the Iberian Peninsula are performed on a PFT basis: the parameter vector used per pixel per PFT originates from the upscaling exercise above.The fraction of each PFT within a grid cell is specified according to the land cover map and remote sensing datasets, for which the carbon fluxes are simulated independently from each other, with no interactions competing for any limiting resources (light, water or nutrients).The ecosystem fluxes per pixel p (e.g.NEP p ) are estimated by integrating all the individual PFT estimates (e.g.NEP p,PFT ) weighted by the fraction of each PFT inside that pixel p (f PFT p ): for each PFT the model is always spun up with a mean yearly dataset until steady state after which a range of nonequilibrium conditions is forced on the more recalcitrant and microbial soil carbon pools.The prescription of different distances to equilibrium is obtained by multiplying these pools by the scalar η after the spin-up and then running the model forward.We create a range of initial conditions, from significant sinks (η 1) to sources (η 1) by changing the η values between 0.01 and 2 (in 0.1 increments).The ensemble of runs obtained allows evaluating the impacts of different distances to equilibrium in the inter-annual variability and temporal trends in net ecosystem fluxes.

Decoupling the drivers effects on ecosystem fluxes from the initial conditions
It is recognized that the steady state estimates of the soil carbon pools are a function of the model parameterization and environmental drivers prescribed for the spin-up (Andrén and Kätterer, 1997).Hence, for the same model parameterization and drivers, any prescribed distance of the soil C pools to equilibrium (η) yields a dynamic recovery response towards equilibrium with the simulation drivers.The recovery from a prescribed η can yield conditions different from the steady state as a response to the simulation's climate and/or phenology time series.It is essential to distinguish or isolate the effects of the drivers from the recovery dynamics on the ecosystem fluxes.Here, we opted for removing the recovery dynamics by performing parallel model runs with constant drivers.In these runs the climate and phenology drivers were identical to the spin-up runs datasets: the mean year of the complete 25 year time series.Further, we prescribed the exact same parameterization and η scalars used in the regular forward model runs.To obtain a climate-phenology, not recovery, driven NEP time series (NEP D , Table 2) we subtracted the fluxes estimated with the constant drivers (NEP K , Table 2) from the NEP time series: With this procedure we can remove the variance and trajectories of NEP resulting from the soil C pools recovery and investigate the responses of fluxes to the climate and phenological drivers.

Sensitivity analysis of net ecosystem fluxes to equilibrium assumptions
The sensitivity of the regional decadal fluxes to the initial conditions is assessed by evaluating the changes in interannual variability (IAV) and the trends in ecosystem fluxes computed assuming different distances to equilibrium.The reference scenario is always the time series of ecosystem fluxes estimated in equilibrium (η = 1, for which NEP is defined as NEP eq ).Inter-annual variability is defined as the variance of the annual ecosystem flux integrals over all the years in the analyzed period: 1982 to 2006.The seasonal cycle is removed from the time series for the estimation of temporal trends on the ecosystem fluxes.The seasonal cycle is estimated on a pixel-by-pixel basis using a local variant of Singular System Analysis, as originally introduced by Yiou et al. (2000), as described by Mahecha et al. (2010).This method allows for detection of a highly phase and amplitude modulated seasonal cycle.Based on the deseasonalized time series, the magnitude of the trends is calculated by the Sen slope (Sen, 1968), which is considered a robust estimator of the magnitude of a monotonic trend (Yue et al., 2002).The significance of these trends in the ecosystem fluxes time series is evaluated through the Mann-Kendall test (Kendall and Griffin, 1975;Mann, 1945).The Mann-Kendall test is a non-parametric rank based method that has been widely used to detect the presence of monotonic trends in environmental variables (e.g., Burn et al., 2004;Hamed and Rao, 1998;Kahya and Kalayci, 2004), given its robustness to outliers and because no assumption on data distribution are required.

Decomposition of ecosystem fluxes
The observation of positive or negative trends in NEP can be due to different processes, since NEP is a balance between net C assimilation (NPP) and emission (R H ). Consequently, a positive trend can be due to increases in NPP and/or decreases in R H , or can be caused by equal signed trends in NPP and R H but where the magnitude of the NPP slope is higher than the R H slope.On the other hand, a negative trend can be caused by the opposite reasons.Here, we propose to decompose the NEP trends into NPP and R H trends by evaluating the latter independently and mapping them in a scatter plot against each other (Fig. 2).
N. Carvalhais et al.: Upscaling ecosystem C fluxes for the Iberian Peninsula Further, the decomposition of NPP and R H trends into their main drivers may clarify the mechanisms behind significant trends in both fluxes.The drivers for modeled NPP are remotely sensed NDVI and climate variables, although ultimately their influences are expressed in terms of f APAR and temperature and water availability stress factors trends: T ε and W ε , respectively.Although the relationships between observed variables and the NPP scalars are non-linear, the scalars themselves all share the same dimensional characteristics -represent fractional properties ranging between 0 and 1 -and yield identical effects on NPP: a 0.1 increase in any of the temperature or water scalars or in f APAR equally yield a 10% increase in NPP (considering the same maximum light use efficiency and solar radiation conditions).In this regard, the isolation of the climatic effects from the f APAR time series is not mechanistically feasible.Hence, the partitioning of NPP drivers is understood as between the changes in phenology (f APAR, which also integrate effects of climate in the canopy's structure) and the effects of climate in the instantaneous light use efficiency (ε).
For R H , the changes in climate drivers may produce trends in the temperature (T s ) and water (W s ) stress scalars but R H is also influenced by substrate availability.The vegetation pools are the main sources of carbon for R H through the transfer of live biomass to detritus via litter fall, wood and root mortality.The carbon transferred to the soil litter pools then cycles through different soil pools mediated by microbial decomposition releasing C through R H .The changes in substrate can be due to changes in inputs from the vegetation pools or due to changes in rates of consumption of the substrate.The latter produces a negative feedback on R H in response to environmental conditions: favorable conditions for R H increase decomposition reducing substrate availability, and vice versa.The connection between these factors hampers distinguishing between trends in soil C availability and environmental conditions.Hence we choose to focus on the carbon available for decomposition from the vegetation, which equals the sum of the vegetation pools magnitude weighted by their respective turnover rates.This approach is consistent with the model formulation where the changes in the vegetation pools are directly proportional to the changes in inputs to the soil.Here, the carbon transfers to the soil pools are driven by constant turnover rates for each vegetation pool.The calculations of trends for the carbon pools were estimated relative to the pools' mean; hence the result is a relative trend, in fractional units, consistent with the units of the environmental scalars.

Model optimization at site level
Throughout eddy covariance sites the parameter optimization yields statistically significant correlations between observations and model simulations (α < 0.0001).The model efficiency (MEF) at site level reflects an overall good agreement between simulations and observations (MEF closer to 1, Table 1) and the median of all site level MEF results is 0.79.Overall, the MEF of DBF, MF and SHR is significantly higher than the MEF of the other PFTs (1-way ANOVA, α < 0.0005).Although the mean model performance is higher in temperate fully humid climates, the difference is not statistically significant from the other climatic regimes (1-way ANOVA, α > 0.25).
The optimized parameters differ among PFTs (Table 3) but the variance explained by this factor is modest compared to the sum of squared errors and the differences between groups are hardly significant (Table 4).The one exception that shows significant differences with PFT is the optimum temperature for NPP (T opt ).In this case, estimated T opt for the grasslands are significant lower than the observed values for other PFTs (α < 0.0001).The removal of grassland sites from the analysis yields non significant differences among groups.Although the optimized parameters show no significant differences among the different climate regimes, this factor often contributes an important part of the explained variability (Table 4).
Table 3. Results of the site level parameter optimization organized by plant functional type (PFT): evergreen broadleaf forests (EBF), deciduous broadleaf forests (DBF), mixed forests (MF), evergreen needleleaf forests (ENF), savannah type ecosystems (SAV), grasslands (GRA), shrublands (SHR) and croplands (CRO).The optimized parameters are: maximum light use efficiency (ε * , [gC MJ −1 APAR]), optimum temperature for photosynthesis (T opt , [ • C]), the sensitivity of light use efficiency to water stress (B wε , [unitless]) and the effect of temperature (Q 10 , [unitless]) and water availability (A ws , [unitless]) on heterotrophic respiration (for a detailed description see Carvalhais et al., 2008).Values in parenthesis report the standard deviation around the mean parameters calculated from the optimized parameters and standard error for all the sites of each PFT.Overall, the optimization results at site level show significant confidence in the CASA model performance for sites representing the PFTs of the IP domain, where the average MEF weighted by PFT fraction is 0.73.

Upscaling parameter vectors for the IP
The map of the distances in climatological and phenological space (CP d ) between individual pixels and the eddy covariance sites shows that in general the IP domain is well represented by the group of eddy covariance sites used in this analysis (Fig. 3).
For the entire Iberia, 94% of the pixels show a weighted CP d between pixel's and site characteristics below 1 -the nominal situation -meaning that, on average, the association between the chosen site for a given pixel is better than just considering the average of the target pixel (Schaefli and Gupta, 2007).The median distance (CP d ) is 0.56, indicating a significant association between the drivers from the eddy covariance sites and from the regionalized datasets.
The NW region shows systematic higher CP d to the eddy covariance sites characteristics (Fig. 4) and emphasizes the low local representativeness for all PFTs in the region.Mixed forest, crops and grasslands each contribute about 20% of the total land cover in this region; hence, substantial improvements in the region's representativeness would be achieved by including in our analysis eddy covariance sites monitoring such PFTs with more similar phenology and climate characteristics.
Clearly, these results are strongly bounded to methodological aspects such as site selection, choice -and sourcesof explanatory variables as well as to the metrics used for distance calculation.Methodological improvements focusing on any of the above aspects may yield different spatial patterns on site representativeness.

Changes in Inter-Annual Variability (IAV)
Inter-annual variability in NEP is higher for the farthest nonequilibrium initial conditions (observed at very low and high η prescriptions, Figs. 5, 4a).The modeling results show regional increases of 47% and 174% in IAV for the lowest and highest η, respectively, considering the average fluxes for the whole domain of the IP (Fig. 5b, inset plot line a).The changes in IAV show significant increases in spatial variability with distance to equilibrium (Fig. 5a).These are driven by the changes in the IAV of R H , which are strongly dominated by the recovery of the carbon pools (Fig. 6a).
The removal of C pools recovery from the prescribed η (Eq.6) significantly reduces the changes in IAV caused by distance to equilibrium (Fig. 5b).The highest changes are still observed for extreme η values although these are modest compared to the previous results, with regional changes of −3% and 7% for η of 0.01 and 2, respectively.Nevertheless, the spatial variability of IAV changes still increases with distance to equilibrium (Fig. 5b).The changes that still exist between the ratio of IAV of NEP D and the IAV of NEP eq for different η values are driven by changes in the IAV of R D H , which respond proportionally to the varying magnitude of soil carbon pools for different η (Fig. 6b).A high η increases the carbon pools which increase the magnitude of the IAV in R D H and its role in the IAV of NEP D .Conversely, low η values reduce the role of R D H in the IAV of NEP D , and increase the role of the IAV in NPP D in the IAV of NEP D .Consequently, the effects on the IAV of NEP D are stronger for larger carbon sources (η 1) and sinks (η 1) (Fig. 5b), but these do not follow the R D H patterns.Overall, the results show that regionally both approaches reveal significant differences between each other (Fig. 5b, inset plot), which originate exclusively from the recovery of carbon pools from the prescribed initial disequilibrium conditions.The similarity between the IAV of the different NEP D suggests quasi-independence from the initial conditions.The slight changes observed in the IAV of NEP D stem not from the recovery from the initial conditions but from the carbon pools sizes.

Temporal trends for the IP
The evaluation of trends in NEP shows that the spatial distribution of its magnitudes (Sen slope) is strongly influenced by the prescribed initial distance to equilibrium (Fig. 7a), as expected.On the other hand, the distribution of NEP D trend magnitudes appears to be η invariant (Fig. 7b) and suggests independence from the pools' initial conditions.The spatial distribution of the significant trends in NEP in equilibrium (NEP eq , Fig. 7c) differs only slightly from the map of mean NEP D trends (Fig. 7d).The Mann-Kendall test yields non significant trends in NEP for 55% of the IP, while about 64% of the IP shows non significant trends for NEP D .The area difference of negative trends in the IP between NEP eq and NEP D is about +3%, while the difference in the area of positive trends is about +6% (Table 5).Overall, both approaches agree almost 90% of the time on the type of trend (negative, positive or non significant) and the spatial NS between NEP eq and NEP D is 0.99.Consequently, these differences are considered minor.These results show significant trends in NEP D driven exclusively by climate and phenology.A direct implication of these results is the ability of the approach to detect climate and phenology induced trends that are independent of the initial carbon pools.

Determinants of temporal trends in the IP
Comparing the NEP D trends (NEP D T ) with trends in NPP D (NPP D T ) and in (RH D T ) (Fig. 8) allows us to identify the reasons behind the temporal behaviour of net ecosystem fluxes.The results reveal that ∼ 97% of the positive trends in NEP D are mostly located in west, northwest and northern regions with positive trends in NPP D and in R D H , but the NPP D T slopes are higher that the slopes in RH D T (Fig. 9; Table 6).These correspond to ∼ 73% of the total significant NEP D trends in the IP.Negative trends in NEP D are mainly located in southern central regions and can be divided into: positive trends in R D H and negative trends in NPP D (covering ∼ 53% of the negative trends in NEP D in the IP); negative or positive trends simultaneously in both fluxes, where the magnitude of the slope of RH D T is higher than the slope of NPP D (covering ∼30% and ∼17%, respectively).These results are only slightly different from the trends resulting from the time series of NEP eq (Table 6).Overall, the trends in the IP for R D H and NPP D are predominantly positive.
The NPP D trends are a result of the trends in stress scalars on light use efficiency driven by temperature and water availability (T ε and W ε , T ε T and W ε T , respectively) and with trends in f APAR (f APAR T ).The trends observed in NPP D (NPP D T ) are mainly driven by trends in f APAR (Fig. 10).The spatial partial correlation between the trends in f APAR and NPP D T is significantly higher than between T ε T and NPP D T , 0.79 compared to 0.10, respectively (Table 7).The magni-tude of the W ε trends was found to be an order of magnitude below f APAR T and T ε T (Fig. 11).However, the spatial partial correlation between Wε T and NPP D T is 0.25: significantly higher than the one found for T ε T .Overall, these results suggest a significant role of f APAR time series in the productivity trends for the Iberian region.
Following the analogous procedure, the trends in R D H are compared to the trends observed in the scalars that translate the effect of temperature and soil moisture in R H (T s and W s , respectively) and in the vegetation C pools available for R H (substrate availability).In the Iberian region the substrate availability has a stronger effect than T s on the trends in R D H (Fig. 12).The same holds true when comparing substrate availability with the effect of soil moisture (Fig. 13).However, the spatial partial correlation between trends in shows for a spatial average of all the Iberian Peninsula the ratio between the IAV of NEP and the IAV of NEP eq (black line a) and the ratio between the IAV of NEP D and the IAV of NEP eq (grey line b).The regional differences in IAV for the Iberian Peninsula between both approaches as a function of distance to equilibrium are significant.).The trends in fluxes in each pixel are considered significant when for all prescribed η s the slopes are significant (α < 0.05), systematically positive or systematically negative and its absolute magnitude at least 1.5 gC m −2 yr −2 .The mapped values report the average of all the significant trends.Table 6.Results for the decomposition of NEP trends into NPP and R H trends for NEP D and NEP eq .The partial frequencies are estimated dividing the occurrence of NPP T and RH T combinations by the total occurrence of positive or negative trends (depending on the NEP T column).Significant frequencies are calculated dividing the occurrence of NPP T and RH T combinations by the total of significant trends for the IP region.The fraction of significant trends over the Iberian Region is identified with IP.
Frequency (%) R D H and its drivers' are not very different between substrate availability (0.25) and temperature (0.23), and are lower for the water stress scalars (0.11).These correlations change significantly when considering only negative trends in the NEP, where the substrate availability explains 86% of the spatial variability, emphasizing the role of substrate availability in the R D H trends (Table 7).

CASA model optimization
Globally, the site level optimization supports significant confidence in CASA's performance.In general, the optimized parameters (Table 3) are within published results (e.g., Kätterer et al., 1998;Kirschbaum, 1995;Ruimy et al., 1994) et al., 1980).Most of these sites are C3, with the exception of PT-Mi2 which is a C3/C4 mix.Site level records show the occurrence of high gross primary production (95 percentile) at temperatures below 10 • C. The coherence between observed climate and the optimized values is supported by previous observations of higher photosynthetic rates at low temperatures in C3 plants native to cool growing seasons (Berry and Björkman, 2008;Bunce, 2008).In Mediterranean climates higher photosynthetic activity of grasses (mostly C3) can be observed during seasonally cooler periods when plants grow due to higher moisture availability (e.g.winter or spring).This is consistent with the comparison between T opt of grasses and the mean monthly temperatures for the Iberian Peninsula, where the confidence bounds of T opt include the monthly temperatures observed during colder months.However, from a model-data integration perspective both model drivers as well as model structure can influence the actual values of parameters.The higher uncertainties in cropland parameters suggest the need to improve/adapt model dynamics in agricultural systems, including the prescription of different root to shoot ratios according to crop type (e.g., Bondeau et al., 2007) and explicit harvest events, with above ground biomass removals (e.g., Hicke and Lobell, 2004;Lobell et al., 2006), as well as different management regimes considering crop rotation, fertilization, irrigation and tillage practices.These dynamics were not included in the current model implementation.Nevertheless, the good model performance as indicated by MEF in croplands (Table 1) suggests that the phenology time series acquired from remote sensing captures most of the variability of the NPP.Further, in site level optimizations, the lack of harvest removal of C from the leaf pools that then feeds the soil pools can be approximated by reductions in η.For the purpose of this study, the absence of such dynamics in model simulations is unlikely to change our results.

Upscaling parameter vectors for the IP
The results reveal systematic lower representativeness in the Northwestern IP region which is generally identified as a high productivity area in regional and global studies (e.g., Jung et al., 2008).Representativeness in this region is low for all PFTs (Fig. 4), which reveals that the under-representation  The current approach identifies the limitations in the sample of eddy covariance sites for representing the regional ecosystems in the IP.Similar approaches would be useful in assisting network design or for selection of possible site locations in future network updates.Certainly these are dependent on the goals of the monitoring and modeling exercises.Here, the main factors explaining the variability of parameters are often the plant function type and the climate regime observed at site level (Table 4).These results lend support to the conceptual approach taken here to upscale the CASA model parameters.However, there is still a significant fraction of unexplained variability in the parameters and our bottom-up approach does not address issues related to the variability of parameters within the same PFT or climate regime.Increasing the number of sites and including other factors in the analysis -such as disturbance or management regimes -are two important issues to consider towards more comprehensive upscaling and modeling exercises.In this regard we recognize the importance of prescribing management (e.g., Bondeau et al., 2007) or disturbance (e.g., van der Werf et al., 2003) regimes, as well as the effects of nutrients dynamics (e.g., Zaehle et al., 2010).

Dynamics of ecosystem fluxes induced by climate and phenology
The effects of different initial conditions on the IAV and temporal trends in NEP were higher for the highest departures from equilibrium, as expected.Since the spinup routines were performed with an average climate and phenology dataset, the modeled recovery from the initial perturbations (η) is expected to lead to similar ecosystem states (carbon pools).The farthest positive departures from equilibrium (NEP 0 since η 1) create an initial sink condition that is then attenuated and consequently leads to the steepest negative slopes of NEP trends.In these cases, carbon is accumulated in the soil pools with time, enhancing R H through increases in substrate availability, and therefore decrementing NEP in time.Oppositely, the farthest negative departures from equilibrium (NEP 0 since η 1) force an initial carbon source by increasing the soil pools.The increase in substrate availability boosts R H in the beginning of the simulation which is reduced in time, generating a positive trend in NEP.Consequently, the highest changes in the IAV of NEP are observed when extreme departures from steady state are prescribed.
The extraction of the carbon pools dynamics (NEP K ) from the NEP time series allows the exploration of the impact of non steady state conditions in the purely climate and phenology driven dynamics of NEP (NEP D ).The significant reduction in IAV changes (Fig. 5) and the distribution of the NEP trends with η (Fig. 7a) emphasize the role of the carbon pools dynamics in the estimation of ecosystem fluxes.These results show that simulated net ecosystem fluxes can strongly diverge depending on initial conditions assumptions.Additionally, the trends in the NEP D time series show a strong insensitivity to η revealing independence from the initial estimates of C pools (Fig. 7b).The results exhibit the potential of NEP D to isolate and diagnose climate and phenology driven trends in ecosystem fluxes.By identifying significant and consistent trends in NEP D we're able to map regions of robust trends in the Iberian region.
The differences between the spatial distribution of trends in NEP D and NEP eq are minor and reveal that here the climate and phenology driven dynamics are quasi-independent from the initial steady state conditions (Fig. 7c, d).Using a mean yearly dataset of drivers during the period 1982-2006 ensures the adjustment of carbon pools to the mean of the drivers of the run.The adjustments follow first order dy-namics that are intrinsic to the model structure (CASA's and many other biogeochemical models).This means that we can isolate the effects of unknown initial conditions by removing NEP K from NEP.The result is the retrieval of a time series of NEP D that is quasi-independent from the initial conditions of carbon pools.In the current experiments these independent trends are analogous to the trends in NEP eq because the dataset that drives NEP K is equal to the spin-up dataset.Consequently, as an alternative to spinning up the models until equilibrium this approach has advantages for exploring the effects of drivers on net ecosystem fluxes independently of the initial conditions of pools.
We should add that ultimately the overall absolute NEP trends are not independent from the initial conditions.Although removing carbon pools dynamics from the ecosystem fluxes may allow independence from equilibrium assumptions, such a procedure does not solve the initial state problem and our ability to quantify temporal trends is still limited.Beside model structure the estimated fluxes depend not only on the model drivers -climate and phenological drivers -but also on the parameters estimated at site level.In this regard, the exhaustive accounting of terrestrial C fluxes implies the explicit integration of disturbances and management regimes that were not simulated -especially agriculture and fire -for the estimation of net biome production (Chapin et al., 2006).In addition, the extension of the current exercise to living pools in non equilibrium conditions can be aided by prospective remotely sensed estimations of above ground biomass.

Decomposition of ecosystem fluxes
The comparison between NEP D trends and the trends in its component fluxes, NPP D  behind the positive and negative trends in NEP D .For most of the Iberian region the positive trends in NEP D are associated with positive trends in NPP D and in R D H , although the slope magnitudes are higher for NPP D .These results are consistent with recent modeling studies (Piao et al., 2009a) and with eddy covariance based studies that advocate the significant role of gross primary production in driving NEP (e.g., Baldocchi, 2008;Reichstein et al., 2007).The positive trends in NPP D are more strongly linked to positive trends in f APAR than to the climate effects on light use efficiency.Since f APAR estimates are based on the NDVI datasets, the modeled positive NPP trends are mostly driven by phenological data rather then climate data.The f APAR time series are estimated from NDVI (following Los et al., 2000), hence the trend results demonstrate the role of NDVI time series in driving the trends in NPP (cf.Jung et al., 2008).
The close association between climate and vegetation implies that the NDVI signal itself contains the effects of climate regimes and patterns in the phenological characteristics of the vegetation (e.g., Myneni et al., 1997).In the CASA model the temperature stress scalars are conceptually associated to adjustments in autotrophic respiration costs while the water stress on canopy productivity can be ascribed to reductions in stomatal conductance (Potter et al., 1993).However, the plants' response mechanisms to environmental stress can yield impacts on APAR, ε or both.In this regard, current studies have highlighted the limitations of using NDVI or f APAR to detect instantaneous light use efficiency changes as a function of water and temperature stress (e.g., Goerner et al., 2009;Grace et al., 2007).It is then implicit in the model structure that the environmental effects on ε and on f APAR are complimentary for the estimation of NPP and act at different time scales.Hence, here, only the effects of temperature and water availability on light used efficiency can be decomposed from the primary productivity signal.
The main mechanism behind positive trends in the net ecosystem fluxes originates from increases in primary production (mostly driven by f APAR) that consequently increase the vegetation carbon pools and the substrate availability for heterotrophic respiration.Here, the impact of substrate availability is in general significantly higher than the effects of temperature or soil moisture on the trends in R D H .In areas of negative NEP D trends the partial correlation between the spatial patterns of trends in R D H and substrate availability is significantly higher than in areas of positive NEP D trends (Table 7).In these cases (NEP D T < 0) the increments in substrate availability are not attributable to increases in NPP, since most of the trends in f APAR (89%) as well as in NPP D (83%) are negative.Here, when RH D T > 0, the positive trends in substrate availability (75%) are mainly associated to positive trends in the root pools (88%).These positive trends in the root pools are contrary to the trends observed in NPP D .This apparent contrary behaviour results from increases in the carbon allocated to the below ground vegetation pools caused by negative trends in the water stress scalars of ε * (88%), increasing water stress.The investment in the root pools is associated with the increasing trends in water stresses and is consistent with the dynamic allocation scheme by Friedlingstein et al. (1999).Due to the higher turnover rates of the fine root pools -compared to the wood pools -the changes in the allocation strategies increase the availability of carbon for decomposition.These observations highlight that the contribution of climate for long term changes in heterotrophic respiration cannot be dissociated from the availability of material for decomposition (Trumbore and Czimczik, 2008).

Conclusions
Our bottom-up approach allowed us to investigate the influences of the initial conditions in modeling regional net ecosystem fluxes.Overall, the site level optimization results showed significant confidence in model performance.According to our approach the set of sites selected is, in general, significantly representative of the climate-phenology conditions observed through different PFTs in the Iberian Peninsula.The resulting spatial patterns highlight the Northwestern region as a location of systematic poorer representativeness ability for the selected climate-phenology descriptors.Consequently, the CASA model reveals itself as a robust diagnostic approach to estimate NEP fluxes and the methodology to upscale parameter vectors allows the identification of less represented regions.The presented bottom-up approach emphasizes the relevance of analogous methods for the design of ecosystem monitoring networks.
The influence of the initial conditions on NEP trends and inter-annual variability is significant, as expected.However, we present a method to distinguish between the model intrinsic dynamics following initialization routines and the flux variability only induced by the driver data.Consequently, we are able to investigate the inter-annual variability and trends in fluxes quasi-independently from the initial conditions.The relevance of such approach emerges from the fact that most of the time the initial conditions of regional or global simulations are unknown.
The results for the Iberian Peninsula show that the trends in the net ecosystem fluxes are strongly linked to trends in primary production.The positive trends in the net ecosystem fluxes are observed in western and northern regions with positive trends in both primary production and heterotrophic respiration fluxes.Here, the magnitudes of the trends are higher for the primary production fluxes.These are generally driven by positive trends in f APAR, that drive positive trends in NPP D , vegetation pools and consequently in the substrate availability for decomposition at the soil level.Further, in regions of negative trends in net ecosystem fluxes located mainly in Southern IP, the positive trends in heterotrophic respiration are associated to positive trends in temperature as well as in substrate availability.The positive trends in substrate availability result from changes in carbon allocation strategies that are driven by positive trends in water stress.Overall, these modeling results suggest a strong link between the component processes of net ecosystem fluxes.Further, the significant role of f APAR in NPP D trends and of substrate availability in R D H trends emphasizes that the underlying mechanisms of trends in net ecosystem fluxes are strongly -but not necessarily linearly -associated with primary production and allocation.
Overall, our results show that isolating the time series of ecosystem fluxes from the initial conditions allows the identification of the effects of driver data on flux trends, which suggest significant advantages for the estimation of the sensitivity of ecosystem fluxes to climate drivers.The alleviation of the initialization uncertainty is significantly relevant to long term modeling studies of net ecosystem fluxes.

Fig. 2 .
Fig. 2. Abacus of NEP T decomposition into NPP and R H trends: NPP T and RH T , respectively.

FrequencyFig. 3 .
Fig. 3. Spatial patterns of the representativeness of the eddy covariance sites for pixels of the Iberian Peninsula.The distance measure (CP d ) of each pixel is computed as the average CP d weighted by the fraction of all the sites used for its parameterization.The right side histogram shows the frequency of the observed CP d for all the pixels in Iberia (Eq.4).

Fig. 4 .
Fig. 4. The classification of the Iberian Peninsula according to the optimized eddy covariance sites per PFT (left column) is based on the selection of the closest site to each pixel.The distance (CP d ) maps per PFT show the minimum distance (CP d ) from the set of sites to the IP region (right column).Overall, the spatial classification shows a predominant selection of southern European and Mediterranean sites (middle column).Systematic higher distances (CP d ) in the NW IP are observed throughout PFT (right column).White regions in the left column maps indicate only a residual presence (< 5%) of the respective PFT.

Fig. 5 .
Fig. 5. Influence of distance to equilibrium (η) in the inter-annual variability (IAV) of net ecosystem fluxes (a) and in the IAV of the de-trended NEP fluxes (removing the sole recovery from the C pools) (b).Each vertical rectangular box represents the spatial distribution of the IAV ratio indicated in the yy-axis for Iberian Peninsula (all pixels).Rectangular boxes are bounded by 25 and 75 percentile bottom and top, respectively, while the horizontal line inside indicates the sample median; dashed lines limited by horizontal bars indicate the extent of the remaining data, excluding outliers; dot signs (•) indicate statistical outliers.The inset in (b)shows for a spatial average of all the Iberian Peninsula the ratio between the IAV of NEP and the IAV of NEP eq (black line a) and the ratio between the IAV of NEP D and the IAV of NEP eq (grey line b).The regional differences in IAV for the Iberian Peninsula between both approaches as a function of distance to equilibrium are significant.

Fig. 6 .
Fig. 6.Influence of distance to equilibrium (η) in the IAV of heterotrophic respiration (R H ) fluxes (a) and in the IAV of the de-trended R H fluxes (R D H ), removing the sole recovery from the C pools) (b).Each vertical rectangular box represents the spatial distribution of the IAV ratio indicated in the yy-axis for Iberian Peninsula (all pixels).Rectangular boxes are bounded by 25 and 75 percentile bottom and top, respectively, while the horizontal line inside indicates the sample median; dashed lines limited by horizontal bars indicate the extent of the remaining data, excluding outliers; dot signs (•) indicate statistical outliers.The inset in (b) shows for a spatial average of all the Iberian Peninsula the ratio between the IAV of R H and the IAV of RH eq (black line a) and the ratio between the IAV of R D H and the IAV of RH eq (grey line b).

Fig. 7 .Fig. 8 .
Fig. 7.The distribution of the NEP trends (Sen slopes, gC m −2 yr −2 ) is strongly dependent on η values (a) while for NEP D its distribution appears invariant with η (b).The spatial pattern of the significant NEP trends for equilibrium conditions (c) only slightly differs from the mean trends for NEP D (d).Per pixel, an NEP D trend is only considered significant when for all η used in the simulation the significance level is systematically lower than 0.05 and the trend's signal is consistent (always positive or always negative).The trend is then calculated as the mean of NEP D time series for all η s.Any slope trend is only considered positive or negative when its absolute value is at least 1.5 gC m −2 yr −2 .

Fig. 9 .
Fig. 9. Decomposition of NEP D trends into R D H and NPP D trends (NEP D T , RH D T and NPP D T , respectively -gC m −2 yr −2 ) (a).The spatial distribution of the underlying NPP D T and RH D T dynamics behind NEP D trends can be mapped by colour-coding each NEP D T according to its position in (a) following the diagram in (c).The colour intensity in (b) is proportional to the magnitude of NEP D T .

Fig. 10 .
Fig. 10.Decomposition of NPP D trends (gC m −2 yr −2 ) into T ε and f APAR trends, T ε T and f APAR T , respectively (fractional units per square meter per year) (a).The spatial distribution of the underlying T ε T and f APAR T dynamics behind NPP D trends can be mapped by colour-coding each NPP D T according to its position in (a) following the diagram in (c).The colour intensity in (b) is proportional to the magnitude of NPP D T .

Table 7 .
Partial correlations between the trends in NPP D (NPP D T ) and R D H (RH D T ) and trends in its drivers: f APAR, temperature (T ε ) and water availability (W ε ) effects on ε * (f APAR T , T ε T and W ε T , respectively) for NPP D T ; and substrate availability (SA), temperature (T s ) and soil moisture (W s ) effects on microbial decomposition (SA T , T s T and W s T , respectively) for RH D

Fig. 11 .Fig. 12 .
Fig. 11.Decomposition of NPP D trends (gC m −2 yr −2 ) into W ε and f APAR trends, W ε T and f APAR T , respectively (fractional units per square meter per year) (a).The spatial distribution of W ε T and f APAR T behind the NPP D trends can be mapped by colour-coding each NPP D T according to its position in (a) following the diagram in (c).The colour intensity in (b) is proportional to the magnitude of NPP D T .

Fig. 13 .
Fig. 13.Decomposition of R D H trends (gC m −2 yr −2 ) into W s and substrate availability (SA) trends, W s T and SA T , respectively (fractional units per square meter per year) (a).The spatial distribution of W s T and SA T behind the R D H trends can be mapped by colour-coding each RH D T according to its position in (a) following the diagram in (c).The colour intensity in (b) is proportional to the magnitude of RH D T .

Table 1 .
Identification of the different sites included in the parameter optimization.The sites are divided by plant functional type (PFT) and the percentage of each site per PFT Kottek et al., 2006)wn; as well as the percentage of each PFT in the Iberian Peninsula (PFT in IP).The mean annual temperature (MAT), total annual precipitation (TAP), incoming solar radiation (Rg), climate classification (KG, according toKottek et al., 2006)and net ecosystem production (NEP) are for the time period (Period).The model efficiency (MEF) is reported as well as the total number of data points per site (total N) and the data points used in the optimization (Filtered N).The different PFTs considered are evergreen broadleaf forests (EBF), deciduous broadleaf forests (DBF), mixed forests (MF), evergreen needleleaf forests (ENF), savannah type ecosystems (SAV), grasslands (GRA), . Monthly climate datasets -air temperature, precipitation and solar radiation -at 10 spatial resolution since 1901 until 2000 are available from the Climate Research Unit of the University of East Anglia (Mitchell et al., 2004).The climate datasets were extended until 2006 using pixel level empirical relationships with coarser climate datasets.For temperature and solar radiation we made use of 0.25 degrees datasets from the Global Land Data Assimilation System

Table 2 .
Acronyms used to identify the different ecosystem flux components and temporal signals.

Table 4 .
Kottek et al., 2006) (%) for each optimized parameter for each factor: PFT -plant functional type; climate classification (KG) (according toKottek et al., 2006); MAT -mean annual temperature; TAP -total annual precipitation; and Rg -incoming solar radiation.The results presented: SS -sum of squares; SST -total sum of squares; α -significance test result for each factor.The KG snow climate classes (D) were integrated within the same group and the warm temperate classes aggregated according to precipitation regimes (second letter), to avoid classes with only one site.We removed SHR from the test, since it only has one site, and grouped the climate classification in one class including all snow climates and the other classes per annual temperature and precipitation regime, to avoid having classes with one site.D.F. stands for degrees of freedom.

Table 5 .
Percentage of positive, negative and non-significant trends in NEP eq and NEP D for the Iberian Peninsula.
although two different situations are worth mentioning: (i) the low T opt values found for grasslands; and (ii) the high uncertainties in croplands ε * and T opt .The low T opt values found for grasses are border line or below the 10 • C to 25 • C range for C3 and 30 • C to 40 • C for C4 plants, although some C3 species are quite active at 5 • C (Breymeyer