Effects of climate variability and functional changes on the interannual variation of the carbon balance in a temperate deciduous forest

. The net ecosystem exchange of CO 2 (NEE) between the atmosphere and a temperate beech forest showed a signiﬁcant interannual variation (IAV) and a decadal trend of increasing carbon uptake (Pilegaard et al., 2011). The objectives of this study were to evaluate to what extent and at which temporal scale, direct climatic variability and changes in ecosystem functional properties regulated the IAV of the carbon balance at this site. Correlation analysis showed that the sensitivity of carbon ﬂuxes to climatic variability was signiﬁcantly higher at shorter than at longer time scales and changed seasonally. Ecosystem response anomalies implied that changes in the distribution of climate anomalies during the vegetation period will have stronger impacts on future ecosystem carbon balances than changes in average climate. We improved a published modelling approach which dis-tinguishes the direct climatic effects from changes in ecosystem functioning (Richardson et al., 2007) by employing the semi empirical model published by Lasslop et al. (2010b). Fitting the model in short moving windows enabled large ﬂexibility to adjust the parameters to the seasonal course of the ecosystem functional state. At the annual time scale as much as 80 % of the IAV in NEE was attributed to the variation in photosynthesis and respiration related model parameters. Our results suggest that the observed decadal NEE trend at the investigated site was dominated by changes in ecosystem functioning. In general this study showed the importance of understanding the mechanisms of ecosystem functional change. Incorporating ecosystem functional change into process based models will reduce the uncertainties in long-term predictions of ecosystem carbon balances in global climate change projections.


Introduction
Terrestrial ecosystems assimilate more than ten times the current annual anthropogenic carbon dioxide (CO 2 ) emission through photosynthesis (Beer et al., 2010;Friedlingstein et al., 2010). At the same time, a similar amount of CO 2 is released back to the atmosphere by respiration from soil microorganisms and plants. The difference between these two opposing fluxes determines the net carbon balance of the ecosystem, which varies across time and space (Luyssaert et al., 2007;Stoy et al., 2009;Yuan et al., 2009). Small perturbations in the climate or ecosystem status may alter the equilibrium between photosynthesis and respiration. Whether the terrestrial ecosystems will act as a sink or a source of CO 2 is important, because approximately 25 % of anthropogenic CO 2 emissions are being re-fixed by the terrestrial ecosystems (Friedlingstein et al., 2010;Houghton, 2007). Therefore, understanding the spatiotemporal variability of the ecosystem carbon balance Published by Copernicus Publications on behalf of the European Geosciences Union. and the mechanisms that control it is crucial for assessing the vulnerability of the terrestrial carbon pools under future changing climate conditions (Reichstein et al., 2007;Heimann and Reichstein, 2008).
One important approach to understand the ecosystem carbon dynamics is to investigate the interannual variation (IAV) of net ecosystem exchange of CO 2 (NEE) with long-term eddy covariance measurements (Baldocchi, 2003). By analyzing the year to year variation in NEE under different climatic conditions, the key factors and processes that determine the ecosystem carbon balance may be identified. The measured NEE is the difference between gross primary production (GPP) and total ecosystem respiration (TER) which are both much larger than the net flux. The responses of GPP and TER to climate are complex. Some processes are direct and instantaneous, for instance the light response of photosynthesis and the temperature effects on the kinetics of photosynthesis (Sage and Kubien, 2007) and respiration (Mahecha et al., 2010). However, there are also indirect responses, especially through changes in phenology , canopy structure (Barr et al., 2004;Ibrom et al., 2006) or physiological acclimation (Luo et al., 2001). Many studies have reported enhanced carbon uptake as warming extended the length of growing seasons (Chen et al., 1999;Black et al., 2000;Suni et al., 2003;Hollinger et al., 2004;Churkina et al., 2005;Penuelas and Filella, 2009;Richardson et al., 2009;Pilegaard et al., 2011). Others show that distribution and intensity of precipitation can also indirectly affect ecosystem carbon balance because the induced water stress could alter the leaf area index (LAI) (Le Dantec et al., 2000;Barr et al., 2007), the carbohydrate reserve status (Sala et al., 2010), the plant allocation pattern (Sack and Grubb, 2002) or the soil microbial community (Sowerby et al., 2005). These indirect responses are often not instantaneous but lagged. Hu et al. (2010) observed that reduced snow cover in the winter led to water stress in the following summer and hence limited photosynthesis in a subalpine forest. Also, climate anomalies, e.g. high temperature in spring, can increase photosynthesis in the following autumn, possibly due to enhanced leaf nitrogen content and canopy photosynthetic capacity as a result of increased nitrogen mineralisation .
Several studies have suggested that to evaluate the climate change impacts on the ecosystem carbon balance, it is important to jointly consider the direct, indirect and lagged responses (Braswell et al., 1997;Delpierre et al., 2009;Dunn et al., 2007). However, to explicitly distinguish between these responses is difficult. Richardson et al. (2010) illustrated with four conceptual models that phenological transitions, which are defined here as an indirect response, can have direct, indirect and lagged impacts on ecosystem productivities. In a simpler approach, the ecosystem responses were classified into direct and biotic responses to environmental forcing . The biotic responses were described by the parameters of an ecosystem model (e.g. maximum canopy photosynthetic capacity and base respiration). Interannual variability in biotic responses was regarded as functional change.
Functional change can be changes in the structure, physiological properties or species composition of an ecosystem and can occur at either short or long time scales. It is a challenge to assess functional change over the whole succession of an ecosystem because most observations and experiments are conducted over relatively short periods (Symstad et al., 2003;Misson et al., 2010). Therefore, any functional change inferred can only represent the change within an ecosystem over the period of observation. Nevertheless, distinguishing the effects of climate variability and functional change on IAV of the carbon balance at a specific stage of an ecosystem is of interest. It allows the evaluation of the necessity of incorporating functional change modules into mechanistic models, which are used to project future ecosystem carbon balances. Hui et al. (2003) used a homogeneity-of-slope model and a stepwise multiple regression approach to assess the IAV of the biotic responses of the Duke Forest, concluding that functional changes accounted for about 10 % of the observed variation in the NEE. Richardson et al. (2007) concluded that it is important to consider the time scale when trying to partition the source of variance in NEE. With the parameters of a modified light response model fitted to NEE in each year, as much as 55 % of the variation could be attributed to the biotic responses at an annual time scale. In contrast, the effect of functional changes were found to be much lower in a peatland (Teklemariam et al., 2010).
Both the abiotic (i.e. direct) and biotic responses to climate variability have seasonal patterns. For instance, the limiting factor for photosynthesis may change at different periods of the year. Also, key parameters (e.g. canopy maximum photosynthetic capacity) can vary seasonally (Hollinger et al., 2004;Wang et al., 2007;Thum et al., 2008). Recently, more studies have concentrated on the analysis of the carbon dynamics in individual seasons (Piao et al., 2008;Wang et al., 2011). Piao et al. (2008) showed that the warming effect in autumn accounts for the annual carbon loss (with a sensitivity of 0.2 Pg C • C −1 ) in northern terrestrial ecosystems. This was supported by a site level study where the autumn temperature dominated the annual carbon balance (Vesala et al., 2010). From these findings we conclude that the IAV of ecosystem carbon balances should not only be analyzed at annual but also at sub-annual time scales (e.g. weekly, monthly or seasonally).
At a temperate beech forest near Sorø, Zealand Denmark, NEE was quasi-continuously measured over the past 13 yr (1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009). The annual NEE (a negative sign corresponding to a net sink) ranged from a source of 32 g C m −2 yr −1 in 1998 to a sink of −344 g C m −2 yr −1 in 2008. A decadal trend of NEE at this site has been reported by Pilegaard et al. (2011), with an average increase in net uptake of 23 g C m −2 yr −2 . The study indicated that an extended carbon uptake period explained 33 % of the trend and Biogeosciences, 9, 13-28, 2012 www.biogeosciences.net/9/13/2012/ hypothesised that an increased maximum photosynthetic capacity during the growing season could be another important factor explaining the trend. The aim of this study was to investigate the IAV of carbon fluxes at Sorø with respect to (1) the seasonal pattern of the ecosystem responses and (2) the source of variability in carbon fluxes. We investigated to what extent and at which temporal scales, climatic variability and changes in ecosystem functional properties determined the IAV of the carbon balance. Empirical models were used as tools to estimate the seasonal and interannual variation of the ecosystem functional properties and to distinguish their impact on the carbon fluxes from the direct climatic effects. An additional methodological focus was to improve published methods with a new approach. We tested the hypothesis that the observed longterm trend of increasing carbon uptake is mostly driven by changes in ecosystem functional properties, despite the significant role of weather variability at short time scales.

Site description
Field measurements were taken at the Euroflux network station Sorø on Zealand, Denmark (55 • 29 N, 11 • 38 E). Mean annual temperature during the measurement period was 8.5 • C and mean annual sum of precipitation was 564 mm. The dominant tree species is European beech (Fagus sylvatica) with approximately 20 % conifers, mainly Norway spruce (Picea abies) and European larch (Larix decidua). In 2010, the stand around the flux tower was 89 yr old, the average tree height was 28 m and the average diameter at breast height was 41 cm. Soils were classified as Alfisols or Mollisols (depending on the base saturation) with 10-40 cm deep organic layers. Leaf area index peaked at 4-5 m −2 m −2 and no significant trend was observed in . Further information on the instrumentation can be found in Pilegaard et al. (2003). The soil and vegetation are described in Ladekarl (2001) and Pilegaard et al. (2001). The fetch and footprint analysis are given in Jensen (2000, 2005), Göckede et al. (2008) and .

Flux measurements and partitioning
Thirteen years (1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)) of half-hourly measurements of NEE and climate data were used in the present analysis. The flux measurements were conducted at 43 m above ground and the data processing followed the standard procedure of Aubinet et al. (2000) with modifications described in . Spectral corrections were applied to the flux data according to Ibrom et al. (2007), using a spectral transfer function approach. The fluxes were filtered for low turbulent mixing at stable stratification when the friction velocity (u * ) was lower than 0.1 m s −1 . The net ecosystem CO 2 exchange was calculated as the sum of CO 2 storage change below the sensor and the measured turbulent CO 2 flux. The gap-filling method was described in Pilegaard et al. (2011). Total ecosystem respiration (a positive sign corresponding to a release of CO 2 to the atmosphere) was estimated based on nighttime data and extrapolated to daytime conditions according to the following equation: where T 0 is the reference soil temperature at 2 cm depth (0 • C), r b is the base respiration at T 0 , T s is the measured soil temperature at 2 cm and Q 10 is the temperature sensitivity parameter and set to a constant value of 2. Base respiration was estimated for every night and Eq. (1) was used to extrapolate the nighttime ecosystem respiration over daytime based on soil temperature measured at daytime. GPP (a negative sign corresponding to a net sink) was subsequently calculated as:

Correlation analysis
To evaluate the ecosystem response to climatic variability, bivariate correlation analysis was first performed between the annual carbon flux integrals (annual sums of NEE, GPP and TER) and the mean annual climate variables, i.e. air temperature T air , global radiation R g , volumetric soil water content SWC in the top soil, and precipitation PPT, using data from 13 yr. In the second step, the same correlation analysis was applied at sub-annual time scale with a 30-day moving window. Pearson's correlation coefficients (n = 13) were calculated between periodical flux integrals (30-day sum of the GPP and TER) and mean periodical climate variables (e.g. mean T air for DOY 1-30) from each of the 13 yr. With the 30-day moving window, the first correlation coefficient calculated represented the period DOY 1-30 and was centred at DOY 15. The second correlation coefficient was calculated by moving the window one day forward, representing the period DOY 2-31. In this way, correlation coefficients were calculated throughout the whole year. The resulting time series of correlation coefficients related the interannual variability of the carbon flux to a potential climatic driver in a certain period of a year. This enabled the analysis of interannual variability at sub-annual time scales.

Model description and parameter estimation
Two statistical modelling approaches were applied in this study to distinguish the impact of climatic variability and ecosystem functional change on carbon fluxes. Firstly, we applied the model (referred to as Model I) and parameter estimation schemes of Richardson et al. (2007). The model consists of a simple Michaelis-Menten light response model (e.g. Hollinger et al., 2004), a Lloyd and Taylor respiration model (1994) and a set of sigmoid scalar functions which are used to calculate potential carbon flux (under optimum environmental conditions) from the actual flux. Parameter estimation was conducted by minimizing the weighted least squares cost function (Richardson and Hollinger, 2005) in two steps. In the first step, 15 parameters were fitted to the entire 13 yr datasets. Then in the second step, 11 of the 15 parameters (10 for five environmental scalar functions and the temperature sensitivity parameter) were fixed and the remaining four parameters were allowed to vary in each year (for details see Richardson et al. 2007). The interannual variation of these four parameters was interpreted as change in the functional properties of the ecosystem. Because the parameters of Model I were estimated at annual or longer time steps, the seasonal variation of the ecosystem functional properties cannot be reproduced. This questions its application to ecosystems with strong seasonality, such as the temperate deciduous forest in this study. Therefore, we implemented a second approach based on a similar model (referred to as Model II) but with a different parameter estimation scheme using short moving time windows (Lasslop et al., 2010b). The model consisted of a rectangular hyperbolic light response model (Falge et al., 2001) and a Lloyd and Taylor (1994) respiration model. The effects of air humidity on photosynthesis was modelled after Körner (1995). The model is described in Eqs. (3) and (4).
where T air is the air temperature, T 0 ( • C) is a constant (−46.02 • C), R g (W m −2 ) is the global radiation; α (µmol CO 2 J −1 ) is the light use efficiency, and β (µmol CO 2 m −2 s −1 ) is the instantaneous maximum canopy photosynthetic capacity (Eq. 4), which represents both the structural and physiological properties (e.g. leaf area index, leaf photosynthesis capacity, C3 or C4 photosynthesis pathway species composition); r b (µmol CO 2 m −2 s −1 ) is the base respiration at reference temperature ( where k (hPa −1 ) is a scaling parameter estimating the effects of VPD on β and VPD 0 is a threshold value set as 10 hPa, above which the stomatal conductances reduce β. The parameters of the model (E 0 , r b , α, β, k) were estimated in two steps: All parameters were allowed to vary throughout the year: E 0 was derived from nighttime NEE every 2 days using a 12 day window. α, β, r b and k were estimated based on daytime NEE every two days, with a 4 day moving window.
In the first step, the parameters k and E 0 were not well constrained. Therefore, we fixed E 0 and k as the median of the previous estimates (E 0 =141.4 • C, k = 0.09 hPa −1 ), after which α, β, and r b were re-estimated based on daytime data every two days, with a 4 day moving window. The temperature sensitivity can be fixed in different ways. Richardson et al. (2007) fixed it by fitting a single set of parameters globally to the data of all years. However, as demonstrated in Reichstein et al. (2005), the estimated temperature sensitivity clearly differs when it was fitted globally or with small time windows. Therefore we decided to use the median of the estimated values from short time windows. Parameter estimation was conducted by minimizing the same weighted least squares cost function as in the first modelling approach (for details see Lasslop et al., 2010b) and the half-hourly fluxes were computed using gap-filled parameter sets (distance weighted average between the two adjacent parameters) and climate data.
Changing model complexity and structure may affect the parameter estimation and thus the subsequent partitioning of variance between the climate and parameter effects (Leuning, 2011). Therefore we performed two additional runs by removing the VPD effect function (Eq. 4) and also imposing a stronger effect (with a lower threshold value VPD 0 ) to test this hypothesis.
In order to compare the seasonal course of the parameters with the structural development of the canopy, the photosynthetically active radiation (PAR) transmittance records (ratio of the PAR measured below and above canopy, details see Pilegaard et al., 2011) were used as a reference, representing the temporal course in the leaf area index.

Distinguishing the direct response from biotic changes
By application of the parameter time series of one year to all the other years and comparing the differences in the modelled fluxes, the effects of direct climate forcing on the carbon fluxes can be investigated, vice versa for the effects of the parameters. After Richardson et al. (2007), model simulations with fixed climate (using climate data of one year for all other years) and fixed parameters (using parameter dataset of one year for all other years) were performed, resulting in a 13 × 13 matrix of model predictions. In each cell of this matrix, the simulated results contained 17 520 half-hourly data points, which were further aggregated by day, week, month, season and year for statistical analysis. Therefore, the datasets in each column represented the modelled fluxes with fixed climate using climate data of a particular year i (F clifix,i ) and the datasets in each row represented the modelled fluxes with fixed parameters using parameter time series of year i (F parfix,i ). The diagonal of the matrix represents the original modelled fluxes (F original ). The differences between F original and F clifix represent the variability that is driven directly by the climate while the differences between F original and F parfix represent the variability that is driven by the changing model parameters, which we interpret as functional change. We applied a sums of squares Biogeosciences, 9, 13-28, 2012 www.biogeosciences.net/9/13/2012/ approach to distinguish between the relative effect of climate and model parameters on the IAV in NEE, GPP and TER at multiple time scales ranging from daily to yearly. The proportion of IAV in the carbon fluxes explained by climate variability (E cli ) and parameter change (E par ) was determined by Eq. (5) and (6), respectively.
where σ is the variance of the interannual variation of the carbon fluxes and 3 Results

Interannual variability in carbon fluxes and climate variables
Variation in NEE and the climate variables were analyzed at both sub-annual (daily values smoothed with a 30-day moving average) and annual time scales (Fig. 1). The variability is displayed as fluctuations (around the mean) relative to the standard deviations in a certain time window over the 13 yr. This presentation enables a standardized comparison of the anomalies in both carbon fluxes and climate variables with their interannual variability, during a particular period of a given year. In 1998, the forest was a source of carbon dioxide (Table 1) where the deviation of NEE from the 13 yr annual mean was about 2 standard deviations (SD) above average (more positive corresponding to less uptake). Following this, a trend of increasing carbon uptake in the forest was observed, manifested in a higher net uptake during the latest years 2008-2009 when the NEE was almost 2 SD below average (Fig. 1). Except for these years, the annual NEE deviations were less than 1 SD. Generally, there was a tendency from above average towards below average except in 2003 and 2006, when the forest appeared to have a reduced rate of carbon uptake. Carbon uptake was less than average in almost all periods of the years from 1997 to 1999, except in autumn 1999 when the uptake was about 1 SD above average. In 2000-2006 the NEE anomaly fluctuated remarkably within each year, ranging from 2 SD above average in winter 2000 to 2 SD below average in autumn 2005. Regardless of the high variability at short time scales, these NEE anomalies tended to compensate each other, resulting in annual fluxes close to the 13 yr average. From 2007 to 2009, NEE was continuously below average with higher carbon uptake except in two short periods of summer 2007 and the beginning of 2008.
The mean annual air temperature of the only source year 1998 was the lowest over the 13 yr period, with more than 1 SD below average. Meanwhile, overcast weather (lowest incoming radiation) and the slightly higher than average precipitation kept the soil water content at 1 SD above average. The opposite conditions prevailed in 2008 compared to 1998, when a higher radiation (almost 2 SD above average) was accompanied with higher air temperature and evaporative demand which significantly reduced the soil water content. Note that the measurement of soil water content in 2008-2009 was slightly different from the previous years because the original TDR sensor (TRIME-EZ, Imko, Ettlingen, Germany) was replaced by a ThetaProbe ML-2x (Delta-T Devices, Cambridge, UK).The new installation was at the same location but in a shallower soil horizon (5 cm instead of 10 cm).

Relationship between the interannual variability in the carbon fluxes and climate variables
The correlation analysis related the interannual variability of the component carbon fluxes to certain climate variables. At the annual time scale, NEE was highly correlated with GPP (r = 0.7, p < 0.01) but not with TER (Table 2). TER was correlated with GPP (r = −0.65, p < 0.05). Apart from soil water content, which was positively correlated with precipitation (0.76, p < 0.01) and negatively correlated with radiation (−0.77, p < 0.01), there were no significant correlations among the other climate variables. Global radiation and soil water content were negatively (−0.73, p < 0.01) and positively (0.62, p < 0.05) correlated with NEE, respectively. Soil water content was low at high R g and correlated with NEE. Surprisingly, none of the climate variables showed a significant correlation with the component fluxes, GPP and TER. Contrary to the analysis at the annual time scale, carbon fluxes were clearly correlated with climatic variables at shorter time scales (Fig. 2). Gross primary production was highly correlated with T air (r < −0.69, p < 0.01) throughout the entire year except during the summer (Fig. 2a). When the soil was usually dry in the summer, it apparently confounded the stimulating effects of temperature on photosynthesis. The correlation coefficient between GPP and T air was lowest (r < −0.8, p < 0.01) in April (leaf flush) and October (leaf fall), indicating that the phenology and leaf development were temperature sensitive. Radiation was also negatively correlated with GPP during spring and autumn. The correlation coefficients increased directly after leaf senescence but decreased during winter, representing the photosynthesis of the 20 % coniferous trees within the footprint. Over a large part of the year, the correlation between GPP and SWC was not significant, possibly because the deep rooting system allowed sufficient water supply even when the water in the top soil was depleted. Nevertheless, a trend of decreasing correlation coefficients during the summer was obvious and the value was lowest in August-September. The Fig. 1. Deviation of NEE, air temperature (T air ), global radiation (R g ), precipitation (PPT) and soil water content (SWC) from the 13 yr annual average (black solid line) and daily average values (grey areas). The daily data were smoothed with a 30-day moving average. The y-axis gives the standard deviation from the mean (e.g. the annual NEE anomaly in 1998 is about 2 SD above the 13 yr average with less uptake of carbon). correlation between T air and TER (Fig. 2b) was positive during the spring, autumn and winter but also turned negative during summer due to the confounding effect of SWC. Similar to the analysis at the annual time scale, TER and GPP were significantly correlated during large parts of spring and summer. Because the climate variables were also inter correlated (Fig. 2c), the attribution of the interannual variability in carbon flux to specific climatic variables was sometimes difficult, especially for T air and R g , which were significantly correlated during spring and summer. Global radiation and SWC were in general negatively correlated and most significantly in June and October. The comparison of the temporal variability of climate variables and carbon fluxes (Fig. 3) during the average and particular years (1998 and 2008, minimum and maximum in annual NEE) illustrates the seasonal pattern of the climate relationship together with the IAV of the carbon fluxes (Fig. 2). Relative to the average, low T air and R g were related with low GPP (Fig. 3a, b and e) in 1998, and high SWC co-occurred with high TER (Fig. 3c and f). In 2008, T air and R g were related with high GPP (Fig. 3a, b and e), and low SWC was only weakly connected to TER ( Fig. 3c and f).

Model predictions and estimated parameter time series
Based on the approach of Richardson et al. (2007), Model I explained on average 77 % of the variance in the annual half-hourly measured (non gap-filled) NEE (Table 1). The root mean squared error (RMSE) was on average 3.92 µmol CO 2 m −2 s −1 and mean absolute error (MAE) was 0.24 µmol CO 2 m −2 s −1 . With the new approach (Model II), the r 2 increased to 85 % and RMSE and MAE was both lower, at 3.13 and 0.08 µmol CO 2 m −2 s −1 , respectively. The model residuals (Fig. 4) showed that Model I tended to be strongly biased across different seasons. The carbon uptake was overestimated during spring and autumn, and underestimated in summer. The predictions based on Model II matched the seasonality of carbon fluxes well; this improvement is particularly important for using this semi-empirical Biogeosciences, 9, 13-28, 2012 www.biogeosciences.net/9/13/2012/ Air temperature, T air ( • C) and global radiation, R g (Wm −2 ) were measured above the canopy. Soil water content, SWC ( %) was measured at 0-10 cm, the system broke down temporarily in 2007-2008. Methods for the gap-filling of the observed NEE, NEE obs (g C m −2 yr −1 ) and modelled NEE (g C m −2 yr −1 ) are described in the text. Model error statistics include, the number of valid observations (N) used for parameter estimation in each year, coefficient of determination (r 2 ), root mean squared error, RMSE ( µmol CO 2 m −2 s −1 ) and mean absolute error, MAE ( µmol CO 2 m −2 s −1 ). model as a tool to distinguish the effect of direct climate variability from the effect of functional change on NEE at subannual time scales. At an annual time scale, both models reproduced the IAV in the measured NEE. The correlation coefficients between gap-filled NEE and modelled NEE of Model I and Model II were 0.85 and 0.87, respectively. Except for 2004, the modelled NEE (based on Model II) was lowest in 2008-2009 (−294 and −348 g C m −2 yr −1 ) and highest in 1998 (82 g C m −2 yr −1 ) which agreed with the gap-filled fluxes ( Table 1). The discrepancy between gapfilled and modelled NEE was large in 2004, as up to 25 % of the data were missing in the growing season. While the measured NEE time series was gap-filled , the corresponding gaps in the simulated NEE were filled using model predictions from distance weighted averages of the two gap-adjacent parameter values and climatic data. These uncertainties contributed to this discrepancy. Interannual variation in the modelled NEE (1SD = 136 g C m −2 yr −1 ) was about 30 % higher than the measured NEE (1SD = 104 g C m −2 yr −1 ).
The estimated parameter time series based on Model II varied between years. During the leafed period, mean light use efficiency (α) was highest in 2000 (0.12 µmol CO 2 J −1 ) and lowest in 2003 and 1997 (0.1 µmol CO 2 J −1 ). The averaged canopy maximum photosynthetic capacity (β) ranged between 35.6 (1997) to 45.5 (2006) µmol CO 2 m −2 s −1 . The annual mean base respiration (r b ) ranged from 5.1 (2005) to 6.2 (1998) µmol CO 2 m −2 s −1 . The seasonal variation in the parameters was stronger than the variation in the annual means. In general, α, β and r b were all significantly higher during the growing season than in winter (Fig. 5) which is consistent with expected effects of phenology, substrate availability and plant respiratory dynamics. Both the parameter α and β were highest in early summer and decreased over the vegetation period. The base respiration ranged from 2.9 to 8.9 µmol CO 2 m −2 s −1 on average, it reached its peak value in July and then reduced until winter, while the SWC reached its lowest value in September and then increased again.
When the VPD effect (Eq. 4) on β was excluded from the model structure, the estimates of β were reduced considerably during the growing season (Fig. 5). This decrease was stronger in July when the relative humidity was generally lower than annual average. Because of parameter correlation, α and r b changed during this period and the canopy photosynthesis at light saturation (P max ) also decreased. When a stronger VPD effect was imposed by reducing the VPD threshold value, the estimates of α, β and r b changed in the opposite direction. The mean absolute error of the model predictions increased from 0.08 to 0.21 and 0.26 µmol CO 2 m −2 s −1 respectively, when the VPD effect  was excluded (underestimation of carbon uptake during summer) or enhanced (overestimated carbon uptake during summer). The comparison between the measured PAR transmission (as a proxy for the seasonal LAI development) and photosynthesis at light saturation (P max ) indicates that the structural change of the forest took place mainly in April and October (Fig. 5d). The sharp decrease in P max marked the start of leaf senescence (mid-September) already one month before the average start time of leaf fall (November). This demonstrated that the aggregated effects of structural and physiological changes on the seasonality of NEE were well represented by our modelling approach.

Sources of interannual variability in the carbon fluxes
One main objective of this work was to separate the influence of functional change from direct effects of climatic variability. This was achieved by analyzing model simulations with Table 2. Bivariate correlation coefficients between annual anomalies in NEE, GPP and TER and mean annual climate variables including air temperature (T air ), global radiation (R g ), soil water content (SWC) and precipitation (PPT). Statistically significant correlations are marked with **(p < 0.01) and *(p < 0.05). Interannual variation in the modelled NEE (Model II) differed when either the parameters or the climate variables were fixed (Fig. 6). When the parameters were fixed, i.e. assuming the ecosystem status during the 13 yr was constant (Fig. 6b), the modelled fluxes NEE parfix ranged between −231 g C m −2 yr −1 in 2009 and −27 g C m −2 yr −1 in 2006. When the climate time series were fixed, NEE clifix ranged between 81 g C m −2 yr −1 in 1998 and −293 g C m −2 yr −1 in 2005 (Fig. 6c). The IAV in the originally modelled fluxes NEE original (Fig. 6a) were better reproduced when the changes in the parameter time series were included. The correlation between the NEE original and NEE clifix was 0.9 (p < 0.01). In contrast, the correlation between the NEE parfix and NEE original was not significant. The differences in these modelled fluxes represent the climate and parameter effects; these effects can be interpreted as the influence of changes in climatic and ecosystem functional properties, respectively.
The analysis of the 13 × 13 matrix of model predictions (Model II), where the parameter and climate time series of each year were applied factorially, attributed 78-82 % of the variance in NEE to the interannual variation in the parameter time series (Fig. 7). Using Model I, we obtained similar results as parameter variability accounted for 80-83 % the IAV in annual NEE. In general, the percentage of the total variance in the carbon flux caused by climatic variability decreased with the level of temporal integration. The estimated impact of climatic variability at shorter time scales using different models did not converge as much as in the annual time scale. The estimated impact of climatic variability, for the IAV in NEE based on Model II at sub-annual scales was clearly lower than those estimated based on Model I (Fig. 7). Because Model I cannot represent the seasonal pattern of parameter change with the same accuracy as Model II, it was expected that the impact of functional change will Biogeosciences, 9, 13-28, 2012 www.biogeosciences.net/9/13/2012/ Fig. 3. The average (black line) and two example years (1998,2008) in terms of climate (a, b, c), carbon flux (d, e, f) and cumulative carbon flux (g, h, i). All data were smoothed with a 14-day moving average. Grey lines indicate the maximum and minimum of the value on the same DOY across the 13 yr. be underestimated. We showed in the previous section that the prediction of Model I was strongly biased throughout the year (Fig. 4) and that the resulting errors were propagated in to the analysis and biased the results. The estimated pa-rameter effect on NEE at annual time scale did not change significantly when the VPD effect in the model was removed (81-85 %) or changed (with a lower VPD 0 , 79-83 %). At other time scales, the difference in the estimated impact of ecosystem functional change caused by changing the VPD function was less than 2 %. GPP was much more sensitive to variation in climate than TER was. At the daily time scale, climate variability accounted for more than 67-71 % of the variation in GPP. The effect of climate variability on TER was very low (19-22 %), even at the daily time scale. This is most likely because the dominant component of TER is soil respiration, a process which is highly dependent on substrate availability and microbial activities (Davidson et al., 2006). This then results in reduced climatic effects on TER, compared to the aboveground processes, such as photosynthesis and phenology. On the other hand, because GPP and TER are only influenced by certain climate variables in the model structure (Eq. 3 and 4), the unaccounted effect of other climate variables on GPP and TER could be propagated into the parameters and thus result in the underestimation of the effect of climate variability. We assessed the influence of temperature on GPP which is not explicitly included in the model by plotting the simulated daily GPP at light saturation (R g = 1000 Wm −2 ) with T air (Fig. 8a) the GPP-related parameters was found when the whole year data was used, however, during the period when the canopy was fully developed and LAI was relatively constant (May-September), no significant relationship between T air and the simulated potential GPP was found. The TER-related parameters were not completely independent from SWC. During the growing season, the predicted TER at 25 • C was still weakly correlated with SWC (Fig. 8b).

Strengths and limitations of the statistical modelling approach
The combination of modelling and data analysis provides insight into the complex interactions between the direct climatic effects and the biotic ecosystem responses, by partitioning their roles in determining the interannual variation of the ecosystem carbon fluxes. In this study, two empirical models were used to this end. The comparison of an existing approach (Model I), and our approach (Model II) showed similar estimates of the impact of functional change at the annual time scale (ca. 80 %). However, our approach, which allowed the parameters to vary seasonally within year, resulted in less seasonal bias (Fig. 4), increased the degree of determination, and thus, increased the confidence in the estimate of the impact of functional change at both the annual and sub-annual time scales.
A major strength of the statistical modelling approach is that it does not require an explicit parameterization of the complex ecosystem process in detail. It has been proposed that in order to realistically distinguish between the effect of climate variability and functional change using a parameter extraction method, a mechanistic model must be used . However, many processes, such as phenology, are considered indirectly by the parameter estimation method. These responses to current and previous climate can be considered to be ecosystem structural change and are captured by this method (Fig. 5). Also, the information contained in the flux data alone is likely insufficient to distinguish the between detailed ecosystem processes (Wang et al., 2009). The risk of over-parameterization is a strong justification to keep the model structure simple. For example, detailed processes such as direct and diffuse radiation  and the separation of autotrophic and heterotrophic respiration were not implemented. The effects of many of those processes are aggregated into the parameters of the empirical model used. A mechanistic model, although desirable for other reasons, might not be an efficient way to achieve the objectives of this study, as it would need much more information to reach the same goal: the realistic representation of the seasonal ecosystem functional state. To address the concerns of the simplicity of our model structure, instead of adding a mechanistic sub-module, we tested the effect of removing the VPD sub-module. The results showed that the partitioning between the climate and parameter effects was not changed. We therefore conclude that the available information on ecosystem functional change within the CO 2 flux data was captured by the flexible modelling approach.
Statistical models have limited potential to elucidate the causes behind the identified functional change. However, this was not the primary aim of this study, as the focus was to distinguish between the direct climatic effects and the changes of ecosystem functioning. The model and fitted parameters cannot be used for predictive purposes or extrapolation. As indicated by Groenendijk et al. (2011), although specific siteyear derived parameters give the best prediction of observed fluxes, they are generally too specific to be used in global studies. However, this modelling approach can be applied to flux data from other sites to obtain site specific parameters and quantify the importance of functional change for the IAV of carbon fluxes in different ecosystems. We advocate using our approach as it offers a high flexibility without requirements for additional information.

Magnitude and uncertainty of the estimated impact of functional change
The most important finding of this study is that the 13 yr trend of increasing carbon uptake was found to be more strongly driven by the aggregated effect of the parameters than the climatic factors. In addition to the changes in the maximum photosynthetic capacity and the carbon uptake period , we found that other photosynthesis and respiration related parameters also changed in different years and indirectly affected the ecosystem carbon balance. The estimated impact of functional change on the IAV of the carbon balance at Sorø was higher than at other sites (Teklemariam et al., 2010;Richardson et al., 2007;Hui et al., 2003), although the methods applied were different. However, using the same model of Richardson et al. (2007), we found also a much higher impact of functional change at Sorø (deciduous) than at Howland Forest (mixed). The results at annual time scale at Sorø did not differ much between Model I and Model II. This is probably due to the fact that the dominance of the functional change over direct climatic effect is so strong at this site that the choice of model is not particularly important. But as our model has not been applied at other sites, this cannot be stated with sufficient certainty. A possible reason for such tendencies could be the type and structure of the ecosystem. Based on a number of published studies, the impact of functional change decreases in the order: deciduous forest (this study), mixed forest , conifer forest (Hui et al., 2003) and peatland (Teklemariam et al., 2010), implying a possible difference in the sensitivity of these ecosystems to environmental change and disturbance. Cross-site studies of IAV in the ecosystem carbon balance also found that deciduous forests tend to be more sensitive than evergreen conifer forests to climate variability (Yuan et al., 2009) and phenological transitions . These differences in the adaptive capacity to changing climate between deciduous and evergreen forests may drive shifts in the composition within mixed forests .
The estimated impact of functional change is not without uncertainties, mainly due to the core assumption that functional change is solely represented by the changes in model parameters. This assumption is challenged by the simplicity of the model structure, e.g. the temperature regulation of GPP and the SWC regulation of TER were not explicitly specified in the model structure. Testing the residual climatic effects on the model predictions showed that, although we did not explicitly account for the temperature effect on GPP in the model, it was at least partly represented by the model, possibly due to the cross correlation between temperature and VPD or temperature and radiation. Therefore the largest part of the correlation between temperature and parameters, when using the whole data set, are derived from functional changes, e.g. canopy development that were either influenced by or coincided with temperature. In contrast, the TER-related parameters were not completely independent from SWC. Approximately 19 % of the parameter effect on TER could be attributed to variation in SWC, which has led to a small overestimation of functional change.

Climate variability and average climate change
The correlation analysis between carbon fluxes and climatic variables revealed that the variation in the mean annual climate could not directly explain the IAV in the annual ecosystem carbon balance. Global radiation and soil water content were significantly correlated with NEE (Table 2). However, contrary to our expectations, neither R g nor SWC were significantly correlated with gross photosynthesis or total ecosystem respiration, i.e. the processes that drive NEE.
Ecosystem respiration was low in years with higher radiation (r = −0.44, p > 0.05), which might have strengthened the correlation between R g and NEE. SWC was low in years with high radiation and was also correlated with NEE. Therefore, despite comparably high correlation coefficients, these correlations do not represent direct cause-effect relationships between mean annual climate variables and NEE. Rather, they are the result of the combination of the various controlling mechanisms in different seasons, as illustrated by the correlation analysis at short time scales.
At shorter time scales, the sensitivity of carbon fluxes to climate variability significantly increased. For more than 80 % of the year, the carbon flux was strongly correlated with at least one of the investigated climate variable. The climate variables, to which the carbon fluxes were most sensitive, differed across seasons (Fig. 2). These variations affected the net ecosystem carbon balance via the different climate sensitivities of GPP and TER. Depending on the sign of the response and the different sensitivity, the climatic impact on NEE is either amplified or attenuated. Analyzing the climate sensitivity of the component fluxes (GPP, TER) contributes more to the understanding of the climate sensitivity of NEE than directly analyzing the climate sensitivity of NEE itself. This significant seasonality of the ecosystem responses to climatic variability demonstrated the importance of the seasonal distribution of the climate anomalies on future ecosystem carbon balances. Our results showed that this beech forest will be sensitive to increases in summer drought. Altered precipitation patterns, i.e. increased rainfall variability rather than changes in annual precipitation are likely to affect ecosystem carbon balances in the future (Knapp et al., 2002). The average climate change is expected to be accompanied by increased variability and weather extremes (Easterling et al., 2000a, b). Climate change projections for Denmark Biogeosciences, 9, 13-28, 2012 www.biogeosciences.net/9/13/2012/ that suggest that a small increase in precipitation will occur, but mainly during winter, while the likelihood of summer droughts will increase (Christensen and Christensen, 2007). Our results suggest that in addition to the changes in average climate, increased climatic variability could alter the ecosystem carbon balance more strongly, as the climate anomalies are projected to take place predominantly during biologically active periods.

Gross primary productivity as a driver of the interannual variability in net ecosystem exchange
Interannual variation in NEE in this beech forest was mostly driven by GPP whereas TER had relatively less influence. This was similar to two deciduous forests in both boreal and temperate zones (Barr et al., 2002). However, the results of these two studies differed from the conclusion of a crosssite synthesis of 15 European forests (Valentini et al., 2000) where respiration was found to determine the spatial variability in ecosystem carbon balance. According to this synthesis, net carbon uptake decreased significantly with increasing latitude, whereas total ecosystem respiration increased and gross photosynthesis tended to be rather independent of latitude. These different findings based on site specific studies and cross-site synthesis indicate that although TER tends to vary significantly and dominate the ecosystem carbon balance over large spatial scales, its influence on the temporal variability in NEE may be much weaker at site level. For deciduous forests at middle and high latitudes, the variability in GPP was much stronger and largely controlled the interannual variation of the ecosystem carbon balance. Total ecosystem respiration was highly correlated with GPP, both at annual and sub-annual time scales. It is still under debate whether this correlation could be artificial because GPP was calculated from TER and NEE (Vickers et al., 2009). Lasslop et al. (2010a), however, showed that only the error in TER that directly propagates into GPP can cause spurious correlation. A large part of the variability of TER was however shown to be independent from GPP, and thus the correlation between GPP and TER was real, rather than spurious. This suggestion is supported by using a TER estimate derived from a fit of a respiration model to nighttime data in combination with a GPP estimate derived using a light response curve fit to daytime data. In this case the errors in one flux component cannot directly propagate to the other and the correlation remains high (Lasslop et al., 2010b). A more likely explanation for the high correlation is the covariance of the main drivers of photosynthesis and respiration, e.g. temperature and radiation. From the physiological perspective the correlation could be an effect of substrate availability on autotrophic (plant activity, i.e. growth respiration and reserve metabolism) and heterotrophic respiration (root exudates and litter production). This interpretation is supported by an increasing number of experimental studies, such as tree girdling studies (Högberg et al., 2001) or direct and quasi parallel measurement of GPP and TER during daytime with ecosystem chambers in shrublands (Larsen et al., 2007).

Implication for mechanistic models
The results of this study suggest that projection of future carbon balance of terrestrial ecosystems could be improved if the biotic responses to climate variability and thus functional change are properly incorporated into mechanistic models. Migliavacca et al. (2011) demonstrated that the spatiotemporal variability in ecosystem respiration can be better modelled when the dynamics in the biotic factors are taken into account. Modelling the IAV of the carbon balance has proven difficult (Siqueira et al., 2006). One important reason is that most parameters are usually kept constant. For example, in global models, parameters are usually static for plant functional types (Krinner, 2005;Sitch et al., 2003). Increasing model complexity by adding processes that enable change in ecosystem state (e.g. nitrogen cycling or dynamics in the microbial community) could possibly improve the situation. However, this could lead to model over-parameterization and increase the demand for validation data, thus limiting model application at large spatial scales. Flux studies should therefore be accompanied by biological process studies and ecosystem structural assessments to overcome the data limitations that result in model over-parameterization. Site level studies clearly show the necessity for further development of these functional change modules. Thus, using an alternative approach by establishing empirical functional relationships between parameters and independent variables and subsequently embedding them in the model could be an intermediate solution. Therefore, combining experimental studies, empirical and mechanistic modelling could elucidate the importance of functional changes when simulating future terrestrial carbon budgets and potentially improve the prognostic capacity of ecosystem models.

Conclusions
An approach to separate the direct climatic effects on the interannual variability of carbon fluxes from ecosystem functional changes  was improved by employing the semi-empirical model of Lasslop et al. (2010b). By fitting the model parameters in small moving windows, the seasonality of the model parameters was accounted for more realistically. This made the representation of the seasonality of ecosystem functioning more flexible, rather than constraining them to prescribed relationships that are constant throughout the year.
Climate variability exerted a strong control on the carbon fluxes at short time scales but this impact weakened as the time integral increased. At longer temporal scales, the effect of ecosystem functional change became progressively larger and appeared to dominate the interannual variability in the www.biogeosciences.net/9/13/2012/ Biogeosciences, 9, 13-28, 2012 ecosystem carbon balance. From the high climate sensitivity at short time scales we conclude that the ecosystem will be sensitive to the seasonal distribution of climate anomalies that are expected to increase in the future. The observed trend of increasing carbon uptake  was found to be driven by ecosystem functional changes rather than direct effects of decadal climatic variability. Such strong effects demonstrate the need to integrate ecosystem functional change into mechanistic models, if they are to realistically predict future ecosystem carbon balances and thus feedbacks to climate change.