This is the accepted version of the following article: Peaucelle, M. et al. Representing explicit budburst and senescence processes for evergreen conifers in global models in Agricultural and forest meteorology (Ed. Elsevier), vol. 266-267 (March 2019), p. 97-108. Which has been published in final form at DOI 10.1016/j.agrformet.2018.12.008 © 2018. This manuscript version is made available under the CC-BY- NC-ND 4.0 license http://creativecommons.org/licenses/by-nc- nd/4.0/ 1 Title: Representing explicit budburst and senescence processes for evergreen conifers in 2 global models 3 4 Running title: A coniferous phenology model for global predictions 5 Authors : Marc Peaucelle1, Philippe Ciais1, Fabienne Maignan1, Manuel Nicolas2, Sébastien 6 Cecchini2, Nicolas Viovy1 1 7 Laboratoire des Sciences du Climat et de l'Environnement, CEA CNRS UVSQ, Gif-sur-Yvette, 8 France 2 9 Office National des Forêts, Département Recherche Développement et Innovation, Fontainebleau, 10 France 11 12 13 Marc Peaucelle: marc.peaucelle@lsce.ipsl.fr (corresponding author, phone:+33 169 088 739 14 /fax:+33 169 087 716 ) 15 Philippe Ciais: philippe.ciais@lsce.ipsl.fr 16 Fabienne Maignan : fabienne.maignan@lsce.ipsl.fr 17 Manuel Nicolas : manuel.nicolas@onf.fr 18 Sébastien Cecchini : sebastien.cecchini@onf.fr 19 Nicolas Viovy: nicolas.viovy@lsce.ipsl.fr 20 21 22 Key words: tree phenology, needle budburst, needle senescence, evergreen needleleaves, 23 ecosystem model, RENECOFOR 24 25 Type of paper: Primary research article 26 27 28 29 30 31 1 32 Abstract: 33 Global ecosystem models lack an explicit representation of budburst and senescence for evergreen 34 conifers despite their primordial role in the carbon cycle. In this study we evaluated eight different 35 budburst models, combining forcing, chilling and photoperiod, for their ability to describe spring 36 budburst, and one model of needle senescence for temperate evergreen coniferous forests. The 37 models’ parameters were optimized against field observations from a national forest monitoring 38 network in France. The best fitting budburst model was determined according to a new metrics 39 which accounts for both temporal and spatial variabilities of budburst events across sites. The best 40 model could reproduce observed budburst dates both at the site scale (±5 days) and at regional scale 41 (±12 days). We also showed that the budburst models parameterized at site scale lose some 42 predictive capability when applied at coarser spatial resolution, e.g., in grid-based simulations. The 43 selected budburst model was then coupled to a senescence function defined from needle 44 survivorship observations in order to describe the full phenology cycle of coniferous forests. 45 Implemented in the process-driven ecosystem model ORCHIDEE, this new conifer phenology 46 module represented accurately the intra and inter-annual dynamics of leaf area index at both the 47 local and regional scales when compared against MODIS remote sensing observations. A sensitivity 48 analysis showed only a small impact of the new budburst model on the timing of the seasonal cycle 49 of photosynthesis (GPP). Yet, due to the faster renewal of needles compared to the standard version 50 of ORCHIDEE, we simulated an increase in the GPP by on average 15% over France, while the 51 simulated needle turnover was doubled. Compared to 1970-2000, projections indicated an 52 advancement of the budburst date of 10.3±2.8 and 12.3±4.1 days in average over the period 2060- 53 2100 with the best forcing and chilling-forcing models respectively. Our study suggests that 54 including an explicit simulation of needle budburst and senescence for evergreen conifers in global 55 terrestrial ecosystem models may significantly impact future projections of carbon budgets. 2 56 1. Introduction 57 58 The phenology of conifers is strongly correlated with local climate (Steiner, 1980; Worrall, 1983; 59 Burr et al., 1989; Leinonen & Hänninen 2002; Hänninen et al., 2007). A number of previous studies 60 concluded that the growing-season length of conifer forests will extend with climate warming and 61 rising CO2 concentration, thus leading to significant modifications of biogeochemical processes 62 being controlled by phenology (Murray et al., 1994; Polgar & Primack, 2011; Gunderson et al., 63 2012; Migliavacca et al., 2012; Richardson et al., 2013). For both deciduous and evergreen species, 64 phenology is commonly divided into three different phases: bud dormancy, bud flush and 65 senescence. The timing of these events partly controls the seasonal cycle of leaf area index and 66 gross primary productivity (GPP; Chen et al., 2016). It also impacts albedo, evapotranspiration, and 67 litter inputs to the soil, of which the latter affects soil respiration (Richardson et al., 2013). On 68 longer time scales, phenology also impacts the competitiveness of a species and its spatial 69 distribution (Baldocchi et al., 2001; Chuine & Beaubien, 2001; Polgar & Primack, 2011). These 70 impacts make it essential to represent phenological events accurately in ecosystem models, both in 71 space and in time, if we seek to improve the simulation of the future role of vegetation in carbon, 72 water and nutrient cycling and its feedbacks on climate. Richardson et al. (2012) pointed out to 73 shortcomings in the representation of phenological processes in global vegetation models. Almost 74 all of the 14 vegetation models with different phenology parameterizations that they compared 75 overestimated the length of the growing season and consequently the GPP for temperate and boreal 76 forests. 77 The state of the understanding of complex molecular pathway processes of dormancy and budburst 78 mechanisms is insufficient to allow a fully mechanistic simulation in global models (Rohde & 79 Bhalerao, 2007; Yakovlev et al., 2008; Rinne et al., 2011; Cooke et al., 2012). Thus, empirical 80 models are used for estimating the response of budburst to temperature. Several conceptual models 81 using temperature to determine the date of budburst have been proposed: they fall into two broad 82 categories. The first category assumes that budburst occurs after a threshold of degree-days has 3 83 been reached during a specific period (forcing). The second class of models assumes that budburst 84 requires both a chilling period during winter followed by a forcing from increasing temperature. We 85 investigated models belonging to both these categories. For conifers, photoperiod, in combination 86 with temperature, has also been proposed as a controlling variable, particularly for boreal regions 87 where the chilling requirement can be quickly reached and photoperiod acts as a safety limitation to 88 prevent a too early budburst and plant exposure to frost (Richardson et al., 2013; Way and 89 Montgomery, 2015). 90 Empirical models for budburst are generally derived from local meteorological data and observed 91 budburst timing, and mostly focus on deciduous species. The direct use of a model calibrated on a 92 site for gridded simulations over a region can be a source of errors, for example because of 93 altitudinal differences not resolved at a given grid horizontal resolution or because the whole range 94 of temperatures was not taken into account during the calibration (Olsson et al., 2013). This 95 problem calls for a multi-site calibration of budburst models with data drawn from a wide area, 96 typical of that used in the grid-based applications of ecosystem models. 97 Regarding the mechanisms involved, past studies did not highlight phenological differences 98 between deciduous and evergreen conifers for budburst. However, compared to deciduous species 99 that shed their leaves in autumn, evergreen conifers keep most of their needles over the year. Needle 100 lifespan can span from 2-3 years (e.g. Pinus Sylvestris) to more than 10 years (e.g. Picea abies) for 101 evergreen conifers. Needle senescence has been less studied (Estrella & Menzel, 2006; Delpierre et 102 al., 2009) than budburst. Some authors observed a peak of senescence during autumn (Sampson et 103 al., 2003; Kivimäenpää & Sutinen, 2007; Wang & Chen, 2012), however integrated over all needle 104 cohorts, needle senescence can be seen as a continuous process in evergreen species. Reich et al. 105 (2014) showed the role of needle longevity (related to nitrogen content) and the impacts of needle 106 senescence on carbon cycling in boreal forest, but very few studies investigated the regulation of 107 needle yellowing and turnover (renewal rate of needles), and none of them proposed any 108 mechanistic model for needle senescence (Muukkonen, 2005; Kayama et al., 2007). 4 109 While current global vegetation models roughly simulate phenological events for deciduous species, 110 these processes are still lacking for evergreen species for which the common approach is to 111 represent phenology implicitly through leaf biomass variations. 112 In this study, we tackle the following objectives: 113 - To calibrate empirical budburst models for temperate needleleaved species in order to 114 reproduce field observations collected in forest monitoring plots for a range of contrasting 115 climate conditions. 116 117 - To evaluate the accuracy of these models when used at low spatial resolution (0.25 and 0.5°) 118 typical of global models and to test the simulation results against independent remote 119 sensing observations. 120 121 - To implement a model for needle senescence on the basis of litterfall observations. 122 123 - To evaluate the potential impact of these model developments in a global vegetation model 124 (ORCHIDEE - Organizing Carbon and Hydrology In Dynamic Ecosystems; Krinner et al., 125 2005) on the representation of forest canopies and the associated carbon balance simulated 126 for temperate needleleaved forests in France. 127 5 128 2. Materials and methods 129 130 We used budburst observations and litterfall samples collected from 1997 to 2011 over the 51 forest 131 plots of the French RENECOFOR (REseau National de suivi à long terme des ECOsystèmes 132 FORestiers) network covering the six main coniferous species in France (number of plots in 133 parentheses): Pseudotsuga mensiezii (Douglas fir; 6), Picea abies (Norway spruce; 11), Pinus nigra 134 (Corsican pine; 2), Pinus pinaster (Maritime pine; 7), Pinus sylvestris (Scots pine; 14) and Abies 135 alba (Silver fir; 11). The parameters of eight different budburst models were calibrated against the 136 RENECOFOR site observations to get the best value for a performance metrics defined specifically 137 to account for both temporal and spatial variabilities. We selected the models that best described 138 budburst for the temperate evergreen needleleaf plant functional type (PFT) as a whole, but also 139 separately for each species. In addition, a new senescence model based on needle age was 140 developed and calibrated based on a literature review of needle survivorship observations. The new 141 phenology module (budburst+senescence) was then incorporated into the process-based model, 142 ORCHIDEE, and evaluated spatially against leaf area index estimated from remote sensing 143 observations. The flow chart of the model calibration and evaluation is given in Fig. 1. 144 145 Fig. 1: Flow chart of the model calibration and evaluation 6 146 147 2.1. Data from the RENECOFOR forest monitoring network 148 149 Visual observations of the timing of budburst were collected in 51 plots located in public forests 150 and stratified according to the major commercial tree species grown in France. Those plots are part 151 of the French national long-term forest monitoring network (RENECOFOR, part of the ICP – 152 International and Co-operative Program - Forest Level II program (http://icp-forests.net/page/level- 153 ii)), which covers a total of 102 permanent plots (51 coniferous) monitored since 1992 (Fig. 2). 154 155 Fig. 2: Distribution throughout France of the 51 RENECOFOR permanent plots dominated by 156 evergreen coniferous tree species. 157 158 Plots are in average 70y old (in 1994) with a range from 23 to 181y old. At each location, 36 trees 159 were chosen for phenological observations (Lebourgeois et al., 2010), and observations were 160 performed at least every week. For budburst, two different dates were measured. The first one is the 161 day of the year when 10% of the trees have open buds for at least 20% of the crown (BD1). The 162 second date corresponds to the day of the year when 90% of the trees have open buds for at least 7 163 20% of the crown (BD9). Observations were performed each year from 1997 to 2011. We 164 approximated a mean bud flushing date (50%; BD5) for budburst model parameterization defined 165 as: (1) 166 The final dataset contains 605 site-years of observations. Litterfall was collected seasonally from 167 1995 to 2007 at all RENECOFOR plots using litter traps. The detailed litter sampling methodology 168 is described in the Supplementary material (Appendix SA). 169 170 2.2. Budburst models 171 172 We tested two types of model for mono-cyclical budburst events, based on a temperature forcing 173 (e.g., degree-days) during spring, or based on “chilling-forcing”, i.e., with a chilling during winter 174 and a forcing period during spring. We used the model M1 as reference (called the static or null 175 model) in which budburst equals the median of the budburst dates observed across all sites (Table 176 1). A model with predictive value (spatially or temporally) should have better performance than 177 M1. 178 Species Mean BD5 (SD) Min BD5 Max BD5 n Abies alba 137.5 (11.5) 109.0 175.5 148 Picea abies 136.8 (11.5) 107.5 169.5 132 Pseudotsuga menziesii 130.3 (10.9) 102.5 151.0 72 Pinus nigra 127.4 (14.7) 106.0 163.5 20 Pinus pinaster 119.6 (18.2) 81.5 164.5 77 Pinus sylvestris 129.5 (13.3) 97.0 175.0 156 179 Table 1: Mean, minimum and maximum observed budburst date (DoY) in the RENECOFOR 180 network over the period 1997-2011. n = number of site-year observations 181 182 Eight models (M2-M9) were taken from the literature and tested (Cannell & Smith, 1983; 183 Hänninen, 1990; Kramer, 1994; Chuine et al., 1998; Chuine, 2000; Harrington et al., 2010; Vitasse 184 et al., 2011). The equations and parameters of each model are detailed in Appendix SB. All models 8 185 are representative of central hypotheses in budburst modeling and all of them except M9 were 186 already applied at regional scales. The models were selected to represent different concepts, but also 187 for their number of parameters to be small enough to allow inclusion in a global vegetation model. 188 The list of models tested has two spring warming forcing models and six chilling-forcing models 189 (Table 2). 190 In addition, we performed sensitivity tests to investigate the role of photoperiod, temperature 191 acclimation of parameters and the use of hourly or daily meteorological data (see Appendix SA for 192 the results of the sensitivity tests). 193 Number of Model Model name and reference parameters (without to, tc and tc,end) Reference model : median value of budburst M1 0 dates observed in the field forcing models Thermal time model: M2 – Kramer, 1994 Spring M2 1 M3 UniForc: M3 – Chuine et al., 1998 2 Sequential model: M4 – Hänninen, 1990; M4 5 Kramer, 1994 Chilling-forcing models Parallel model M5 – Cannell & Smith, 1983; M5 6 Kramer, 1994 M6 Unified model M6 – Chuine, 2000 5 M7 Alternating model M7 – Kramer, 1994; 2 Logistic alternating model M8 – Vitasse et al., M8 3 2011 M9 Harrington model M9 – Harrington et al. 2010 4 The daily forcing temperature calculated in Photoperiod models M2-M8 is weighted by the day length +1 following Blümel & Chmielewski (2012) 194 195 Table 2: Names and references of the eight budburst thermal models optimized in this study and the 196 corresponding number of optimized parameters when the starting date for chilling (t c), forcing (t o) 197 and the ending date for chilling (t cend) are fixed. See Appendix SB for a full description of the model 198 parameters and their equations. Note that critical temperature thresholds are not optimized in this 199 study but are estimated from the observed budburst dates. For example in M2 for which classical 200 optimization studies optimized both the base temperature (Tb) and the critical forcing threshold 201 (Fcrit) for budburst, here only Tb is optimized, Fcrit being the median Fcrit simulated at each 202 observed budburst date with M2. 9 203 204 2.3. Budburst model selection criterion 205 206 The purpose of this study being to calibrate a budburst model for regional simulations, we needed a 207 specific performance metric that characterizes the ability of this model to capture both spatial and 208 temporal budburst gradients. Most studies have used root mean square error (RMSE, Eqn 2) or 209 linear regression between simulated and observed budburst dates to select the best phenological 210 model (Vitasse et al., 2011; Olsson et al., 2013). We argue that two performance criteria are 211 desirable: a representation of both spatial and temporal extremes of budburst across a region, and a 212 representation with minimal systematic spatial or temporal bias. A single metric is insufficient to 213 account for these two criteria. To select the best set of parameters for each model, we thus propose a 214 combined metric accounting for three key aspects of the model: the ability to reproduce extreme 215 values, the average bias to observations and the effect of outliers. This new criterion is the 216 Euclidean distance to optima (DIST, Eqn 3) of six weighted different metrics across sites and years. 217 This new metrics maximizes model performances by catching both spatial and temporal variabilities 218 of budburst and by reducing the chance to converge to local minimum during calibration. This 219 combined metric has two components of the RMSE to limit outliers, the squared bias (BS, Eqn 4) 220 and the squared difference between standard deviation (SDSD, Eqn 5; Kobayashi & Salam, 2000) 221 normalized by the mean natural difference observed between sites SDSD nat and BSnat; the Spearman 222 rank correlation coefficient (R s), the linear regression slope (lms; Eqn 6) and its associated 223 coefficient of determination (lmR2) between observed and modeled budburst, which captures 224 extreme values; finally the average bias (AB, Eqn 7), and temporal and spatial biases between 225 modeled and observed budburst dates. For the purpose of large scale simulations, more weight was 226 given to capture bias and extremes than outliers (Eqn 3), with the sum of weighting factors equaling 227 1. With this metrics, the best performance is achieved when DIST=0. The different components of 228 DIST are: 10 (2) (3) (4) (5) (6) (7) 229 230 with x and y the observed and simulated budburst date, n the number of observations t, and SD the 231 standard deviation. 232 SDSDnat and BSnat were defined as the SDSD between sites in similar conditions and represent 233 variability in the observed data. Pairs of sites being within 0.5° (~55 km) maximum distance, with a 234 mean annual temperature difference less than 0.5 °C, are considered as having “similar conditions”. 235 Here, species similarity was not specified as a “similar condition”, in order to have enough 236 observations. Thus, we hypothesized that SDSDnat and BSnat are caused by species differences but 237 also non-resolved biotic or edaphic factors (local adaptation, age, soil effect, etc...), and define the 238 smallest value that an optimized model should approach when considering all conifers species. The 239 RMSEnat was estimated from nine sites with “similar conditions” to 7.6±3.5 days, SDSD nat to 240 3.16±5.09 days² and BSnat to 41.76±38.46 days². 241 The optimization of the parameters of models based on cross-site spatial variability only (DIST-S 242 averaging budburst years across site) resulted in a different best model to the one based on the 243 temporal variability only (DIST-T averaging budburst dates across all the sites each year) or 244 considering both the spatial and temporal variability (DIST-ST). Thus, we optimized below the 245 parameters of the eight budburst models described in Table 2 and we selected as “best predictive 11 246 model” the one corresponding to the minimum value of the DIST-ST metrics. 247 248 2.4. Budburst model optimization 249 250 The parameters of each model were optimized to minimize the value of DIST-ST against a subset of 251 the RENECOFOR observations (optimization dataset) consisting of 455 sites-years randomly 252 selected from the full dataset, with at least one observation per site. The remaining 150 observations 253 were used as cross-validation data. Note that models have different numbers of parameters (Table 254 2), i.e., different degrees of freedom. In addition to the cross-validation, the overall model accuracy 255 was assessed by coupling DIST-ST results to the Akaike’s information criterion corrected for 256 sample size (AICc) in order to select the best predictive (DIST-ST) and parsimonious (AICc) 257 model. Parsimonious models were selected by calculating the AICc difference (ΔAIC) between 258 AICc and the minimal AICc obtained among all models. Thus, the higher is ΔAIC, the less 259 parsimonious is a model. Models with ΔAIC higher than 10 were excluded (Burnham & Anderson, 260 2003). 261 We optimized models with a generalized simulated annealing algorithm (R package genSA; 262 Chuine, 2000; Xiang et al., 2013) considering parameters 1) per species, 2) grouping pines (Pinus 263 pinaster, Pinus sylvestris and Pinus nigra) versus firs and spruces (grouping Abies alba, Picea 264 abies and Pseudotsuga mensiezii, hereafter ‘fir’) into two groups according to Peaucelle et al. 265 (2016), and 3) pooling all conifer species together. Models were also fitted site by site to assess 266 possible emerging relationships between local parameter values and environmental conditions 267 (Appendix SA). For chilling-forcing models (M6-M9), an exponential relationship between chilling 268 units and forcing units is commonly used to estimate budburst. However, this exponential 269 relationship is potentially an artifact (Chuine, 2000) and is not observed for all species. We thus 270 decided to compare relationships fitted by exponential or by linear functions (Appendices SA & 271 SB). Different optimizations were performed by fixing or by optimizing the starting date of forcing 272 (to; 1st of January) and chilling (tc; 1st of November) in order to assess the models’ robustness with 12 273 fewer degrees of freedom. For M6 to M8 we assessed the impact of optimizing the end date for 274 chilling accumulation (tc,end), thus representing the fulfillment of the chilling requirement, otherwise 275 chilling is summed until budburst (Vitasse et al. 2011). 276 277 2.5. Meteorological dataset 278 279 We used the SAFRAN (Système d’Analyse Fournissant des Renseignements Adaptés à la 280 Nivologie) meteorological data (Vidal et al., 2010) for model optimization and for ORCHIDEE 281 site-scale simulations. This dataset produced by the Centre National de Recherches 282 Météorologiques (CNRM) provides hourly weather data over France at a spatial resolution of 0.07° 283 (8 km). At 0.07° resolution, each of the 51 coniferous forest sites is located in an independent grid 284 cell. Once the best set of parameters was retrieved for each model, we assessed the effect of the 285 spatial scale of climate data by applying the same models with SAFRAN data aggregated at 0.25° 286 (~28 km) and 0.5° (~55 km) resolution respectively. At 0.25° resolution, 49 grid cells contained at 287 least one site, and at 0.5°, 43 grid cells contained at least one site. All temperatures were corrected 288 in a simple way for local altitude following Eqn 8 (U.S. Standard Atmosphere, 1976, Olsson & 289 Jönsson, 2015): (8) 290 where Tobs (°C) is the mean observed temperature of the site, Tsaf (°C) the mean temperature of the 291 site from SAFRAN dataset, Asaf (km) the mean altitude of the SAFRAN cell and Aobs (km) the 292 altitude of the site. 293 294 2.6. Senescence model 295 296 The senescence model is not a stand-alone model as is the case for the model of budburst dates, but 297 rather a modification of the ORCHIDEE original phenology described in Krinner et al. (2005). 298 Compared to budburst models that are functions of environmental conditions, the senescence model 13 299 is based on needle age. The original version of ORCHIDEE includes two types of senescence for 300 needles. Firstly, a base rate of leaf mortality is applied each day (Krinner et al., 2005). It represents 301 the probability for needles to fall independently of needle age or meteorological conditions. 302 Secondly, senescence is triggered when needle age (calculated for four cohorts as in Section A1 of 303 Krinner et al., 2005) reaches a pre-defined longevity parameter for each PFT. As no phenological 304 process is explicitly defined in the default model, needle age is implicitly estimated from needle 305 biomass with the assumption that newly assimilated biomass through photosynthesis is used to 306 create new needles at the beginning of the year. 307 We did not find any suitable needle senescence model for coniferous species in the literature. We 308 thus decided to fit a senescence function against field observations of needle survival probability 309 from different studies (all studies and species are listed in Table SC1, Appendix SC). We retrieved 310 45 needle survivorship curves (determines the probability of needles to survive (0-1) over time 311 according to their age) from the literature and used these to calibrate a logistic function given by: (9) 312 Where S(t) is the survivorship probability of a needle, t the needle age (days), and µ (days-1) and λ 313 (days) parameters to be fitted on literature observations. 314 315 2.7. Modification of the ORCHIDEE model 316 The inclusion of an explicit phenology for evergreen conifers in ORCHIDEE (Krinner et al. 2005) 317 needed the modification of the original model. The needle maximum age parameter for evergreen 318 conifers, fixed at 910 days (average lifespan of Pinus needles) in the standard version of 319 ORCHIDEE was modified to depend on species. This maximum age can vary considerably, going 320 from 2 years in pine species to more than 10 years for Abies alba and Picea abies (Peaucelle et al., 321 2016). In ORCHIDEE, the Vcmax parameter (maximal rate of the RUBISCO carboxylation activity 322 in µmol m-2 s-1) increases with needle age, reaching a maximum value when the relative age of the 323 needle (the ratio of the needle age to its maximum) is 0.03 and then linearly decreasing to its 14 324 minimum value when the relative age reaches 0.5. This function describing the evolution of Vcmax 325 roughly represents species with short-lived needles such as pines (Niinemets, 2002). Observations 326 show a rapid decrease of Vcmax after only 1 year even in high needle-longevity species (Porté & 327 Loustau, 1998; Niinemets, 2002; Warren, 2006). We thus adapted this relationship prescribing Vcmax 328 to reach its maximum 3 months after formation and then starting to decrease linearly after one year, 329 until reaching 0.5 Vcmax at the maximal needle age of the species. In following simulations, we used 330 two different needle maximum age, 1275 and 2340 days for pines and spruces/firs species, 331 respectively (Peaucelle et al., 2016). 332 Given the senescence function in Eqn 9, we also modified the way carbon is distributed in the 333 crown by representing each cohort of needle (i.e. groups of needles developed the same year). All 334 the biomass gained during the current year is placed in the youngest needle cohort. Other cohorts do 335 not receive new biomass, but lose needles according to the senescence function from Eqn 9. To 336 exclude simulations where budburst never occurs, we imposed a maximum needle onset at day 182 337 (1 July), the latest observed budburst in our dataset. 338 339 2.8. Validation of the phenology models 340 341 Simulated budburst date was evaluated against observed data from the optimization and the cross- 342 validation dataset. We also investigated the ability of each model to predict spatial and temporal 343 variations in budburst across sites by comparing median modeled and observed budburst dates at 344 each site, and by looking at the interannual variability in the timing of budburst at each site. 345 Because of the high cross-site variability, an evaluation of the senescence model could not be 346 performed against litterfall observations, which also depend on non-modeled factors such as stand 347 health, stand age, density, species composition and management events. Thus, observed and 348 modeled litterfall were simply compared for information, and we validated the senescence model 349 through the indirect comparison with satellite-derived leaf area index (Wang et al., 2004). We 350 compared LAI simulated at a spatial resolution of 0.07° with ORCHIDEE forced by SAFRAN 15 351 against MODIS MCD15A3 LAI (1 km, 4-day frequency; Myneni et al., 2002). This local 352 comparison between modeled and satellite LAI was performed on a few grid cells where the 353 coverage of two representative coniferous species (Picea abies and Pinus sylvestris) exceeds 80% at 354 1 km, based on the European tree species map of Brus et al. (2012). The correlation coefficient 355 between modeled and estimated LAI was used to assess modeled LAI seasonality. 356 We also performed grid-based simulations for coniferous forests in Europe at a 0.25° spatial 357 resolution to compare simulated and satellite LAI at a larger scale over the period 2000-2007. Initial 358 conditions and forcing data used for simulations are detailed in the Supplementary Material 359 (Appendix SA). 360 Finally, we ran the model over France at a resolution of 0.07° for the period 1970-2100 to assess 361 patterns in budburst timing in future decades. The climate forcing (daily data) was from the A2 362 scenario of ARPEGE v4 model downscaled and bias corrected by Pagé et al. (2008). In these 363 simulations of the future, land cover was imposed from the IGBP map (Loveland & Belward, 1997) 364 and soil depth and texture (used to derive wilting points and field capacities and thence to give plant 365 water stress) from the FAO dataset (Vetter et al., 2008). 366 367 3. Results 368 3.1. RENECOFOR budburst and litterfall observations 369 On average 50% of buds flushed for the sites of Fig. 2 within a range of 14 days around the 12 May 370 (day 132; Table 1). The earliest budburst was recorded on 22 March 2007 for Pinus pinaster in “Les 371 Landes” forest, in southwest France. The latest budburst dates were observed in mountainous 372 regions for Abies alba and Pinus sylvestris, 24 June 1999 and 24 June 2008, respectively. On 373 average, coniferous stands dropped 2336.5 kg ha−1 yr−1 (dry matter) of their needles as litter from 374 1997 to 2007, which represents 65% of the total annual litterfall for all compartments on average 375 (all stands, only considering the dominant species; 49% considering secondary species) and more 376 than 70% for Abies alba, Picea abies and Pinus pinaster (Table 3). We could see large differences 16 377 in the mass of needles lost per year among species. Species with minimum losses were Abies alba 378 (mean over all sites was 1892.0 kg ha−1 yr−1) and Pinus sylvestris (1859.5 kg ha−1 yr−1) whereas the 379 maximum was observed for Pinus pinaster (3175.4 kg ha−1 yr−1). 380 Species Leaves CVL Prop. Total CVT n kg ha-1 yr-1 % % kg ha-1 yr-1 % A.alb 1892.2 37.2 58.0 3263.2 33.0 145 P.abi 2692.6 30.6 72.8 3698.6 30.0 121 P.men 2036.7 31.8 84.5 2411.3 30.4 66 P.neg 2447.2 24.3 51.6 4740.4 18.4 21 P.pin 3175.4 30.3 76.1 4172.4 28.0 81 P.syl 1859.5 34.6 49.6 3750.4 30.0 171 All 2336.5 20.4 49.2 3580.4 21.4 605 381 382 Table 3: Litterfall mass (kg dry matter ha -1 yr-1) measured in the RENECOFOR network over the 383 period 1997-2007 for each compartment. CV corresponds to the coefficient of variation for each 384 compartment (leaves, and total) and Prop. corresponds to the proportion of the compartment 385 compared to the total litterfall. P.men=Pseudotsuga menziesii, P.abi=Picea abies, P.neg=Pinus 386 negra, P.pin=Pinus pinaster, P.syl=Pinus sylvestris, A.alb=Abies alba. 387 388 3.2. Budburst models comparison and selection of a best model 389 Best models (parsimonious and predictive) retained for each species are listed in Table 4 and the 390 corresponding model parameters are given in Table SC2 (Appendix SC). Figure SC1 shows the 391 DIST-ST evaluation metric after parameter optimization for each model forced by daily 392 temperatures. For all species together, the best model (DIST-ST criterion) is the simple spring 393 forcing model M3 (DIST-ST=0.25, RMSE=12.5 days, Fig. 3) with a starting date fixed to 1 394 January. The most parsimonious is the chilling-forcing model M7 (DIST-ST=0.39). With M3, both 395 spatial (DIST-S=0.32, RMSE=10.7 days) and temporal (DIST-T=0.17, RMSE=7.3 days) variability 396 was well reproduced for the validation dataset. The DIST-ST values obtained with models M2, M6, 397 M7, M8 and M9 are close to this best model, with DIST-ST of 0.35, 0.26, 0.39, 0.4 and 0.29, 398 respectively, but only models M7, M8 and M9 are considered parsimonious according to ΔAIC. 399 The sequential (M4) and parallel models (M5) could not reproduce observations properly (DIST- 400 ST=1.58 to 2.56, RMSE=19.0 to 28.0 days). 401 17 402 Species Model DIST-ST All M3 0.25 Firs M7 0.47 Pines M7 0.46 Abies alba M3 0.56 Picea abies M7 1.25 Pseudotsuga mensiezii M9 0.32 Pinus pinaster M3 1.15 Pinus sylvestris M7 0.45 403 404 Table 4: Best models retained for each species according to both the predictive power considering 405 spatial and temporal variability (DIST-ST) and the parsimony (with the lowest number of 406 parameters). Models were optimized against daily temperatures. See Appendix SB for a detailed 407 description of each model. 408 409 For groups of species, in the case of firs, the best model is the chilling-forcing model M7 (DIST- 410 ST=0.47, RMSE=11.9 days). The temporal variability (DIST-T=0.45) representation was 411 equivalent to the spatial variability for M7 for firs (DIST-S=0.48). Models M2, M3, M6, M8 and 412 M9 have performances close to M7 but only M3, M7 and M8 have a ΔAIC<10. For pines, all 413 optimized models produced better DIST-ST values than the null model M1, but higher RMSE. The 414 best model for pines is again the chilling-forcing model M7 (DIST-ST=0.46, RMSE=16.3 days), 415 while the best parsimonious is M2 (DIST-ST=0.63). Model M7 better represented the spatial 416 pattern of observed budburst (DIST-S=0.20, RMSE=10.2 days) than for temporal variability 417 (DIST_T=0.62, RMSE=13.5 days). 418 For individual species, the most parsimonious model is also M7. In the case of Abies alba the best 419 DIST-ST is obtained with M6 (DIST-ST=0.48) while the best parsimonious model is the model M3 420 (DIST-ST=0.56). For Picea abies, none of the models reproduced accurately the observations, M7 421 having the best score and parsimony with DIST-ST=1.25 (RMSE=21.0 days). For Pseudotsuga 422 mensiezii the best model was M9 (DIST-ST=0.32, RMSE=10.7 days), which concords with the 423 results of Harrington et al. (2010) for this species. However, for both Pseudotsuga mensiezii and 18 424 Pinus pinaster the most parsimonious model led to high DIST-ST values, with 2.19 and 2.64 425 respectively. For Pinus pinaster, the best DIST-ST was obtained with M8 (DIST-ST=0.33, 426 RMSE=14.0 days), but was not selected as parsimonious. The best model was thus M3 with DIST- 427 ST=1.15. Finally, the best DIST-ST for Pinus sylvestris was obtained with M4 (DIST-ST=0.38, 428 RMSE=12.0 days) while the most parsimonious model is again M7. 429 Lower performances of the optimization at the species level compared to groups of species can be 430 explained by the smaller training dataset available for parameter calibration. We argue that models 431 calibrated with all species or groups of species should thus be more robust and more suitable for 432 large scale simulations. 433 434 19 435 Fig. 3: Results for each site and each year for the best model with all species together (model M3). 436 Representation of the validation dataset (150 obs.). (a) and (b) correspond to mean dates by site 437 (Error bars correspond to the inter-annual variability), (c) and (d) correspond to mean dates by year 438 (all sites together, error bars correspond to the inter-site variability). The y-axis of (b) represents 439 each site for which the dominant species is represented by symbols listed in (a). For (b) and (d), 440 black dots correspond to mean observations, red diamonds correspond to mean modeled budburst 441 dates. 442 443 Results with or without optimizing the starting date for temperature accumulation (to, tc) have 444 similar model performance (Fig. SC1). However, we can see different performances of the same 445 model depending on the species. For example, the optimization of model M6 with fixed tc led to 446 better DIST-ST for Picea abies and worse DIST-ST for Pseudotsuga mensiezii compared to M6 447 with optimized tc. Thus, we preferentially selected models with a fixed starting date for large scale 448 simulations when optimization results were equivalent. The same conclusion applies to 449 optimizations with varying tcend (end of chilling accumulation). 450 451 For all species, pines and firs, we found quiet similar performances for both forcing and chilling- 452 forcing models. Note that chilling-forcing models may still be more physiologically realistic for 453 future predictions where warmer winters may exacerbate the effects of incomplete fulfillment of 454 chilling, or for applications in cold regions where chilling should be more important than in France 455 and western Europe. We selected M3 (forcing) and M7 (chilling-forcing) for inclusion in 456 ORCHIDEE. 457 458 3.3. Model performances from site-scale to grid-based resolution 459 460 We checked for model robustness at lower spatial resolutions, representative of typical forcing data 461 for global vegetation models. At 0.25° resolution, most of DIST-ST values were higher than at a 462 resolution of 0.07° and the best models differed. For all species, the best DIST-ST increased from 463 0.25 (M3-to) to 0.43 (M8+to+tc). But the best parsimonious model was still M7 with a DIST-ST 464 value of 0.49. Some models were no longer able to work correctly and DIST-ST values diverged 20 465 (DIST-ST= 132.2 for Pinus sylvestris-M8). At 0.5° the most parsimonious model remained M7 but 466 DIST-ST increased from 0.49 to 0.77. If we compare the evolution for the best models M3 and M7 467 at 0.07° and 0.5°, we can observe that the result is more degraded for M3 (DIST-ST increased from 468 0.25 to 1.77) than for M7 (DIST-ST from 0.39 to 0.77), but is still a much better performance than 469 the null model M1 with DIST-ST=9.17. 470 At a lower spatial resolution, some models could not be used. This was the case for models with a 471 fixed threshold for chilling accumulation (M4, M5 and M6). By averaging temperatures, the critical 472 threshold for chilling accumulation of these models was never reached at some sites and 473 consequently forcing temperatures could never accumulate. In the rare cases when the model 474 succeeded in estimating a budburst date, we could see that the performance was lower than the null 475 model M1 (DIST>10). The implications of using models derived at the site scale for low-resolution 476 prediction are further addressed in the discussion section. 477 478 3.4. Senescence model parameters 479 480 The minimum and maximum needle lifespan in Eq. 9 retrieved from literature studies were 4 and 15 481 years, respectively. Independently of environmental factors, species or tree health, needle 482 survivorship follows almost the same pattern in each study: the needle biomass turnover is 483 relatively low during the first years of the needle life and then rapidly increases over time. The 484 logistic relationship we fitted on those data was strongly correlated to the needle lifespan of the tree 485 (Eqn 10, Eqn 11) with R2=0.93 and 0.94 for parameters µ and λ, respectively (Fig. SC3, SC4; 486 Appendix SC). (10) (11) 487 with Needleage the maximal needle age. 488 The strong correlation between senescence parameters and the maximal needle age thus allows the 489 use of one unique relationship for all species. 21 490 491 3.5. Comparison against satellite data 492 493 We first compared LAI simulated at site scale for two representative species (Picea abies (Fig. 4a) 494 and Pinus sylvestris (Fig. 4b)) with MODIS satellite observations over the MODIS pixel of 1 km 495 containing each site. Simulations were performed here with budburst results obtained for the pines 496 and fir/spruces groups. For both models M3 (best forcing model) and M7 (best chilling-forcing 497 model), results are equivalent. The amplitude of the LAI cycle with the improved phenology 498 (1.4±0.1 for Pinus sylvestris; 1.2±0.1 for Picea abies - model M7) was closer to the amplitude of 499 observed LAI (1.7±0.6; 1.9±0.7) than with the standard version of the model which does not have 500 an explicit needle budburst equation (0.4±0.1; 0.4±0.1, respectively). For the 15 sites compared 501 with MODIS LAI for each species, the mean correlation coefficient between modeled and estimated 502 LAI improved from 0.45 ±0.2 to 0.77±0.1 for Picea abies and from 0.47±0.2 to 0.71±0.1 for Pinus 503 sylvestris. 22 504 505 Fig. 4: Comparison between observed (MODIS) and simulated (model M7) LAI dynamics with 506 ORCHIDEE for a) a Picea abies stand (lat=50.16º, long=5.46º) and b) a Pinus sylvestris stand 507 (lat=49.25º, long=8.06º). All data are centered on the average observed (or simulated) LAI value 508 (2003-2008). Gray and black lines represent the observed MODIS LAI and the moving average 509 over a 30-days window, respectively. The orange dotted line represents the simulated LAI with the 510 standard version of ORCHIDEE without phenological processes. The green dashed line represents 511 the simulated LAI with ORCHIDEE including budburst and senescence processes. 512 513 Figure 5 shows the correlation coefficient between satellite-observed and simulated LAI at the 514 European scale with the best forcing model (M3 optimized with all species) and the senescence 515 model. The spatial correlation between modeled and satellite LAI improved by 0.24 (from 23 516 R=0.48±0.3 to 0.72±0.2) over Europe even if the calibration was performed only over France. We 517 could observe the same improvement with the best parsimonious chilling-forcing model M7 518 (R=0.69±2). Moreover, we observed that the modeled budburst with M3 reached the imposed limit 519 of budburst date in very high latitudes and altitudes. On the contrary, the chilling-forcing model M7 520 was better able to predict the LAI seasonality at high altitudes and high latitude, because it never 521 reached the imposed budburst date in these areas. 522 523 Fig. 5: Spatial representation of the correlation coefficient between MODIS estimated LAI and 524 ORCHIDEE simulated LAI at a 0.25° spatial resolution averaged over the period 2000-2007: a) in 525 the standard configuration without explicit phenological processes, or b) with the budburst model 526 M3 and the senescence model. Map c) gives the absolute difference (b - a). Only the dominant 527 species used in this study are represented (Abies alba, Picea abies, Pinus pinaster, Pinus sylvestris, 528 Pseudotsuga mensiezii). The correlation is calculated for pixels with a minimal coniferous fraction 529 cover of 20%. 530 531 532 533 3.6. Budburst evolution from 1970 to 2100 534 535 Figure 6 shows the simulated evolution of the mean budburst date over France from 1970 to 2100 536 with the best model M3 and the most parsimonious forcing-chilling model M7. For the two models, 537 we simulated an earlier needle unfolding over time (mean slope=-0.126±0.01 days yr−1, R2 538 =0.57±0.06). Compared to 1970-2000, projections indicated an advancement of the budburst date of 539 coniferous species of 10.3±2.8 and 12.3±4.1 days in average over the period 2060-2100 with M3 540 and M7, respectively. However, model M7 exhibited higher variability in the prediction of needle 541 unfolding over time and an earlier budburst date on average compared to M3 (4.8±0.8 days) over 24 542 the period 1970-2100. 543 544 Fig. 6 : Temporal evolution of the mean budburst dates (DOY) simulated over France from 1970 to 545 2100 for models M3 and M7 (all species together). For each line, the corresponding colored area 546 indicates one s.d. either side of the mean. 547 548 549 3.7. Impact of the new conifer phenology model on GPP 550 We assessed the impact of the new phenology model on the simulated GPP. We found in our 551 simulations an increase of GPP by 15±1% when compared to the standard version of ORCHIDEE. 552 For example, GPP increased from an average 5.5±0.2 g C m−2 d−1 to 6.3±0.2 g C m−2 d−1 for spruce 553 stands over France during the historical period (2000-2010). GPP increase was mainly induced by a 554 modification of the canopy composition. Due to the production of new needles each year and 555 senescence removing old needles, needle cohorts forming the canopy are younger (see Materials 556 and methods) in the new model (845.3±55.5 days for spruces versus 1740±18.8 days in the standard 557 version of ORCHIDEE). New needles having higher photosynthetic rates, this results into a higher 558 simulated GPP. Because old needles are already present at the beginning of the growing season (i.e. 559 photosynthetic activity starts as soon as climate conditions are favorable) and because new needles 25 560 reach their maximal activity only 3 months after unfolding, simulated GPP was not sensitive to 561 changes in budburst date. By imposing the needle onset (from day 90 to 160), we observed a mean 562 GPP difference of 0.002 g C m−2 d−1, which sums up to only 0.8±0.1 g C m−2 over the whole year 563 for each day difference in budburst. 564 With the new senescence model instead of the standard parameterization, we found an increase of 565 the needle turnover from 0.16±0.002 to 0.43±0.025 g C m−2 d−1 for spruce. Consequently, litterfall 566 for Picea abies stands without explicit senescence has an average of 1155.0±12.6 kg ha−1 yr−1 over 567 France while we simulated losses of 3146.5±182.7 kg ha−1 yr−1 with explicit senescence, which is 568 closer to observations for this species (2692.6 kg ha−1 yr−1); this is a relative difference of -57% and 569 +17%, respectively, between observed and simulated litterfalls. 570 571 4. Discussion 572 573 4.1. Uncertainties in model validation 574 575 Phenology is a central function in stands of conifers. In this study we optimized different 576 phenological models against in situ budburst observations for six coniferous species. We showed 577 that most models reproduced budburst dates with a precision of ±12 days across France. However 578 we highlighted the need to consider both spatial and temporal variability when calibrating a 579 phenological model. Figure SC5 (Appendix SC) illustrates the differences in results for the two 580 spring forcing models M2 and M3 considering either the simple RMSE or our more comprehensive 581 DIST metrics to select the best model after calibration. With the model M2 the best optimization 582 based on RMSE does not allow to simulate budburst occurring before day 130 (Fig. SC5a, b), while 583 with the model M3 (Fig. SC5c, d) the best optimization based on RMSE resulted in a high 584 variability of simulated budburst. Even when considering both temporal and spatial variability of 585 budburst with the DIST metrics, our results highlighted that, while some models managed to 586 reproduce the observed budbursts (for instance M7 for firs, section 3.2), the same model calibrated 26 587 on different dataset generally reproduce better the spatial variability than the temporal variability of 588 budburst. The DIST metrics proposed in this study is a first attempt to take into account both 589 temporal and spatial variability of budburst. The combination of multiple metrics is promising to 590 improve model calibration but more investigation is needed to improve this metrics, like the number 591 of component we need to consider and their weight. 592 This implementation of better phenology models for evergreen conifers in the global model 593 ORCHIDEE shows an improved ability to reproduce the seasonal LAI dynamics observed from 594 MODIS both at the site and regional scale, despite representativeness differences between one site 595 and a MODIS pixel, and the contribution of understory vegetation to the MODIS signals. Previous 596 studies have shown a good correspondence between field-measured LAI and MODIS products 597 (Jensen et al., 2011; Rautiainen et al., 2012). However, needle senescence is locally dependent of 598 stand health, age and disturbances. We indirectly validated our senescence model against MODIS 599 LAI data, despite uncertainties of this product. Comparisons with litterfall observations further 600 allowed to show that our senescence model represent an improvement. By comparing site 601 measurements of carbon fluxes and MODIS products, Verma et al. (2014) highlighted the 602 uncertainties linked to the heterogeneity of the vegetation at a larger spatial resolution. In this study 603 we visually checked for canopy openness at the site scale, however LAI from MODIS products 604 integrates contributions from both the dominant species and the understory vegetation, which will 605 result in a bias when comparing PFT specific simulations and observed amplitude and dynamics of 606 the LAI (Wang et al., 2004; Jensen et al., 2011; Rautiainen et al., 2011; Rautiainen & Heiskanen, 607 2013). Moreover, LAI is related to stand age and health (Pokornỳ et al., 2013), which were not 608 taken into account in our study. 609 The senescence model proposed in this study defines a continuous process over the needle lifespan. 610 Moreover, we fixed the length of the period of needle accretion in ORCHIDEE, which can result in 611 biases in LAI at the beginning of the growing season. We argue that an optimization of the 612 senescence model parameters against carbon fluxes and remote sensing observations could further 27 613 improve the senescence model. 614 615 4.2. Relevance of site-calibrated models for gridded simulations 616 617 All tested models in this study could be optimized with good accuracy at the site scale 618 (RMSE=6.0±3.4 days for M3, 4.3±2.7 days for M7) or across sites (RMSE=12.5 days with all 619 sites). We also show that forcing models performed equally compared to chilling-forcing models at 620 the regional scale. Olsson & Jönsson (2014) indicated that simple models with few parameters are 621 more accurate over larger regions in general. In our study, most phenological models were more 622 efficient than the null model (M1), i.e., setting a fixed date equal to the median observed one, even 623 at a degraded spatial resolution. However, some models were not able to reproduce budburst events 624 when calibrated using large-scale temperature forcing data, especially models M4 and M5. 625 Modelers should be cautious when applying empirical models fitted at site scale for large-scale 626 predictions for two aspects: 1) If a model developed for site scale studies does not work at larger 627 resolutions, it means either that the model is not generic enough: the sites were too specific or some 628 processes are missing, like adaptation or acclimation for example; 2) The spatial aggregation of 629 temperatures can smooth and modify the response of chilling and forcing. Modelers have to check 630 that models calibrated on sites, and thus dependent of site conditions, are still able to reproduce 631 average responses of budburst globally and not only in limited environmental conditions. Thus the 632 model validity should be assessed at different spatial and temporal scales. The metrics developed in 633 our study lowers the weight of outliers in the calibration and thus limits this effect. 634 Here, the best forcing model M3 has more degraded performance than the best chilling-forcing 635 model M7 when applied at 0.25° and 0.5°, even if M7 exhibits a lower performance than M3 when 636 calibrated at 0.07°. This result suggests that the model M7 may be more suitable than M3 to be used 637 in a global model such as ORCHIDEE. 638 We also found that the best model calibrated with in situ observations is not necessarily the best 639 model to reproduce the seasonality of the satellite LAI cycle. Our results suggest that looking at the 28 640 whole phenological cycle, in addition to in situ observations, could be a better way of 641 discriminating between budburst models intended to be used in global vegetation models. In a 642 recent study, Gamon et al. (2016) demonstrated a method to track photosynthetic phenology in 643 evergreen conifers using a remotely sensed reflectance chlorophyll/carotenoid index (CCI). This 644 new indicator could allow an accurate calibration of phenological processes for conifers. However, 645 because of the land cover heterogeneity over Europe, a calibration of phenological models based 646 only on satellite observations would also benefit in the future from a very high temporal and spatial 647 resolution dataset (Delegido et al., 2011; Verrelst et al., 2012; Klosterman et al. 2014). 648 649 4.3. Impact of phenology for large-scale simulations of GPP and LAI 650 651 The new phenological processes incorporated into ORCHIDEE led to a better representation of the 652 seasonal cycle of LAI both at site and regional scale (Figs 5, 6). Simulated LAI was however not 653 sensitive to the timing of needle onset, while the production of new needles and the use of an 654 explicit representation of needle senescence had a strong impact on simulated variables. In our 655 simulations, all the new needle biomass was allocated to younger needles with high photosynthetic 656 efficiency, thus leading to simulated GPP being higher than in the original model. The explicit 657 representation of the senescence also led to a higher needle turnover, and litterfall. A higher 658 litterfall rate will strongly impact soil carbon pools and heterotrophic respiration. Here, the 659 estimated amount of needle-fall was in the range of values observed in French stands, while the 660 standard version of ORCHIDEE underestimated the amount of litter. However, the lack of 661 information about the living needle biomass at each site did not allow an accurate comparison of 662 simulated and observed litterfall. Needle-fall is closely related to stand age, stand health, climate 663 and disturbances (Balster & Marshall, 2000; Choi et al., 2006; Reich et al., 2014). The relation 664 defined in this study does not take into account all the factors influencing needle-fall and more 665 investigation is needed. 666 The strong correlation between the senescence parameters and the maximal age of needles (Eq. 7, 29 667 8) makes our proposed senescence model relevant for a large variety of evergreen coniferous 668 species and may be generalized to other evergreen species in other biomes. With the generalization 669 of models with varying traits (Pavlick et al., 2012, Verheijen et al., 2015), we argue that our 670 senescence model could be easily implemented with trade-offs concerning the maximal age of 671 needles, as for example relationships between the needle lifespan and mean annual temperatures, 672 recently implemented in the CABLE model by Reich et al. (2014). As for budburst, we can expect a 673 change in the senescence rate with global changes. The inclusions of trade-offs between needle 674 longevity and climate in the ORCHIDEE model as it was done by Reich et al. (2014) will be a first 675 step in understanding the impact of such changes on the carbon balance of forest ecosystems. 676 For evergreen species, a significant amount of leaves/needles is already present at the beginning of 677 spring. The presence of old needles thus allows the recovery of carbon and water fluxes when 678 temperatures become favorable — this explains the low sensitivity of simulated fluxes to budburst 679 date. However, even if the needle onset had few impacts in this study, the implementation of an 680 explicit budburst model was shown to play a key role in other mechanisms such as ozone sensitivity 681 of needles (Watanabe et al. 2010; Verbeke, 2015), frost risk (Hänninen, 2006; Man et al., 2015), 682 biogenic emissions (Richardson et al., 2013) or vegetation dynamics (Lu et al., 2016). The 683 mechanisms presented in this study could be extended to improve the representation of other 684 evergreen species in global models. 685 5. Conclusion 686 Phenology plays a central role in bio-geochemical cycles in conifers stands. In this study we 687 optimized different phenological models against budburst observations of six conifers species. We 688 show that all models managed to reproduce needle emergence at the site scale with good accuracy 689 (±5d). At the national scale, most models reproduce budburst dates with a precision of 12 days. 690 However, we highlight the need to consider both spatial and temporal variability when calibrating 691 phenological models. 692 Different budburst models performed equally independently of their complexity and the process 30 693 they include (forcing, chilling, photoperiod) and models calibrated at a fine spatial resolution were 694 not able to predict budburst dates when applied at coarse resolutions typical of the grid of global 695 models. This suggests that common models developed for site scale experiments might be 696 inadequate for large scale simulations. This first attempt in implementing an explicit phenological 697 model for evergreen conifer PFTs for large scale simulations managed to reproduce the observed 698 LAI dynamics both at the site and regional scale. A sensitivity analysis highlighted that the new 699 phenology module has a significant impact on the simulated carbon fluxes. We showed that needle 700 onset will be ±11 days earlier in 2060-2100 compared to 1970-2000 and more analysis are needed 701 to quantify the effect of evergreen conifers phenology on the projected carbon budget. The findings 702 will help future research to better improve current and future predictions of carbon, water, nutrient 703 and heat cycles using ecosystem model. 704 705 Acknowledgment 706 This work was granted access to the HPC resources of TGCC under the allocation 2015-6328 made 707 by GENCI. The authors would like to acknowledge the financial support from the European 708 Research Council Synergy grant ERC-SyG-2013-610028 IMBALANCE-P. Many thanks go to all 709 the people who contributed to setting up and maintaining the RENECOFOR network, to collecting 710 data and to ensuring its quality. RENECOFOR is part of the ICP Forests monitoring program and 711 has been continuously supported by French public funds (Office National des Forêts, Ministry of 712 Agriculture, Ministry of Ecology, ADEME) and by the European Union from 1991 until 2006. 713 31 714 References 715 Atmosphere, U.S., 1976. US standard atmosphere. National Oceanic and Atmospheric 716 Administration. 717 Baldocchi, D., Falge, E., Gu, L., Olson, R., Hollinger, D., Running, S., Anthoni, P., Bernhofer, C., 718 Davis, K., Evans, R., Fuentes, J., Goldstein, A., Katul, G., Law, B., Lee, X., Malhi, Y., Meyers, T., 719 Munger, W., Oechel, W., Paw, K.T., Pilegaard, K., Schmid, H.P., Valentini, R., Verma, S., Vesala, 720 T., Wilson, K., Wofsy, S., 2001. FLUXNET: A New Tool to Study the Temporal and Spatial 721 Variability of Ecosystem–Scale Carbon Dioxide, Water Vapor, and Energy Flux Densities. Bull. 722 Am. Meteorol. Soc. 82, 2415–2434. doi:10.1175/1520-0477 723 Balster, N.J., Marshall, J.D., 2000. Decreased Needle Longevity of Fertilized Douglas-Fir and 724 Grand Fir in the Northern Rockies. Tree Physiol. 20, 1191–1197. doi:10.1093/treephys/20.17.1191 725 Blümel, K., Chmielewski, F.-M., 2012. Shortcomings of classical phenological forcing models and 726 a way to overcome them. Agric. For. Meteorol. 164, 10–19. doi:10.1016/j.agrformet.2012.05.001 727 Brus, D.J., Hengeveld, G.M., Walvoort, D.J.J., Goedhart, P.W., Heidema, A.H., Nabuurs, G.J., 728 Gunia, K., 2012. Statistical mapping of tree species over Europe. Eur. J. For. Res. 131, 145–157. 729 doi:10.1007/s10342-011-0513-5 730 Burnham, K.P., Anderson, D.R., 2003. Model Selection and Multimodel Inference: A Practical 731 Information-Theoretic Approach. Springer Science & Business Media. 732 Burr, K.E., Tinus, R.W., Wallner, S.J., King, R.M., 1989. Relationships among cold hardiness, root 733 growth potential and bud dormancy in three conifers. Tree Physiol. 734 Cannell, M.G.R., Smith, R.I., 1983. Thermal Time, Chill Days and Prediction of Budburst in Picea 735 sitchensis. J. Appl. Ecol. 20, 951–963. doi:10.2307/2403139 736 Chen, M., Melaas, E.K., Gray, J.M., Friedl, M.A., Richardson, A.D., 2016. A new seasonal- 737 deciduous spring phenology submodel in the Community Land Model 4.5: impacts on carbon and 738 water cycling under future climate scenarios. Glob. Change Biol. 22, 3675–3688. 739 doi:10.1111/gcb.13326 740 Choi, D.S., Kayama, M., Jin, H.O., Lee, C.H., Izuta, T., Koike, T., 2006. Growth and 741 photosynthetic responses of two pine species (Pinus koraiensis and Pinus rigida) in a polluted 742 industrial region in Korea. Environ. Pollut. 139, 421–432. doi:10.1016/j.envpol.2005.06.006 743 Chuine, I., 2000. A unified model for budburst of trees. J. Theor. Biol. 207, 337–347. 744 doi:10.1006/jtbi.2000.2178 745 Chuine, I., Beaubien, E.G., 2001. Phenology is a major determinant of tree species range. Ecol. 746 Lett. 4, 500–510. doi:10.1046/j.1461-0248.2001.00261.x 747 Chuine, I., Cour, P., Rousseau, D.D., 1998. Fitting models predicting dates of flowering of 748 temperate-zone trees using simulated annealing. Plant Cell Environ. 21, 455–466. 749 doi:10.1046/j.1365-3040.1998.00299.x 750 Cooke, J.E.K., Eriksson, M.E., Junttila, O., 2012. The dynamic nature of bud dormancy in trees: 751 environmental control and molecular mechanisms. Plant Cell Environ. 35, 1707–1728. 752 doi:10.1111/j.1365-3040.2012.02552.x 32 753 Delegido, J., Verrelst, J., Alonso, L., Moreno, J., 2011. Evaluation of Sentinel-2 Red-Edge Bands 754 for Empirical Estimation of Green LAI and Chlorophyll Content. Sensors 11, 7063–7081. 755 doi:10.3390/s110707063 756 Delpierre, N., Dufrêne, E., Soudani, K., Ulrich, E., Cecchini, S., Boé, J., François, C., 2009. 757 Modelling interannual and spatial variability of leaf senescence for three deciduous tree species in 758 France. Agric. For. Meteorol. 149, 938–948. doi:10.1016/j.agrformet.2008.11.014 759 Estrella, N., Menzel, A., 2006. Responses of leaf colouring in four deciduous tree species to climate 760 and weather in Germany. Clim. Res. 32, 253. 761 Gamon, J.A., Huemmrich, K.F., Wong, C.Y.S., Ensminger, I., Garrity, S., Hollinger, D.Y., 762 Noormets, A., Peñuelas, J., 2016. A remotely sensed pigment index reveals photosynthetic 763 phenology in evergreen conifers. Proc. Natl. Acad. Sci. 113, 13087–13092. 764 doi:10.1073/pnas.1606162113 765 Gunderson, C.A., Edwards, N.T., Walker, A.V., O’Hara, K.H., Campion, C.M., Hanson, P.J., 2012. 766 Forest phenology and a warmer climate – growing season extension in relation to climatic 767 provenance. Glob. Change Biol. 18, 2008–2025. doi:10.1111/j.1365-2486.2011.02632.x 768 Hänninen, H., 2006. Climate warming and the risk of frost damage to boreal forest trees: 769 identification of critical ecophysiological traits. Tree Physiol. 26, 889–898. 770 Hänninen, H., 1990. Modelling bud dormancy release in trees from cool and temperate regions. 771 Acta For. Fenn. 213, 1–47. 772 Hänninen, H., Slaney, M., Linder, S., 2007. Dormancy release of Norway spruce under climatic 773 warming: testing ecophysiological models of bud burst with a whole-tree chamber experiment. Tree 774 Physiol. 27, 291–300. doi:10.1093/treephys/27.2.291 775 Harrington, C.A., Gould, P.J., St.Clair, J.B., 2010. Modeling the effects of winter environment on 776 dormancy release of Douglas-fir. For. Ecol. Manag., doi:10.1016/j.foreco.2009.06.018 777 Jensen, J.L.R., Humes, K.S., Hudak, A.T., Vierling, L.A., Delmelle, E., 2011. Evaluation of the 778 MODIS LAI product using independent lidar-derived LAI: A case study in mixed conifer forest. 779 Remote Sens. Environ. 115, 3625–3639. doi:10.1016/j.rse.2011.08.023 780 Kayama, M., Kitaoka, S., Wang, W., Choi, D., Koike, T., 2007. Needle longevity, photosynthetic 781 rate and nitrogen concentration of eight spruce taxa planted in northern Japan. Tree Physiol. 27, 782 1585–1593. 783 Kivimäenpää, M., Sutinen, S., 2007. Microscopic structure of Scots pine (Pinus sylvestris (L.)) 784 needles during ageing and autumnal senescence. Trees 21, 645–659. doi:10.1007/s00468-007-0157- 785 8 786 Klosterman, S.T., Hufkens, K., Gray, J.M., Melaas, E., Sonnentag, O., Lavine, I., Mitchell, L., 787 Norman, R., Friedl, M.A., Richardson, A.D., 2014. Evaluating remote sensing of deciduous forest 788 phenology at multiple spatial scales using PhenoCam imagery. Biogeosciences 11, 4305–4320. 789 doi:10.5194/bg-11-4305-2014 790 Kobayashi, K., Salam, M.U., 2000. Comparing simulated and measured values using mean squared 791 deviation and its components. Agron. J. 92, 345–352. 792 Kramer, K., 1994. Selecting a Model to Predict the Onset of Growth of Fagus sylvatica. J. Appl. 793 Ecol. 31, 172–181. doi:10.2307/2404609 33 794 Krinner, G., Viovy, N., Noblet-Ducoudré, N. de, Ogée, J., Polcher, J., Friedlingstein, P., Ciais, P., 795 Sitch, S., Prentice, I.C., 2005. A dynamic global vegetation model for studies of the coupled 796 atmosphere-biosphere system. Glob. Biogeochem. Cycles 19, 33. doi:200510.1029/2003GB002199 797 Lebourgeois, F., Pierrat, J.-C., Perez, V., Piedallu, C., Cecchini, S., Ulrich, E., 2010. Simulating 798 phenological shifts in French temperate forests under two climatic change scenarios and four 799 driving global circulation models. Int. J. Biometeorol. 54, 563–581. doi:10.1007/s00484-010-0305- 800 5 801 Leinonen, I., Hänninen, H., 2002. Adaptation of the timing of bud burst of Norway spruce to 802 temperate and boreal climates. Silva Fenn. 803 Loveland, T.R., Belward, A.S., 1997. The IGBP-DIS global 1km land cover data set, DISCover: 804 First results. Int. J. Remote Sens. 18, 3289–3295. doi:10.1080/014311697217099 805 Lu, X., Wang, Y.-P., Wright, I.J., Reich, P.B., Shi, Z., Dai, Y., 2016. Incorporation of plant traits in 806 a land surface model helps explain the global biogeographical distribution of major forest functional 807 types. Glob. Ecol. Biogeogr. n/a-n/a. doi:10.1111/geb.12535 808 Man, R., Colombo, S., Lu, P., Dang, Q.-L., 2015. Effects of winter warming on cold hardiness and 809 spring budbreak of four boreal conifers. Botany 94, 117–126. doi:10.1139/cjb-2015-0181 810 Migliavacca, M., Sonnentag, O., Keenan, T.F., Cescatti, A., O’Keefe, J., Richardson, A.D., 2012. 811 On the uncertainty of phenological responses to climate change, and implications for a terrestrial 812 biosphere model. Biogeosciences 9, 2063–2083. doi:10.5194/bg-9-2063-2012 813 Murray, M.B., Smith, R.I., Leith, I.D., Fowler, D., Lee, H.S.J., Friend, A.D., Jarvis, P.G., 1994. 814 Effects of elevated CO2, nutrition and climatic warming on bud phenology in Sitka spruce (Picea 815 sitchensis) and their impact on the risk of frost damage. Tree Physiol. 14, 691–706. 816 doi:10.1093/treephys/14.7-8-9.691 817 Muukkonen, P., Lehtonen, A., 2004. Needle and branch biomass turnover rates of Norway spruce 818 (Picea abies). Can. J. For. Res. 34, 2517–2527. doi:10.1139/x04-133 819 Myneni, R.B., Hoffman, S., Knyazikhin, Y., Privette, J.L., Glassy, J., Tian, Y., Wang, Y., Song, X., 820 Zhang, Y., Smith, G.R., Lotsch, A., Friedl, M., Morisette, J.T., Votava, P., Nemani, R.R., Running, 821 S.W., 2002. Global products of vegetation leaf area and fraction absorbed PAR from year one of 822 MODIS data. Remote Sens. Environ., 83, 214–231. doi:10.1016/S0034-4257(02)00074-3 823 Niinemets, Ü., 2002. Stomatal Conductance Alone Does Not Explain the Decline in Foliar 824 Photosynthetic Rates with Increasing Tree Age and Size in Picea Abies and Pinus Sylvestris. Tree 825 Physiol. 22, 515–535. doi:10.1093/treephys/22.8.515 826 Olsson, C., Bolmgren, K., Lindström, J., Jönsson, A.M., 2013. Performance of tree phenology 827 models along a bioclimatic gradient in Sweden. Ecol. Model. 266, 103–117. 828 doi:10.1016/j.ecolmodel.2013.06.026 829 Olsson, C., Jönsson, A.M., 2015. Budburst model performance: The effect of the spatial resolution 830 of temperature data sets. Agric. For. Meteorol. 200, 302–312. doi:10.1016/j.agrformet.2014.10.003 831 Olsson, C., Jönsson, A.M., 2014. Process-based models not always better than empirical models for 832 simulating budburst of Norway spruce and birch in Europe. Glob. Change Biol. 20, 3492–3507. 833 doi:10.1111/gcb.12593 34 834 Pagé, C., Terray, L., Boé, J., 2008. Projections climatiques à échelle fine sur la France pour le 835 21ème siècle: les scénarii SCRATCH08. Clim. Model. Glob. Change CERFACS. 836 Pavlick, R., Drewry, D., Bohn, K., Reu, B., Kleidon, A., 2012. The Jena Diversity-Dynamic Global 837 Vegetation Model (JeDi-DGVM): a diverse approach to representing terrestrial biogeography and 838 biogeochemistry based on plant functional trade-offs. Biogeosciences Discuss 9, 4627–4726. 839 Peaucelle, M., Bellassen, V., Ciais, P., Peñuelas, J., Viovy, N., 2016. A new approach to optimal 840 discretization of plant functional types in a process-based ecosystem model with forest 841 management: a case study for temperate conifers. Glob. Ecol. Biogeogr. doi:10.1111/geb.12557 842 Pokornỳ, R., Stojnič, S., others, 2013. Leaf area index of Norway spruce stands in relation to age 843 and defoliation. Beskydy 5, 173–180. 844 Polgar, C.A., Primack, R.B., 2011. Leaf-out phenology of temperate woody plants: from trees to 845 ecosystems. New Phytol. 191, 926–941. doi:10.1111/j.1469-8137.2011.03803.x 846 Porté, A., Loustau, D., 1998. Variability of the Photosynthetic Characteristics of Mature Needles 847 Within the Crown of a 25-Year-Old Pinus Pinaster. Tree Physiol. 18, 223–232. 848 doi:10.1093/treephys/18.4.223 849 Rautiainen, M., Heiskanen, J., 2013. Seasonal contribution of understory vegetation to the 850 reflectance of a boreal landscape at different spatial scales. IEEE Geosci. Remote Sens. Lett. 10, 851 923–927. 852 Rautiainen, M., Mõttus, M., Heiskanen, J., Akujärvi, A., Majasalmi, T., Stenberg, P., 2011. 853 Seasonal reflectance dynamics of common understory types in a northern European boreal forest. 854 Remote Sens. Environ. 115, 3020–3028. doi:10.1016/j.rse.2011.06.005 855 Reich, P.B., Rich, R.L., Lu, X., Wang, Y.-P., Oleksyn, J., 2014. Biogeographic variation in 856 evergreen conifer needle longevity and impacts on boreal forest carbon cycle projections. Proc. 857 Natl. Acad. Sci. 111, 13703–13708. doi:10.1073/pnas.1216054110 858 Richardson, A.D., Anderson, R.S., Arain, M.A., Barr, A.G., Bohrer, G., Chen, G., Chen, J.M., 859 Ciais, P., Davis, K.J., Desai, A.R., Dietze, M.C., Dragoni, D., Garrity, S.R., Gough, C.M., Grant, 860 R., Hollinger, D.Y., Margolis, H.A., McCaughey, H., Migliavacca, M., Monson, R.K., Munger, 861 J.W., Poulter, B., Raczka, B.M., Ricciuto, D.M., Sahoo, A.K., Schaefer, K., Tian, H., Vargas, R., 862 Verbeeck, H., Xiao, J., Xue, Y., 2012. Terrestrial biosphere models need better representation of 863 vegetation phenology: results from the North American Carbon Program Site Synthesis. Glob. 864 Change Biol. 18, 566–584. doi:10.1111/j.1365-2486.2011.02562.x 865 Richardson, A.D., Keenan, T.F., Migliavacca, M., Ryu, Y., Sonnentag, O., Toomey, M., 2013. 866 Climate change, phenology, and phenological control of vegetation feedbacks to the climate system. 867 Agric. For. Meteorol. 169, 156–173. doi:10.1016/j.agrformet.2012.09.012 868 Rinne, P.L.H., Welling, A., Vahala, J., Ripel, L., Ruonala, R., Kangasjärvi, J., Schoot, C. van der, 869 2011. Chilling of Dormant Buds Hyperinduces FLOWERING LOCUS T and Recruits GA- 870 Inducible 1,3-β-Glucanases to Reopen Signal Conduits and Release Dormancy in Populus. Plant 871 Cell 23, 130–146. doi:10.1105/tpc.110.081307 872 Rohde, A., Bhalerao, R.P., 2007. Plant dormancy in the perennial context. Trends Plant Sci. 12, 873 217–223. doi:10.1016/j.tplants.2007.03.012 35 874 Sampson, D.A., Albaugh, T.J., Johnsen, K.H., Allen, H.L., Zarnoch, S.J., 2003. Monthly leaf area 875 index estimates from point-in-time measurements and needle phenology for Pinus taeda. Can. J. 876 For. Res. 33, 2477–2490. doi:10.1139/x03-166 877 Steiner, K.C., 1980. Patterns of variation in bud-burst timing among populations in several Pinus 878 species. Silvae Genet. 879 Verbeke, T., 2015. Développement et quantification des impacts de l’ozone sur la biosphère 880 continentale dans un modèle global de végétation. PhD thesis. Saint Quentin en yveline, Saint 881 Quentin en yveline, France. 882 Verheijen, L.M., Aerts, R., Brovkin, V., Cavender-Bares, J., Cornelissen, J.H.C., Kattge, J., van 883 Bodegom, P.M., 2015. Inclusion of ecologically based trait variation in plant functional types 884 reduces the projected land carbon sink in an earth system model. Glob. Change Biol. 21, 3074– 885 3086. doi:10.1111/gcb.12871 886 Verma, M., Friedl, M., Richardson, A., Kiely, G., Cescatti, A., Law, B., Wohlfahrt, G., Gielen, B., 887 Roupsard, O., Moors, E., others, 2013. Remote sensing of annual terrestrial gross primary 888 productivity from MODIS: an assessment using the FLUXNET La Thuile dataset. Biogeosciences 889 BG 10, 11627–11669. 890 Verrelst, J., Muñoz, J., Alonso, L., Delegido, J., Rivera, J.P., Camps-Valls, G., Moreno, J., 2012. 891 Machine learning regression algorithms for biophysical parameter retrieval: Opportunities for 892 Sentinel-2 and -3. Remote Sens. Environ. 118, 127–139. doi:10.1016/j.rse.2011.11.002 893 Vetter, M., Churkina, G., Jung, M., Reichstein, M., Zaehle, S., Bondeau, A., Chen, Y., Ciais, P., 894 Feser, F., Freibauer, A., Geyer, R., Jones, C., Papale, D., Tenhunen, J., Tomelleri, E., Trusilova, K., 895 Viovy, N., Heimann, M., 2008. Analyzing the causes and spatial pattern of the European 2003 896 carbon flux anomaly using seven models. Biogeosciences 5, 561–583. 897 Vidal, J.-P., Martin, E., Franchistéguy, L., Baillon, M., Soubeyroux, J.-M., 2010. A 50-year high- 898 resolution atmospheric reanalysis over France with the Safran system. Int. J. Climatol. 30, 1627– 899 1644. doi:10.1002/joc.2003 900 Vitasse, Y., François, C., Delpierre, N., Dufrêne, E., Kremer, A., Chuine, I., Delzon, S., 2011. 901 Assessing the effects of climate change on the phenology of European temperate trees. Agric. For. 902 Meteorol. 151, 969–980. doi:10.1016/j.agrformet.2011.03.003 903 Wang, R., Chen, J.M., 2012. Seasonal leaf area index variations derived from needle growth and 904 fall measurements in two eastern white pine (Pinus Strobes L.) stands, in: 2012 IEEE 4th 905 International Symposium on Plant Growth Modeling, Simulation, Visualization and Applications. 906 pp. 413–417. doi:10.1109/PMA.2012.6524866 907 Wang, Y., Woodcock, C.E., Buermann, W., Stenberg, P., Voipio, P., Smolander, H., Häme, T., 908 Tian, Y., Hu, J., Knyazikhin, Y., Myneni, R.B., 2004. Evaluation of the MODIS LAI algorithm at a 909 coniferous forest site in Finland. Remote Sens. Environ. 91, 114–127. 910 doi:10.1016/j.rse.2004.02.007 911 Warren, C., 2006. Why does photosynthesis decrease with needle age in Pinus pinaster? Trees - 912 Struct. Funct. 20, 157–164. doi:10.1007/s00468-005-0021-7 913 Watanabe, M., Matsuo, N., Yamaguchi, M., Matsumura, H., Kohno, Y., Izuta, T., 2010. Risk 914 assessment of ozone impact on the carbon absorption of Japanese representative conifers. Eur. J. 915 For. Res. 129, 421–430. doi:10.1007/s10342-009-0316-0 36 916 Way, D. A., & Montgomery, R. A. (2015). Photoperiod constraints on tree phenology, performance 917 and migration in a warming world. Plant, Cell & Environment, 38(9), 1725-1736. 918 919 Worrall, J., 1983. Temperature–bud-burst relationship in amabilis and subalpine fir provenance 920 tests replicated at different elevations. Silvae Genet. 921 Xiang, Y., Gubian, S., Suomela, B., Hoeng, J., 2013. Generalized simulated annealing for global 922 optimization: the GenSA Package. R J. 5. 923 Yakovlev, I.A., Asante, D.K.A., Fossdal, C.G., Partanen, J., Junttila, O., Johnsen, O., 2008. 924 Dehydrins expression related to timing of bud burst in Norway spruce. Planta 228, 459–472. 925 doi:10.1007/s00425-008-0750-0 926 927 928 Supporting information captions 929 Appendix SA: Litter sampling methodology, initial conditions and forcing data used for 930 ORCHIDEE simulations, sensitivity analysis 931 Appendix SB: Description, equations and parameters for the eight budburst models 932 Appendix SC: Supplementary figures and tables 37 38