Overestimation of the effect of climatic warming on spring phenology due to misrepresentation of chilling

Phenological and climatic data

We used data from the Pan European Phenology Project (PEP725)²⁹, an open-access database with long-term plant phenological observations across 25 European countries (http://www.pep725.eu/). The regional/national network partners of PEP725 are following a consistent guideline for phenological observations³⁰ and prepare the data for submission to the PEP725 database curators²⁹. We selected 30 species for which sufficient observational data were available: 21 deciduous broadleaved trees or shrubs, 6 herbaceous perennials, 2 evergreen coniferous trees, and 1 deciduous coniferous tree (Supplementary Table 3). Particularly, our data set included one fruit tree (Prunus avium) and one nut tree (Corylus avellana) since some of the chilling models are specifically developed for fruit and nut trees. A total of 2,493,644 individual records from 15,533 phenological stations were used. The stations were mainly distributed in moderate climates in Central Europe (Supplementary Fig. 3). Four spring events based on the BBCH code were investigated: BBCH 10, 11, 60, and 69, representing first leaves separated, first leaves unfolded, first flowers open, and end of flowering, respectively³¹.

We used the E-OBS v19.0eHOM data set³² with a spatial resolution of 0.1 × 0.1° for 1950–2018 for calculating CA and HR of the in situ phenological records. This data set is provided by the European Climate Assessment & Data set project and includes homogenized series of daily mean, minimum, and maximum temperatures. We also use the daily maximum and minimum temperature data from the GHCN data set³³ to assess the scale effect. The GHCN data set contains station-based measurements from over 90,000 land-based stations worldwide, but only parts of PEP725 stations match with the GHCN stations.

For future climatic data (2019–2099), we used daily minimum and maximum temperatures simulated by the HADGEM2-ES model (with a spatial resolution of 0.5 × 0.5°) under two climatic scenarios (RCP 4.5 and RCP 8.5). These data have been bias-corrected by applying the method used in the Inter-Sectoral Impact Model Intercomparison Project (ISIMIP)³⁴, which were available on the ISIMIP server (https://esg.pik-potsdam.de/projects/isimip2b/).

Chilling models

We used 12 chilling models to measure the amount of chilling. One type of chilling model is based on several specific temperature thresholds. The most commonly used model, developed in the 1930s and 1940s for peach³⁵, calculates the number of hours or days with temperatures <7.2 or 7 °C. Another commonly used upper-temperature threshold for chilling is 5 °C^11,18. Some studies, however, have suggested that subfreezing temperatures were effective^36,37, so we also tested chilling models using −10 °C as the lower limit. Many studies have also assumed that freezing temperatures did not contribute to winter CA and only included temperatures >0 °C for calculating CA^11,38. Models C₁–C₆ were developed based on various combinations of the upper and lower temperature limits. The rate of chilling was 1 for daily temperatures <5 °C for Model C₁ (Eq. (1)), between −10 and 5 °C for Model C₂ (Eq. (2)), between 0 and 5 °C for Model C₃ (Eq. (3)), <7 °C for Model C₄ (Eq. (4)), between −10 and 7 °C for Model C₅ (Eq. (5)), and between 0 and 7 °C for Model C₆ (Eq. (6)). The equations for Models C₁–C₆ are as follows:

where CU_i is the rate of chilling for Model C_i, and T is the daily mean temperature (°C).

Model C₇ is also known as the Utah Model³⁹, which assigned different weights to different ranges of temperatures and was first used to measure the chilling requirements of peach (Eq. (7)). The Utah Model was modified to produce Model C₈ (Eq. (8)), which removed the negative contributions of warm temperatures to accumulated chilling²⁰.

where CU_i is the rate of chilling for Model C_i, and T is the daily mean temperature (°C).

Model C₉ is a dynamic model developed for peach in Israel and South Africa^40,41 and now adopted for apricot cultivars⁴². The most important characteristic of Model C₉ was that a previous intermediate product affected the rate of chilling in the current hour or day. We did not provide equations for Model C₉ for simplicity (see the equation in Luedeling et al.²⁰).

Harrington et al.⁴³ summarized published results for chilling units and constructed a chilling function based on a three-parameter Weibull distribution, coded as Model C₁₀ (Eq. (9)). Model C₁₁ has a triangular form, which was fitted by Hänninen⁴⁴ using previous experimental results for Finnish birch seedlings (Eq. (10)). Zhang et al.⁴⁵ recently fitted observational data to the triangular model for 24 plant species and found that a mean optimal chilling temperature of 0.2 °C and an upper limit of the chilling temperature of 6.9 °C were most effective. Model C₁₂, therefore, uses the triangular form with parameters of 0.2 and 6.9 °C (Eq. (11)).

where CU_i is the rate of chilling for Model C_i, and T is the daily mean temperature in °C.

Forcing models

Forcing models were used to measure HR for the spring events of plants. The GDD model is the most commonly used forcing model, which assumes that the rate of forcing is linearly correlated with temperature if the temperature is above a particular threshold. We mainly used Model F₁ (Eq. (12)), which adopts a temperature threshold of 0 °C^46,47,48, for examining the relationship between CA and HR:

where FU₁ is the rate of forcing for Model F₁, and T is the daily mean temperature (°C).

We also validated the chilling models by correlating them with HR based on seven other forcing models to assess the impact of the choice of forcing model on the results. Model F₂ (Eq. (13)) is also a GDD model but has a temperature threshold of 5 °C^12,18:

where FU₂ is the rate of forcing for Model F₂, and T is the daily mean temperature (°C).

Piao et al.⁴⁸ found that leaf onset in the Northern Hemisphere was triggered more by daytime than nighttime temperature. They thus proposed a GDD model using maximum instead of mean temperature. Models F₃ (Eq. (14)) and F₄ (Eq. (15)) are based on maximum temperature with thresholds of 0 and 5 °C, respectively.

where FU_i is the rate of forcing for Model F_i, and T_max is the daily maximum temperature (°C).

A recent experiment demonstrated that the impact of daytime temperature on leaf unfolding for temperate trees was approximately threefold higher than the impact of nighttime temperature⁴⁹. Model F₅ (Eq. 16) thus uses two parameters (0.75 and 0.25) to weigh the impact of daytime and nighttime temperatures on HR.

where FU₅ is the rate of forcing for Model F₅. T_max and T_min are the daily maximum and minimum temperatures (°C), respectively.

Many studies have suggested that the rate of forcing followed a logistic function of temperature^44,50. Model F₆ (Eq. (17)) uses a logistic function proposed by Hänninen⁴⁴, and Model F₇ (Eq. (18)) uses another logistic function proposed by Harrington et al.⁴³.

where FU_i is the rate of forcing for Model F_i, and T is the daily mean temperature (°C).

Model F₈ is a growing degree hour (GDH) model, where species have an optimum temperature for growth and where temperatures above or below that optimum have a smaller impact⁵¹. Model F₈ (Eq. (19)) was first designed for calculating HR at hourly intervals, but we applied it at a daily interval. The stress factor in the original GDH model was ignored, because we assumed that the plants were not under other stresses.

where FU₈ is the rate of forcing for Model F₈, T is the daily mean temperature (°C), T_L = 4, T_u = 25, and T_c = 36.

Analysis

We assessed the ability of each chilling model to represent long-term trends in the chilling conditions by calculating CA using each chilling model for each station for 1951–2018. CA was calculated as the sum of CU_i from 1 November in the previous year to 30 April. The trend of CA at each station was visualized as the slope of the linear regression of CA against year. We also calculated Pearson’s r between each pair of chilling models for each station to determine if the chilling models were interrelated.

We calculated HR and CA of spring events for each species, station, and year to determine if the chilling models are consistent with the physiological assumption that the reduction in chilling would increase HR (Fig. 1). HR was calculated as the sum of FU_i from 1 January to the date of onset of spring events using Model F₁, and the performances of the other forcing models (F₂–F₈) were also tested. We also compared 1 January with the other two starting dates of temperature accumulation (15 January and 1 February) to test any potential difference causing by the date when temperature accumulation begins.

CA was calculated as the sum of CU_i from 1 November in the previous year to the date of onset of spring events. We chose 1 November as the start date for CA because the endodormancy of temperate trees began around 1 November⁵². We only tested the linear relationship because the data were better fitted by the linear regression than the exponential model (Fig. 3), even though CA was linearly or nonlinearly negatively correlated with HR¹⁷. Pearson’s r between CA and HR for all records was calculated for each species, with a significantly negative Pearson’s r (p < 0.05) indicating that the chilling model met the physiological assumption. We also analyzed the relationship between CA and HR at stations with at least 15-year records to determine if the results were robust at the station level.

The above analysis is based on the E-OBS data set. Given the large grid size of the E-OBS data set (about 10 km) and the non-uniform cover of plants and temperature in the grid cells, we assessed the scale effect by comparing the HR and CA based on E-OBS data set with GHCN data set. We only retained the phenological station where a corresponding meteorological station (distance and altitude difference should be less than 5 km and 100 m, respectively) existed in the GHCN data set and compared the HR, CA, and their relationship.

We developed an empirical model based on the linear regression function between CA and HR (e.g., the linear fitted line in Fig. 3) to simulate the past and future spring phenological change. We calculated CA (from the previous 1 November to the current date) and heat accumulation (from 1 January to the current date) for each species in each year using a daily step. HR for the current date was calculated using the predefined linear regression function between HR and CA. The day when heat accumulation began to be larger than HR was determined as the date of onset of spring events. Compared to our empirical models, in current terrestrial biosphere models, the simulation of leaf onset is usually only based on GDDs⁵³ (e.g., Model F₁ and F₂ in this study), while only one model (ORCHIDEE) consider the effect of chilling⁵⁴ (Model C₁ in this study). Thus, at least currently, the CA–HR mechanism has not been well represented in ecosystem models. Our modeling efforts could timely provide the basis for a better representation of vegetation phenology in terrestrial biosphere models.

We predicted the annual spring phenological change (2019–2099) for each species using the above process-based models under RCP 4.5 and RCP 8.5. To produce a consistent phenological time series from past to future, we simulated the past spring phenological change from 1980 to 2018 by using the E-OBS v19.0eHOM data set, which was resampled to the same spatial resolution (0.5 × 0.5°) with the climatic data projected by HADGEM2-ES. First, we assessed whether different chilling models could reproduce spatial gradients of spring events across warm and cold regions in Europe. For each species, the simulated mean date (1980–2018) was correlated to the observed mean dates across grids with the observation data. Second, we assessed whether different chilling models could reproduce temporal trends of spring events in Europe. For each grid with at least 15-year observation data, the simulated and observed trends were estimated as the slope of the linear regression of spring phenology against year. At last, we compared the simulated trends estimated by each chilling model with the observed trends.

Reporting summary

Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Source: Ecology - nature.com