Jianhong Lin
email: jianhong.lin@universite-paris-saclay.fr

Daniel Berveiller

Christophe François

Heikki Hänninen

Alexandre Morfin

Gaëlle Vincent

Rui Zhang

Cyrille Rathgeber

Nicolas Delpierre
email: nicolas.delpierre@universite-paris-saclay.fr

A model of the within-population variability of budburst in forest trees

Keywords: budburst variability, model, temperate trees, climate warming, budburst duration, population

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Phenology, as the study of recurrent biological events such as budburst in spring, has attracted increasing attention due to climate warming [START_REF] Piao | Plant phenology and global climate change: Current progresses and challenges[END_REF]. The timing of leaf phenology in spring is a major indicator of climate warming [START_REF] Parmesan | A globally coherent fingerprint of climate change impacts across natural systems[END_REF] and is mainly modulated by temperature [START_REF] Menzel | European phenological response to climate change matches the warming pattern[END_REF][START_REF] Zhang | Spring phenology in subtropical trees: Developing process-based models on an experimental basis[END_REF][START_REF] Zhang | The chilling requirement of subtropical trees is fulfilled by high temperatures: A generalized hypothesis for tree endodormancy release and a method for testing it[END_REF][START_REF] Chen | Spring phenology at different altitudes is becoming more uniform under global warming in Europe[END_REF]Vitasse et al., 2009a) and photoperiod [START_REF] Delpierre | Temperate and boreal forest tree phenology: from organ-scale processes to terrestrial ecosystem models[END_REF][START_REF] Fu | Daylength helps temperate deciduous trees to leaf-out at the optimal time[END_REF][START_REF] Vitasse | What role for photoperiod in the bud burst phenology of European beech[END_REF][START_REF] Meng | Photoperiod decelerates the advance of spring phenology of six deciduous tree species under climate warming[END_REF]. In the northern hemisphere, it is well established that spring phenological events have been advanced by climate warming [START_REF] Walther | Ecological responses to recent climate change[END_REF][START_REF] Menzel | European phenological response to climate change matches the warming pattern[END_REF], although this advancement is currently slowing down [START_REF] Fu | Declining global warming effects on the phenology of spring leaf unfolding[END_REF][START_REF] Chen | Long-term changes in the impacts of global warming on leaf phenology of four temperate tree species[END_REF]. To date, massive efforts have been made to study the little attention to date [START_REF] Scotti | Fifty years of genetic studies: what to make of the large amounts of variation found within populations?[END_REF][START_REF] Delpierre | Tree phenological ranks repeat from year to year and correlate with growth in temperate deciduous forests[END_REF], which is in line with the general focus of ecological studies on average traits [START_REF] Violle | The return of the variance: intraspecific variability in community ecology[END_REF]. This is intriguing, since the within-population (i.e., tree-to-tree) variability of phenological events is vast and can even be equivalent to that observed among populations [START_REF] Delpierre | Tree phenological ranks repeat from year to year and correlate with growth in temperate deciduous forests[END_REF]Vitasse et al., 2009a;[START_REF] Rathgeber | Cambial activity related to tree size in a mature silver-fir plantation[END_REF]. It typically takes 1 to 4 weeks from the first to the last tree to burst buds in a population [START_REF] Denechere | The withinpopulation variability of leaf spring and autumn phenology is influenced by temperature in temperate deciduous trees[END_REF], with an average of 19 days [START_REF] Delpierre | Tree phenological ranks repeat from year to year and correlate with growth in temperate deciduous forests[END_REF]. Furthermore, the duration from the first to last tree to burst buds in a given population varies annually [START_REF] Denechere | The withinpopulation variability of leaf spring and autumn phenology is influenced by temperature in temperate deciduous trees[END_REF].

The large within-population variability (WPV) of budburst observed in natural tree populations is considered to result from their exposure to a large range of fluctuating environmental (e.g., frost) and biotic (e.g., herbivores and pathogens) selection pressures, which alternatively favor trees that burst buds early or late [START_REF] Alberto | Adaptive responses for seed and leaf phenology in natural populations of sessile oak along an altitudinal gradient[END_REF]. From an evolutionary point of view, this phenotypic diversity has an adaptive value at the population scale, because the environment is likely to change across the lifetime of trees [START_REF] Petit | Some evolutionary consequences of being a tree[END_REF][START_REF] Morente-Lopez | Linking ecological niche models and common garden experiments to predict phenotypic differentiation in stressful environments: Assessing the adaptive value of marginal populations in an alpine plant[END_REF][START_REF] Blanquart | A practical guide to measuring local adaptation[END_REF]. For instance, if a local climate becomes suitable in early spring under climate warming, trees that burst buds early will benefit from an extended growing season, thus maximizing their carbon assimilation and possibly their biomass production [START_REF] Zohner | Leaf-out in northern ecotypes of wide-ranging trees requires less spring warming, enhancing the risk of spring frost damage at cold range limits[END_REF][START_REF] Delpierre | Modelling interannual and spatial variability of leaf senescence for three deciduous tree species in France[END_REF][START_REF] Richardson | Influence of spring and autumn phenological transitions on forest ecosystem productivity[END_REF], which will allow them to gradually occupy a dominant position in the population. Moreover, early budburst enables trees to escape pathogens (e.g., for oak, see [START_REF] Dantec | Escape of spring frost and disease through phenological variations in oak populations along elevation gradients[END_REF]. On the contrary, if freezing events occur frequently in early spring with the advance of budburst, late trees can grow better by avoiding freezing injury [START_REF] Delpierre | Tree phenological ranks repeat from year to year and correlate with growth in temperate deciduous forests[END_REF][START_REF] Zohner | Leaf-out in northern ecotypes of wide-ranging trees requires less spring warming, enhancing the risk of spring frost damage at cold range limits[END_REF][START_REF] Puchalka | Did the late spring frost in 2007 and 2011 affect tree-ring width and earlywood vessel size in Pedunculate oak (Quercus robur) in northern Poland?[END_REF]. Moreover, the WPV also affects interactions with competing plants and herbivores [START_REF] Hart | How variation between individuals affects species coexistence[END_REF][START_REF] Renner | Climate Change and Phenological Mismatch in Trophic Interactions Among Plants, Insects, and Vertebrates[END_REF].

The WPV of budburst is probably underpinned by genetic diversity, as evidenced by the variability of phenological traits among individual trees that experience similar environmental conditions [START_REF] Bontemps | In situ marker-based assessment of leaf trait evolutionary potential in a marginal European beech population[END_REF][START_REF] Delpierre | Tree phenological ranks repeat from year to year and correlate with growth in temperate deciduous forests[END_REF]. This genetic determinism is further reflected in the year-to-year repeatability of the phenological ranking of individuals within tree populations [START_REF] Delpierre | Tree phenological ranks repeat from year to year and correlate with growth in temperate deciduous forests[END_REF]. In addition to this genetic determinism, the WPV is also likely influenced by micro-environmental variations such as the unbalanced distribution of soil-water content within populations, edaphic conditions, or microtopography [START_REF] Delpierre | Tree phenological ranks repeat from year to year and correlate with growth in temperate deciduous forests[END_REF][START_REF] Denechere | The withinpopulation variability of leaf spring and autumn phenology is influenced by temperature in temperate deciduous trees[END_REF][START_REF] Scotti | Fifty years of genetic studies: what to make of the large amounts of variation found within populations?[END_REF].

To the best of our knowledge, the question of whether and to what extent would the WPV of budburst be modified in the current context of climate change has not been addressed so far. Quantifying WPV as the duration (in days) from the first to the last tree to burst buds in one population, we identify three alternative hypotheses for the modification of WPV with climate change: (i) the duration of budburst remains unchanged because all trees in the population advance to the same extent; (ii) the duration of budburst decreases because of the increasing warming rate during the budburst period [START_REF] Malyshev | Inter-Individual Budburst Variation in Fagus sylvatica Is Driven by Warming Rate[END_REF]; (iii) the duration of budburst increases because of insufficient chilling accumulation, as hypothesized previously from experimental studies [START_REF] Zohner | Global warming reduces leaf-out and flowering synchrony among individuals[END_REF][START_REF] Zhang | The chilling requirement of subtropical trees is fulfilled by high temperatures: A generalized hypothesis for tree endodormancy release and a method for testing it[END_REF].

Phenological research has made extensive use of modeling to study the response of the spatiotemporal variability of budburst to climate warming [START_REF] Zhang | Spring phenology in subtropical trees: Developing process-based models on an experimental basis[END_REF][START_REF] Meng | Photoperiod decelerates the advance of spring phenology of six deciduous tree species under climate warming[END_REF][START_REF] Delpierre | Modelling interannual and spatial variability of leaf senescence for three deciduous tree species in France[END_REF][START_REF] Chuine | Process-Based Models of Phenology for Plants and Animals[END_REF].

The models postulate that temperature and photoperiod are the main environmental cues that trigger budburst in boreal and temperate [START_REF] Delpierre | Modelling interannual and spatial variability of leaf senescence for three deciduous tree species in France[END_REF][START_REF] Kramer | A Modeling Analysis of the Effects of Climatic Warming on the Probability of Spring Frost Damage To Tree Species in the Netherlands and Germany[END_REF][START_REF] Hänninen | A framework for modelling the annual cycle of trees in boreal and temperate regions[END_REF], subtropical [START_REF] Zhang | Spring phenology in subtropical trees: Developing process-based models on an experimental basis[END_REF] et al., 2019), and tropical trees [START_REF] Chen | Delayed response of spring phenology to global warming in subtropics and tropics[END_REF]. In these models, temperature plays a dual role. Endodormancy is released by chilling in late autumn or winter, with high temperatures allowing for ontogenetic growth during the ecodormancy stage [START_REF] Hänninen | Boreal and temperate trees in a changing climate: Modelling the ecophysiology of seasonality[END_REF][START_REF] Jewaria | A hundred years after: endodormancy and the chilling requirement in subtropical trees[END_REF]. Meanwhile, there is an interaction between these two stages in the models, namely, ontogenetic growth is influenced by dormancy release [START_REF] Hänninen | Boreal and temperate trees in a changing climate: Modelling the ecophysiology of seasonality[END_REF][START_REF] Hänninen | A framework for modelling the annual cycle of trees in boreal and temperate regions[END_REF][START_REF] Vegis | Dormancy in Higher Plants[END_REF]. [START_REF] Lundell | Beyond rest and quiescence (endodormancy and ecodormancy): A novel model for quantifying plant-environment interaction in bud dormancy release[END_REF] further proved that this interaction can be affected by prevailing temperatures.

One important point is that these models do not pay attention to the WPV of phenological traits. They have been parameterized and applied to predict the mean or median date of budburst in a given tree population [START_REF] Lundell | Beyond rest and quiescence (endodormancy and ecodormancy): A novel model for quantifying plant-environment interaction in bud dormancy release[END_REF][START_REF] Kramer | A Modeling Analysis of the Effects of Climatic Warming on the Probability of Spring Frost Damage To Tree Species in the Netherlands and Germany[END_REF][START_REF] Zhang | Spring phenology in subtropical trees: Developing process-based models on an experimental basis[END_REF]. In other words, these models simulate the timing of budburst as a discrete event in the population without considering the WPV of leaf phenology. To the best of our knowledge, only two studies to date, notably [START_REF] Rousi | Temperature sum accumulation effects on within-population variation and long-term trends in date of bud burst of European white birch (Betula pendula)[END_REF] in Birch (Betula pendula) and [START_REF] Langvall | Frost damage to planted Norway spruce seedlings -influence of site preparation and seedling type[END_REF] in Norway spruce (Picea abies (L.) Karst.), have attempted to establish a link between WPV and environmental conditions through the temperature sum required for the opening of buds at the scale of individual trees. At the scale of tree populations, a distribution of temperature sums to budburst was also used in the so-called physio-demo-genetic (PDG) models [START_REF] Kramer | Bridging the gap between ecophysiological and genetic knowledge to assess the adaptive potential of European beech[END_REF][START_REF] Oddou-Muratorio | Simulating local adaptation to climate of forest trees with a Physio-Demo-Genetics model[END_REF] to simulate the adaptive potential of tree populations. However, a systematic model for the WPV of budburst is still lacking.

Here we developed a model that simulates the WPV of budburst in temperate deciduous trees. We calibrated and validated the model over an extensive budburst dataset acquired from five tree populations at the individual tree scale over 23 years (representing 48,442 observations). Specially, we aim to 1) develop the WPV model and validate its ability for predicting the progress of budburst in tree populations, 2) use the model to in a retrospective simulation exercise testing whether the duration of budburst period in the population changed with climate warming in the recent decades.

Materials and Methods

Study sites

We used budburst data collected from two forests located near Paris (France): Barbeau (48.476° N, 2.780° E, 95 m asl) and Orsay (48.705° N, 2.167° E, 105 m asl). At these sites, the progress of budburst was observed at the individual scale in populations of three major temperate deciduous tree species, namely, hornbeam (Carpinus betulus L.), oak (Quercus petraea (Matt.) Liebl) and chestnut (Castanea sativa Mill.). Hornbeam is an early leafing tree species, chestnut is a late species while oak is intermediate. Hornbeam and oak are present in both forests, while chestnut is present in Orsay only (Table 1). For each species, we focused on healthy and dominant trees, except for hornbeam (an understory species). We collected budburst observations from 2000 to 2022, which yielded a dataset comprising 5 populations and 103 population-years. In each population, we observed between 28 and 309 individuals (mean 90) (Table 1).

Phenology dataset

A team of eight local observers (including most of the authors of this paper) conducted the observations of developing buds in the tree crowns throughout spring. The observers used binoculars and occasionally received training in order to reduce observer bias [START_REF] Liu | Higher sample sizes and observer inter-calibration are needed for reliable scoring of leaf phenology in trees[END_REF]. The interval between phenological observations was of 4 days on average (from 2 to 7 days). A tree was considered to have burst its buds when at least 50% of the buds in the upper third of the crown presented leaves that extended beyond the tip of the scales, which corresponded to stage BBCH 9 [START_REF] Meier | Growth stages of mono-and dicotyledonous plants[END_REF].

At each observation date, we calculated the percentage of trees that had reached budburst in the tree population, dividing the number of trees at BBCH 9 by the total number of trees observed on that date and multiplying the result by 100.

Temperature data

We obtained the mean daily temperature data from the meteorological station nearest to the study sites (Table 1).

However, there were missing values in the temperature data collected from the stations, especially before 1970. To fill these gaps and predict the missing data in order simulate budburst in previous years, we used the SAFRAN reanalysis data (grid-resolution of 8*8 km²) [START_REF] Vidal | A 50-year high-resolution atmospheric reanalysis over France with the Safran system[END_REF], which we de-biased by establishing a linear regression between the local and corresponding SAFRAN temperature data from September of previous year to June.

Model description

We introduce a novel model, named the within-population variability (WPV) model, which was constructed to predict the progress of budburst in tree populations (i.e., percentage of trees having burst buds at a given date in a tree population). We hypothesized that the difference between individuals in the population was reflected in the difference of the forcing accumulation requirement (F*).

We built the WPV model by modifying a state-of-the-art process-based model that simulated a discrete budburst event (i.e., budburst of an individual plant or mean budburst date in a tree population) [START_REF] Lundell | Beyond rest and quiescence (endodormancy and ecodormancy): A novel model for quantifying plant-environment interaction in bud dormancy release[END_REF]. In short, the model represents the release of endodormancy through the accumulation of "chilling" temperatures and simulates the ontogenetic growth of buds through the accumulation of "forcing" temperatures. One particularity of the model is that ontogenetic growth is regulated by the state of rest break and the prevailing temperature [START_REF] Lundell | Beyond rest and quiescence (endodormancy and ecodormancy): A novel model for quantifying plant-environment interaction in bud dormancy release[END_REF][START_REF] Hänninen | Modelling bud dormancy release in trees from cool and temperate regions[END_REF][START_REF] Hänninen | A framework for modelling the annual cycle of trees in boreal and temperate regions[END_REF][START_REF] Vegis | Dormancy in Higher Plants[END_REF]. The ontogenetic competence, Co (a dimensionless [0, 1] multiplier), is applied to represent this regulation [START_REF] Lundell | Beyond rest and quiescence (endodormancy and ecodormancy): A novel model for quantifying plant-environment interaction in bud dormancy release[END_REF][START_REF] Hänninen | A framework for modelling the annual cycle of trees in boreal and temperate regions[END_REF][START_REF] Hänninen | Boreal and temperate trees in a changing climate: Modelling the ecophysiology of seasonality[END_REF]. In the model, budburst is considered to occur at date t when a given sum of the forcing temperature is reached such that 𝐹(𝑡) ≥ 𝐹 * .

In the WPV model, we assumed that F* followed a normal distribution at the level of the tree population (see Fig. S1 for a flow chart of the model). At each date t, the model simulates the proportion of the population (BP, for budburst percent) that has fulfilled the forcing accumulation requirement:

𝐹 * = (𝜇, 𝜎 2 ) eq.1 𝐵𝑃(𝑡) = 0.5 × (1 + erf ( 𝐹(𝑡)-𝜇 (𝜎×𝑠𝑞𝑟𝑡(2))
)) × 100 eq.2

Where F(t) is the forcing degree-day accumulation reached on day t, 𝜇 is the mean of normal distribution, 𝜎 is the standard deviation of normal distribution, and erf is the Gaussian error function.

The forcing accumulation F(t) is calculated as the integral of a "forcing rate" as follows:

𝐹(𝑡) = ∑ 𝑅𝑓 𝑎𝑐𝑡 𝑡 𝑑=270 eq.3
Where d is the start date of forcing accumulation (d = DoY 270 in the previous year). In this model, the stage of dormancy release and the stage of ontogenetic growth can occur simultaneously (i.e., the model belongs to the "parallel" model category) [START_REF] Hänninen | Boreal and temperate trees in a changing climate: Modelling the ecophysiology of seasonality[END_REF][START_REF] Chuine | Process-Based Models of Phenology for Plants and Animals[END_REF]. However, the forcing rate Rfact, which is the actual rate of ontogenetic growth, is affected by both temperatures and ontogenetic competence (Co). It is calculated as follows:

𝑅𝑓 𝑎𝑐𝑡 (𝑡) = 𝑅𝑓(𝑡) * 𝐶𝑜(𝑡) eq. 4

Where Rf(t) is the potential rate of ontogenetic growth at time t, and Co is the ontogenetic competence at time t; these two variables are calculated as follows:

𝑅𝑓(𝑡) = { 0, 𝑇(𝑡) < 𝑇 𝑏 𝑇(𝑡) -𝑇 𝑏 , 𝑇(𝑡) ≥ 𝑇 𝑏 eq. 5
Where Tb is the temperature threshold (°C) above which forcing accumulation occurs.

The ontogenetic competence Co varies over time and is simulated as:

𝐶𝑜(𝑡) = max (0; min (1; 𝑔 × 𝑇(𝑡) + ℎ + 𝑆𝑟(𝑡) 100 * (1 -ℎ))) eq.6
Where Co(t) is the ontogenetic competence at time t in the range [0, 1], which modulates the effect of the state of rest break on the rate of ontogenetic growth (see Fig. S2). When Co=0, ontogenetic growth is stopped. The ability of ontogenetic growth is restored between Co=0 and Co=1 with rest breaking. Finally, g and h are parameters [START_REF] Lundell | Beyond rest and quiescence (endodormancy and ecodormancy): A novel model for quantifying plant-environment interaction in bud dormancy release[END_REF] Where Ccri is the chilling requirement for rest completion, and Ctot is the actual accumulation of chilling temperature, quantified as the number of chilling units (in chill units C.U.) and calculated from DoY=270 of the previous year up to time t as follows:

𝐶 𝑡𝑜𝑡 = ∑ 𝑅𝑐 𝑡 𝑑=270
eq.8

Where the daily rate of chilling accumulation (Rc) is calculated as follows:

𝑅𝑐 = { 1, 𝑇(𝑡) < 𝑇 𝑐 0, 𝑇(𝑡) ≥ 𝑇 𝑐 eq.9
Where Tc is the temperature threshold (°C) below which chilling accumulation occurs.

Parameter estimation

We calibrated the model using budburst data obtained during the period 2000-2016 in Orsay (all three species: hornbeam, oak, chestnut) and then validated it using data from 2017-2022 in Orsay (three species) and from 2000-2022 in Barbeau (two species: hornbeam and oak). The model was therefore calibrated over 17 years for the three species (Orsay populations, representing 52, 71 and 50 observation dates for hornbeam, oak and chestnut, respectively)

and validated over 29 site-years for hornbeam and oak (representing 89 and 114 observation dates, resp.), and 6 years (29 observation dates) for chestnut. A previous study (Vitasse et al., 2009b) provided evidence of similar apparent phenological responses to temperature among populations of the same species located as far as 650 km apart, which also suggests the low differentiation of phenological traits across populations. Orsay and Barbeau populations are separated by a distance of 50 km and experience a similar climate. This is why we used the Barbeau data as a validation counterpart to the Orsay data used for calibration. The model predicts the percentage of budburst in the population (from 0% to 100% budburst) along with the corresponding date. Thus, we calculated the root mean square error (RMSE) over two dimensions (Fig. S3). First, we calculated RMSE over the percentage of budburst in the tree population (i.e., comparing the difference between the observed and predicted budburst percent occurring on the same day of the year, DoY).

𝑅𝑀𝑆𝐸 𝐵𝑃 = √ ∑ (√𝑛𝑢𝑚×(𝐵𝑃 𝑜𝑏𝑠,𝑖 -𝐵𝑃 𝑝𝑟𝑒𝑑,𝑖 ) 2 ) 𝑛 𝑖=1 ∑ √𝑛𝑢𝑚 𝑛 𝑖=1
eq.10

Where RMSEBP is the root mean square error for budburst percent (expressed in percent), num is the number of trees observed on a given day of the year, BPobs,i is the observed percentage of budburst of datum i, BPpred,i is the predicted percentage of budburst of same datum, and n is the total number of data (e.g., n=50 in a hypothetic case where the percentage of budburst has been observed five times per year on average over 10 years in a given population). We used √𝑛𝑢𝑚 as a weight in the calculation of squared errors to compensate for the fact that a very large number of trees (i.e., >300 trees) were observed at some dates: these observations are more representative of the actual percentage of budburst in the population (as compared to observations established for a smaller number of trees), although they also tend to overrepresent them in the calculation of errors. We then calculated the RMSE of dates (i.e., comparing the difference, in number of days, between the observations and predictions for the same percentage of budburst; Fig. S3).

𝑅𝑀𝑆𝐸 𝐷𝑜𝑌 = √ ∑ (√𝑛𝑢𝑚×(𝐷𝑜𝑌 𝑜𝑏𝑠,𝑖 -𝐷𝑜𝑌 𝑝𝑟𝑒𝑑,𝑖 ) 2 ) 𝑛 𝑖=1 ∑ √𝑛𝑢𝑚 𝑛 𝑖=1
eq.11

Where RMSEDoY is the root mean square error for the budburst date (in days), num is the number of trees observed, DoYobs, i is the observed date of budburst of datum i (e.g., the date when we observed 24% budburst for the population of interest in a given year), DoYpred,i is the predicted date of budburst of the same datum (e.g., the date when the model predicted 24% budburst in the same tree population and year), and n is the total number of data.

Finally, we calculated the total RMSE as follows:

𝑅𝑀𝑆𝐸 𝑡𝑜𝑡 = 𝑅𝑀𝑆𝐸 𝐵𝑃 + 𝑅𝑀𝑆𝐸 𝐷𝑜𝑌 eq.12

In the calibration stage, we determined the best parameter set as the one that minimized RMSEtot.

In addition to RMSE, we also used mean bias error to evaluate the model forecast accuracy (in terms of budburst percentage or DoY), which is calculated as follows:

𝑚𝑒𝑎𝑛 𝑏𝑖𝑎𝑠 = 1 𝑁 ∑ (𝑜𝑏𝑠 𝑖 -𝑝𝑟𝑒𝑑 𝑖 ) 𝑁 𝑖=1
Where obsi and predi are the i-th observation and prediction, respectively, N is the number of observations.

Evaluating the modelled F* distributions

To validate the modelled F* distribution, we simulated the distribution of the forcing accumulation at the date of each BP observation. Because there are different observed BP in each year. We binned the observed BP data into 11 groups (i.e., BP0, BP10, BP20…BP100, e.g., we regard the data between BP5 (date at which 5% trees burst buds) to BP15

(date at which 15% trees burst buds) as group "BP10"; note that BP0 refers to dates at which 5% or less trees have burst buds, and BP100 refers to dates at which 95% or more trees have burst buds). Then we used a sigmoid function to simulate the relation between BP and averaged corresponding forcing accumulation across all the years. We also calculated their first derivatives (i.e., the increasing of BP per unit of forcing accumulation). Moreover, we calculated the distribution of observed BP across all the years.

Evaluating the response of the within-population variability of budburst to climate warming

We used our model to predict budburst in the past (1961-2022) using historical daily mean temperature data and gapfilled data using debiased SAFRAN reanalysis of temperatures (see above).

As explained earlier, our model simulates the percentage of budburst in a tree population at a given date. To evaluate the response of the WPV of budburst to climate warming, we focused on the particular dates at which 20% and 80%

of trees in a given population had reached budburst (termed BP20 and BP80, respectively) and the duration between these two dates (DurBB = BP80-BP20), which we consider to represent the variability of budburst within the population for a given year. BP20 represents the "beginning" of budburst in the tree population, whereas BP80 represents its "end." We chose these quantiles instead of more extreme quantiles of distribution (e.g., 5% and 95%), because they are well represented in our dataset (Fig. 1), thus implying higher model accuracy. For sake of model evaluation, we calculated the DurBB in observed phenology data. Specifically, we selected years which had records before BP20 and after BP80. Then the date of BP20 or BP80 were calculated by using the nearest two data (one is below BP20 or BP80, another is above BP20 or BP80) through interpolation (e.g., 15 % budburst percent is on DoY 80 and 25 % budburst percent is on DoY 84. We can obtain the date of BP20 by interpolation, that is DoY 82).

Statistical analyses

For each population, we quantified by linear regression the sensitivity of budburst date (BP20 and BP80) and the DurBB to time (days year -1 ) and to Jan-May temperature (days o C -1 ). Analysis of Variance (ANOVA) was used to analysis whether the significance of the regression slopes (=0.05). All simulations and statistical analyses were carried out with R statistical software v.4.0.3 (R Development Core Team, 2020). Orsay were very close, with average differences of 2 and 1 days (Table S1). The duration of budburst in the population (DurBB) (i.e., time interval, in days, during which the proportion of trees having reached budburst increases from 20%

to 80%) differs for each species depending on the site and year, with a mean of 8 days over the whole dataset and ranging from 3 days for hornbeam at Orsay in 2018 and 2021 to 21 days for oak at Orsay in 2012 (Fig. 1).

Model performance

For all the populations considered here, the WPV model predicted with good accuracy the progress of budburst in tree populations during spring as well as the interannual variability of budburst (Fig. 2, Fig. 3; see Fig. S4 for a comparison of observed and simulated time series). The model predicted the percentage of budburst in tree populations with an error (RMSEBP) of 16% ± 0.2% for the calibration dataset and 20% ± 2.1% for the validation dataset. This corresponded to prediction errors for the date of budburst (RMSEDOY) of 3.9 ± 0.6 days for the calibration dataset and 5.6 ± 0.3 days for the validation dataset. This compared well to the time resolution of the phenological observations (3-7 days). The mean bias was within 1 day (Fig. 3). F* had a highest mean and standard deviation for oak compared with hornbeam and chestnut (Fig. 4, Table 2). The distributions of F* compared well to the actual distribution of forcing accumulation established from observations (Fig.

5b, e, h), validating the choice of the normal distribution. However, the modelled distribution did not overlap exactly the distribution established from observed data, because the distribution of observations along the BP scale was uneven (Fig. 5c,f,i). The temperature threshold for chilling accumulation (Tc) ranged from 9.7 o C for chestnut to 10.5 o C for hornbeam and oak (Table 2). The temperature threshold for forcing accumulation (Tf) ranged from 3.9 o C for hornbeam to 7.7 o C for chestnut (Table 2, Fig. S2). In all species, buds could not begin ontogenetic growth until the accumulation of chilling to a certain extent (i.e., parameter h was negative for all populations, Table 2). We found that the threshold of chilling accumulation necessary for the onset of forcing accumulation (i.e., value of Sr(t) from which Co becomes positive) was very high for early species and decreased for late species (e.g., value of h increased approximately from -0.98 in hornbeam to 0 in chestnut; see Table 2 and Fig. S2). Prevailing temperatures could compensate for the lack of chilling accumulation (positive parameter g; Table 2) in hornbeam and oak, but not in chestnut (g=0).

Retrospective analysis for within-population variability of budburst

Over the past six decades , spring average temperature increased by +1.9 o C in Orsay and +1.4°C in Barbeau (Fig. S5). Over this time period, our retrospective simulations suggest that the beginning (20%, BP20) and end (80%, BP80) of budburst in tree populations has advanced significantly for all the species (Fig. 6), with respectively 1.7 ± 0.6 days decade -1 (mean ± SD across species) and 1.8 ± 0.4 days decade -1 and apparent temperature sensitivities of 5.8 ± 0.4 days o C -1 and 5.6 ± 0.4 days o C -1 . These similar trends regarding the beginning and end of budburst result in an unchanged duration of the budburst period (DurBB in the considered populations over the past 62 years (no trend in DurBB is significantly different from zero in Fig. 7, P>0.05). Notably, the interannual variability of DurBB was large (Fig. 6), and fairly simulated by our model (RMSE of 3.4 ± 1.8 days).

Discussion

To the best of our knowledge, this paper presents the first model simulating the within-population variability of budburst in tree populations. An added value of this model is that it can simulate the duration of budburst in tree populations. The central hypothesis of the model is that F*, the amount of accumulated forcing temperature required for trees to budburst, follows a normal distribution in tree populations. The ability of the model to simulate the dynamics of budburst over the calibration and validation data, as well as the good agreement between the observed and the simulated F* distributions (Fig. 5), lend support to this hypothesis for all the species and populations considered.

Our model yielded RMSE for the validation data (5.4 to 5.9 days), which are close to the temporal resolution of the https://doi.org/10.5194/egusphere-2023-1075 Preprint. Discussion started: 6 July 2023 c Author(s) 2023. CC BY 4.0 License.

spring phenology observation (from 2-7 days) and similar to the typical prediction accuracy of models simulating discrete (i.e., population average) budburst dates (e.g., [START_REF] Basler | Evaluating phenological models for the prediction of leaf-out dates in six temperate tree species across central Europe[END_REF].

The variability in the timing of budburst among individuals in tree populations is considered to be mainly determined by genetic diversity [START_REF] Bontemps | In situ marker-based assessment of leaf trait evolutionary potential in a marginal European beech population[END_REF][START_REF] Delpierre | Tree phenological ranks repeat from year to year and correlate with growth in temperate deciduous forests[END_REF][START_REF] Jarvinen | DNA sequence variation in BpMADS2 gene in two populations of Betula pendula[END_REF][START_REF] Rousi | Temperature sum accumulation effects on within-population variation and long-term trends in date of bud burst of European white birch (Betula pendula)[END_REF][START_REF] Rusanen | Genetic structure of Acer platanoides and Betula pendula in northern Europe[END_REF] followed by the influence of the microenvironment [START_REF] Delpierre | Tree phenological ranks repeat from year to year and correlate with growth in temperate deciduous forests[END_REF][START_REF] Rousi | Temperature sum accumulation effects on within-population variation and long-term trends in date of bud burst of European white birch (Betula pendula)[END_REF]. The phenological ranking of individuals is largely conserved in tree populations [START_REF] Delpierre | Tree phenological ranks repeat from year to year and correlate with growth in temperate deciduous forests[END_REF], leading to the identification of "early," "intermediate," and "late" trees [START_REF] Malyshev | Inter-Individual Budburst Variation in Fagus sylvatica Is Driven by Warming Rate[END_REF]. Further, the distribution of budburst categories is not uniform in natural tree populations, with numerous "intermediate" individuals and comparatively fewer "early" and "late" trees [START_REF] Malyshev | Inter-Individual Budburst Variation in Fagus sylvatica Is Driven by Warming Rate[END_REF][START_REF] Chesnoiu | Bud burst and flowering phenology in a mixed oak forest from Eastern Romania[END_REF][START_REF] Zohner | Global warming reduces leaf-out and flowering synchrony among individuals[END_REF][START_REF] Caradonna | Shifts in flowering phenology reshape a subalpine plant community[END_REF], which lends further support to a unimodal distribution such as the normal law. Our model reproduces this phenomenon, with categories of "early," "intermediate," and "late" trees corresponding to increasing values of F*. This core assumption of the model is supported by previous empirical studies, which observe that the variability of F* could represent the variability of budburst among trees [START_REF] Langvall | Frost damage to planted Norway spruce seedlings -influence of site preparation and seedling type[END_REF][START_REF] Rousi | Temperature sum accumulation effects on within-population variation and long-term trends in date of bud burst of European white birch (Betula pendula)[END_REF]. Nevertheless, we could have chosen to assign the variance among individuals to one or several other parameters of the model, related to the fact that genetic variations may affect any of the plant traits determining the modelled parameters. For instance, [START_REF] Gauzere | Sensitivity analysis of tree phenology models reveals increasing sensitivity of their predictions to winter chilling temperature and photoperiod with warming climate[END_REF] found that the temperature yielding mid-forcing during ecodormancy (T50) was more sensitive than F* in the UniChill model, which suggests that this parameter is another good candidate for identifying the phenological behavior of individual trees in a population. Thus, we constructed a model assuming that the threshold for forcing temperature (Tb, i.e., parameter of our model analogous to T50) followed a normal distribution, whereas F* was fitted as a constant parameter for the population. This model fitted the data less effectively in both the calibration and validation steps (see Fig. S6 and S7 compared with Fig. 2 and 3), which further supports our decision to assign the among-individual variance to F*. Questions remain regarding the actual shape of the F* distribution. Indeed, natural selection can lead to traits that are not normally distributed [START_REF] Caradonna | Shifts in flowering phenology reshape a subalpine plant community[END_REF], and uneven distribution of observations may contribute to the non-perfect overlapping of observed and simulated F* distributions (Fig. 5). However, earlier results [START_REF] Vallet | Modélisation de la dynamique intra-populationnelle du débourrement en Ile-de-France[END_REF] showed that the form of the distribution had little influence on the prediction accuracy.

We built the WPV model based on a two-phase parallel model framework, which describes the cumulative effect of chilling and forcing temperatures on the endodormancy and ecodormancy phases, respectively [START_REF] Hänninen | Boreal and temperate trees in a changing climate: Modelling the ecophysiology of seasonality[END_REF][START_REF] Hänninen | A framework for modelling the annual cycle of trees in boreal and temperate regions[END_REF][START_REF] Lundell | Beyond rest and quiescence (endodormancy and ecodormancy): A novel model for quantifying plant-environment interaction in bud dormancy release[END_REF][START_REF] Chuine | Process-Based Models of Phenology for Plants and Animals[END_REF]. This model structure is in line with our current understanding of the physiological and molecular basis of dormancy in which the dynamics of the dormancy mechanism are emphasized as opposed to a strict classification between the dormancy stages [START_REF] Lundell | Beyond rest and quiescence (endodormancy and ecodormancy): A novel model for quantifying plant-environment interaction in bud dormancy release[END_REF][START_REF] Cooke | The dynamic nature of bud dormancy in trees: environmental control and molecular mechanisms[END_REF]. In this study, the threshold of chilling accumulation is up to 10.5 o C for oak and hornbeam.

It is consistent with the experimental results in [START_REF] Baumgarten | Chilled to be forced: the best dose to wake up buds from winter dormancy[END_REF] which challenge the common assumption that optimal chilling temperatures range ca. 4-6°C, showing 10 o C is also effective for chilling accumulation in six dominant temperate European tree species including oak. Furthermore, the model uses the concept of ontogenetic competence (Co) to simulate the process of regulation for the rate of ontogenetic growth by the state of rest break, a phenomenon that has found support in phenological experiments [START_REF] Lundell | Beyond rest and quiescence (endodormancy and ecodormancy): A novel model for quantifying plant-environment interaction in bud dormancy release[END_REF][START_REF] Zhang | Spring phenology in subtropical trees: Developing process-based models on an experimental basis[END_REF] that in the investigated species, Co is 0 until dormancy is released to a certain extent (Fig. S2), that is, ontogenetic growth cannot start before a certain amount of chilling accumulation has been reached, which is consistent with previous findings [START_REF] Lundell | Beyond rest and quiescence (endodormancy and ecodormancy): A novel model for quantifying plant-environment interaction in bud dormancy release[END_REF][START_REF] Zhang | Spring phenology in subtropical trees: Developing process-based models on an experimental basis[END_REF]. According to the calibrated parameter values, ontogenetic competence is also influenced by the prevailing temperature, although the effect is minimal. Indeed, parameter g, which is related to the effect of the prevailing temperature, ranges from 0 to 0.0080 (Table 2), which is comparable to values found in a previous study [START_REF] Lundell | Beyond rest and quiescence (endodormancy and ecodormancy): A novel model for quantifying plant-environment interaction in bud dormancy release[END_REF]. To some extent in this model, one consequence is that the effect of the prevailing temperature can compensate for the deficiency of chilling accumulation.

Beyond introducing a model to describe the WPV of budburst in tree populations, our study aimed to quantify the response of the duration of budburst (DurBB) to climate warming. We used temperature data to simulate the occurrence of 20% (BP20) and 80% (BP80) budburst, and DurBB over the past decades. Our results suggest that the start and end of budburst in tree populations have advanced over the past 62 years with climate warming (Fig. 6), which is consistent with previous results showing advances in the population average dates of budburst [START_REF] Wenden | Shifts in the temperature-sensitive periods for spring phenology in European beech and pedunculate oak clones across latitudes and over recent decades[END_REF][START_REF] Menzel | European phenological response to climate change matches the warming pattern[END_REF][START_REF] Fu | Declining global warming effects on the phenology of spring leaf unfolding[END_REF]. In addition, our model simulates sensitivities of budburst to time and temperature that are comparable to values reported earlier (Vitasse et al., 2009b, see Table S2). Our results point to significant sensitivities to both time and temperature for oak as well as significant sensitivity to temperature for hornbeam, which is consistent with the results of Vitasse et al. (2009b).

Our retrospective simulations suggest that there was not trend in the duration of budburst in tree populations, DurBB, over the past 62 years (Fig. 7), in spite of climate warming (Fig. S5). Since both BP20 and BP80 advanced at a similar rate, DurBB did not evolve over time over the 1961-2022 period. Interestingly, the analysis of temperature data revealed no significant warming in the period of time from BP20 to BP80 over the past decades (P>0.05, Fig. S8). This could explain why DurBB (time interval between BP20 and BP80) did not change over time, in spite of the strong trends in both BP20 and BP80, caused by climate warming. Moreover, our study sites are located in the temperate zone, at the heart (for oak and hornbeam) and at the north (chestnut) of our study species distribution areas [START_REF] Caudullo | Chorological maps for the main European woody species[END_REF]. At those sites, trees can accumulate enough chilling, or at least, chilling accumulation is not a limitation for ontogenetic growth in nature so far, meaning that budburst is still advancing [START_REF] Wenden | Shifts in the temperature-sensitive periods for spring phenology in European beech and pedunculate oak clones across latitudes and over recent decades[END_REF][START_REF] Piao | Plant phenology and global climate change: Current progresses and challenges[END_REF].

Thus, the phenomenon by which DurBB increased with insufficient chilling accumulation in a given population (see [START_REF] Zhang | The chilling requirement of subtropical trees is fulfilled by high temperatures: A generalized hypothesis for tree endodormancy release and a method for testing it[END_REF], their Fig. 2, 3, 4 for evidence in subtropical trees) did not appear in our retrospective simulations.

However, we can infer that if chilling accumulation can't be fulfilled under future, continuous climate warming, it will take more time to fulfill the forcing requirement for late trees with a high forcing requirement, leading to the prolonging of DurBB. A longer duration of budburst would increase the possibility of damage (i.e., freezing, insect damage).

These results suggest that the WPV of budburst should be given greater attention, because the longer duration of budburst may be an important factor in the future when researchers project the damage in forests or determine the best strategy for forest management.

We acknowledge that the projections of the WPV of budburst produced by the model are uncertain, first and foremost because the parameter values were inferred from observation data collected in natural conditions as opposed to controlled experiments [START_REF] Hanninen | Experiments Are Necessary in Process-Based Tree Phenology Modelling[END_REF]. Another cause of uncertainty is the ability of the phenological response of plants to acclimatize to the changing climate [START_REF] Bennie | Predicting spatial and temporal patterns of bud-burst and spring frost risk in north-west Europe: the implications of local adaptation to climate[END_REF]. Under the hypothesis of plant acclimatization, the parameters of the WPV model could have changed over the past decades, and would further change with ongoing climate warming. Consequently, related experiments are urgently needed to improve our understanding of the WPV of budburst to infer more reliable parameters and analyze the behavior of phenology models in different climates [START_REF] Hanninen | Experiments Are Necessary in Process-Based Tree Phenology Modelling[END_REF]. However, because our model addresses for the first time explicitly the within-population variation of the physiological traits affecting phenology, it can contribute as a framework for future experimental studies. In our study, we only considered the effect of temperature on budburst. However, other environmental factors may also affect budburst (e.g., photoperiod). Previous studies showed that photoperiod is expected to modulate the timing of budburst in late-successional species such as oak and chestnut, but not in early-successional species such as hornbeam [START_REF] Basler | Photoperiod sensitivity of bud burst in 14 temperate forest tree species[END_REF], but see a counter-example on oak in [START_REF] Malyshev | Temporal photoperiod sensitivity and forcing requirements for budburst in temperate tree seedlings[END_REF]. Moreover, photoperiod may have a more complex interaction mechanism with temperature in terms of regulating the time of budburst [START_REF] Meng | Photoperiod decelerates the advance of spring phenology of six deciduous tree species under climate warming[END_REF]. We envision that improved versions of the WPV of budburst could be proposed based on a more comprehensive understanding of the potential mechanism between phenology and environmental factors in the future.

Conclusion

In conclusion, our work presents a novel model, simulating the continuity of budburst in tree populations in spring.

This phenological model can be adapted to the study of other stages of the tree phenological cycle, which are all of continuous nature in tree populations (e.g., leaf senescence, wood formation etc.). We found budburst was advanced in the past 62 years due to climate warming. However, the duration of budburst period of population was not affected by increasing temperature. This is the first model simulating the within population variability of budburst in the population. It provides a basis for implementation of a module in models directly interested in the within-population variability of phenological and other functional traits (e.g., physio-demo-genetic models). It can also be used as a standalone, to study the dynamics of phenological traits from the scale of individuals to the population and community in the context of climate change.

  https://doi.org/10.5194/egusphere-2023-1075 Preprint. Discussion started: 6 July 2023 c Author(s) 2023. CC BY 4.0 License.

  https://doi.org/10.5194/egusphere-2023-1075 Preprint. Discussion started: 6 July 2023 c Author(s) 2023. CC BY 4.0 License.

  https://doi.org/10.5194/egusphere-2023-1075 Preprint. Discussion started: 6 July 2023 c Author(s) 2023. CC BY 4.0 License.

  https://doi.org/10.5194/egusphere-2023-1075 Preprint. Discussion started: 6 July 2023 c Author(s) 2023. CC BY 4.0 License.

  https://doi.org/10.5194/egusphere-2023-1075 Preprint. Discussion started: 6 July 2023 c Author(s) 2023. CC BY 4.0 License.

Fig. 1 .

 1 Fig. 1. Observed percentage of budburst in five tree populations during the period 2000-2022. The size of the points is scaled

Figure 1

 1 Figure1shows the observed percentages of budburst in the five tree populations monitored from 2000 to 2022. These

  https://doi.org/10.5194/egusphere-2023-1075 Preprint. Discussion started: 6 July 2023 c Author(s) 2023. CC BY 4.0 License.

Fig. 2 .

 2 Fig. 2. Evaluation of the within-population variability (WPV) model predicting the budburst percentage over calibration (red points) and validation (blue points) data. The points of circle are observed in Orsay and of triangle are observed in Barbeau. The points establish the correspondence between the observed and predicted percentage of budburst on an observation day in the population of interest. The one-to-one relation is shown as the black line. RMSE which is root mean square error for the budburst percentage and mean bias are shown. There are 52, 71 and 50 points (i.e., observation dates) for calibration and 89, 114, 29 points for validation for hornbeam, oak and chestnut, respectively.

Table 2 .

 2 Parameter values of the WPV model for three populations. μ (°C-days) and σ (°C-days) are the mean and standard deviation of the distribution of F*, respectively (Eqn. 1). Tb and Tc (°C) are the threshold temperatures for the accumulation of forcing and chilling temperatures, respectively (Eqns. 5 and 9). g ( o C -1 ) and h (dimensionless) are the parameters determining the interactive effect of the state of rest break and the prevailing air temperature on the ontogenetic competence (Eqn. 6). Ccri (number of days) is the chilling requirement of rest completion.

Fig. 4 .Fig. 5 .

 45 Fig. 4. Normal distribution of the forcing requirement (F*) for three tree species.

Fig. 6 .Fig. 7 .

 67 Fig. 6. Simulated occurrence of the beginning (20%, BP20 in blue) and end (80%, BP80 in red) of budburst using the WPV model for three tree species during the period 1961-2022. The fitted lines highlight the trends over the past 62 years. Text in blue (red) shows the sensitivity of BP20 (BP80) to time and mean spring temperature (from January to May), respectively.

  https://doi.org/10.5194/egusphere-2023-1075 Preprint. Discussion started: 6 July 2023 c Author(s) 2023. CC BY 4.0 License.

Table 1 . Description of the phenological and meteorological datasets.

 1

								Number of
	Phenology Site	Coordinate	Meteorological station	Coordinate Species	Number of year	Number of data	trees (min / max /	Observation years
								average)
					Quercus	23	153	29/190/85	2000-2022
	Orsay	48.705° N, 2.165° E	Gometz-le-Châtel	48.677° N, 2.136° E	Carpinus Castanea	20 21	124 112	29/146/50 29/192/80	2002-2006, 2008-2022 2000-2007, 2010-2022
	Barbeau	48.476° N, 2.780° E	Châtelet-en-Brie	48491° N, 2.802° E	Quercus Carpinus	20 19	87 64	29/309/154 28/241/114	2003-2022 2004-2022

Acknowledgements

We acknowledge Eric Dufrêne for setting up the phenological surveys. We are also grateful to Eric Dufrêne and Jean-Yves Pontailler for their invaluable contributions regarding the collection of phenological data. This work was supported by the China Scholarship Council (202008330320).

Code and data availability

The related phenology data and R code for the phenological model are openly accessible under https://doi.org/10.5281/zenodo.7962840 and https://doi.org/10.5281/zenodo.7188160, respectively.

273

As mentioned earlier, we assumed that the forcing requirement (F*) followed a normal distribution. The calibration 274 procedure yielded a set of distribution curves that differed across species (Fig. 4). We observed that the distribution of

Authors' contributions

ND and JL designed the research. ND, JL, AM, GV, DB collected phenological data. JL and ND performed the research. JL wrote the manuscript with substantial inputs from all co-authors.

Competing interests

The authors declare that they have no conflict of interest.