The study area covered 30-90° N, with a spatial resolution of 0.5° × 0.5° . Information on vegetation cover was obtained from the land surface data of the Community Land Model, version 4.5 (CLM4.5). We focused on the land area dominated by plant functional types (PFTs) with summer-green phenology, including broadleaf or needleleaf deciduous tree (temperate or boreal), broadleaf deciduous shrub (temperate or boreal), and C3 grass (arctic or non-arctic). The number of samples was 12 079, which was about 60.5% of the natural vegetation coverage.

The data on leaf onset used in this research were derived from a long-term Global Land Surface Satellite (GLASS) dataset (Liang et al., 2013). The GLASS products include leaf area index (LAI), broadband albedo, broadband emissivity, downward shortwave radiation, and photosynthetically active radiation. The LAI product is derived from Advanced Very High Resolution Radiometer (AVHRR) and Moderate Resolution Imaging Spectroradiometer (MODIS) ranging from 1982 to 2012. A modified piecewise logistic fitting method was utilized to derive leaf phenology from LAI (Zhang et al., 2003), and zero point of the forth order derivative of the fitted logistic models was used to determine the four transition dates, including green-up, maturity, senescence, and dormancy, for each phenology cycle (in these phenological data, there were up to two phenological cycles for each year).

Since most natural vegetation has only one phenology cycle, only green-up in the first phenology cycle from 1982 to 2000 was used as the leaf onset date in this research. We also calculated its multi-year average to obtain the climatology of leaf onset, which we named Phe8200. The spatial distribution of Phe8200 is shown in Fig. 1a.

Figure 1

	Figure Option ViewDownloadNew Window
	Figure 1 (a) Spatial distribution of leaf onset dates for summer-green plants in Global Land Surface Satellite (GLASS) dataset (b) Fitted Tb; (c) Simulated leaf onset date in Experiment 1 (EX1); (d) The difference between model and observation.

We utilized a global land surface air temperature (SAT) dataset produced by Wang and Zeng (2013), which was developed from four reanalysis products: Modern-Era Retrospective Analysis for Research and Applications (MERRA), for 1979-2009; 40-yr European Centre for Medium-Range Weather Forecasts (ECMWF) Reanalysis (ERA-40), for 1958-2001; ECMWF Interim Reanalysis (ERA-Interim), for 1979-2009; and National Centers for Environmental Prediction-National Center for Atmospheric Research (NCEP-NCAR) reanalysis, for 1948-2009. Climatic Research Unit (CRU) monthly-mean maximum and minimum data were applied for bias correction. Wang and Zeng (2013) showed that this dataset performed well, in comparison with in-situ hourly measurements, over six sites and with a regional daily SAT dataset over Europe; hence, it could provide better estimation of daily mean temperature.

The spatial resolution of the original dataset was 0.5° × 0.5° , and the temporal resolution was one hour. We aggregated the data into daily temperature. The climatology from 1982 to 2000 was calculated and named Tem8200.

The TM is a mechanistic model developed by Cannell and Smith (1983). It is based on known or assumed cause-effect relationships between biological processes and environmental driving factors (Chuine et al., 2003). The TM is described by GDD, which is defined as

(1)

where i is the day of year, T_i is the temperature of the ith day, and T_b is the baseline temperature. In some applications, e.g., CLM-DGVM (dynamic global vegetation model) 3.0 (Levis et al., 2004), GDD is reset to 0 when T_i< T_b, an approach used to avoid the influence of cold injury. Plants start leaf shooting until GDD exceeds a preset critical value (V_c).

Due to the lack of large-scale phenological observations during the early development of the TM, T_b was preset as a global constant or a parameter depending only on species, biomes or PFTs. Although in different applications T_b may have different preset values, they have no spatial or temporal variation.

The TM has been introduced into DGVM; for instance, CLM-DGVM 3.0 (Levis et al., 2004), CLM-CNDV 4.5 (Bonan et al., 2013), SEIB (spatially explicit individual-based dynamic global vegetation model)-DGVM (Sato et al., 2007), and IAP-DGVM (Zeng et al., 2014). Such simplification, however, has caused large errors when applied at global or large scales (Richardson et al., 2012). For example, as GLASS phenology data were analyzed in the present study, the temperature on leaf onset dates in northern Canada, Russia, and Europe was below 0° C, much lower than in temperate regions. Hence, if T_b was set to 0° C for the whole study area, it was not surprising that the model predicted a systematic delay in leaf onset, increasing with latitude (see section 4.2 and Fig. 2b for more information).

Figure 2

	Figure Option ViewDownloadNew Window
	Figure 2 (a) The zonal mean of model Tb, shaded area is between the 25th and 75th percentile of model Tb; (b) The error of the zonal mean in the control Experiment (EX0) (dash line) and EX1 (solid line).

Regarding the dramatic geographical range of the distribution of summer-green vegetation, it is reasonable to assume that both T_b and V_c vary with location. Actually, the value of T_b has been assigned within a wide range in previous studies. For instance, Hunter and Lechowicz (1992) selected -5° C, 5° C, and 15° C for 26 species of native hardwood trees in Ohio, and Botta et al. (2000) chose -5° C, 0° C, and 5° C for different boreal/temperate biomes. In this study, we developed an MTM in which T_b was estimated from global phenological observations and local climate conditions.

Technically speaking, the spatial distribution of T_b can be inferred from satellite-observed global phenological data (and hence is called “ observed” T_b). By applying GLASS phenological data and the temperature data of Wang and Zeng (2013), and setting V_c to be 100, the value of T_b can be calculated using Eq. (1), by employing the method of dichotomies to search between the temperatures on 1 January and leaf onset date, and controlling the error of GDD below 0.5. However, such a method might not be of practical use in climate models and DGVMs, because the errors in observations (especially in phenological data) may be amplified in plant phenology. Furthermore, it could not produce the historical variation of T_b, and hence may be inapplicable to studies of paleo- or future climate.

An alternative approach is to relate T_b with climate factors that are available in climate models. Our study indicated that there was a significant linear relationship between T_b and spring temperature g(T) (the average temperature from March to May), and the regression function can be written as:

, (2)

where α and β are fitted parameters, and ε is the residual error.

To obtain a robust result, the climatological data, i.e., Phe8200 and Tem8200, were used for the regression. We obtained α =0.72 and β =-4.50, while the residual standard error of the fitting was 0.0027, and the coefficient of determination R² was 0.84. In the following analysis, we denote T_b simulated by Eq. (2) as T_{b_mod}. The average of T_{b_mod} throughout the study area was -8.1° C, ranging between -19.9° C and 11.5° C, with standard deviation of 5.7° C— roughly the same as the “ observed T_b” . Figure 1b shows the geographical structure ofT_{b_mod}, and Fig. 2a shows that its zonal average decreased with latitude.

Three experiments were performed for model verification. In Experiment 1 (EX1) and Experiment 2 (EX2), T_{b_mod} was applied. While EX1 used Tem8200 to predict the climatology of leaf onset date and was compared with Phe8200, EX2 used temperature in each year from 1982 to 2000 to show the interannual variability of predicted phenology. The control experiment, EX0, had T_{b_mod} ? 0, which represented a standard TM scheme (the same as in IAP-DGVM and CLM-DGVM 3.0), and also Tem8200 to predict a climatology for comparison with EX1. In all experiments, the GDD requirement of V_c was 100.

For each experiment, the mean error (ME) and mean absolute error (MAE) of the simulated leaf onset date were calculated as follows:

(3)

where n is the number of samples, x_i is prediction, and y_i is observation.

In order to evaluate the interannual variability, the error in the predicted climatology of leaf onset, i.e., EX1, should be removed from EX2. Thus, the anomaly of model error (denoted as EA) was calculated as follows:

EA_j=ME_j-ME_cli, (4)

where j is the index of year, ME_j is the mean model error in year j in EX2, and ME_cli is the mean model error for Phe8200 in EX1. To further analyze EA, the climate interannual variability, i.e., the anomaly of spring temperature (denoted as STA), was also introduced as:

SAT_j=ST_j-ST_cli , (5)

where ST_j is the mean spring temperature in year j, and ST_cli is the mean spring temperature for the climatology of 1982-2000.

Figure 1c shows the spatial pattern of leaf onset simulated in EX1. The ME of EX1 was 1.11, and the MAE was 6.8 (days). The 5th, 25th, 75th, and 95th percentile values of the model error were -12.4, -3.1, 6.9, and 13.7 (days), respectively. The MTM captured the overall spatial pattern of leaf onset as in the observation (Fig. 1a), matched well in Russia, but the predicted leaf onset dates were later than observed in most of North America and Nordic Europe, and earlier in most of northern China (Fig. 1d). Figure 2 shows that the error of the zonal mean in EX1 was close to zero between 40° N and 75° N.

On the contrary, much larger errors existed in EX0, in which the ME was 31 and the MAE was 32.4 (days). The bias of the zonal mean was larger than 30 days between 50° N and 80° N (Fig. 2b), consistent with the region where T_b was much lower than 0° C (Fig. 2a). This shows the considerable advantage of a spatially variable T_b over a constant T_b in predicting leaf onset date.

Figure 3 shows that model error varied for different years in EX2. For all 19 years, the average ME of EX2 was 1.87 and the average MAE was 10.6. For each year, the ME was between -3.1 and 7.16 and MAE was between 8.23 and 12.9. Sixty-one percent of model error fell in to the interval between -10 and 10 days, and 87.2% was between -20 and 20 days.

Figure 3

	Figure Option ViewDownloadNew Window
	Figure 3 Box plot of model error from 1982 to 2000 in Experiment 2 (EX2). The bottom and top of the box are the first and third quartiles, the band inside the box is the median, and the ends of the whiskers are the 5th and 95th percentiles.

Figure 4 shows that the anomaly of model error (EA) and the corresponding anomaly of spring temperature (STA) were to a certain extent correlated. The correlation coefficient between EA and STA was -0.623, significant at 0.01 level. Generally, EA was positive when spring temperature was cooler, and negative when spring temperature was warmer than climatology, i.e., the model overestimated the phenological shift induced by the temperature anomaly. This implies that, while plants may adapt to the local climate during their evolutionary history, they are less capable of adjusting to interannual climate variability or disturbance.

Figure 4

	Figure Option ViewDownloadNew Window
	Figure 4 Anomaly of model error (denoted as EA, dash line) and anomaly of spring temperature (denoted as STA, solid line) from 1982 to 2000 in EX2.