1. Introduction
Data from the Chinese Seventh National Forest Resources Inventory show that the poplar-forested area has reached 10.10 × 106 hm2 in China. The poplar plantation area was 7.57 × 106 hm2, accounting for 18.9% of the artificial arbor area, the highest in the world [1] . Populus simonii Carr., one of the main poplar tree species, has been widely planted as protection forest with windbreak, sand-fixa- tion, and solid-soil functions, and as timber forests in Northeast and Northwest China. Timber production is strictly controlled by plant phenology, which plays a major role in the carbon, water, and energy exchanges in the terrestrial ecosystem [2] . The time of leave sprouting of deciduous trees is controlled by the air temperature [3] ; therefore, climate warming has affected the deciduous plant leafing phenology significantly. It is important to understand the control mechanism and improve the accuracy of phenology simulation for P. simonii and to improve the simulation precision of the land-surface process models [4] in order to determine accurately the timber production of P. simonii forests.
Various plant phenology models have been developed, for example based on different control factors. These can be divided into the temperature model, hydrothermal model [5] , light and temperature model [6] , and the water model [7] [8] . Depending on the construction principles, such models can be regard as statistical [9] , mechanism [10] [11] , and the theory models [12] . It has been shown that the leafing phenology of deciduous tree species is correlated highly with the air temperature in high latitudes or the boreal climatic zone [10] . Considering the release process of dormant of deciduous plants, some models consider the chilling process a positive effect. Various studies have indicated the exponential relationship between the accumulation of the chilling temperature and the accumulation of the forcing temperature [13] .
Climate warming has led to leafing advances for P. simonii [14] . In view of climate change, it is important to determine the quantitative relationship between the leafing of P. simonii and climate warming in order to evaluate quantitatively the production potential of the plantations. In recent years, P. simonii leafing phenology has been studied continuously and widely, and the effects of climate warming on the phenology have been reported widely [14] [15] . However, relatively few studies have reported on the quantitative simulation of the poplar phenology.
The main objectives in this study are to: 1) quantitatively analyze the influence of climate warming on P. simonii leaf sprouting, 2) evaluate a P. simonii leaf- sprouting phenology model, and 3) quantitatively simulate the distribution characteristics of P. simonii leaf sprouting in 2015 in China. Such analyses are expected to provide a theoretical basis for improving the evaluation accuracy for P. simonii productivity in China.
2. Materials and Methods
2.1. Study Region
The study region encompassed 10 agricultural meteorological experimental stations in the Heilongjiang Province, located between 121˚W - 135˚W and 43˚N - 53˚N (Figure 1), with the northern most of China. The region has a temperate continental monsoon climate, with the mean annual temperature ranging from −5˚C to 5˚C, and the mean annual precipitation ranging from 400 to 650 mm.
2.2. Phenology Data
Plant phenology observation by the China Meteorological Administration started in agricultural experimental stations in 1980. The observers recorded the dates of the main phenophase events of plants (including leaf unfolding). The date of leaf unfolding was expressed by day of the year (DOY). During this study, ten datasets were collected, from which the data on the leaf unfolding of the individual trees were extracted. The leaf unfolding time was defined as the first flat leaf appearing from the buds of trees with simple leaves. The final available dataset has 186 site-years in 10 sites (Table 1).
Figure 1. The location of stations in the study area.
Table 1. Observation sites and data sets in Heilongjiang Province.
2.3. Meteorological Data
The daily average air temperature was used to calibrate and validate the plant phenology models. The air temperature data were collected from 10 agricultural meteorological experimental stations where leaf unfolding was being observed. To evaluate P. simonii leaf sprouting, we obtained the daily average temperatures, recorded by the 832 observation stations in 2015 in China, from the CIMSS (China Integrated Meteorological Information Service System).
2.4. Phenology Model
In this study, we selected two models that simply consider thermal forcing and three models that combine thermal forcing with chilling.
2.4.1. Two Thermal Forcing Models
Accumulated temperature model (AT)
The AT regards the thermal forcing process as a linear increase.
(1)
(2)
where y is the data on leaf unfolding, xt is the daily average temperature, Sf is the state of forcing, is the forcing rate, Tb is the base temperature, and t0 is the start date of the thermal forcing calculation.
Spring warming model (SW)
In the SW, the thermal forcing accumulation is expressed as a logistic function [16] .
(3)
(4)
where va, vb, and vc are optimized to the least root mean square error of the model.
2.4.2. Three Models Combining Thermal Forcing with Chilling
Sequential model (SM)
In the SM, thermal forcing accumulation starts after the accumulated chilling unit has reached a critical threshold [17] .
(5)
(6)
(7)
(8)
where t0 is the start date of the chilling calculation, t1 is the start date of the thermal forcing calculation, and To, Tl, and Th are the optimal, lowest, and highest temperature of the chilling rate, respectively.
Parallel model (PM)
In this model, the chilling and forcing accumulations are computed simultaneously. In the accumulation process, thermal forcing is adjusted with a chilling unit [18] .
(9)
(10)
(11)
(12)
(13)
where Km is the minimum potential of un chilled buds to respond to the forcing temperature, and is the critical value of the state of chilling.
Alternating model (AM)
In this model, chilling and forcing in turn accumulate from an initial start date, relative to a single base temperature [19] .
(14)
(15)
(16)
(17)
(18)
where a, b, c, and d are optimized to the least root mean square error of the model.
In order to test whether reducing the model parameters would result in poor model precision, two methods were considered for the AT, SW, SM, and PM models. One method was to set fixed constants for some parameters, and in the other method, all parameters were obtained optimally.
Half of the observation data were used to calibrate the models, whereas the other half of the data set was used to evaluate the performance of the models. The parameters of each model were determined and optimized by using the simulated annealing method. To select the best model, the corrected Akaike information criterion (AICc) was analyzed [20] .
where l is maximized log-likelihood, K is the number of estimable parameters. n is the sample size.
The validation of the models was carried out with the root mean square error (RMSE), mean bias error (MBE), and correlation coefficient (CC).
3. Results
3.1. Leaf Unfolding and Air Temperature
A significantly negative correlation was indicated between P. simonii leaf unfolding and the air temperature. A rising air temperature during the early stage of leaf unfolding (e.g., average air temperature in April) significantly advanced the P. simonii leaf unfolding. A rise of 1˚C advanced leaf unfolding by 2.47 days (Figure 2). This result is consistent with the results of various other reports; for example, an amplitude variation of 1 - 4 days/˚C indicated in Inner Mongolia [21] . A difference of 1˚C from the mean temperature led to a change of 2.8 days in the spring phenology date for deciduous forests, and a change of 4.8 days for evergreen forests [22] . The temperature sensitivity of leaf unfolding of P. simonii is remarkable. A significant linear correlation between the air temperature and leaf unfolding suggests that climate warming has substantial potential to affect leaf unfolding of P. simonii.
3.2. Model Calibration and Validation
A summary of the best fit for all the species-specific models and the relative performance of the models are provided in Table 2. As regards the simplicity and accuracy of the models, the AICc values of SW2 are the smallest and the SW2 model is therefore regarded as superior for simulating the leaf unfolding of P.
Figure 2. Relationship between leaf unfolding of Populus simonii and the average air temperature in April.
simonii. The AICc value increased with the parameter numbers. As regards the forcing temperature model, using a fixed value for some parameters could be a better option. For example, in the AT and SW models, a lower AICc value was obtained by using the constant value, referred to in the documents, for some parameters. In respect of the models concluding chilling and forcing temperature, the lowest AICc value appears in the SM2 model. Accordingly, increasing the parameters could improve the accuracy of the model, but would also increase the complexity. As January is the coldest month in most regions of China, the T0 is set to January 1, and calculations of chilling and forcing started on this date (Table 2). In respect of the independent data of 97 site-years for the calibration of all the models, the absolute error of 91.6% of the calibration data is within five days. The absolute error of 76.7% of the calibration data is less than three days. Therefore, the simulation accuracy of the model that considers only the temperature to simulate the leaf unfolding of P. simonii is acceptable (Figure 3).
In comparing the RMSE, CC, and MBE of the nine models, the precision of the PM1 model was found the highest; however, the complexity of the model
Table 2. Parameter values of all models for leaf unfolding of Populus simonii.
Figure 3. Distribution of average value of absolute error for nine models in predicting leaf unfolding of Populus simonii, using independent data.
was ignored. The relevant RMSE, CC, and MBE average values were 2.88, 0.78, and 2.5 days, respectively. The maximum value of CC was 0.91 in the PM1 and the minimum value was 0.49 in the SW1 models. The average CC value range in all the models was 0.71 to 0.79. The minimum value of 1.73 days for RMSE appeared in PM1, whereas 4.96 days was the maximum value. The range of the average RMSE values was between 2.88 days in PM1 and 3.70 days in SM1. The minimum MBE value, only 1.50 days, simultaneously appeared in PM1, PM2, SM2, and AT1. The average MBE range was 1.5 to 4.0 days (Figure 4). The leaf unfolding of P. simonii can be simulated by using the PM1 model. The determine coefficient of PM1 reaches 70% when leaf unfolding is predicted with independent data. Accordingly, 70% of the variation in the leaf unfolding of P. simonii can be explained with the PM1 model (Figure 5).
3.3. Model Application
In 2015, we simulated the spatial distribution pattern of P. simonii leaf unfolding by employing nine models (Figure 6). Leaf unfolding occurs from DOY (day of the year) 100 to 120 in the main distribution region of Populus in China, namely, the SongNen, SanJiang, and SongLiao plains and the Inner Mongolia plateau
Figure 4. Comparison of correlation coefficient, root mean square error (RMSE), and mean bias error (MBE) of nine models for P. simonii leaf unfolding. The top and bottom of the box are maximum and average values, and the lowest point of the line indicates the lowest value.
Figure 5. Comparison of predicted and observed values for leaf unfolding of Populus simonii in PM1, using independent data.
[23] . In the Haihe River plain, the coast of the Bohai area, and the loess plateau, the range of leaf unfolding is from DOY 80 to 100. In the HuangHuai and the WeiHe river basins, it is between DOY 60 and 80, and in the JiangHuai River basin, between DOY 40 and 60. In the northwest of China, P. simonii leaf unfolding occurs between DOY 80 and 140. In the south of XinJiang, leaf unfolding occurs the earliest, from DOY 80 to 100. In the HeXi corridor, the north of XinJiang, and the Ili River valley of XinJiang, leaf unfolding occurs between DOY 100 and 120. Leaf unfolding takes place the latest in the QingHai plateau, namely, from DOY 120 to 140.
In addition to the simulation by the SM1 and SM2 models, the P. simonii pattern of leaf unfolding simulated by the other models was consistent. The earliest unfolding occurs in the south of Yunnan Guangxi and Guangdong and is less than 20 DOY. However, in the SM1 and SM2 models, the earliest leaf unfolding appears in the Sichuan basin and Jiangxi, Fujian, and others, within DOY 20 to 40. However, in the south of Yunnan, Guangxi, and Guangdong, leaf unfolding of P. simonii cannot be successful because the chilling threshold cannot be reached (Figure 5).
4. Discussion
Changes in phenology significantly affect the carbon balance of terrestrial ecosystems [22] . Plant leafing is more sensitive to climate change than plant yellowing [24] ; however, recent literature indicates that temperate deciduous forest species have shown a wide range of phonological responses to climate changes over the last half-century [19] . As P. simonii is an important forestation tree species, it is planted widely in China. Forestation success is dependent on climate adaptability and sensitivity. Populus simonii is a deciduous tree species, and most of the research results indicate that leaf unfolding is significantly associated with air temperature. Therefore, climate warming could improve the productivity of P. simonii forests.
Phenological models have been used widely to simulate plant leaf unfolding [11] [25] . An interesting result of the modeling is that, generally, the same level of accuracy was demonstrated by all the models that considered only the effect of temperature on leaf unfolding. Increasing the number of parameters does not improve the model accuracy significantly. On the contrary, such increase could result in model complications and could increase the AICc value. In this instance, the results obtained by the simple model were superior, as has been indicated in many other studies [19] [24] . The AT and SW methodologies maintain the simplicity of the model. In addition, these models have the ability to explain most of the influence of the temperature on the phenology. As regards the plant growth mechanism, employing the constant mentioned in the literature for some parameters not only simplified these models but also reduced the AICc value. Simulation precision can be maintained with this model; moreover, the model is not only physiologically significant but also simple and suitable for application in the simulation of plant leaf unfolding.
Populus simonii leaf unfolding occurs when a critical state of forcing temperature is reached in the AT and SW models. The chilling and forcing unit of temperature was considered in SM, PM, and AM, with these models being fitted and their parameters estimated. Obviously, the parameters of these models have physical dimensions that in principle could be measured directly instead of being estimated by a fitting process. However, this is rarely possible [26] . When the chilling of temperature is considered, if only by considering its auxiliary function, e.g. in PM and AM, the plant leaf unfolding under climate warming is consistent with the AT and SM models. Conversely, when considering the crucial role of chilling against the background of climate warming, predicting the value of P. simonii leaf unfolding by SM is clearly different from that with AT, SW, PM, and AM. In the southern parts of Yunnan, Guangdong, and Guangxi provinces, leaf unfolding could not take place because the critical threshold chilling value could not be reached. Accordingly, the SM model suggests that P. simonii is not suitable for planting in warmer regions. Chilling will often be insufficient to break bud dormancy may become a cruicial factor limiting leafing [27] .
Suppose Sc = 0 in warmer areas, PM is changed into SW multiplied by Km and AM is changed into AT1. Therefore, the predicted value by PM2 is less than that by SW2 for the difference. The predicted value by AM is approximately the same as that by SW2 because the in AT1 is slightly equal to the a value plus the c value in the AM parameter. When Sc = 0 in SM, the calculation of forcing could not start. According to the AT and SW, P. simonii will change from deciduous trees into an evergreen tree species when it is hot enough. As indicated by SM, P. simonii is not suitable for growing in tropical regions, and the characteristic of shedding leaves should be maintained. Regardless of how it came about, the PM and AM models have lost the significance of the chilling effect on leaf unfolding in the tropical regions. Carter et al. (2017) [28] elaborated that greater amounts of warming units accumulated prior to leaf emergence during the extreme warm year relative to non-extreme years. Which is consistent with the principle of PM and AM. Special experimental validation is required to determine which of SW and AT is reasonable or SM is correct. However, in terms of the distribution area of poplar trees in China, the southern parts of the country is not a typical growing region, and poplar trees are not suitable for planting in the south of Guangdong and Guangxi provinces [23] . Kin view of the above, the SM model is more reasonable for simulating P. simonii leaf unfolding.
5. Conclusion
Analysis of the relationship between P. simonii leaf unfolding and the air temperature has indicated that the leaf unfolding is sensitive to the air temperature. Climate warming is expected to advance the leaf unfolding of P. simonii, lengthen the growing season, and increase the productivity. Model analysis shows that a simple model is superior to a complicated model. Without due consideration being given to the influence of the parameters, the accuracy of the PM model was indicated superior. The phenological model based on air temperature could be superior at simulating P. simonii leaf unfolding. If the chilling threshold needs to be met (e.g., the SM model), the southern areas of Yunnan, Guangdong, and Guangxi would be unsuitable; however, other areas can be used to plant P. simonii.
Acknowledgements
The authors gratefully acknowledge the CIMISS for supplying the meteorology data. This research was supported by grants from the National Basic Research Program of China (973 Program) (2012CB416906).