# Characteristics of Spring Phenological Changes in China over the Past 50 Years

Authors: Junhu Dai; Huanjiong Wang; Quansheng Ge
Journal: Advances in Meteorology
 (2014)

---

## Abstract

In order to understand past plant phenological responses to climate change in China (1963–2009), we conducted trends analysis of spring phenophases based on observation data at 33 sites from the Chinese Phenological Observation Network (CPON). The phenological data on first leaf date (FLD) and first flowering date (FFD) for five broad-leaved woody plants from 1963 to 2009 were analyzed. Since most phenological time series are discontinuous because of observation interruptions at certain period, we first interpolated phenological time series by using the optimal model between the spring warming (SW) model and the UniChill model to form continuous time series. Subsequently, by using regression analysis, we found that the spring phenophases of woody plants in China advanced at a mean rate of 0.18 days/year over the past 50 years. Changes of spring phenophases exhibited strong regional difference. The linear trends in spring phenophases were −0.18, −0.28, −0.21, −0.04, and −0.14 days/year for the Northeast China Plain, the North China Plain, the Middle-Lower Yangtze Plain, the Yunnan-Guizhou Plateau, and South China, respectively. The spatial differences in phenological trends can be attributed to regional climate change patterns in China.


---

## Body

## 1. Introduction

Plant phenology, which is the study of seasonal plant development events and their relationship to environmental factors [1], has attracted much attention in the context of climate change [2, 3]. Plant phenophases can be directly affected by the interannual variations of climate factors, such as temperature, light, and moisture [4]. Also, phenology can in turn affect climate [5, 6]. For example, a longer presence of green cover in large areas should alter physical processes such as albedo, latent and sensible heat, and turbulence [5]. Moreover, a longer growing season can influence ecosystem productivity and vegetation-atmosphere CO2 exchange [7, 8]. Thus, the study of past phenological changes is beneficial for assessing the impacts of climate change [9].In recent years, pronounced phenological shifts have been detected on all the continents of the world based on satellite reflectance data [10–12] or ground observation data [13–17]. For example, using the Advanced Very High Resolution Radiometer (AVHRR) Normalized Difference Vegetation Index (NDVI) dataset, Stöckli and Vidale [10] found a prolonged growing season in Europe over the past two decades. Similar results were also found in North America and China based on the same satellite dataset [11, 12]. Regarding ground-based phenological change, Menzel et al. [14] did a systematic assessment of European phenological responses to climate change and found the spring and summer phenophases such as timing of leaf unfolding and flowering had advanced by a mean trend of 2.5 days/decade from 1971 to 2000. Such earlier spring phenophases in recent decades were also found in North America [18]. In East Asia, several studies showed clear phenological responses to climate change in Japan and South Korea based on ground observations [13, 19, 20]. In China, although some work has been done to assess phenological changes in several locations [16, 21–24], systematic studies of phenological shifts at a national scale over a long time period are still lacking.In China, phenological observations have been conducted by the Chinese Phenological Observation Network (CPON) in 1963. Since that time, phenophases of typical plant species at each site of the CPON network were monitored by relevant local organizations, such as botanical gardens, weather stations, research institutions, universities, and middle schools. There was no observation, however, carried out in certain periods. Given that discontinuous time series of phenophases affects the assessment of phenological trends, we firstly interpolated the missing phenological data of five widespread woody plants using phenological models. Subsequently, we systematically studied the characteristics of spatiotemporal changes in spring phenology so as to evaluate the specific impacts of climate change on biological systems.


## 2. Materials and Methods

### 2.1. Phenological and Meteorological Data

Phenological data used in this study were derived from CPON. Spring phenophases of five deciduous woody plants, includingFraxinus chinensis,Ailanthus altissima,Melia azedarach,Paulownia tomentosa, andUlmus pumila, were investigated. We only used time series with more than 11 years of data from 1963 to 2009. As a result, a total of 33 sites have been included (Figure 1; Table 1). We took first flowering date (FFD) and first leaf date (FLD) as being representative of spring phenophases. According to the uniform observation criteria and guidelines of CPON [25], FFD and FLD are defined as the date when a fixed individual formed its first full flower and first full leaf, respectively. Overall, a total of 118 phenological time series (a phase of a specific species at certain site is considered as one time series) are analyzed (Table 1). As shown in Figure 2, there are a lot of missing data in each year for all phenophases. Especially during the period of 1969–1972, no observation data was made due to the social upheavals of the Cultural Revolution. Similar problems with data collection also existed in the period of 1997–2002 due to funding shortages. The available data (a total of 2316 observations) could only account for 41.8% of the full dataset (118 time series × 47 years), so these missing data would significantly affect the results of trend estimates. Thus, we interpolated the missing data using phenological models. Although these interpolated phenological data may introduce considerable bias into the results, these uncertainties could be quantified in our analyses (see Section 2.2).Table 1
Species selected and corresponding observation sites in the study. Site numbers correspond to site numbers shown in Figure1.
Species
Numbers of sites
Distribution sites (site number)
Fraxinus chinensis
3
9, 15, 28
Ailanthus altissima
7
8, 9, 12, 14, 15, 17, 21
Melia azedarach
16
14, 15, 17, 18, 23, 24, 24, 25, 26, 29, 30, 31, 33
Paulownia tomentosa
15
9, 10, 13, 15, 18, 19, 20, 21, 24, 25, 26, 28, 29, 30, 32
Ulmus pumila
18
1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 14, 15, 16, 17, 18, 28Figure 1
Illustration of phenological observation sites and spring phenological trends at each site. Site numbers are shown.Figure 2
Number of phenological observations (including first leaf date and first flowering date) for five plant species from 1963 to 2009.Meteorological data were derived from the Chinese Meteorological Administration (http://cdc.cma.gov.cn/) and included daily mean temperatures from 33 meteorological stations (Figure 1). Most of these meteorological stations are relatively close to the corresponding phenological sites (usually less than 5 km), though for site 8 the corresponding phenological site is about 30 km away.


### 2.2. Interpolation of Discontinuous Series Using Phenological Models

The spring phenophase of a number of trees has been modeled successfully using accumulated forcing units, which are often calculated by the accumulated degree days above a threshold temperature, regardless of the presence of the additional constraint of a chilling requirement [26]. Since climates over the study area vary in type from site to site, the most applicable models are bound to be different in different places. Thus, we developed two types of models for each phenological time series and then chose the most accurate one to interpolate the time series. The first model used is the spring warming (SW) model, which is a model without the constraint of a chilling requirement [27]. Regarding the models with the constraint of a chilling requirement, we chose the UniChill model developed by Chuine [28].


#### 2.2.1. SW Model

In the SW model, plant development occurs in response to aggregated heat sums, or heating degree days (HDD), measured as the sum ofR
f
(
x
t
) (daily mean temperature x
t above a base temperature T
b), starting at a day of year t
0 (1). The HDD achieve the critical forcing temperature (F
*) at time t
y, which represent the FFD or FLD that occurred. The equations for SW model are
(1)
HDD
=
∑
t
0
t
y
R
f
(
x
t
)
=
F
*
R
f
(
x
t
)
=
{
0
if
x
t
≤
T
b
x
t
-
T
b
if
x
t
>
T
b
.


#### 2.2.2. UniChill Model

The UniChill model divided the process of bud development into two phases: dormancy and quiescence. Dormancy refers to the physiological state during which development and cell growth are prevented by internal factors despite favorable external conditions [29]. Plants enter the phase of quiescence (during which development and cell growth are triggered by warm temperature) when dormancy is broken. The equations are as follows:
(2)
S
c
=
∑
t
0
t
1
R
c
(
x
t
)
=
C
*
(3)
S
f
=
∑
t
1
t
b
R
f
(
x
t
)
=
F
*
(4)
R
c
(
x
t
)
=
1
1
+
e
a
(
x
t
-
c
)
2
+
b
(
x
t
-
c
)
(5)
R
f
(
x
t
)
=
1
1
+
e
d
(
x
t
-
e
)
,
with t
0 set at September 1, and seven species-specific parameters (a, b, c, d, e, C
*, F
*) fitted to phenological observations. Parameters a, b, and c define the response function to temperature R
c
(
x
t
), also called “chilling units”, which models the effect of “cold” temperatures in breaking dormancy (4). The accumulative sum of chilling units is called the state of chilling (S
c). Parameter C
* is the S
c threshold at which bud dormancy is broken (2). Parameters d and e define the response function to temperature (R
f
(
x
t
), also called “forcing units” which conditions the effect of “warm” temperatures during quiescence (5). Forcing units are accumulated as soon as C
* is reached at t
1. The model predicts that phenophase (T
b) occurs when the state of forcing (S
f) reaches a particular critical value F
* (3).The above two models were both fitted for respective species at each site using daily mean temperatures and the FLD or FFD time series through the least square method. The functionf
(
x
)
=
∑
i
[
r
i
(
P
)
]
2 is minimized in the parameter space P, where r
i
(
P
) is the residual:
(6)
r
i
(
P
)
=
d
i
(
P
)
-
d
i
obs
with d
i
(
P
) and d
i
obs being the simulated date and the observed date in the year i, respectively. We only used odd years of each time series to fit the parameters. The internal validity of each series was measured by the percentage variance explained by the model (R
2) and the root mean square error (RMSE) between the observed dates and the simulated dates.The external validity was also measured by theR
2 and RMSE between remaining observed dates in even years and corresponding simulated dates. The external validity measures the goodness of simulation for the years not used to fit the parameters, so the uncertainties of models could be represented by the RMSE of external validity. Between two models we chose the optimal one (with fewest uncertainties) to interpolate the missing data. In addition, in order to minimize the errors, only the time series with associated uncertainties less than seven days (RMSE of external validity <7) were retained for further analysis.


### 2.3. Estimation and Comparison of the Phenological Trends

The temporal trends of time series could be calculated as the slope coefficient of a linear regression model with phenophases as the dependent variable and years as the independent variable. Subsequently, we calculated the mean trend of all the phenological time series for each site and investigated the spatial difference of phenological changes in China. Finally, based on the mean phenological series for all five species across the 33 sites, we applied the Mann-Kendall trend test method for detecting monotonic trends in phenological time series [30]. In addition, we used a moving-trend analysis method to investigate the temporal evolution of phenological trends [31], that is, calculating the regression slope for each 31-year period by moving the center year with a step length of one year. For highlighting the decadal variation of spring phenophases, the anomalies (relative to the mean over the 1963–1990 period) of spring phenophases in each decade (1960s, 1970s, 1980s, 1990s, and 2000s) were calculated.


## 2.1. Phenological and Meteorological Data

Phenological data used in this study were derived from CPON. Spring phenophases of five deciduous woody plants, includingFraxinus chinensis,Ailanthus altissima,Melia azedarach,Paulownia tomentosa, andUlmus pumila, were investigated. We only used time series with more than 11 years of data from 1963 to 2009. As a result, a total of 33 sites have been included (Figure 1; Table 1). We took first flowering date (FFD) and first leaf date (FLD) as being representative of spring phenophases. According to the uniform observation criteria and guidelines of CPON [25], FFD and FLD are defined as the date when a fixed individual formed its first full flower and first full leaf, respectively. Overall, a total of 118 phenological time series (a phase of a specific species at certain site is considered as one time series) are analyzed (Table 1). As shown in Figure 2, there are a lot of missing data in each year for all phenophases. Especially during the period of 1969–1972, no observation data was made due to the social upheavals of the Cultural Revolution. Similar problems with data collection also existed in the period of 1997–2002 due to funding shortages. The available data (a total of 2316 observations) could only account for 41.8% of the full dataset (118 time series × 47 years), so these missing data would significantly affect the results of trend estimates. Thus, we interpolated the missing data using phenological models. Although these interpolated phenological data may introduce considerable bias into the results, these uncertainties could be quantified in our analyses (see Section 2.2).Table 1
Species selected and corresponding observation sites in the study. Site numbers correspond to site numbers shown in Figure1.
Species
Numbers of sites
Distribution sites (site number)
Fraxinus chinensis
3
9, 15, 28
Ailanthus altissima
7
8, 9, 12, 14, 15, 17, 21
Melia azedarach
16
14, 15, 17, 18, 23, 24, 24, 25, 26, 29, 30, 31, 33
Paulownia tomentosa
15
9, 10, 13, 15, 18, 19, 20, 21, 24, 25, 26, 28, 29, 30, 32
Ulmus pumila
18
1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 14, 15, 16, 17, 18, 28Figure 1
Illustration of phenological observation sites and spring phenological trends at each site. Site numbers are shown.Figure 2
Number of phenological observations (including first leaf date and first flowering date) for five plant species from 1963 to 2009.Meteorological data were derived from the Chinese Meteorological Administration (http://cdc.cma.gov.cn/) and included daily mean temperatures from 33 meteorological stations (Figure 1). Most of these meteorological stations are relatively close to the corresponding phenological sites (usually less than 5 km), though for site 8 the corresponding phenological site is about 30 km away.


## 2.2. Interpolation of Discontinuous Series Using Phenological Models

The spring phenophase of a number of trees has been modeled successfully using accumulated forcing units, which are often calculated by the accumulated degree days above a threshold temperature, regardless of the presence of the additional constraint of a chilling requirement [26]. Since climates over the study area vary in type from site to site, the most applicable models are bound to be different in different places. Thus, we developed two types of models for each phenological time series and then chose the most accurate one to interpolate the time series. The first model used is the spring warming (SW) model, which is a model without the constraint of a chilling requirement [27]. Regarding the models with the constraint of a chilling requirement, we chose the UniChill model developed by Chuine [28].


### 2.2.1. SW Model

In the SW model, plant development occurs in response to aggregated heat sums, or heating degree days (HDD), measured as the sum ofR
f
(
x
t
) (daily mean temperature x
t above a base temperature T
b), starting at a day of year t
0 (1). The HDD achieve the critical forcing temperature (F
*) at time t
y, which represent the FFD or FLD that occurred. The equations for SW model are
(1)
HDD
=
∑
t
0
t
y
R
f
(
x
t
)
=
F
*
R
f
(
x
t
)
=
{
0
if
x
t
≤
T
b
x
t
-
T
b
if
x
t
>
T
b
.


### 2.2.2. UniChill Model

The UniChill model divided the process of bud development into two phases: dormancy and quiescence. Dormancy refers to the physiological state during which development and cell growth are prevented by internal factors despite favorable external conditions [29]. Plants enter the phase of quiescence (during which development and cell growth are triggered by warm temperature) when dormancy is broken. The equations are as follows:
(2)
S
c
=
∑
t
0
t
1
R
c
(
x
t
)
=
C
*
(3)
S
f
=
∑
t
1
t
b
R
f
(
x
t
)
=
F
*
(4)
R
c
(
x
t
)
=
1
1
+
e
a
(
x
t
-
c
)
2
+
b
(
x
t
-
c
)
(5)
R
f
(
x
t
)
=
1
1
+
e
d
(
x
t
-
e
)
,
with t
0 set at September 1, and seven species-specific parameters (a, b, c, d, e, C
*, F
*) fitted to phenological observations. Parameters a, b, and c define the response function to temperature R
c
(
x
t
), also called “chilling units”, which models the effect of “cold” temperatures in breaking dormancy (4). The accumulative sum of chilling units is called the state of chilling (S
c). Parameter C
* is the S
c threshold at which bud dormancy is broken (2). Parameters d and e define the response function to temperature (R
f
(
x
t
), also called “forcing units” which conditions the effect of “warm” temperatures during quiescence (5). Forcing units are accumulated as soon as C
* is reached at t
1. The model predicts that phenophase (T
b) occurs when the state of forcing (S
f) reaches a particular critical value F
* (3).The above two models were both fitted for respective species at each site using daily mean temperatures and the FLD or FFD time series through the least square method. The functionf
(
x
)
=
∑
i
[
r
i
(
P
)
]
2 is minimized in the parameter space P, where r
i
(
P
) is the residual:
(6)
r
i
(
P
)
=
d
i
(
P
)
-
d
i
obs
with d
i
(
P
) and d
i
obs being the simulated date and the observed date in the year i, respectively. We only used odd years of each time series to fit the parameters. The internal validity of each series was measured by the percentage variance explained by the model (R
2) and the root mean square error (RMSE) between the observed dates and the simulated dates.The external validity was also measured by theR
2 and RMSE between remaining observed dates in even years and corresponding simulated dates. The external validity measures the goodness of simulation for the years not used to fit the parameters, so the uncertainties of models could be represented by the RMSE of external validity. Between two models we chose the optimal one (with fewest uncertainties) to interpolate the missing data. In addition, in order to minimize the errors, only the time series with associated uncertainties less than seven days (RMSE of external validity <7) were retained for further analysis.


## 2.2.1. SW Model

In the SW model, plant development occurs in response to aggregated heat sums, or heating degree days (HDD), measured as the sum ofR
f
(
x
t
) (daily mean temperature x
t above a base temperature T
b), starting at a day of year t
0 (1). The HDD achieve the critical forcing temperature (F
*) at time t
y, which represent the FFD or FLD that occurred. The equations for SW model are
(1)
HDD
=
∑
t
0
t
y
R
f
(
x
t
)
=
F
*
R
f
(
x
t
)
=
{
0
if
x
t
≤
T
b
x
t
-
T
b
if
x
t
>
T
b
.


## 2.2.2. UniChill Model

The UniChill model divided the process of bud development into two phases: dormancy and quiescence. Dormancy refers to the physiological state during which development and cell growth are prevented by internal factors despite favorable external conditions [29]. Plants enter the phase of quiescence (during which development and cell growth are triggered by warm temperature) when dormancy is broken. The equations are as follows:
(2)
S
c
=
∑
t
0
t
1
R
c
(
x
t
)
=
C
*
(3)
S
f
=
∑
t
1
t
b
R
f
(
x
t
)
=
F
*
(4)
R
c
(
x
t
)
=
1
1
+
e
a
(
x
t
-
c
)
2
+
b
(
x
t
-
c
)
(5)
R
f
(
x
t
)
=
1
1
+
e
d
(
x
t
-
e
)
,
with t
0 set at September 1, and seven species-specific parameters (a, b, c, d, e, C
*, F
*) fitted to phenological observations. Parameters a, b, and c define the response function to temperature R
c
(
x
t
), also called “chilling units”, which models the effect of “cold” temperatures in breaking dormancy (4). The accumulative sum of chilling units is called the state of chilling (S
c). Parameter C
* is the S
c threshold at which bud dormancy is broken (2). Parameters d and e define the response function to temperature (R
f
(
x
t
), also called “forcing units” which conditions the effect of “warm” temperatures during quiescence (5). Forcing units are accumulated as soon as C
* is reached at t
1. The model predicts that phenophase (T
b) occurs when the state of forcing (S
f) reaches a particular critical value F
* (3).The above two models were both fitted for respective species at each site using daily mean temperatures and the FLD or FFD time series through the least square method. The functionf
(
x
)
=
∑
i
[
r
i
(
P
)
]
2 is minimized in the parameter space P, where r
i
(
P
) is the residual:
(6)
r
i
(
P
)
=
d
i
(
P
)
-
d
i
obs
with d
i
(
P
) and d
i
obs being the simulated date and the observed date in the year i, respectively. We only used odd years of each time series to fit the parameters. The internal validity of each series was measured by the percentage variance explained by the model (R
2) and the root mean square error (RMSE) between the observed dates and the simulated dates.The external validity was also measured by theR
2 and RMSE between remaining observed dates in even years and corresponding simulated dates. The external validity measures the goodness of simulation for the years not used to fit the parameters, so the uncertainties of models could be represented by the RMSE of external validity. Between two models we chose the optimal one (with fewest uncertainties) to interpolate the missing data. In addition, in order to minimize the errors, only the time series with associated uncertainties less than seven days (RMSE of external validity <7) were retained for further analysis.


## 2.3. Estimation and Comparison of the Phenological Trends

The temporal trends of time series could be calculated as the slope coefficient of a linear regression model with phenophases as the dependent variable and years as the independent variable. Subsequently, we calculated the mean trend of all the phenological time series for each site and investigated the spatial difference of phenological changes in China. Finally, based on the mean phenological series for all five species across the 33 sites, we applied the Mann-Kendall trend test method for detecting monotonic trends in phenological time series [30]. In addition, we used a moving-trend analysis method to investigate the temporal evolution of phenological trends [31], that is, calculating the regression slope for each 31-year period by moving the center year with a step length of one year. For highlighting the decadal variation of spring phenophases, the anomalies (relative to the mean over the 1963–1990 period) of spring phenophases in each decade (1960s, 1970s, 1980s, 1990s, and 2000s) were calculated.


## 3. Results

### 3.1. Validity of the Models

The SW model performed better (the RMSE of external validity is less) than the UniChill model in 67 time series, while the UniChill model showed better simulation in the other 51 time series. For each time series of spring phenophases, we used the optimal model to interpolate the missing data. In order to minimize the uncertainty of the interpolation as far as possible, we only retained the time series with models uncertainties less than 7 days. As a result, 92 of 118 time series were chosen for further analysis. The results of model validity for these series are shown in Table2. The average R
2 of internal validity for each phenophase ranged from 0.41 to 0.85 with an overall mean of 0.67, while the overall mean of R
2 for external validity was 0.62 (0.31–0.77). Accordingly, the overall mean of RMSE for internal validity and external validity was 3.7 days and 4.2 days, respectively. The variances of phenophases (standard deviation) were often more than 3 times stronger than the RMSE (Table 2); therefore, the error introduced by the models was acceptable.Table 2
Internal and external validity of optimal model for 92 time series of spring phenophases. Each column shows the mean ± standard deviation.
Species
Phase
N
DOY
Internal validity
External validity
R
2
RMSE
R
2
RMSE
Fraxinus chinensis
FFD
2
101.9 ± 12.2
0.41 ± 0.18
3.3 ± 0.3
0.31 ± 0.18
4.6 ± 1.6
FLD
2
100.3 ± 13.9
0.78 ± 0.29
2.4 ± 1.5
0.75 ± 0.28
3.3 ± 1.0
Ailanthus altissima
FFD
5
140.1 ± 9.9
0.69 ± 0.18
2.4 ± 1.2
0.42 ± 0.21
3.6 ± 0.5
FLD
7
105.9 ± 10.0
0.69 ± 0.15
3.3 ± 1.0
0.77 ± 0.12
3.4 ± 1.5
Melia azedarach
FFD
16
111.2 ± 19.2
0.67 ± 0.28
3.6 ± 1.8
0.66 ± 0.22
4.1 ± 1.6
FLD
9
89.3 ± 19.7
0.67 ± 0.28
3.9 ± 1.4
0.55 ± 0.28
4.7 ± 1.4
Paulownia tomentosa
FFD
9
89.3 ± 16.7
0.85 ± 0.17
2.3 ± 0.9
0.69 ± 0.18
3.4 ± 1.3
FLD
12
94.8 ± 15.2
0.74 ± 0.17
4.4 ± 2.3
0.59 ± 0.21
4.6 ± 1.1
Ulmus pumila
FFD
13
82.0 ± 23.8
0.63 ± 0.34
4.2 ± 2.1
0.60 ± 0.17
5.2 ± 1.4
FLD
17
103.4 ± 19.1
0.59 ± 0.21
4.1 ± 1.1
0.63 ± 0.23
4.1 ± 1.2
Overall
92
N/A
0.67 ± 0.24
3.7 ± 1.7
0.62 ± 0.22
4.2 ± 1.4


### 3.2. Spatial Patterns of Spring Phenophases Change

The frequency distribution of the phenological trends is summarized in Figure3. 88 of 92 spring phenological series showed earlier trends (62 series reached a 0.05 level of significance). Therefore, the earlier trends of spring phenophases from 1963 to 2009 are notable in China. The temporal trends of spring phenophases ranged from −0.35 to −0.03 days/year among different sites (Figure 1). The trends of spring phenophases in the main areas of China are as follows:(1)
The spring phenophases in the Northeast China Plain, represented by sites 1–6 and 8, showed a strongly earlier trend of 0.18 days/year;
(2)
the spring phenophases in the North China Plain, represented by sites 9 and 11–13, exhibited a very marked trend of −0.28 days/year, which was the strongest change in China;
(3)
the spring phenophases in the Middle-Lower Yangtze Plain, represented by sites 17–19, 21–23, and 26–28, advanced by a mean trend of 0.21 days/year. In the southern part of this area, however, the trends in spring phenophases were relatively weak (e.g., trends at sites 26 and 27 are only −0.07 and −0.14 days/year, resp.);
(4)
mean spring phenophases trend in the area of the Yunnan-Guizhou Plateau (represented by site 28) and the Sichuan Basin (represented by sites 20 and 24) was only −0.04 days/year, which was the weakest change in China;
(5)
the spring phenophases in the South China area (represented by sites 29–32) advanced by 0.14 days/year.Figure 3
Frequency distributions of trends (1963–2009) in spring phenophases. Solid values in columns represent significant trends (P
<
0.05).Furthermore, the spring phenophases at other sites also showed consistent shifts (Figure1). For example, the spring phenophases at site 7 (Hohhot), the only site located in the temperate grassland region, showed a trend of −0.22 days/year, while the spring phenophases at site 10 (Minqin), located in the temperate desert region, showed an advancing trend of 0.18 days/year. The trend at site 14 (Luoyang), located in the Yi-Luo River Basin, demonstrated a very weak trend (−0.07 days/year). The mean trend of −0.28 days/year at site 15 (Xi’an), located in the Weihe River Plain, was comparable with the trend in the North China Plain area. Site 33 (Mengla), the only one located in the tropical monsoon forest area, shows an obvious trend of −0.23 days/year towards earlier spring. Overall, the spring phenophases for most areas over China had advanced very significantly, but the strength of the advance varied among different regions.


### 3.3. Temporal Evolutions of Spring Phenophases

The annual changes of mean spring phenophases are shown in Figure4(a). We found that the spring phenophases have large interannual variation with about 18 days’ amplitude. In the context of the last 47 years, 1980 marks the latest spring phenophases and 2002 marks the earliest spring phenophases. The spring phenophase in 1980 was 5.0 days later than the 1963–1990 average, while the spring phenophase in 2002 was 13.2 days earlier than the 1963–1990 average. For the overall period of 1963–2009, the linear trend of spring phenophases was −0.18 days/year (P
<
0.001), suggesting a significant advance of 8.3 days. The Mann-Kendall trend test also indicates a significantly earlier trend in the mean spring phenophases time series (Z
=
-
3.72; P
<
0.001).The change of spring phenophases from 1963 to 2009 (a) and the 31-year moving trends of the spring phenophases (b) in China. The statistically significant trends (P
<
0.05) are marked with red circles.
(a)
(b)The 31-year moving linear analysis of the spring phenological series indicated that the time periods would affect the estimation of phenological trends (Figure4(b)). The spring phenophases showed a pronounced advancing trend of 0.30 days/year for the period of 1979–2009. The 31-year periods with the center year from 1989 to 1993 also showed significant advancing trends of around 0.23 days/year, though they were less than the trend over the 1979–2009 period (Figure 4(b)). In the 31-year periods with the center year before 1989, the trends of spring phenophases were insignificant. In general, the trends in the recent 31 years showed an unusual advancing trend that surpasses all previously observed trends in the past 31-year periods before 1979.Furthermore, in terms of the temporal evolution of spring phenophases, a very apparent decadal change was detected (Figure5). The spring phenophases in the 1990s and 2000s occurred 2.53 and 6.93 days earlier, respectively, than the 1963–1990 average. However, in the first three decades (1960s–1980s) the spring phenophases were close to the 1963–1990 average. Therefore, the advances in spring phenology can be said to have become evident after the 1980s and then strengthened in the 2000s.Figure 5
The spring phenology anomaly (relative to the 1963–1990 mean) in 1960s (1963–1970), 1970s (1971–1980), 1980s (1981–1990), 1990s (1991–2000), and 2000s (2001–2009).


## 3.1. Validity of the Models

The SW model performed better (the RMSE of external validity is less) than the UniChill model in 67 time series, while the UniChill model showed better simulation in the other 51 time series. For each time series of spring phenophases, we used the optimal model to interpolate the missing data. In order to minimize the uncertainty of the interpolation as far as possible, we only retained the time series with models uncertainties less than 7 days. As a result, 92 of 118 time series were chosen for further analysis. The results of model validity for these series are shown in Table2. The average R
2 of internal validity for each phenophase ranged from 0.41 to 0.85 with an overall mean of 0.67, while the overall mean of R
2 for external validity was 0.62 (0.31–0.77). Accordingly, the overall mean of RMSE for internal validity and external validity was 3.7 days and 4.2 days, respectively. The variances of phenophases (standard deviation) were often more than 3 times stronger than the RMSE (Table 2); therefore, the error introduced by the models was acceptable.Table 2
Internal and external validity of optimal model for 92 time series of spring phenophases. Each column shows the mean ± standard deviation.
Species
Phase
N
DOY
Internal validity
External validity
R
2
RMSE
R
2
RMSE
Fraxinus chinensis
FFD
2
101.9 ± 12.2
0.41 ± 0.18
3.3 ± 0.3
0.31 ± 0.18
4.6 ± 1.6
FLD
2
100.3 ± 13.9
0.78 ± 0.29
2.4 ± 1.5
0.75 ± 0.28
3.3 ± 1.0
Ailanthus altissima
FFD
5
140.1 ± 9.9
0.69 ± 0.18
2.4 ± 1.2
0.42 ± 0.21
3.6 ± 0.5
FLD
7
105.9 ± 10.0
0.69 ± 0.15
3.3 ± 1.0
0.77 ± 0.12
3.4 ± 1.5
Melia azedarach
FFD
16
111.2 ± 19.2
0.67 ± 0.28
3.6 ± 1.8
0.66 ± 0.22
4.1 ± 1.6
FLD
9
89.3 ± 19.7
0.67 ± 0.28
3.9 ± 1.4
0.55 ± 0.28
4.7 ± 1.4
Paulownia tomentosa
FFD
9
89.3 ± 16.7
0.85 ± 0.17
2.3 ± 0.9
0.69 ± 0.18
3.4 ± 1.3
FLD
12
94.8 ± 15.2
0.74 ± 0.17
4.4 ± 2.3
0.59 ± 0.21
4.6 ± 1.1
Ulmus pumila
FFD
13
82.0 ± 23.8
0.63 ± 0.34
4.2 ± 2.1
0.60 ± 0.17
5.2 ± 1.4
FLD
17
103.4 ± 19.1
0.59 ± 0.21
4.1 ± 1.1
0.63 ± 0.23
4.1 ± 1.2
Overall
92
N/A
0.67 ± 0.24
3.7 ± 1.7
0.62 ± 0.22
4.2 ± 1.4


## 3.2. Spatial Patterns of Spring Phenophases Change

The frequency distribution of the phenological trends is summarized in Figure3. 88 of 92 spring phenological series showed earlier trends (62 series reached a 0.05 level of significance). Therefore, the earlier trends of spring phenophases from 1963 to 2009 are notable in China. The temporal trends of spring phenophases ranged from −0.35 to −0.03 days/year among different sites (Figure 1). The trends of spring phenophases in the main areas of China are as follows:(1)
The spring phenophases in the Northeast China Plain, represented by sites 1–6 and 8, showed a strongly earlier trend of 0.18 days/year;
(2)
the spring phenophases in the North China Plain, represented by sites 9 and 11–13, exhibited a very marked trend of −0.28 days/year, which was the strongest change in China;
(3)
the spring phenophases in the Middle-Lower Yangtze Plain, represented by sites 17–19, 21–23, and 26–28, advanced by a mean trend of 0.21 days/year. In the southern part of this area, however, the trends in spring phenophases were relatively weak (e.g., trends at sites 26 and 27 are only −0.07 and −0.14 days/year, resp.);
(4)
mean spring phenophases trend in the area of the Yunnan-Guizhou Plateau (represented by site 28) and the Sichuan Basin (represented by sites 20 and 24) was only −0.04 days/year, which was the weakest change in China;
(5)
the spring phenophases in the South China area (represented by sites 29–32) advanced by 0.14 days/year.Figure 3
Frequency distributions of trends (1963–2009) in spring phenophases. Solid values in columns represent significant trends (P
<
0.05).Furthermore, the spring phenophases at other sites also showed consistent shifts (Figure1). For example, the spring phenophases at site 7 (Hohhot), the only site located in the temperate grassland region, showed a trend of −0.22 days/year, while the spring phenophases at site 10 (Minqin), located in the temperate desert region, showed an advancing trend of 0.18 days/year. The trend at site 14 (Luoyang), located in the Yi-Luo River Basin, demonstrated a very weak trend (−0.07 days/year). The mean trend of −0.28 days/year at site 15 (Xi’an), located in the Weihe River Plain, was comparable with the trend in the North China Plain area. Site 33 (Mengla), the only one located in the tropical monsoon forest area, shows an obvious trend of −0.23 days/year towards earlier spring. Overall, the spring phenophases for most areas over China had advanced very significantly, but the strength of the advance varied among different regions.


## 3.3. Temporal Evolutions of Spring Phenophases

The annual changes of mean spring phenophases are shown in Figure4(a). We found that the spring phenophases have large interannual variation with about 18 days’ amplitude. In the context of the last 47 years, 1980 marks the latest spring phenophases and 2002 marks the earliest spring phenophases. The spring phenophase in 1980 was 5.0 days later than the 1963–1990 average, while the spring phenophase in 2002 was 13.2 days earlier than the 1963–1990 average. For the overall period of 1963–2009, the linear trend of spring phenophases was −0.18 days/year (P
<
0.001), suggesting a significant advance of 8.3 days. The Mann-Kendall trend test also indicates a significantly earlier trend in the mean spring phenophases time series (Z
=
-
3.72; P
<
0.001).The change of spring phenophases from 1963 to 2009 (a) and the 31-year moving trends of the spring phenophases (b) in China. The statistically significant trends (P
<
0.05) are marked with red circles.
(a)
(b)The 31-year moving linear analysis of the spring phenological series indicated that the time periods would affect the estimation of phenological trends (Figure4(b)). The spring phenophases showed a pronounced advancing trend of 0.30 days/year for the period of 1979–2009. The 31-year periods with the center year from 1989 to 1993 also showed significant advancing trends of around 0.23 days/year, though they were less than the trend over the 1979–2009 period (Figure 4(b)). In the 31-year periods with the center year before 1989, the trends of spring phenophases were insignificant. In general, the trends in the recent 31 years showed an unusual advancing trend that surpasses all previously observed trends in the past 31-year periods before 1979.Furthermore, in terms of the temporal evolution of spring phenophases, a very apparent decadal change was detected (Figure5). The spring phenophases in the 1990s and 2000s occurred 2.53 and 6.93 days earlier, respectively, than the 1963–1990 average. However, in the first three decades (1960s–1980s) the spring phenophases were close to the 1963–1990 average. Therefore, the advances in spring phenology can be said to have become evident after the 1980s and then strengthened in the 2000s.Figure 5
The spring phenology anomaly (relative to the 1963–1990 mean) in 1960s (1963–1970), 1970s (1971–1980), 1980s (1981–1990), 1990s (1991–2000), and 2000s (2001–2009).


## 4. Discussion

The significant advance in FFD and FLD of the five tree species observed in this study is consistent with spring phenophase changes in other parts of the world (Table3). In Europe, an enormous systematic phenological network dataset of more than 100 000 observational series of 542 plants indicates an advance of 0.25 days/year in spring/summer phase for the period of 1971–2000 [14]. In this study, when restricting the time period to 1971–2000, the spring phenophase trend (−0.11 days/decade) is shown to be weaker than that observed in Europe (Table 3). Compared with the Northeastern US, however, the observed trend in the flowering time of lilac from 1965 to 2001 is close to our estimates [32]. This evidence indicates that the onset of spring across the Northern Hemisphere has appeared earlier over the past several decades. In the Southern Hemisphere, Australia has also experienced a warming climate in recent decades, leading to earlier wine grape maturity dates with a trend of −0.8 days/year (1985–2009) [33], which is about two times a stronger advance than indicated by spring phenophases in China (−0.43 days/year).Table 3
Comparison of trends in spring phenophases between this study and other studies based on ground phenological observations.
Regions
Objects
Periods
Trends (days/year)
Trends in this study
Source
Europe
Spring/summer events of 542 plant
1971–2000
−0.25
−0.11
[14]
Northeastern USA
The flowering time of lilac
1965–2001
−0.09
−0.12
[32]
Australia
Wine grape maturity date
1985–2009
−0.80
−0.43
[33]
China
The leaf unfolding time of 13 plant species
1982–2006
−0.41
−0.36
[34]In addition, another study from China suggested that the spring started 0.41 days earlier per year on average from 1982 to 2006 [34]. The result is stronger than our estimates (−0.33 days/year, Table 3) over the same period. The possible reason is that both the number of plant species and phenological observation sites involved are different from those relied on in this study. In [34], the authors discussed the FLD of 13 plant species at 20 observation sites, while our study consists of FLD and FBD for five plant species at 33 sites. Therefore, differences in species and phases selection, as well as distribution of phenological observation sites, could affect the results of trend estimates.The advance of spring phenophases in China has obvious regional differences. In general, the trends in northern China are stronger than those in southern China (Figure1). As a previous study suggested, the trend of annual temperature increase (1961–2000) was about 0.2-0.3°C/decade in northern China and less than 0.1°C/decade in southern China [35]. Therefore, in China, the phenological response to climate change matches the observed warming pattern. However, the distribution of phenological observation sites involved in this study is uneven, so these sites have limited spatial representativeness. In future, it will be necessary to enlarge CPON and enlist volunteers to help acquire more phenological data, which could minimize the effect of site distribution in estimating phenological change.Apart from the temperature in spring and the chilling temperature in winter, the photoperiod and evaporative demand also play an important role in regulating the growth and development of plants in temperate regions [36, 37]. Decreasing day lengths are reliable cues of the impending end of the growing season and winter onset for many temperate biomes, while increasing day lengths indicate the arrival of spring [37]. Because photosynthesis and growth are likely to be significantly limited if the vapor pressure deficit (VPD) is at a high value [38], the distribution of vegetation with different phenological patterns is very sensitive to seasonal changes in VPD. Among such influencing factors, the SW model used in this study only considers the effects of spring temperature, but the UniChill model also considers chilling requirements. Thus, the impacts of other potential factors on the model’s predictions need to be further studied.


## 5. Conclusions

Based on observational data from CPON, this study investigated the changes in FFD and FLD of five woody plants (includingFraxinus chinensis,Ailanthus altissima,Melia azedarach,Paulownia tomentosa, andUlmus pumila) in China. The results show that the spring phenophases in China became remarkably earlier at a mean rate of 0.18 days/year over the period from 1963 to 2009. The spring phenophases were stable from the 1960s to the 1980s but advanced by 2.5 days during the 1990s and 6.9 days during the 2000s (compared to the 1963–1990 mean). In addition, the changes of spring phenophases showed noticeable regional difference. The magnitudes of advance in the North China Plain, the Middle-Lower Yangtze Plain, and the Northeast China Plain are the strongest, while the magnitudes of advances in the Yunnan-Guizhou Plateau, the Sichuan Basin, and South China are weaker. In general, the change of spring phenophases in China matches the warming pattern observed over the past 50 years.



---

*Source: 843568-2014-04-07.xml*