A Temperature-Precipitation Based Leafing Model and Its Application in Northeast China

Plant phenology models, especially leafing models, play critical roles in evaluating the impact of climate change on the primary production of temperate plants. Existing models based on temperature alone could not accurately simulate plant leafing in arid and semi-arid regions. The objective of the present study was to test the suitability of the existing temperature-based leafing models in arid and semi-arid regions, and to develop a temperature-precipitation based leafing model (TP), based on the long-term (i.e., 12–27 years) ground leafing observation data and meteorological data in Northeast China. The better simulation of leafing for all the plant species in Northeast China was given by TP with the fixed starting date (TPn) than with the parameterized starting date (TPm), which gave the smallest average root mean square error (RMSE) of 4.21 days. Tree leafing models were validated with independent data, and the coefficient of determination (R 2) was greater than 0.60 in 75% of the estimates by TP and the spring warming model (SW) with the fixed starting date. The average RMSE of herb leafing simulated by TPn was 5.03 days, much lower than other models (>9.51 days), while the average R 2 of TPn and TPm were 0.68 and 0.57, respectively, much higher than the other models (<0.22). It indicates that TPn is a universal model and more suitable for simulating leafing of trees and herbs than the prior models. Furthermore, water is an important factor determining herb leafing in arid and semi-arid temperate regions.


Introduction
Plant leafing and yellowing stages both play critical roles in accurately estimating carbon and water flux exchanges between the land and atmosphere [1,2] and the changes in land surface characteristics [3,4]. Moreover, plant leafing is more sensitive to climate change than plant yellowing [5,6]. The ability to precisely predict plant leafing is crucial to modeling the impacts of climate change on plant primary productivity.
Currently, there are many phenological models to predict the changes in plant spring phenology, including bud, leafing, and flowering stages [7][8][9][10][11][12]. The simplest spring phenological models consider only temperature, as exemplified by the cumulated temperature model [13,14]. More complex models based on the intensive study of plant physiology incorporate the dual roles of temperature (i.e., chilling and forcing); such models include the sequential model [9,12,15], the parallel model [9,10,16], and the alternating model [7,17]. The most complex phenology models consider the impact of day length as well as temperature, for example, the light and temperature phenology model [10]. Overall, the current temperature-based spring phenology models are mainly used for tree species in temperate humid and semihumid areas; the efficacy of these models for simulating spring phenology for trees and herbs in temperate arid and semi-arid regions is untested.
Precipitation is also a key determinant of plant leafing, especially in semi-arid and arid area, however, the impact of precipitation on spring phenology has seldom been considered in plant leafing models. Yuan et al. [18] initially explored the leafing responses of dominant herbs (Leymus chinensis and Stipa grandis) to soil moisture in Inner Mongolia and developed a leafing model for L. chinensis and S. grandis based on the effects of temperature and soil moisture. However, the model was established only for two herbs and soil moisture is rarely measured. This model, moreover, was not validated by other external data [18].
Northeast China, located at high latitude in the northern hemisphere, is highly sensitive to climate change and has experienced the increase in air temperature twice that of the global average [19]. Furthermore, precipitation varies significantly in Northeast China. There exists a latitude-based thermal gradient, including warm (south), moderate, and cold (north) temperate regions. From east to west, there are various precipitation zones, including humid, semi-humid and semi-arid regions. Thus, precise modeling of tree and herb leafing in Northeast China is rather complicated, yet it is crucial to simulate the environmental consequences of climate change in this region.
The present study is based on long-term (i.e., 12-27 years) ground observations of leafing and simple meteorological data (i.e., daily mean temperature and precipitation). The leafing response of 13 plants dominant in Northeast China (including eight trees and five herbs) to hydrothermal factors would be evaluated and simulated. Our main objectives were to determine: (1) whether the existing temperature-based models accurately simulate tree and herb leafing in arid and semi-arid regions; and (2) whether a temperature-precipitation based leafing model can more precisely simulate plant leafing under different water and heat conditions in arid and semi-arid regions.

Study site and plant species
Northeast China consists of Heilongjiang, Jilin, Liaoning provinces and four leagues of Inner Mongolia (Fig. 1). The region has a continental monsoon climate. Mean annual air temperature is 4.5uC, with an average temperature of 218.2uC in the coldest month (January), and an average of 22.4uC in the warmest month (July). Annual average precipitation is 514 mm, 77% of which falls from May to August. The region has a minimum annual precipitation of 245 mm in the west, a typical semi-arid area, and a maximum annual precipitation of 1079 mm in the east.
In the east of Northeast China, Great Xing'an Mountains is the cold temperate coniferous forest zone, dominated by typical tree species Larix dahurica. Changbai mountains, typical temperate coniferous and broadleaved mixed forest zone, dominated by Pinus koraiensis. The west of Northeast China is the typical zone of temperate steppe with herbs such as Leymus chinensis, Stipa krylovii, Agropyron cristatum and S. baicalensis [20]. In Northeast China, Populus simonii is the main afforestation tree species, and Salix matsudana, Armeniaca vulgaris, and Ulmus pumila are the common garden species. These plant species are strictly controlled by hydrothermal conditions. Thus, thirteen dominant plant species are selected in the present study, including eight trees (Salix matsudana, Armeniaca vulgaris, Ulmus pumila, Populus simonii, Syringa oblate, Pinus koraiensis, Larix dahurica, and Picea koraiensis) and five herbs (Leymus chinensis, Stipa krylovii, S. baicalensis, Elymus dahuricus, and Agropyron cristatum) ( Table 1).

Phenological data collection
The leafing data of dominant plants were collected from nine Agricultural Meteorological Experiment Stations, China Meteorological Administration, located in Northeast China. Plant leafing status was observed daily; plants were considered to have leafed if: (1) the first flat leaf had appeared from the buds of trees with simple leaves; (2) young leaves had emerged from the leaf sheaths of conifers; (3) one or two leaflets of compound leaves had unfolded; or (4) old exposed leaves of over-wintering herbs had turned from yellow to green, and the first leaf of herbs had emerged above the ground [21].

Meteorological data collection
Meteorological data, including daily mean temperature, daily precipitation, daily minimum temperature and relative humidity, were collected from nine Agriculture Meteorological Stations where plant leafing was observed.
Water plays a critical role in regulating plant phenology in arid and semi-arid areas [18,24,25], and has been included in some phenology models. Examples include cumulative precipitation in the current year [26], vapor pressure deficit in the growth season index (GSI) [25], and soil moisture [18]. However, the combination of hydrological and thermal factors has not been considered. Usually, when both hydrological and thermal conditions reach certain thresholds, plants begin leafing. Previous studies of plant leafing phenology [18,24] have identified the accumulated precipitation in the previous year and current year as an important hydrological factor and the accumulated temperature in the current year as an important thermal factor affecting plant leafing. Thus, a new plant leafing model (so-called TP) based on the effects of both temperature and precipitation could be expressed as: where k 1 , k 2 , P crit , T b , and F * are parameters obtained through optimization. k 1 , k 2 are the efficiency of precipitation in the previous and current (prior to leafing) year in affecting leafing in the current year. P b is the annual precipitation in the previous year. y is the day of plant leafing in the current year. R i is the daily precipitation in the current year (mm). P crit is the water threshold (mm). T i is the average daily temperature (uC) in the current year. T b is the base temperature (uC). F * is the temperature sum critical threshold (uC). Models should always be validated with independent data not used to construct the model itself [27]. In this study, the phenology model parameters were estimated using leafing data from odd-numbered years (12 years), and the simulation accuracy was tested with the independent even-year data (12 years).

Parameter estimation
Phenological data were converted to Julian day (DOY). Model parameters were estimated using the least root mean square error (RMSE) method: where d i (x) is the predicted date of plant leafing in the ith year and d iobs is the observed value of plant leafing in the ith year. The model was evaluated by F-tests. The optimized parameters of the model were determined using the simulated annealing method [8].
The average RMSE, coefficient of determination (R 2 ) were used to validate the model in the present study.
Using odd-year data, the parameter sets and three statistic variables (RMSE, R 2 and F) were given in Tables 3, 4, 5 for five models (i.e., SW, SM, PM, AM, and TP) ( Table 2). Both fixed and parameterized starting date were considered for SW, TP and AM, i.e., the fixed starting date was set on 1 January for SW and TP and 1 September for AM. The starting date (t 0 ), the minimum (T low ) and maximum (T high ) values of chill temperature, and the parameters values (va, vb, vc) of response curves for forcing temperature were fixed (e.g., t 0 , T low , T high , va, vb and vc were set to September 1,23.4,10.4,28.4,20.185 and 218.4, respectively) and parameterized for SM and PM, based on the reference and parameterization in the present study.

Model fitting
The model comparison would be evaluated by RMSE, coefficients of determination (R 2 ) and sum of residual squares (e.g., F-test). The average RMSE for all models fitting ranged from 2.00 to 4.18. The average RMSEs of PM with all parameters estimated (PMm) was the smallest, and the average RMSEs of SM with all parameters estimated (SMm) was close behind (Tables 3, 4, 5). However, the result of models fitting did not represent an accurate model test. The model validation would be assessed using independent data.
The average RMSEs of temperature -precipitation based leafing model (TP) with fixed and parameterized starting date were 3.22 and 2.69 days, respectively. k 1 (average = 0.05) of TP was smaller than k 2 (average = 0.45) ( Table 3). The hydrological conditions for plant leafing varied with available water and climate in different regions. Thus, model coefficients comprehensively reflected the hydrological requirements for plant leafing. The mean hydrological requirement for woody plants was more than herbs in Northeast China.

Model validation
We validated these models using plant leafing data of evennumbered years in Northeast China. TP with fixed starting date (TPn) had the smallest average RMSE (4.21 days). The RMSEs of TPn for all plant species ranged from 2.22 to 6.44 days. The average RMSE of TP with parameterized starting date (TPm) was 5.12 days. The average RMSE of SW with fixed starting date (SWn) was smaller than SW with parameterized starting date (SWm) (i.e., 6.09 days ,6.95 days). The average RMSEs of the other models were more than 10 days (Fig. 2, Fig. 3). For tree leafing, R 2 was greater than 0.60 for 75% of the estimates by SWn and TPn, but less than 50% for the others. Both SWn and TPn had RMSE less than five days in 87.5% of the estimates, while they met this criterion less than 75% of the time (Table 6). Based on R 2 and RMSE, SWn and TPn were the best models for simulating tree leafing in Northeast China.
When we validated the models using even-year data for herb leafing, TPn model yielded RMSE of 4-6 days (average 5.03 days) and R 2 of 0.4-0.82 (average 0.675); followed by TPm with average Table 2. Equations of the phenology models compared in the present study.

Model type Equation of models
Photo §39600s iGSI~iT min |iVPD|iPhoto y: date of leafing; x t : daily mean air temperature in degrees Celsius; R f (x t ): forcing rate function; R c (x t ): chilling rate function; S f : state of forcing; S c : state of chilling; K m : minimum potential of unchilled buds to respond to forcing temperature; C*: critical value of state of chilling for the transition from rest to quiescence; F*: temperature sum critical threshold; t 0 : starting day of the heat sum calculation or date of onset of rest; t 1 : date of onset of quiescence; T b : base temperature; T o : optimal temperature of the rate of chilling; T low : the lowest temperature of the rate of chilling; T high : the highest temperature of the rate of chilling; a, b, va, vb and vc: constants; iT min : daily indicator for minimum temperature; T min : observed daily minimum temperature in degrees Celsius; iVPD: daily indicator for vapor pressure deficit; VPD: observed daily vapor pressure deficit in Pascals; iPhoto: daily photoperiod indicator; Photo: daily photoperiod in seconds; iGSI: daily growing season index. doi:10.1371/journal.pone.0033192.t002 RMSE and R 2 of 6.21 days and 0.572, respectively. For the other models, the minimum and average RMSE were more than 4.50 and 9.51 days, respectively (Fig. 2), and the maximum and average R 2 were less than 0.62 and 0.22, respectively (Table 6). Consequently, TPn was the best model for simulating herb leafing in the present study.

Discussion
Previous works have suggested that R 2 [28] or RMSE [29] are better than other criteria for comparing different nonlinear models. Generally, the smaller the RMSE value, the larger the R 2 value. However, the opposite case (i.e., a smaller RMSE was associated with a smaller R 2 ) occasionally appeared in the model calibration of SW and TP for leafing of Salix matsudana (Table 5). SWn could simulate the leafing of Salix matsudana precisely (R 2 = 0.67, P,0.001) ( Table 6). This result supports the opinion that RMSE is a more reliable measure of fit than R 2 for nonlinear regression [29]. The opposite case may be due to the noise of parameterized data for S. matsudana. Furthermore, the RMSE of plant leafing in Northeast China given by TPn ranged from 2.22 to 6.44 days, and the predicted DOY values for leafing were close to the observed values (Fig. 4). It indicated that TP can precisely simulate the leafing stages of both woody plants and herbs in Northeast China.
The model accuracy for tree leafing at high latitudes could not necessarily be improved with more complex models, consistent with the result of Hannerz [30]. In the present study, the most complicated model, PM, could give better simulation with the least precise, while a much simpler model, SW, with more accurate. SM, AM, and PM consider all the factors in the chilling process, whereas SW and TP do not, i.e., the former models include more information on temperature change through time during parameterization process. However, these three models (SM, AM, and Table 3. Parameter values of spring warming model (SW), temperature-precipitation based leafing model (TP), and alternating model (AM) for plant leafing in Northeast China.     PM) performed poorly when validated using actual observed data at different times, as seen with Ulmus pumila (Table 6). SW and TP without chilling requirement considering the starting date can reflect the effect of temperature beyond base temperature at different periods on plant leafing. Models with both fixed starting date (SWn and TPn) on 1 January and parameterized starting day (SWm and TPm) were analyzed in the present study. We found that after the validation by independent data, SWn (average RMSE of 6.05 days) could give better leafing simulation than SWm (average RMSE of 6.95 days), and TPn (average RMSE of 3.59 days) was better than TPm (average RMSE of 4.31 days). This result indicated that the temperature beyond base temperature at most periods could affect plant leafing. Generally, the models with parameterized starting date can consider only the temperature beyond base temperature after starting date, but the temperature before starting date might have important effect on plant leafing. In Northeast China, the coldest month is January with mean air temperature of 218.2uC. Compared with SWm and TPm, SWn and TPn can give better simulation in Northeast China because the effective temperature during longer period can be considered. According to the observation data, tree leafing tended to advance with the extent of 0.23 days yr 21 during 1980-2005, and was significantly negatively correlated with temperature in February, March and April. Furthermore, the effect of average temperature in April and February on plant leafing was the largest (2.35 days uC 21 ) and the smallest (1.18 days uC 21 ), respectively [31]. Therefore, regarding to climate warming, the main driving factor of plant leafing should be temperature instead of light [32], and SWn and TPn could give better simulation of plant leafing.  The simulation precision of phenology models with the consideration of chilling is critically influenced by the starting date. In Europe, phenology models are set to start on 1 September [23] or 1 November [7,10]. In this study, SMn, PMn and AMn were set to start on 1 September and the starting dates of SMm, PMm and AMm were parameterized using odd-year data. There are larger RMSEs for those models from independent data, and the reasons are attributed to the starting date: (1) the starting date is set early, resulting in untimely meeting the chilling and forcing thermal requirements. For example, the RMSE of SM with starting date on 1 September and 16 October for Agropyron cristatum were 58.9 and 11.57 days, respectively (Table 4, Fig. 3); and (2) the starting date is set late, leading that chilling thermal requirement can not meet. e.g., the RMSE of SM with starting date on 1 September was smaller (3.87 days) than 13 November (92.33 days) for Larix dahurica (Table 3, Fig. 3).  The fixed parameters (t 0 , T low , T high , va, vb and vc) in SM and PM, from the result of Europe [9,11,12,16], were widely used in plant phenology models [10,21,23]. Furthermore, these parameters were estimated using local observing data [31]. In the model fitting process of the present study, all thresholds and parameters of SM and PM were estimated. We found that the models with many parameters could be fitted well (Tables 3, 4, 5), but could not give accurate simulation when tested with the independent data (Fig. 2). This finding was consistent with the result from Linkosalo et al. [33], and it might be because the models were overparameterized and able to adapt to noise in addition to the phenological phenomenon itself [33]. The parameter estimate lies on the parameter space boundary [34].
Variation in the base temperature has no significant effect on the precision of spring phenology models [14,30,35]. Generally, the heat unit total depends on the threshold used [35], as with Leymus chinensis in this study (Tables 3,4,5). Therefore, changes in the base temperature induce different thresholds for the accumulated temperature, resulting in no significant variation in model accuracy. The threshold measure is a mathematical construct which may or may not be related to the physiological threshold [34]. Physiological parameters can be estimated from simulation experiments, but can not be obtained from the process of parameter optimization. This is because the optimization process is mostly dependent on the precision of observed field data, sample number, and local climate conditions. The biological interpretation of model parameters should not be considered as absolute [34,36]. The base leafing temperature of the same plant can be different simply due to different models, as seen with all plant species in this study (Tables 3,4,5). The base temperature in the growth season index (GSI) is fixed, and the minimum temperature is derived from experimental data [36,37]. In the present study, the leafing dates of woody plants actually changed can be, to a certain extent, explained by GSI using the fixed parameter. However, considering the threshold of 0.5 from the original model [25], the predicted dates for plant leafing in Northeast China was earlier than the observed values. Thus, the original parameter threshold for GSI was too small for the present study, and the optimal threshold varied with different species.
Previous studies indicated that the phenology of woody plants in temperate regions can be accurately predicted by a temperaturebased model (e.g., SW). For example, SW is considered to be the best model to accurately simulate the bud development of Picea abies [30] and has been validated [23]. Chilling has been introduced into some temperature-based models (e.g., AM) to improve accuracy. For example, AM is much more suitable for simulating budburst of Picea abies [7]. Furthermore, there is a correlation between chilling and forcing, i.e., forcing takes the effect after the chilling has finished [38].
In this study, plant leafing in Northeast China was simulated using four temperature-based and two temperature-precipitation based phenology models. When validated with independent data, SW and TP could give best simulation of the woody plant leafing. The effect of temperature in TP was the same with SW, and its accuracy was consistent with SW in moist conditions. The effect of precipitation in TP does not change the model manner of temperature, therefore (1) TP with fixed starting date (TPn) could be used to simulate the leafing affected by hydrological factors. For example, leafing of Leymus chinensis and Stipa krylovii in Xilinhot, Inner Mongolia was delayed 27 and 22 days due to water stress in some years; in other locations, the average delay time for Leymus chinensis, Stipa baicalensis, Elymus dahuricus and Agropyron cristatum reached 7, 5, 14 and 16 days, respectively; and (2) the leafing of woody plants in Northeast China was mainly driven by thermal conditions, and hydrological conditions were not limiting factors. For example, the average RMSE of SWn was much close to TPn for tree leafing (Fig. 2). This finding was consistent with other studies in which the precision simulating the leafing of broadleaved deciduous plants could not be substantially improved by adding precipitation into the model [26]. Kramer et al. [39] also believed that, in temperate zone, water availability mainly affects leaf area index and has little effect on leaf phenology. However, different species have different sensitivities to water conditions, and the leafing of some trees may be affected by hydrological conditions. Thus, TP is suitable for a wider range of plants because it considers the effects of both temperature and precipitation on plant leafing. TP was far superior to other five models in estimating herb leafing in Northeast China, and all simulated dates were very close to the observed dates (Fig. 4B). Other five models can not accurately simulate herb leafing. In addition to temperature, water availability was shown to be an important controlling factor for the phenology of herbs [18]. For example, the annual precipitation of 5 years (1984, 1989, 2000, 2001, and 2002) was less than 210 mm in Xilinhot, and was negatively correlated with herb leafing in next year (R 2 .0.9; Fig. 5). Therefore, the reduction of previous annual precipitation leads to the delay of herb leafing. It is reasonable to consider the effect of previous annual precipitation in TP. Because herbs generally have shallow roots, they tend to be strongly affected by hydrological conditions. Herb leafing in response to environmental conditions can be estimated using hydrological factors incorporated into the TP model. Although GSI can model the effects of hydrological factors with vapor pressure deficit (VPD), the validation of the model was rather poor in the present study, due to the lack of available VPD value or the lower sensitivity of VPD to hydrological conditions (Fig. 4). In addition, two instances can indicate that TP is superior over SW in more  arid areas: (1) the RMSE of TPn was smaller than that of SWn for Leymus chinensis and Stipa krylovii in arid Xilinhot, and (2) the predicted leafings from SWn and TPn with the same T b and F* were very close in moist years, and much closer to the observed value from TPn in drought years than SWn (data not shown).
Phenology model parameters can be obtained experimentally, but most phenology models use parameters optimized based on long-term observed data, e.g., the four temperature-based models (SW, SM, PM, and AM) and the temperature-precipitation based leafing model (TP). Different parameters are selected for different plant species in various regions. Therefore, sufficient data should be used to ensure the effectiveness and reliability of the model parameters. In this study, leafing simulations of Pinus koraiensis and Picea koraiensis were poor because the parameters were optimized based on only six years' observed data ( Table 6, Fig. 2). There is a strong possibility that more errors occur when the model is based on less data. However, the parameters of leafing phenology models for the other six woody plants were optimized using 12 years' observed data, and all models accurately simulated plant leafing. Overall, TP will be more suitable and reliable for modeling both woody and herbaceous plant leafing given climate changes (especially variation in hydrological conditions), while other leafing models that do not consider water will be less applicable in semiarid and arid areas.

Author Contributions
Conceived and designed the experiments: GZ RL. Performed the experiments: RL. Analyzed the data: RL. Contributed reagents/materials/analysis tools: GZ RL. Wrote the paper: RL. Contributed to the revision of the manuscript: GZ.