Interactive comment on “ On the uncertainty of phenological responses to climate change and its implication for terrestrial biosphere models ”

The authors investigated the uncertainties in model predictions related to phenology. Bud burst observations of the Harvard site were used to determine model parameters and their uncertainties of 12 different phenology models. These models and parameters were then used to predict the bud burst at the Harvard site in the next century for two different scenario’s. In addition the BEPS model is used to investigate the influence of phenology on the uncertainty of photosynthesis and evapotranspiration fluxes.


Introduction
Phenology is the study of the timing of recurrent biological events and the causes of their temporal change in response to biotic and abiotic forces (Lieth and Radford, 1971).As variability in the timing of phenology is tightly coupled to variability in climate, phenology can be considered as an important indicator of climate change (e.g.IPCC, 2007;Menzel et al., 2006).Numerous studies have documented the impacts of climate change on plant and tree phenology.For instance, the Fourth Assessment Report (AR4) of the IPCC reported an overall trend towards earlier spring phenological events (e.g.bud-burst, leaf unfolding, flowering and pollen release) between 2 and 5 day decade −1 .Menzel et al. (2006) estimated an average advance of spring phenology in Europe of 2.5 day decade −1 while Schwartz et al. (2006) showed a similar earlier bud-burst of 1.1 day decade −1 across most temperate Northern Hemisphere land regions over the 1955-2002 period with different dynamics among major continental areas.Jeong et al. (2011)  regional scale and suggested a reduction of the rate of advancement of the start of the season in the period 2000-2008 (0.2 days) compared to the period 1982-1999 (5.2 days).
Although these studies highlighted that spring phenology has responded to recent climate change, large uncertainties remain as to how phenology will respond to projected future climate change.There are several conflicting reports in the literature about the relative roles of different environmental factors such as chilling (i.e.amount of accumulated cool temperature exposure required for breaking dormancy) and photoperiod (daylight duration control of the induction or release from dormancy) in controlling tree phenology.K örner and Basler (2010) have argued that because of photoperiodic constraints, observed effects of temperature on spring life-cycle events cannot be extrapolated to future temperature conditions without simultaneously accounting for photoperiod.Other studies (e.g.Chuine et al., 2010;Morin et al., 2009;Vitasse et al., 2011) suggest that photoperiod has not been shown to be more dominant than temperature when predicting bud-burst or flowering.These latter studies argue that temperature has a more dominant role than photoperiod during both the endodormancy phase (i.e. a non-growing -resting, quiescent or inactive-phase caused by conditions or factor within a plant or seed itself) and the ecodormancy phase (i.e. the cessation of growth induced by environmental factors) that controls spring phenology.Another approach (Schleip et al., 2008) suggests the importance of weighting temperature forcing, in a sensitive time span, to determine the effective temperature which controls each phenophase.These different hypotheses lead to different possible structures for phenology models, but the associated uncertainties in forecast responses of phenology to climate change have yet to be quantified.
Uncertainty in model projections can be classified in three categories: uncertainty model assumptions and formulations, with different processes described by different models (i.e.accumulation of degree-days, chilling requirement, photoperiod limitation); model driver uncertainty is due to uncertainty in future climate scenarios used for ecological forecasts (Cook et al., 2010).Model-data fusion (e.g.Keenan et al., 2011;Wang et al., 2009) provides a promising approach for assessing uncertainties in ecological forecasting.Also referred to as model-data integration or data assimilation, model-data fusion relies on the combination of models with observational constrains through an optimization approach (e.g.simulating annealing, quasi-newton methods etc.).In this way, model parameters, model states and their respective uncertainties can be estimated, conditional on the data (here: consistent phenological observations by a human observer).With a model-data fusion approach it is possible to objectively incorporate data, along with associated uncertainties, allowing for a full characterization of posterior distributions of model parameters.In this way, confidence estimates of model projections can be obtained, both for current climate conditions and for future climate change scenarios.
Phenology, and in particular bud-burst, controls numerous land surface feedbacks to the climate systems and ecological interactions through the regulation of exchanges of carbon, water and energy between the biosphere and atmosphere (e.g.Richardson et al., 2009Richardson et al., , 2010;;Baldocchi, 2008;Morisette et al., 2009;Fitzjarrald et al., 2001).Uncertainty in the prediction of spring phenology can therefore feed-forward to generate uncertainty in estimates of carbon and water cycling from terrestrial biosphere models.Several studies have shown the sensitivity of different biogeochemical and terrestrial biosphere models to bud-burst and other phenological transitions (e.g.Randerson et al., 2009;Levis and Bonan, 2004;White et al., 2000;Migliavacca et al., 2009).Furthermore, a recent multi-model synthesis study has shown that spring phenology is poorly simulated by different terrestrial biosphere models, resulting in large biases in model estimates of carbon cycling, in particular for deciduous broadleaf forests (Richardson et al., 2012).Hence, modeled future carbon, water and energy fluxes, as well as many biosphere-climate interactions projected by terrestrial biosphere models, Introduction

Conclusions References
Tables Figures

Back Close
Full Here we present a phenological forecasting study using phenological data for 11 different North American tree species observed at Harvard Forest over the last 20 yr (Richardson and O'Keefe, 2009).We advance beyond the recent work by Morin et al. (2009), by characterizing the sources of uncertainty in bud-burst forecasts for the 21st century, and what these uncertainties mean for modeling forest-atmosphere fluxes of carbon and water.
We combine phenological observations collected at the Harvard Forest with 12 different phenological models using a model-data fusion approach.With this analysis we characterize the uncertainty of model parameters and model structures.We then project model estimates of phenology forward, along with the associated parameter uncertainties, using statistically downscaled climate projections (Delworth et al., 2006;Hayhoe et al., 2007) for two different IPCC climate change scenarios (A1fi, or high CO 2 emissions scenario, and B1 or low CO 2 emissions scenario).This allows us to explore how the uncertainty characterized using current phenological observations is propagated in the future and to quantify how uncertainty in model parameters and model structure interacts with uncertainty in climate scenarios.
Then, we analyze the impact of the uncertainty of future bud-burst in a widely used terrestrial biosphere model (Boreal Ecosystems Productivity Simulator BEPS, Ju et al. (2006)).We evaluate the differences between gross primary productivity (GPP) and evapotranspiration (ET) as simulated by BEPS with the native phenological model and forced by the bud-burst forecasts obtained with the best model formulation selected according to data.Finally, we evaluate the sensitivity of GPP and ET to different levels of uncertainty in bud-burst (±10, ±1 days).Introduction

Conclusions References
Tables Figures

Back Close
Full 2 Materials and methods

Site description and phenological observations
The Harvard Forest (42.54In the present analysis we focus on bud-burst dates from 1990 to 2011.We define bud-burst as the date when 50 % of all buds on an individual tree had recognizable leaves emerging (Richardson et al., 2009).Our analysis uses temperature and photoperiod as drivers of phenology.Mean daily air temperatures were computed from the maximum and minimum daily temperatures recorded for the period of study at the Shaler (1964( -2002( ) and Fisher (2001( -2011) ) meteorological stations.Photoperiod was computed as described by a standard equation based on latitude and day of the year (Monteith and Unsworth, 1990).

Phenological models
A large number of different models exist for the simulation of bud-burst for different species (e.g.Chuine et al., 1999;Schaber and Badeck, 2003;Morin et al., 2009).The application of different models against different datasets of the same species, however, Introduction

Conclusions References
Tables Figures

Back Close
Full gave contrasting results about which modeling approach is best (Hunter and Lechowicz, 1992;Chuine et al., 1998Chuine et al., , 1999)).The models provide a context for interpreting observed interannual and inter-specific variability in phenology and to assess which model structure is best supported by the available data.
The models used in this study are largely based on those presented by Chuine et al. (1999) and updated in Richardson and O'Keefe (2009).We define two main categories of models (Table 2).The model categories differ in their assumptions of how warm and cold temperatures control developmental processes (Fig. 1).In the spring warming models, temperatures above a base temperature accumulate until a threshold (in degree-days) is reached, thus triggering bud-burst.In the chilling models, cold weather also plays a role.In the sequential chilling model, a chilling threshold must be reached before warming is effective; in the alternating and parallel chilling models, an increase in the amount of chilling experienced reduces the amount of warming that is required.For all models, the rates of forcing (or Spring Warming) and chilling are calculated based on the threshold approach (CF1) and on the Sarvas' function (CF2) (Sarvas, 1972;Chuine et al., 1999).Either in Spring Warming or Chilling models, photoperiod can be a factor by controlling the point in time (i.e. a day-length threshold) at which chilling and forcing begin to have an effect.
For example, parameter t 2 in Table 2 controls the date at which forcing and/or chilling begins to be accumulated in some models.We also compared versions of the model with no photoperiod control; in this instance (models denoted by the suffix t 0 ; see Table 2) parameter t 2 is fixed to 1 January, which is the value for onset of degree-day accumulation most commonly used in other studies.Introduction

Conclusions References
Tables Figures

Back Close
Full

Model parameters and uncertainty estimations
Model parameters (listed in Table 2) were estimated separately for each species and model using bud-burst observations from Harvard Forest as constrains.This allowed for the characterization of species-specific biological responses to environmental cues.
Model optimization and uncertainty analysis was performed using a model-data fusion platform based on simulated annealing-type routines and Monte Carlo techniques (Metropolis et al., 1953) as described by Richardson et al. (2010).The cost function selected for this purpose was the sum of squared error between observed and modeled data.Once the best parameter set was identified, the parameter space was further explored until 1000 parameter sets that gave statistically equivalent fits to the data were accepted.A specific parameter set was accepted if passed a χ 2 test (at 95 % confidence) for acceptance/rejection after a normalization of the variance based on the minimum of the cost function (Franks et al., 1999).The resulting posterior distributions defined the parameter space within which approximately equally good agreement between data and model simulations can be obtained.Uncertainty estimates for model parameters and model predictions were thus provided directly by the model-data fusion framework, conditional on the data and the cost function.
By running an ensemble of models, with parameter sets selected from the posterior distributions, we can characterize the uncertainty in model predictions (both under current and future climate scenarios) that is attributed to parameter uncertainty.

Model selection and evaluation
We used the Akaike Information Criterion (AIC), a method based on information theory (Akaike, 1973;Anderson et al., 2000) for model selection purposes.AIC is a measure of the trade-off between the goodness-of-fit (model explanatory power) and model Introduction

Conclusions References
Tables Figures

Back Close
Full complexity (number of parameters).The small sample corrected criterion, AICc, developed for cases where the number of samples is small (Burnham and Anderson, 2002), is calculated as in Eq. ( 1): where n is the number of samples (i.e.observation years), p is the number of model parameters and σ 2 is the residual sum of square divided by n.The competing model formulations proposed in Table 2 can be ranked from the best (i.e.lowest AICc) to the less supported by data (i.e.highest AICc).To rank models for each species, the absolute difference in AICc (∆AICc) between alternatives models and the best model have been used (Burnham and Anderson, 2002): if the difference, ∆AICc, is zero or very small, both models are essentially equally likely to be the best model.If ∆AICc is around 2.0, then the model with the lower AICc is almost three times more likely to be best.If, however, ∆AICc is around 6.0, then the model with the lower AICc is about 20 times more likely to be best (Table 3).Finally, for each species, we tested the best model selected by AIC using phenological observations from 2010 and 2011, which had not been used for model calibration.This evaluation was conducted computing the R 2 , the slope of the linear regression analysis (observed vs modeled) as well as the root mean square error (RMSE) between observed and modeled bud-burst dates (Pineiro et al., 2008).

Phenological forecast and propagation of model uncertainty
With the aim of assessing the potential effects of climate change on phenology, we ran the models from 1960 to 2099 using climate projections for Harvard forest.These were generated by Hayhoe et al. (2007) using the NOAA Geophysical Fluid Dynamics Laboratory (GFDL) CM2 global coupled climate model (Delworth et al., 2006), statistically downscaled to approximately 10 km spatial resolution.The CM2 model was run using two scenarios of CO 2 and other greenhouse gas emissions: the IPCC Special Report Introduction

Conclusions References
Tables Figures

Back Close
Full on Emission Scenarios (SRES) higher (A1fi) and lower (B1) scenarios (Nakicenovic et al., 2000).Compared to a 1960-1990 baseline of 7.1 • C mean annual temperature and 1100 mm annual precipitation, corresponding values at the end of simulation (mean 2070-2099) are 12.0 • C and 1270 mm for the A1fi scenario and 9.5 • C and 1240 mm for the B1 scenario (Fig. 2a).
Uncertainty was propagated by running the models forward with the ensemble of parameter sets that passed the χ −2 test at 95 % confidence (Sect.2.3.1),yielding for each year a range of model predictions (Fig. 2b).
The uncertainty in model parameters was analyzed by evaluating the propagation of the uncertainty for each scenario and model structure independently.The ratio of the mean width of the 95 % confidence interval, computed for the last and the first decade of the projections, was used to quantify the degree to which the parameter uncertainty changes between current and future climate conditions.
The uncertainty related to the model structure was characterized by analyzing the average smoothed projected trend across model formulations and climatic scenarios.
The smoothed bud-burst projection was extracted from the time series by using a local polynomial regression fitting (Cleveland and Devlin, 1988).The interannual variability was computed as the standard deviation of the residual between the forecasted budburst and its smoothed time-series.
The uncertainty related to future climate forcing was analyzed by computing, for the best model selected for each species, the trends in advancing bud-burst under the two different climatic scenarios.Trends were characterized by using the Sen's slope estimator and the nonparametric Mann-Kendall test (Helsel and Hirsch, 2002).

Description of BEPS
The Boreal Ecosystem Productivity Simulator (BEPS) was originally developed to simulate the carbon and water fluxes of the Canadian landmass at daily time steps in a Introduction

Conclusions References
Tables Figures

Back Close
Full remote sensing framework (Liu et al., 1997(Liu et al., , 2002)).Derivatives of this original version have been developed and tested at various boreal forest and peatland ecosystems (e.g.Govind et al., 2009Govind et al., , 2011;;Sonnentag et al., 2008;Parton et al., 2000;Ju and Chen, 2005;Ju et al., 2006).Spring phenological events such as bud-burst are modeled as functions of air temperature (Ju et al., 2006), with a model structure classifiable as Spring Warming without photoperiod limitation.After bud-burst (accumulated growingdegree days above 5 • C reach 75 • C), leaf area index keeps increasing linearly up to growing-degree days of 500 • C when prescribed maximum growing season leaf area index for understory and overstory are reached.
For this study, we used the version of BEPS described by Ju et al. (2006), parameterizing the soil-vegetation continuum described by BEPS to consist of five soil and two vegetation layers with site-specific information for Harvard Forest from the literature (Urbanski et al., 2007) and the Harvard Forest data archive.Spanning the period 1960-2099, the model was driven by half-hourly meteorological forcing data including incoming shortwave radiation, air temperature, relative humidity, wind speed and precipitation.
For the future climate projection, we used downscaled data (Hayhoe et al., 2007) from the regionalized projection of the GFDL-CM global coupled climate-land model Delworth et al. (2006) driven with scenario A1fi IPCC (2007).For each run, the model was initialized following the procedure outlined in Ju et al. (2006).-Run 1: BEPS with bud-burst simulated by the native phenological routine (Reference Run);

Modeling strategy
-Run 2: BEPS forced by bud-burst dates simulated by the best model formulation selected according to the AICc as described in Sect.2.3.1 and the optimized parameters set; -Run 3: BEPS forced by 4 different bud-burst dates.The bud-burst forcing time series has been computed by adding to the bud-burst simulated with the Reference Run ±1 and ±10 days.These runs are hereafter referred as BEPS +1 ; BEPS −1 ; BEPS +10 and BEPS −10 .
The differences in cumulated annual GPP and ET between Run 1 and Run 2 allow for the quantification of uncertainty in carbon and water fluxes associated with the BEPS native bud-burst sub-model.The separate runs of Run 3 allow for the characterization of the sensitivity of carbon and water fluxes in BEPS to variations (shift) in bud-burst dates.In other words, after the characterization of the uncertainty around individual future bud-burst dates, we look at the effect of constant "extra" days in spring along the simulation period in terms of carbon and water fluxes as described by BEPS.Hereafter we referred to Run 3 as the "Sensitivity Runs".

Evaluation of phenological models
The AICc values, reported in Table 3, show that models belonging to the class Spring Warming with photoperiod control are overall better supported by the data than competing model structures.In particular, the simple Spring Warming models (SW-CF2 plus SW-CF1) are more often (for 7 species) selected as best models (Tables 1 and 3).
For several species, however, AICc gave support for Chilling models.Among Chilling class models, Alternating (Alt) is the class that is more often selected as the best (3 Introduction

Conclusions References
Tables Figures

Back Close
Full times and 5 times with ∆AICc < 2), while Sequential (Seq) models are selected once.Spring Warming models without photoperiod control are never selected as the best model.Moreover, Table 3 shows that the model SW-CF2 is the model with the lowest AICc across species and then the most generalizable model, followed by Alt-CF1.The time series of bud-burst dates (1960-2099) simulated for red oak with the best model selected (SW-CF2) is shown in Fig. 2b while simulations with different model structures are reported in Fig. 2c.
The R 2 and the slope of the linear regression between observed and predicted budburst dates with the model selected as best for the period 1990-2009 (i.e. years used for models calibration) and for each species are reported in Fig. 3a and b, respectively.
Figure 3c shows that these calibrated models were able to successfully predict 2010 and 2011 bud-burst dates (R 2 = 0.79; RMSE = 4.3 day).

Uncertainty of model parameters
The impact of the uncertainty of model parameters has been evaluated by running the models forward with the 1000 realizations of model parameter sets accepted by the χ 2 test.The uncertainty in individulal years results in an uncertainty in the projected budburst trends.A summary of the resulting uncertainty in the magnitude of the projected bud-burst trends simulated for each species with the best model selected as in Table 1 is represented in the violin plot in Fig. 4. The uncertainty in future trends varied across species and is larger for simulations conducted under the scenario A1fi.The average uncertainty, computed as the average of the 95 % CI, is 2.4 day century −1 (ranging from ±0.7 day century −1 for ACRU to ±4.1 day century −1 for BELE) for scenario B1 and 4.5 day century −1 (ranging from ±0.7 day century −1 for ACRU to ±9.2 day century −1 for BELE) for A1fi.
The average ratio of the mean width of the 95 % confidence interval, computed for the last and the first decade of the projections with the best models is 1.2 (1sd: ± Introduction

Conclusions References
Tables Figures

Back Close
Full 2) for both scenarios.For many models and many species, uncertainty at the end of the simulation is similar to what it is with the present climate.For models without photoperiod limitations (more often for SW-CF1t0 and SW-CF2t0) and some Chilling models (Par and Par2 models), the uncertainty at the end of simulation doubles (ratio between 2.0 and 5.3) and, in particular under scenario A1fi, there is a large increase in uncertainty.

Uncertainty of models structure
The smoothed trends computed for a selection of models, averaged across species, for the two scenarios are reported in Fig. 5.The average of the model ensembles (grey line in Fig. 5) and the average computed with the Akaike's weights (black line in Fig. 5) according to Turkheimer et al. (2003) are also shown.Weighted averages and the standard deviation of trends computed across the model structures are reported in Table 1.The average standard deviation of the projected trends is a measure of the uncertainty in trends related to model structure and varies between ±3.6 day century −1 for B1 and ±7.7 day century −1 for A1fi.
Models more frequently selected as the best (i.e.Spring Warming limited by photoperiod and Alternating models), predict a response to future warming in the middle of the model ensemble.Spring Warming models without photoperiod limitations showed a strong trend toward early bud-burst, particularly relevant under scenario A1fi, that affect also the arithmetic average of multi-model ensemble (grey line in Fig. 5).By 2099, differences in forecast bud-burst date across model structures could reach about 10 days for the A1fi scenario.Chilling models and Spring Warming models with photoperiodic limitation show less steep projected trends.The trends simulated with these two classes of models level-off around 2060 for scenario B1 while for scenario A1fi a slow reduction of the slope in time is observed.For scenario B1, the leveling off reflects the projected decrease of the rate of temperature warming (Fig. 2a) while for scenario A1fi the interaction between model structure, parameters and warm temperatures (see Discussion) contributed to reduce the rate of bud-burst advancing.The Introduction

Conclusions References
Tables Figures

Back Close
Full future sensitivity of bud-burst to temperature is variable across model structures as depicted in the time series in Fig. 6.We observe a stronger sensitivity for Spring Warming models without photoperiod limitation (e.g.SW-CF2t0: −5.3 day • C −1 for A1fi, −5.0 day • C −1 for B1) compared with Spring Warming models limited by photoperiod (SW-CF1: −2.5 day ) and Chilling models (e.g.Alt, −2.2 day • C −1 for both scenarios).For these models (often selected as the best) and for both scenarios, sensitivity to temperature is stable or decrease in time and spans from −2.4 to −2.9 day • C −1 in the 1960s to values ranging from −2.3 to −2.4 day • C −1 by the end of simulation.Under the scenario B1, the average reduction of the bud-burst sensitivity to temperature for the best models selected is about 8.3 % compared to the sensitivity at the beginning of the simulation.Models without photoperiod or chilling limitation show an increase of the sensitivity of bud-burst to temperature by the end of simulation reaching about −6.0 day • C −1 while Parallel models showed a quite variable and higher sensitivity to temperature.
The interannual variability predicted by all models and by all species is reported in Fig. 7, where the blue histograms represent the year-to-year variations under climate scenario B1 and the red histograms under the scenario A1fi.For Chilling models and Spring Warming models without photoperiod limitations, the future interannual variability predicted is larger than the one predicted by Spring Warming models with photoperiod limitation and Alternating models.Within different species we do not observe large differences, except for black cherry.Interannual variability under the scenario A1fi is slightly larger as shown by the differences in red and blue histograms in Fig. 7.However, differences within the 2 scenarios are not statistically significant.

Uncertainty of model drivers
The uncertainty of model drivers is related to the uncertainty in future climate scenarios used for ecological forecasts and here we analyzed the differences in phenological forecasts obtained using the scenarios A1fi and B1.Introduction

Conclusions References
Tables Figures

Back Close
Full The uncertainty of trends for each species was assessed by computing the Sen's slope estimator for the 1000 projected bud-burst time series and reported for the best model for each species (Fig. 4).Projected trends computed for the best model selected for each species for the two scenarios are reported in Table 1.A larger uncertainty in bud-burst trends was evident for scenario A1fi compared to scenario B1.

Uncertainty in Gross Primary Productivity and Evapotranspiration simulated with BEPS
A summary of the differences in ET and GPP simulated with the runs conducted with BEPS are reported in Table 4. On average the BEPS native phenology model simulates an earlier bud-burst of about 17 days (up to a maximum of 59 days) compared to that simulated with the best model formulation selected for the Harvard forest (i.e.

SW-CF2
).An example of the time series of bud-burst simulated with BEPS, SW-CF2 and SW-CF2t0 are shown in Fig. 2c.Differences between the Reference Run and BEPS SW−CF2 represent the impact of uncertainty in model structure on carbon and water fluxes modeled.We observe an average overestimation of the Reference Run of 136.0 (±59.2, 1 sd) gCm −2 day −1 yr for GPP (±9.6 % of the annual GPP) and 1.38 (±0.77, 1 sd) cm yr −1 for ET.The correlation between the differences in fluxes (∆GPP and ∆ET) and in bud-burst (∆BB) computed between the Reference Run and BEPS SW−CF2 is stronger for GPP (r = −0.85p < 0.001, slope = −6.49gCm −2 day −1 ) than for ET (r = −0.55,p < 0.001, slope = −0.055cm day −1 ). Figure 8 shows that ∆ET and ∆BB are strongly negatively correlated (r = −0.82,p < 0.001) when the bud-burst occur late (bud-burst > 115) while are poorly correlated (r = −0.38,p < 0.01) when the bud-burst occur early (bud-burst < 115).This highlights that year-to-year variations in GPP are well correlated to year-to-year variations of bud-burst dates, while year-to-year 895 Introduction

Conclusions References
Tables Figures

Back Close
Full variations of ET depend on variations in bud-burst in years with late bud-burst but not in years with early bud-burst, suggesting the role of others meteorological factors controlling ET late in the season when the bud-burst occurs early.
The results of the sensitivity runs (i.e.differences between the Reference Run and BEPS −1 , BEPS +1 , BEPS −10 and BEPS +10 ) are reported in Table 4 and represent the sensitivity of annual GPP and ET to the uncertainty in bud-burst dates due to model parameters uncertainty (reported as sensitivity of fluxes as consequence of ±1 day of bud-burst uncertainty).The time series of the residuals (∆GPP, ∆ET) to the Reference Run are reported in Fig. 9.We observe an average increase in annual GPP of 4.3 % (67.3 ± 16.4 (1 sd) gCm −2 ) for 10 days of advance in bud-burst and a decrease −5.1 % (−80.5 ± 20.3 (1 sd) gCm −2 ) for 10 days of delay in bud-burst.The sensitivity runs confirm that ET respond less to variations of bud-burst than GPP: ET varies of 0.92 % (0.42 ± 0.25 (1 sd) cm) and −1.6 % (−0.76 ± 0.38 (1 sd) cm) for variation in bud-burst of −10 and +10 days, respectively.The average sensitivity of annual GPP to ±1 day of variation of bud-burst simulated by BEPS is 8.15 gCm −2 day −1 , while for ET is 0.07 cm day −1 .The sensitivity of ET to ±10 days of variations in bud-burst is asymmetric (green and grey boxes in Fig. 9) and larger for +10 days of bud-burst variations.Moreover, the year-to-year variations of ∆ET is larger compared to ∆GPP.
A schematic representation of the within-year sensitivity of ∆GPP and ∆ET to variations in ∆BB is reported in Fig. 10a, b and confirms that the sensitivity of GPP is similar both for early and late bud-burst, while the sensitivity of ET is less pronounced for the run BEPS −10 and in particular for the years with earlier bud-burst (Blue line in Fig. 10b).Within-years sensitivity of the average soil moisture in summer to ∆BB and the within-years sensitivity of the percentage of variations of summer soil moisture to ∆BB are shown in Fig. 10c, d.The average soil moisture in summer is lower for years with earlier bud-burst and the sensitivity of summertime soil moisture to variations in bud-burst dates is lower for seasons with earlier bud-burst.Introduction

Conclusions References
Tables Figures

Back Close
Full Table 5 summarizes the correlation and the slopes of the linear regression computed between ∆GPP and ∆ET simulated with the different sensitivity runs and different meteorological and environmental conditions.Year-to-year variations in ∆GPP are correlated with the bud-burst dates, spring and mean annual temperature and the sensitivity (i.e.slope of the linear regression analysis) is stronger for the run BEPS −10 .The interannual variability of ∆ET is poorly correlated with environmental conditions except for temperature and snow depth for the run BEPS −10 .The interannual variability of summer ∆ET is more correlated with the main meteorological variables (Table 5) compared to annual ∆ET.

Sources of uncertainty of phenological forecasts
Our analysis has characterized different sources of uncertainty in phenological forecasts.Here we have explicitly accounted for: uncertainty in estimated model parameters (parameter uncertainty), uncertainty due to model structure (model uncertainty) and uncertainty due to the forcing drivers (driver uncertainty).
Our analysis suggests that, once model parameters are optimized according to the data, parameter uncertainty is the smallest (average 95 % CI: ±2.4 day century −1 for scenario B1 and ±4.5 day century −1 for A1fi) and the most straightforward to quantify by using a model-data fusion framework (model-data-optimization routine), whereas driver uncertainty is the largest (up to 8.4 day century −1 of difference in trends computed with the two scenarios and the best models), and the most complicated to quantify given the uncertainty in the future climate.Model uncertainty is also large and comparable with the driver uncertainty: the predicted bud-burst trends, as well as the shape of the smoothed projections (Fig. 5), showed a large variability across models (1sd: ±3.6 day century −1 for B1 and ±7.7 day century −1 for A1fi).Introduction

Conclusions References
Tables Figures

Back Close
Full With the best model selected for each species and its optimized parameters set, an average trend towards earlier spring by 4.4 day (scenario B1) and 8.8 day (scenario A1fi) over the coming century is predicted.
Although the parameter uncertainty is smaller than other sources of uncertainty, we observe that in some species (such as yellow birch, white oak, red oak and sugar maple), forecasted bud-burst trends vary greatly depending on the parameters set used in some models.
Forecast uncertainty due to parameter uncertainty scales with time although does not increase for all the model structures.Increasing uncertainty by the end of simulation is model-dependent and it is expected mainly for those model structures less supported by data (e.g.Spring Warming without photoperiodic limitation and Parallel models).Whether these models are used for long-term forecast of spring phenology, the nonstationary of uncertainty in time should be carefully considered.
The above considerations suggest an interaction between the parameter uncertainty and the driver uncertainty.More in detail, the parameter uncertainty scales depending on the climate scenario used, with higher uncertainty for the high CO 2 emissions scenario (Fig. 4).The expansion of the propagated uncertainty under A1fi scenario necessitates efforts toward the reduction of the uncertainty in model parameters and, more important, in the driver uncertainty (i.e.climate scenario).
The driver uncertainty is related to uncertainty in future climate scenarios used to run the phenological models.The driver uncertainty is enhanced for models less supported by data (e.g. up to 2 weeks for Spring Warming class models without photoperiod limitation and for some versions of Parallel and Sequential models).The mean temperature increases expected with the scenario A1fi is about 3 times larger than the one expected for B1.However, some models predict an adjustment of the apparent sensitivity of bud-burst to temperature over the next century (Fig. 6).In general, there is a reduction in the amount that bud-burst advances for a given increase in temperature.We note that this should not be seen as an adaptation of trees to climate change.Rather, we believe this reflects photoperiod and/or chilling constrains (according to Introduction

Conclusions References
Tables Figures

Back Close
Full the species), which effectively limit the degree to which the bud-burst can advance in spring.
Here, we focused on the uncertainty related to the socio-economic scenario used (A1fi vs B1).However, regarding climate data, other uncertainties can also be explored.Cook et al. (2010), showed the impact of the temporal and spatial resolution of climate scenario in phenological modeling, as well as the uncertainty related to the use of different general circulation models (GCM).This uncertainty could be assessed by reporting the ensemble spread obtained by running one phenological model with climate drivers from different GCMs.
The model uncertainty is comparable to driver uncertainty and both the predicted bud-burst trends, and the shape of the smoothed projections, show a large variability across models (Fig. 5).In this regard, the importance of photoperiod limitation is highlighted by (1) the large difference (up to 10 days by the end of our simulations) between Spring Warming models with and without photoperiod limitation, and (2) by the similarity of the sensitivity of bud-burst to temperature for models that feature photoperiod limitation.
For some Chilling class models (e.g.Parallel Models) with specific set of parameters we observe that in some years the predicted bud-burst dates were highly uncertain and unrealistic (e.g.no bud-burst predicted).This effect, already noted in an earlier study (Morin et al., 2009), might be related to interactions between parameter uncertainty and model structural uncertainty.Thus, in anomalous years (i.e.winter very warm), for certain extreme parameter sets, the winter chilling requirements were not fulfilled and models failed to predict bud-burst in spring (i.e.no bud-burst predicted).This occurred more often for scenario A1fi than scenario B1, highlighting the need to reduce this source of uncertainty by either reformulating the model structure or better constraining model parameters by increasing the number of available observational datasets.
To summarize, certain model structures and parameterizations resulted in predicted bud-burst dates that are not sensitive to warm temperatures during winter and early spring, which limits the degree to which modeled bud-burst can advance under future Introduction

Conclusions References
Tables Figures

Back Close
Full climate scenarios.Up to a certain level of projected warming in climate, different model structures predict a different absolute value of bud-burst but a similar behavior of the shape and trend of the projections.In these condition, spring temperature is the main driver, while photoperiod and chilling are less limiting.For high CO 2 emissions scenarios, instead, photoperiod and chilling become important as additional limiting factors leading to a different shape of projected spring phenology.Our model selection (Spring Warming with photoperiodic limitation have been selected more frequently than Chilling models) supports the position of K örner and Basler (2010) in the current debate about the interactions between warm temperatures in spring, cold temperatures during dormancy and photoperiod and the degree at which they will control phenology under climate warming scenario.However, it should be considered that, although this study is conducted over a relatively long time series, it includes only one site, which may limit our ability to effectively constrain the photoperiod parameter.Moreover, one of the major limitations related to currently available observational datasets is that these might be too short to span a wide range of climatic conditions needed to constrain more complex models (Richardson et al., 2012).For this reason exploiting, in a model-data fusion framework, data from warming experiments or common garden experiments (e.g.Morin et al., 2010) as well as new long-term dataset of consistent phenological observations from across a wide geographic range and climatic conditions, would of great value for constraining model parameterization (i.e. to reduce parameters uncertainty) and for testing and developing phenology models (i.e. to reduce model uncertainty).Finally, the uncertainty described above and multi-model trends reported in Fig. 5 suggest that when constructing a multi-model ensemble, not all model structures should be weighted equally.More weight should be given to those models that are better supported by the data, and less weight to those that are not supported by the data.Akaike weights Turkheimer et al. (2003) can be used for this.Introduction

Conclusions References
Tables Figures

Back Close
Full

Impacts of uncertainty in bud-burst forecasts on carbon and water cycle modelling
The uncertainties of phenological forecasts might have important implications for the land surface energy budget and for carbon and water cycling modeling because of the sensitivity of biogeochemical models to bud-burst dates (Levis and Bonan, 2004;White et al., 2000).The prognostic routine in commonly used terrestrial biosphere models (e.g.Thornton et al., 2002;White et al., 1997;Lawrence et al., 2011;Ju et al., 2006) is often based on the heat sums or growing degree-day assumption (i.e.roughly similar to models SW-CF1 t0 and SW2-CF2 t0) and model parameters are not optimized for the application on different species.Model structures relying only on the degree-days description at the Harvard Forest (as in other studies) are the least supported by data being the Spring Warming models with photoperiodic control or some Chilling class models selected as best models.This might lead to a poor description of spring phenology either in terms of bud-burst or seasonal development of leaf area index (Richardson et al., 2012;Lawrence et al., 2011).Hence, the seasonality of carbon, water and energy fluxes projected by these models are likely incorrect as well (Richardson et al., 2012).
In our modeling exercise with BEPS we explored the effects on water and carbon cycling due first to the uncertainty in phenological model routines not optimized embedded in the model (i.e.differences between the Reference Run and BEPS SW−CF2 ) and second to the uncertainty in model parameters (i.e.sensitivity runs: BEPS −10 , BEPS −1 , BEPS +1 , BEPS +10 ).
Model uncertainty (i.e.differences in GPP and ET obtained with the Reference Run and the BEPS SW−CF2 ) lead to biases of 9.6 % for annual GPP and 2.9 % for annual ET indicating the importance of improving phenological routines.
Beside the uncertainty related to model structure, phenological forecasts for future climate scenarios are uncertain because of the uncertainty in phenological model parameters that can be considered, conservatively, as maximum ±10 days.The impact Introduction

Conclusions References
Tables Figures

Back Close
Full of this uncertainty on terrestrial biosphere models has been explored with the sensitivity analysis and showed (Table 4) that a variation of ±10 and ±1 day in bud-burst dates lead to a variation of about ±5 % on annual GPP (8.15 gCm −2 day −1 , i.e. sensitivity of annual GPP to 1 day of variation in bud-burst) and about ±2 % on annual ET (0.07 cm −1 day −1 , i.e. sensitivity of annual ET to 1 day of variation in bud-burst).
About the sensitivity of GPP to bud-burst uncertainty, values found here are comparable with recent empirical estimates of the sensitivity of gross productivity to phenology: Richardson et al. (2010) found a relationships between interannual phenological anomalies in spring onset and interannual CO 2 flux anomalies ranging from 1.20 (se: ±6.20) to 20.2 (se: ±2.90) gCm −2 day −1 (Supplementary Table 5 in Richardson et al., 2010); Richardson et al. (2009a), in an empirical studies in two different North American forests (Harvard Forest and Howland) indicated that a one-day advance in spring onset date was associated with an increase in GPP of 12.30 (se: ±2.50) gCm −2 day −1 .
Although the impact on GPP of uncertainty in model parameters is lower than uncertainty due to phenological model structure, this is of the same magnitude (or slightly lower) of the sensitivity of GPP to anomalies in bud-burst reported in empirical studies.This emphasizes the importance of correctly modeling phenological transition dates in order to successfully predict annual CO 2 uptake, as well as other biosphereatmosphere interactions and climate system feedbacks.
As shown in Fig. 10a, b the sensitivity of GPP to bud-burst is rather constant while the sensitivity of ET to variations in bud-burst decreases substantially for the runs BEPS −10 , in particular in years with earlier bud-burst.We observe that an earlier budburst of 10 days in spring it is not always associated to an increase of ET given the reduction of soil moisture and its sensitivity to variations in bud-burst (Fig. 10b, c, d).
The effect is a reduction of ET (and of its sensitivity) and might be associated to a late summer water deficit due to a depletion of the soil water reservoir (Fig. 10d) due to the long growing season.As shown in Fig. 10c, according to the BEPS structure, the earlier is the bud-burst the lower is the average soil moisture in summer.Moreover, the sensitivity of summer soil moisture to variations of ±10 days in bud-burst (Fig. 10d) is Introduction

Conclusions References
Tables Figures

Back Close
Full very little in years with early bud-burst, increasing in years with late bud-burst.A similar behavior has been found in boreal forest stands (Kljun et al., 2006) and in broadleaved forest in Central Europe (Leuzinger et al., 2005) where an earlier spring phenology increased transpiration rates, leaving less moisture in the soil in summer.However, in these studies was reported a reduction of summer (and annual) productivity which is not observed with BEPS.The carbon cycle is therefore more affected by, and sensitive to, uncertainty in phenology than the water cycle (Figs.8, 9 and 10).Regarding water cycling, Lawrence and Slingo (2004) and White et al. (1999), showed that the three primary components of evapotranspiration, namely, transpiration, soil evaporation and canopy evaporation respond in a opposite direction to changes in phenology.For instance, a certain amount of soil evaporation can occur independently of the seasonal development of leaves in years with late bud-burst.Or, while bud-burst and spring leaf development enhances transpiration, it may limit surface soil evaporation because of the lower amount of transmitted solar radiation below the canopy.Moreover, in years with early bud-burst, the enhancement of evapotranspiration in spring might lead to lower soil moisture in summer limiting the increase of transpiration in summer and then buffering the sensitivity of total ET to bud-burst.However, these considerations are based on modeling studies and empirical studies would be useful to better understand the sensitivity of ET to spring phenology.
Although the uncertainty forced in BEPS was kept constant (i.e.±10 days), a large interannual variability of the sensitivity of carbon and water fluxes to bud-burst has been observed (Fig. 10).The interannual variability of GPP depends on environmental conditions such as springtime and annual temperature and snow-depth, which is related to water availability during the spring growth.The same amount of uncertainty (e.g. 10 days) might have different effects in terms of carbon and water cycle depending on the timing of bud-burst.In models with a tight link between water and carbon cycle the earlier bud-burst and its associated uncertainty can lead to an increase of the probability of summer water deficit and a consequently reduction of carbon uptake.Thus,

Conclusions References
Tables Figures

Back Close
Full the description of the feedbacks between spring phenology and carbon cycle under a warming scenario might be hampered.

Conclusions
The evaluation of different models with a model-data fusion approach and the long-term Harvard Forest phenology record indicate support for Spring Warming models with photoperiod limitation and to a less extent Chilling models with alternating description of chilling accumulation.Models without the explicit description of the photoperiod limitation are never supported by data.We characterize and quantify the three different sources of uncertainty in phenological forecasts: parameter uncertainty, model uncertainty, and driver uncertainty.Parameter uncertainty is the smallest, and most readily quantified whereas driver uncertainty is larger and the most unquantifiable.Model structure uncertainty is comparable with driver uncertainty, but can be quantified by combining models with different complexity and observational data.The uncertainty in model structure affects the sensitivity of bud-burst to temperature simulated by phenology models.
Phenology regulates many ecosystem feedbacks to climate, and here we quantify the impact of uncertainties in bud-burst forecasts to simulated carbon and water cycling.We conclude that the structure and the optimization of parameters of the phenological routine is a relevant source of uncertainty for the simulation conducted by terrestrial biosphere models.This highlights the importance of developing better phenological routines.The carbon cycle is more sensitive to uncertainty in phenology than the water cycle, and the impact of uncertainty in phenological model parameters on carbon cycling is of the same magnitude of the sensitivity of productivity to year-toyear variations of bud-burst.For evapotranspiration, the impact of this uncertainty is more relevant in years with late bud-burst for evapotranspiration given the asymmetric response of water cycling to variations of bud-burst observed with our model.

Conclusions References
Tables Figures

Back Close
Full To reduce model uncertainties it is necessary to continuously strive to re-evaluate model predictions against new observational datasets.Of particular value would be data from diverse populations (i.e.including genetic variability) growing under as divergent environmental conditions as possible.To push models to their limits, these data should include variations in both temperature (data from transplant or warming experiments would be particularly valuable) and, to the fullest extent possible, photoperiod.These data will permit the development of better phenological routines, and the reduction of uncertainties in forecasts of ecosystem responses to climate change that are mediated by phenology.

Conclusions References
Tables Figures

Back Close
Full  Full  Full  Full   Model representation of the seasonal cycle of terrestrial ecosystems, from senescence through dormancy (a period of rest, followed by quiescence) and then active growth.In temperate and boreal systems, the transition through dormancy to active growth has been described using a variety of approaches.Threshold points t 0 and t 1 may be triggered either by photoperiod or temperature (i.e., various forms of chilling, e.g.sequential, alternating or parallel, as indicated).At threshold point t2, bud-burst is triggered when accumulated forcing reaches a critical state.The manner in which chilling and forcing accumulates varies among models.
Chilling models describe phenology as a combination of temperature forcing, photoperiod limitation and chilling limitation; Forcing models describe phenology as function of temperature forcing and photoperiod limitation (Tab. 2 for a list of model).Modified from Chuine (2000). 34 Fig. 1.Model representation of the seasonal cycle of terrestrial ecosystems, from senescence through dormancy (a period of rest, followed by quiescence) and then active growth.In temperate and boreal systems, the transition through dormancy to active growth has been described using a variety of approaches.Threshold points t 0 and t 1 may be triggered either by photoperiod or temperature (i.e.various forms of chilling, e.g.sequential, alternating or parallel, as indicated).At threshold point t 2 , bud-burst is triggered when accumulated forcing reaches a critical state.The manner in which chilling and forcing accumulates varies among models.Chilling models describe phenology as a combination of temperature forcing, photoperiod limitation and chilling limitation; Spring Warming models describe phenology as function of temperature forcing and photoperiod limitation (Table 2 for a list of model).Modified from Chuine (2000).Introduction

Conclusions References
Tables Figures

Back Close
Full  ) and slope of the linear regression between observed and predicted bud-burst with the best model selected (Table 1) for each year (a), for each species (b).Vertical lines represent the average bud-burst day of the year in each year (a) and for each species (b).Predicted versus Observed bud-burst dates for the years 2010-2011 simulated for the 12 species and the best model selected as described in Table 1 (c).Error bars for the predictions represent the propagated uncertainty of model parameters while for the observed dates is the minimum and maximum bud-burst dates observed for each individual.
36 Fig. 3. Determination coefficient (R 2 ) and slope of the linear regression between observed and predicted bud-burst with the best model selected (Table 1) for each year (a), for each species (b).Vertical lines represent the average bud-burst day of the year in each year (a) and for each species (b).Predicted versus Observed bud-burst dates for the years 2010-2011 simulated for the 12 species and the best model selected as described in Table 1 ( BB projections were modeled by using the best model selected and the parameters ensemble (1000 parameter set) for each species (Table 1).The white circle represent the mean of the distribution.Red violins represent the predictions under the Scenario A1fi, while blue violins represent the predictions under the Scenario B1.
37 Fig. 4. Violin-plot of the trends computed with the Sen's slope estimator for each species.Bud-burst dates (BB) BB projections were modeled by using the best model selected and the parameters ensemble (1000 parameter set) for each species (Table 1).The white circle represent the mean of the distribution.Red violins represent the predictions under the Scenario A1fi, while blue violins represent the predictions under the Scenario B1.Introduction

Conclusions References
Tables Figures

Back Close
Full   Full reported several trends of start of season at global and Figures Back Close Full Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | due to (1) model parameters; (2) model structure; and (3) model drivers (i.e.uncertainty of future climate).The evaluation of phenological model parameter uncertainty is necessary in order to estimate the uncertainty in performance of a particular model structure and for parameter optimization; model structural uncertainty stems from different Figures Back Close Full Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | might be subject to uncertainty due to the uncertain representation of phenological responses to climate change.
Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Several model runs have been conducted over the period 1960-2099 with BEPS and forced by different bud-burst forecasts.Due to the computational demand of running BEPS with the full posterior distribution of uncertainty from each model, we designed specific experiments to test both the BEPS native phenology model, and the sensitivity of carbon and water cycles estimated by BEPS to errors in simulated bud-burst dates.The different runs are described below: Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | 0.
Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | H., Vargas, R., Verbeeck, H., Xiao, J., and Xue, Y.: Terrestrial biosphere models need better representation of vegetation phenology: Results from the North American Carbon Program-Site Synthesis, Glob.Change Biol., doi:10.1111/j.1365-2486.2011.02562.x, in press, 2012.883, 900, 901 Sarvas, R.: Investigations on the annual cycle of development on forest trees active period, Discussion Paper | Discussion Paper | Discussion Paper | White, M. A., Thornton, P. E., and Running, S. W.: A continental phenology model for monitoring vegetation responses to interannual climatic variability, Global Biogeochem.Cy., 11, 217-234, 1997.901 White, M. A., Running, S. W., and Thornton P. E.: The impact of growing-season length variability on carbon assimilation and evapotranspiration over 88 years in the eastern US deciduous Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Fig.1.Model representation of the seasonal cycle of terrestrial ecosystems, from senescence through dormancy (a period of rest, followed by quiescence) and then active growth.In temperate and boreal systems, the transition through dormancy to active growth has been described using a variety of approaches.Threshold points t 0 and t 1 may be triggered either by photoperiod or temperature (i.e., various forms of chilling, e.g.sequential, alternating or parallel, as indicated).At threshold point t2, bud-burst is triggered when accumulated forcing reaches a critical state.The manner in which chilling and forcing accumulates varies among models.Chilling models describe phenology as a combination of temperature forcing, photoperiod limitation and chilling limitation; Forcing models describe phenology as function of temperature forcing and photoperiod limitation (Tab. 2 for a list of model).Modified fromChuine (2000).
Fig. 2. a) Time-series of mean annual temperature projections (1960-2099) under the Scenario A1 (orange line) and Scenario B1 (blue line); b) Time series of bud-burst date (BB) modeled for red oak (Quercus rubra) with the best model selected (Tab 1).Blue line represents the BB projected under the Scenario b1, orange line represents the BB projected under Scenario A1fi.Grey area represents the uncertainty for the Scenario A1, while light purple area represents the uncertainty for Scenario B1; c) Time series of bud-burst date (BB) modeled for red oak with SW-CF2 (best model), SW-CF2t0 (as SW-CF2 but without photoperiod limitation) and with the internal phenological routine in BEPS (green line).Grey crosses represent the phenological observations collected at the Harvard Forest.Models are forced by temperatures from the scenario A1.

Fig. 2 .Fig. 3 .
Fig. 2. (a) Time-series of mean annual temperature projections (1960-2099) under the Scenario A1fi (orange line) and Scenario B1 (blue line); (b) time series of bud-burst dates (BB) modeled for red oak (Quercus rubra) with the best model selected (Table1).Blue line represents the BB projected under the Scenario B1, orange line represents the BB projected under Scenario A1fi.Grey area represents the uncertainty for the Scenario A1fi, while light purple area represents the uncertainty for Scenario B1; (c) time series of bud-burst date (BB) modeled for red oak with SW-CF2 (best model), SW-CF2t0 (as SW-CF2 but without photoperiod limitation) and with the internal phenological routine in BEPS (green line).Grey crosses represent the phenological observations collected at the Harvard Forest.Models in Fig.2care forced by temperatures from the scenario A1fi.

Fig. 4 .
Fig.3.Determination coefficient (R 2 ) and slope of the linear regression between observed and predicted bud-burst with the best model selected (Table1) for each year (a), for each species (b).Vertical lines represent the average bud-burst day of the year in each year (a) and for each species (b).Predicted versus Observed bud-burst dates for the years 2010-2011 simulated for the 12 species and the best model selected as described in Table1(c).Error bars for the predictions represent the propagated uncertainty of model parameters while for the observed dates is the minimum and maximum bud-burst dates observed for each individual.

Fig. 5 .Fig. 5 .
Fig. 5. Smoothed bud-burst (BB) projected time series under the Scenario a1 (a) and the Scenario b1 (b).Different lines represents different models listed in Tab 1. Orange lines represent Spring Warming models with photoperiod limitations; green lines represent Spring Warming models without photoperiod limitations (starting dates fixed at 1 s t January); blue lines represent some chilling class models with (Alternating, Alt) and without photoperiod limitation (Alt-CF1t0 and Sequential, Seq-CF1).Black line represent the Akaike weighted average of model ensembles wile grey line represent the average.

Fig. 6 .Fig. 7 .
Fig. 6.Future sensitivity of bud-burst to temperature (d BB/d T ) as projected by different phenological models under the Scenario A1fi (a) and the Scenario B1 (b).Different lines represents some of the most representative models listed in Table2.

Fig. 8 .Fig. 9 .
Fig.8.Relationships between differences in BB (∆BB) simulated with the best model (SW-CF2 as described in Table1) and the internal phenological routine of BEPS and the differences in annual ET (a) and GPP (b) simulated with Reference Run and BEPS SW−CF2 (Table4).In (a) circles are represented with different colors according to the BB date simulated with the internal routine of BEPS.

Fig. 10 .
Fig. 10.(a) Within-years sensitivity of GPP (∆GPP) to variations of BB (∆BB); (b) within-years sensitivity of ∆ET to ∆BB; (c) within-years sensitivity of the average soil moisture in summer to ∆BB; (d) within-years sensitivity of the percentage of variations of summer soil moisture to ∆BB.Different lines represent the within-years sensitivity.Years are grouped according to the distribution of the day of the year (DOY) of bud-burst (i.e. from years with early bud-burst to years with late bud-burst as reported in the legend).The purple line represent the within-years sensitivity computed using all data.
W, el.220 410 m a.s.l.) site used in this study is located in central Massachusetts, about 100 km west of Boston, USA.The climate is classified as humid-continental, with a mean July temperature of 20• C and mean January temperature of 7 • C. Mean annual precipitation is 1100 mm, and is distributed evenly across the seasons.Forests are dominated by transition hardwoods: red oak (Quer- • N, 72.18 • cus rubra), red maple (Acer rubrum) black oak (Quercus velutina), white oak (Quercus alba) and yellow birch (Betula alleghaniensis).Conifers include eastern hemlock (Tsuga canadensis), red pine (Pinus resinosa) and white pine (Pinus strobus).Since 1990, springtime phenology observations have been made at 37 day intervals.Leaf development was monitored on three or more individuals (a total of 39 permanently marked trees or shrubs) of 11 woody species (Table1).Phenological observations used here are available online (http://harvardforest.fas.harvard.edu).

Table 1 .
List of the 11 species used in the analysis, species number (Spec), Species identifier (Species ID), species Latin and common name, average bud-burst date (BB), best model according to data (Best Model), determination coefficient of fitting (R 2 ), root mean square error (RMSE, day), Sen's slope trend estimated for Scenario A1fi and Scenario B1 with the best model, with the Akaike Weigthed average (Average wAICc) and the standard deviation.All trends are expressed in day century −1

Table 3 .
∆AICC values for a range of different models (see text and Tables1, 2for additional information) fit to Harvard Forest bud-burst data.Species ID are as given in Table1.The best model, based on Akaikes Information Criterion corrected for small samples (AICc) has ∆AICc = 0 and is indicated by bold type.

Table 4 .
Summary of BEPS simulation protocol.Bud-burst forcing represents the bud burst dates used for model runs.GPP is the mean Gross Primary Productivity, ET is the mean evapotranspiration while GPP and ET sensitivity represent the average difference to the Reference Run for annual GPP and annual ET expressed as variations of gCm −2 and cm of water per 1 day of variation of bud-burst.The Reference Run is the run conducted with BEPS and with bud-burst simulated with its internal phenological routine.

Table 5 .
Slope of the linear regression analysis computed between ∆GPP, ∆ET and summertime ∆ET (from June to August) and different meteorological variables.In parentheses the Pearson's correlation coefficient is reported.The differences of GPP (∆GPP) and ET (∆ET) are computed as the differences in annual GPP [gC m −2 yr −1 ] and ET [cm yr −1 ] between the runs BEPS −10 , BEPS −1 , BEPS +1 , BEPS +10 and the Reference Run (Table3).Tann is the mean annual temperature, Tspring is the springtime mean temperature (i.e. from March to May), Tsoil is the mean soil temperature, Prec [mm] is the annual cumulated precipitation and Snow [cm] is the average snow depth.Bold characters represent statistical significant correlations (p < 0.001).

Table 2 (
Histogram of the interannual variability predicted by all model structures reported in Table2(left panel) and for all the species reported in Table1(right panel).Red histograms represent the interannual variability under Scenario A1fi, blue histograms represent the interannual variability under Scenario B1.The standard deviation (1 sd) of the interannual variability is also reported.Introduction left panel) and for all the species reported in Table1 (right panel).Red histograms represent the interannual variability under Scenario A1fi, blue histograms represent the interannual variability under Scenario B1.The standard deviation (1 sd) of the interannual variability is also reported 40 Fig. 7.