The process-based spring phenology model for pecan leaf-out was developed by applying a modular model structure consisting of three sub-models (Supplementary Figure 1; Hänninen, 1990, 2016; Hänninen and Kramer, 2007; for recent applications of the approach, refer to Lundell et al., 2020; Zhang et al., 2021c). Sub-models I and II correspond to the various chilling and forcing accumulation models (Kramer, 1994; Chuine et al., 1998), respectively, but in order to emphasize the physiological phenomena modeled, we use the biological names of the rate and state variables here (Hänninen, 1995, 2016; Hänninen et al., 2019; Supplementary Figure 1).

Sub-model I addresses the effects of the chilling temperatures on the rate of a rest break, R_t(t) (Supplementary Figure 1). The state of rest break, S_r(t), is obtained by integrating R_r(t) with respect to time from the onset of rest, t₀, to time instant t. R_t(t) is scaled on a percentage scale, so that rest completion is predicted to take place when S_t(t) = 100%. Sub-model II addresses the ontogenetic development, i.e., the microscopic anatomical changes within the bud that lead to the visible phenological event, such as leafing out. In particular, Sub-model II addresses the effects of high forcing temperatures on the potential rate of ontogenetic development, R_o,pot(t) (Supplementary Figure 1). The potential rate indicates the rate after rest completion when the rate is no longer restricted by the rest status of the bud.

The restrictions caused by the rest status are addressed by sub-model III, where the value of a theoretical variable, ontogenetic competence, C_o(t), is calculated as a function of the state of rest break, S_r(t). C_o is a [0,1]-multiplier mediating the effects of rest status to the (realized) rate of ontogenetic development, R_o(t) (Supplementary Figure 1). With C_o(t) = 0, no ontogenetic development takes place regardless of the ambient temperature. At rest completion with R_t(t) = 100%, the value ontogenetic competence is C_o(t) = 1 by definition. Otherwise, the relationship between C_o(t) and S_r(t) can take different forms, making it possible to address contrasting models, such as the parallel model and the sequential model, within the same framework (Hänninen, 1990, 2016). Finally, the state of ontogenetic development, S_o(t), is obtained by integrating R_r(t) with respect to time from the onset of rest, t₀, to time instant t. Leaf-out is predicted to occur when S_o(t) = 100%.

In the autumn experiment, the differences in BB% between the treatments were analyzed by means of logistic regression with a binary response (bud burst: no/yes) using the type (natural vs. controlled in +6°C) and duration of chilling as the explaining factors. The differences in DBB between the treatments were analyzed by means of a two-way ANOVA using the same explaining factors as with BB%. A two-way ANOVA was also applied to analyze the differences in DBB between the treatments in the spring experiment. The explaining factors were the time of transfer from natural conditions to forcing conditions and the temperature under the forcing conditions.

With both natural and controlled chilling, the bud burst percentage, BB%, increased significantly from 10% at zero days of chilling to 100% in the fourth transfer after 42 days of chilling (Figures 1A,B and Table 1). The days to bud burst, DBB, decreased significantly with prolonged chilling, both natural and controlled (Table 1). For the controlled chilling, the decrease was clearly exponential, showing leveling off with long durations of chilling (Figure 1C), whereas for natural chilling, the decrease in DBB with prolonged chilling was only slightly exponential, with no clear leveling off (Figure 1D). In the spring experiment, DBB decreased significantly with increased forcing temperature (Figure 2 and Table 1). Furthermore, DBB also decreased significantly from the first transfer on February 13 to the last one on March 13 (Figure 2 and Table 1).

For constructing sub-model I, the value of Δt = 64 days was obtained in the autumn experiment for the duration of chilling required for rest break in the controlled chilling at constant +6°C (Figure 1C). On an hourly basis, this implies an empirical rate of rest break R_r,obs = 100%/(64 days × 24 h day^–1) = .065% h^–1 at +6°C (Equation 1; Figure 3A and Table 2).

The upper threshold, T_upp, of the rest-breaking temperature range in sub-model I was estimated by fitting the overall model (Supplementary Figures 1, 3) to the leaf-out observations in the WTC experiment (Figure 4). Simultaneously, the value of parameter t₀, starting day of rest period, was estimated by carrying out the calculations for estimating T_upp with three different values of t₀. The values of RMSE obtained for the three t₀ values November 1, November 15, and November 23 were 2.11, 1.67, and 2.24 days, respectively, implying the value of 15 November for t₀ (Table 2). Using that value of t₀, the RMSE between the predicted and the observed timing of leaf-out was minimized with T_upp = +13°C (Figure 4); so, that value was adopted as the estimate of T_upp (Table 2).

Sub-model II addresses the air temperature response of the potential rate of ontogenetic development, R_o,pot. This implies that the experiment for determining sub-model II should be started as close as possible to the estimated time of rest completion under natural conditions (refer to “Constructing the sub-models on the basis of experimental data” in the section “Materials and Methods”). In the autumn experiment, rest completion under natural conditions was estimated to take place at Δt = 90 days after the onset of the experiment on November 21, 2017 (Figure 1D), i.e., February 19, 2018. This means that the first (February 13), second (February 27), and third transfers (March 13) in the spring experiment (Figure 2) were carried out 6 days before, 8 days after, and 22 days after the estimated rest completion, respectively. Accordingly, the results of the first transfer in the spring experiment (Figure 2) were used in the construction of sub-model II (Figure 3B and Table 2).

The empirical value of ontogenetic competence, C_o, was zero in the first transfer with S_r = 0% (Figure 3C). This was because of the low value of BB% in the first transfer (Figure 1A; refer to Equation 9). Starting from the second transfer, the empirical value of C_o increased with increasing S_r (= with prolonged chilling). Accordingly, a piece-wise linear relationship with a discontinuity at S_r,min = 14% was determined for sub-model III (Figure 3C and Table 2).

According to the model prediction for the dormancy dynamics for the period of autumn 2017 to spring 2018, rest break started immediately at the beginning of the simulation in November 2017, as shown by the increasing value of the state of rest break, S_r, at that time (Supplementary Figure 4B). This was because the air temperature (Supplementary Figure 4A) at the time fluctuated mainly in the range causing rest break (Figure 3A). The simulated rest break progressed at an almost constant rate, leading to predicted rest completion on February 1, 2019, 20 days before the independently observed rest completion on February 21 (Supplementary Figure 4B), estimated on the basis of the autumn experiment results (Figure 1D). No ontogenetic development was predicted for the beginning of the simulation, as shown by the values of S_o = 0% at that time (Supplementary Figure 4B). That was because with S_r < 14%, the ontogenetic competence stayed at C_o = 0 (Figure 3C; ontogenetic competence not shown in Supplementary Figure 4B). With further rest break beyond S_r = 14%, ontogenetic development also started, but the S_o increased slowly at first because the ontogenetic competence was still relatively low and the air temperatures (Supplementary Figure 4A) were generally too low to promote ontogenetic development (refer to Figure 3B). Toward spring, the predicted rate of ontogenetic development increased, and leaf-out was predicted to occur on April 2, 2018, 1 day after leaf-out was independently observed under natural conditions (Supplementary Figure 4B).

The results presented in Supplementary Figure 4 provide an independent test of our model with observational data. However, as only 1 year was included in the test, the test is relatively weak, but it serves to illustrate the dormancy dynamics predicted by our process-based leaf-out model. In the main test, the model simulated the timing of leaf-out in the nine treatments of the WTC experiment with good accuracy (RMSE = 1.67 days, Figure 5).

For both RCP4.5 and RCP8.5, our model projected an overall advancing trend in pecan seedling leaf-out in all three locations addressed (Figure 6). There were, however, considerable differences between the locations and the warming scenarios. For the two northern locations Hangzhou and Hefei, an overall advancement at rates of about −1 and −1.5 days per decade, respectively, was projected under RCP4.5 (Figures 6A,C), and the rate was about double under RCP8.5 (Figures 6B,D). In the southern location Nanping, the advancement rate was about −1.5 days per decade under RCP4.5 and slightly smaller (less negative) under RCP8.5 (Figures 6E,F). In all three locations, there was more year-to-year variation in the timing of leaf-out under RCP8.5 than under RCP4.5 (Figure 6). In the southern location Nanping, the large year-to-year variation led to a delay in leaf-out in many years during the first decades of the 21st century despite the overall advancement trend (Figure 6F).

Our model projected no cold damage for Hangzhou under RCP4.5 (cold damage year percentage CDY% = 0), and under RCP8.5, projected cold damage was rare (CDY% = 3.8) (Figure 7A). Projected cold damage was also quite rare in Hefei, with the CDY% values at 5.1 and 6.4 for RCP4.5 and RCP8.5, respectively (Figure 7B). For Nanping, much more frequent cold damage was projected, and contrary to the two northern locations, there was more projected cold damage under RCP4.5 (CDY% = 16.7) than under RCP8.5 (CDY% = 9) (Figure 7C). No clear differences among the three locations were found in the annual indices of cold damage (incidence and severity), with the following two exceptions: first, in the simulation for Nanping under RCP4.5, there was one exceptional year with an exceptionally high incidence (4 days) of cold damage (Figure 7C). This is in line with the projection of cold damage being common in Nanping (CDY% = 16.7) under RCP4.5. Second, even though cold damage was projected for only 3 years in Hangzhou under RCP8.5 (CDY% = 3.8), it was more severe in Hangzhou than in the other two locations in one of those years (Figure 7).

In the sensitivity analysis, only a minor effect of the parameter T_upp on the CDY% values was found for Hangzhou and Hefei, indicating that our results for these two locations are robust (Supplementary Table 1). For Nanping, however, low values of T_upp (+11 and +12°C) decreased the CDY%, and high values of T_upp (+14 and +15°C) increased it in comparison with the CDY% values obtained in the main simulation with the estimated value of T_upp = +13°C (Supplementary Table 1).

Sparks (1993) demonstrated experimentally the existence of rest period (endodormancy) and chilling requirement for pecan trees growing over a large geographical range in the United States. Our experimental findings confirm that seedlings of the subtropical pecan grown in southeastern China also evince these traits. This is in line with recent studies, which show that several subtropical tree species have a rest period and a chilling requirement (Du et al., 2019; Song et al., 2020; Pan et al., 2021; Zhang et al., 2021a, b). The novelty and scarcity of these findings in the case of subtropical trees are in striking contrast with the case of boreal and temperate trees, where the existence of these traits has already been known for a 100 years (Coville, 1920; Jewaria et al., 2021).

We introduced a process-based spring leaf-out model for the subtropical first-year pecan seedlings we examined in this study. Our main emphasis was on the biological realism of the models (Levins, 1966, 1968; Charrier et al., 2015; Hänninen, 2016), which is why the model development was based on a hypothetico-deductive (HDM) approach, which makes use of experiments explicitly designed for examining the responses modeled (sub-models I–III). The approach is regarded as an HDM one because the model development is based on the implications of the models for the occurrence and timing of the readily observable leaf-out under experimental conditions rather than direct observations and measurements of the modelled physiological and morphological phenomena of rest break and ontogenetic development. Accordingly, all the data points in Figure 3 are determined by analyzing the leaf-out data gathered under the experimental conditions.

This approach, referred to as the ecophysiological approach by Hänninen (2016), has been taken in process-based tree phenology modeling only rarely (Sarvas, 1972, 1974; Campbell and Sugano, 1975, 1979; Caffarra et al., 2011; Zhang et al., 2021c). This is obviously because the experiments are time-consuming and laborious. Most earlier studies with process-based tree phenology modeling have been based on the technique of inverse modeling, where the model parameters are estimated by fitting the model to observational long-term phenological and air temperature records (Kramer, 1994; Chuine et al., 1998; Basler, 2016; this is a phenological approach in the terminology of Hänninen, 2016).

This is unfortunate, for it has been known since the pioneering study of Hunter and Lechowicz (1992) that inverse modeling with observational data runs a high risk of producing biologically unrealistic process-based models of tree phenology. To take just two different examples of the problems related to inverse modeling (for a review of the problems, refer to Hänninen et al., 2019), first, by applying inverse modeling to the growth onset of Scots pine (Pinus sylvestris) saplings in a warming experiment, Hänninen (1995) found that a model based on photoperiod outcompeted all models in which chilling causes rest break. However, in an explicit experiment where photoperiod was controlled, he found that photoperiod does not affect the rest break of Scots pine saplings at all. The cause of the problem was that in the warming experiment the natural photoperiod correlated with chilling accumulation, which actually caused the rest break. Second, using an inverse modeling approach, Luedeling et al. (2021) recently found that temperatures of up to +30°C cause rest breaks in apples (Malus domestica). As discussed by the authors themselves, this result is most probably an artifact caused by the lack of temperatures above +10°C in the observational data set.

Despite the experimental approach taken, our model still needs further testing and development. This is especially true in the case of sub-model I, where we had just one empirical data point for the model. This is where our study deviated most from the study of Zhang et al. (2021c), which we otherwise followed closely in building the model on the basis of experimental data. We solved the problem following Lundell et al. (2020), who also had just one controlled chilling temperature in their experiment: (1) the lower threshold of the rest-breaking temperature range was set at −3.4°C (Sarvas, 1974), (2) the upper threshold T_upp was estimated using an inverse modeling approach of fitting the model to the experimental results, and (3) a plateau response was assumed between the two threshold temperatures. In this way, we obtained a reasonable fit of the overall model to the experimental results (Figure 5). It should also be noted that in the mainline inverse modeling approach that uses observational phenological records (Kramer, 1994; Chuine et al., 1998; Basler, 2016), there are no experimentally determined data points available for estimating the value of any parameter of any sub-model.

Our estimated value of T_upp = +13°C is higher than the upper thresholds usually applied in process-based models for boreal and temperate trees (Sarvas, 1974; Cannell and Smith, 1983; Wang et al., 2020; Xu et al., 2020). However, the value is well in line with the experimental results available for subtropical trees. Zhang et al. (2021a, c) found that in two subtropical tree species the rest-breaking effect of +15°C was equal to that of +10 and +5°C, suggesting that in subtropical trees T_upp can even be well above +15°C (refer to also Jewaria et al., 2021). Further experimental studies with several controlled chilling temperatures and chilling durations (Baumgarten et al., 2021; Zhang et al., 2021a, c) are needed to test and improve the sub-model I developed for pecan seedlings in this study.

The starting date of the rest period (endodormancy) used in this study (t₀ = November 15) also involves a great deal of uncertainty. In that respect, our model is no exception among the published process-based tree phenology models. To the best of our knowledge, the value of t₀ is not based on explicit ecophysiological data in any of them (but refer to Fuchigami et al., 1982 and the related discussion in Hänninen, 2016, p. 119–120). Arbitrary starting dates, such as September 1 (Hänninen, 1990, 1991; Chuine et al., 1998, 2016) and November 1 (Cannell and Smith, 1983; Kramer, 1994), have been used with boreal and temperate trees earlier on. In subtropical trees, bud set typically occurs from September to October (Zhang et al., 2021a), which is why we used November 1 as the first candidate when estimating the value of t₀.

For our experiments, we used seedlings raised from seeds collected in the seminatural pecan forests of our study area in the Hangzhou region. It is well-known that seminatural pecan trees are of subtropical origin, but no exact information about their genetic background is available (Zhang et al., 2015). In pecan nut prediction, bred cultivars are normally used, which implies that the generality of our model should also be tested with bred cultivars in future studies.

Depending on the geographical location, the warming scenario, and the year of simulation, all three of our hypotheses gained support from some part of our simulation results. This even held for the delay hypothesis, for in many years of the first decades of the 21st century, leaf-out was delayed in the simulations for the southern location Nanping under the RCP8.5 warming scenario, despite the overall advancement trend found for Nanping under RCP8.5. For most years, however, advancement of leaf-out was projected for all the three locations and under both climatic scenarios addressed. In comparison with earlier studies with temperate trees (Kramer, 1994; Murray et al., 1994; Wang et al., 2020), the advancement rates found in this study were relatively low. This suggests that reduced chilling restricts the advancement of spring phenology more under the warm subtropical conditions than under the cooler temperate conditions.

A clear although not very rapid advancement of leaf-out was projected for the two northern locations, Hangzhou and Hefei. This held for both RCP4.5 and RCP8.5, and the rate of advancement under RCP8.5 was about double the rate under RCP4.5. These results show that in Hangzhou and Hefei, the effect of increased forcing overrode the effect of decreased chilling on the projected timing of leaf-out. For these two locations, cold damage was projected quite infrequently: once every 26 years (CDY% = 3.8%) at most in Hangzhou and once every 16 years (CDY% = 6.4%) in Hefei. In all, despite some cold damage projected, the prolonged growing season hypothesis was mainly supported by our findings for Hangzhou and Hefei.

The results projected for the southern location Nanping deviated from those found for the two northern locations in many ways. First, the advancement rate of leaf-out in Nanping was almost identical for the two warming scenarios. This showed that the effect of reduced chilling was stronger effect under the warming of the warmer climate in Nanping than under the warming of the cooler climates of the two northern locations and that Nanping is ultimately receiving less chilling with warming, which is impacting leaf-out timing. Second, as the projected advancement rates for Nanping were quite similar under the two climatic scenarios, more cold damage was projected under RCP4.5 than under the warmer scenario of RCP8.5. In Nanping, the value of CDY% was 16.7 for RCP4.5 and 9 for RCP8.5, implying that cold damage would occur once every 6 (RCP4.5) or once every 11 (RCP8.5) years. However, the results of the sensitivity analysis show a higher uncertainty for Nanping than for the other two locations. For Nanping, much less cold damage was projected for the low values (+11 and +12°C) of T_upp than for the estimated value of it (T_upp = +13°C) (Supplementary Table 1). This finding emphasizes the need for further experimental studies to improve sub-model I (refer to above).

In Nanping, CDY% was exceptionally higher under the cooler scenario RCP4.5 than under the warmer scenario RCP8.5. This was caused by the similar timing of leaf-out projected in Nanping for the two climatic scenarios. In other words, if leaf-out occurs simultaneously in a cooler and a warmer climate, then probably there is more cold damage in the cooler than in the warmer climate. For Nanping, then, a paradoxical conclusion arises: warming as such causes cold damage by false springs, but rapid warming causes such damage less than slow warming because with rapid warming the climate gets so warm after some decades that the potentially damaging temperatures below +5°C do not occur around the time of leaf-out any longer. This is why the projected damage in Nanping under RCP8.5 concentrates on the early decades, whereas under RCP4.5 cold damage is also projected for the late years. A similar phenomenon was actually found for Hefei; however, the results for Hefei were different in other respects (refer to the previous paragraph). This finding suggests that in the southern part of the subtropical zone, where the climate is relatively warm already, rapid warming creates a time window of high risks of false springs for a few decades. After the time window closes, the risk of cold damage decreases.

The modeling approach adopted in this study has been applied earlier for assessing the effects of climatic warming on spring phenology and the risk of false springs in various boreal and temperate trees (Cannell, 1985; Hänninen, 1991; Kramer, 1994; Linkosalo et al., 2000). The approach holds great potential for assessing the suitability of perennial horticultural crops for cultivation under climatic warming at different geographical locations. For instance, under climatic warming, many subtropical crops may find new cultivation areas under temperate conditions. The present modeling approach is well-suited for assessing the suitability of subtropical crops to such transfers to temperate conditions.