Trend Evolution of Vegetation Phenology in China during the Period of 1981–2016

: The trend of vegetation phenology dynamics is of crucial importance for understanding vegetation growth and its responses to climate change. However, it remains unclear how the trends of vegetation phenology changed over the past decades. By analyzing phenology data including start (SOS), end (EOS) and length (LOS) of growth season with the Ensemble empirical mode decomposition (EEMD), we revealed the trend evolution of vegetation phenology in China during 1981-2016. Our study suggests that:


Introduction
Phenology refers to the study of annually or periodically recurring patterns of growth and development of plants and animal behavior [1,2]. Vegetation phenology can effectually reflect the complicated interactions among atmospheric, biospheric, and pedospheric biogeochemical properties, and thus the Intergovernmental Panel on Climate Change (IPCC) pointed out that phenology is one of the most direct and effective indicators of climate change [2,3]. Since the 20th century, especially in the past 20 years, the average surface temperature of the world has risen remarkably, and vegetation phenology has undergone tremendous changes as a response [4]. It may make a profound impact on all living creatures as plants are primary ecosystem producers [5]. Hence, monitoring variation of vegetation phenology is critical for understanding responses of ecosystems to climate change at different time scales [6][7][8].
Ground-based observation, a traditional method which can be dated back to 17th century, is widely used in phenology studies and can provide first-hand accurate data [9,10]. However, limited by the highly uneven spatial distribution of ground phenology observations and ambiguous seasonal and inter-annual variation in some regions, it is nearly impossible to detect phenology dynamics at regional or even global scale [2,5]. During the past several decades, the rapid improvement of satellite remote sensing technique makes it become an important tool to track vegetation phenology dynamics over large areas [11][12][13]. Satellite-derived phenology is defined as land surface phenology (LSP), which is commonly determined from vegetation indices (VI) such as normalized difference vegetation index (NDVI), the enhanced vegetation index (EVI) and two-band enhanced vegetation index (EVI2) [8,[13][14][15][16]. There are three commonly used phenological events: the start (SOS), end (EOS), and length (LOS) of growing season [17][18][19]. Correctly characterizing vegetation phenology dynamics and the inherent shifts are the foundation for understanding the climate change impacts on vegetation and ecosystems [2]. Previous studies have studied phenology dynamics during different periods and in different regions. A significant advanced trend in SOS occurred at the Northern Hemisphere, including 7.9 days per decade SOS earlier in temperate China from 1982-1999 [15] and 2.75 days per decade SOS earlier in China by meta-analyses [20]. Meanwhile, a significant delayed trend in EOS occurred, including 0.2 days per decade EOS later in Europe from 1971-2000 [21], 2.6 days per decade EOS later in China from 1981-2011 [22] and 2.4-3.6 days per decade EOS later in USA from 1982-2008 [23]. These studies showed that the magnitude of delayed trends in EOS was much weaker than that in SOS. However, most analyses consider only linear trends by assuming that the vegetation phenology will keep its change rate unchanged throughout the study period. Accounting for the nonlinear response of vegetation to climate changes (i.e., temperature), phenology trends should thus show nonlinearity rather than linearity [24]. Therefore, linear methods can hardly monitor trend evolution (i.e., changes of trends with time) and reveal the inherent trend shifts. Though methods like piecewise linear regression models have been proposed to gain trend shifts of phenology [25,26], these methods are sensitive to short-term fluctuations, which makes it hard to identify trend shifts for long time series [27]. These methods have a predetermined assumption that the change rates may abruptly alter at the turning points, but do not change with time neither before nor after the turning points. It means that these methods may be able to reveal the trends shifts, but unable to reveal the trend evolution. This assumption is true in several ecosystems such as vegetation changes before and after disturbances [28]. However, gradual climate change and land degradation are more likely to induce continuous changes in vegetation growth [29]. Therefore, phenology trends should have gradual temporal variation rather than abrupt change.
A potentially effective approach for extracting a trend from phenology data is ensemble empirical mode decomposition (EEMD) [30]. EEMD is an empirical but clearly efficient and adaptive time-frequency analysis method suitable for processing nonlinear and nonstationary time series [31,32]. EEMD decomposes complex time series into a small and finite number of simple intrinsic mode functions (IMFs) with different frequencies and a residual (trend). This analysis method is quite different from curve fitting methods, in which the functions are predetermined. The trend, containing at most one extremum and varying with time, does not follow a predetermined functional form. Moreover, this trend does not contain shorter timescales and also has low sensitivity to the extension of time domain [33], which means that even if lengthening the time domain, the physical interpretations will not change significantly. These attributes ensure that EEMD trends can robustly uncover more reliable hidden information on the nonlinear and nonstationary time series. EEMD has been proved to be effective in analyzing the trend evolution in vegetation growth [34][35][36][37][38][39].
Previous studies paid much attention to the spatial distribution of phenology trends, but little is known about trend evolution. In this study, using EEMD and Mann-Kendall test, we comprehensively analyzed the spatiotemporal variation of phenology dynamics in China during the period of 1981-2016. We aimed to address the following questions: (1) trend evolution of SOS, EOS and LOS at the national scale;(2) the trends in SOS, EOS and LOS in China and their spatial patterns; and (3) how the trend evolution changed in different latitudes, longitudes and land cover types, and discuss the potential factors.
We interpolated the missing value by using cubic spline interpolation for the grids whose number of missing value was less than 7 (~20% of the length of time series) and excluded the grids whose number of missing value was more than 7 [42]. This is because that the EEMD method has low sensitivity to the addition (extension) of new data [33]. Moreover, we also excluded the grids which had the same value during the whole period because they were unable to be analyzed in EEMD.

Land Cover Data
In this study, we used the MODIS MCD12C1 product, a yearly dataset providing the dominant land cover class and the sub-grid frequency distribution of land cover classes of each grid at a spatial resolution of 0.05 • from 2001-2016 [43].
The 17-class International Geosphere-Biosphere Programme classification (IGBP) was adopted [44]. To reduce the possible influence of classification errors and land cover change, we also conducted land-cover-specific analyses [45]. In the land-cover-specific analyses, only the grids with the same dominant land cover class during 2001-2016 were kept back. The 'closed shrublands', 'permanent wetlands' and 'open shrublands' were excluded because they contained too few grids. The 'water bodies', 'snow and ice' and 'barren or sparsely vegetated' were also excluded because they had no or little vegetation cover. Finally, the land cover data used in this study contain eleven land cover classes (Figure 1

Methods of Analysis
All the analyses were performed using MATLAB (R2016b).

EEMD
The EEMD method decomposes nonlinear and nonstationary time series data into n intrinsic mode functions (IMF; imfi, i=1, 2, …, n) and a secular trend. The decomposition process 'sifting' can be described as follows [30][31][32][33][34]: Step 1: Add a Gaussian white noise series w1(t) to the original data x(t). The amplitude of the Gaussian white noise series was set to 0.2 times standard deviation of the original data.
Step 2: Form the upper and lower envelope curves of the time series data x1(t) by connecting local maxima and minima with cubic splines respectively; then the time series data x1(t) minus the mean value m1(t) of the upper and lower envelope curves.
Step 3: Determine whether f1(t) satisfies the given criterion (close enough to zero at anywhere). If yes, stop sifting. Otherwise, take f1(t) as a new time series data and repeat step 2. In this way, obtain the first IMF: imf1(t).
Step 4: Obtain the remainder R1(t) by subtracting imf1(t) from x1(t). If R1(t) still contains oscillatory component, repeat 2 and 3 but with R1(t) being the new time series data.

Methods of Analysis
All the analyses were performed using MATLAB (R2016b).

EEMD
The EEMD method decomposes nonlinear and nonstationary time series data into n intrinsic mode functions (IMF; imf i , i = 1, 2, . . . , n) and a secular trend. The decomposition process 'sifting' can be described as follows [30][31][32][33][34]: Step 1: Add a Gaussian white noise series w 1 (t) to the original data x(t). The amplitude of the Gaussian white noise series was set to 0.2 times standard deviation of the original data.
Step 2: Form the upper and lower envelope curves of the time series data x 1 (t) by connecting local maxima and minima with cubic splines respectively; then the time series data x 1 (t) minus the mean value m 1 (t) of the upper and lower envelope curves.
Step 3: Determine whether f 1 (t) satisfies the given criterion (close enough to zero at anywhere). If yes, stop sifting. Otherwise, take f 1 (t) as a new time series data and repeat step 2. In this way, obtain the first IMF: imf 1 (t).
Step 4: Obtain the remainder R 1 (t) by subtracting imf 1 (t) from x 1 (t). If R 1 (t) still contains oscillatory component, repeat 2 and 3 but with R 1 (t) being the new time series data.
Remote Sens. 2020, 12, 572 Thus, x 1 (t) is decomposed into n IMFs with decreasing frequencies and a residual trend that is monotonic or has at most one extremum.
Step 5: Repeat steps 1-4 for l times (l was set to 100 in this study) with different Gaussian white noise series added each time. Obtain the ensemble means of the corresponding IMFs and trends of the decompositions as the final results.
A MATLAB EEMD package with the stoppage criterion and end treatment is downloadable at http://rcada.ncu.edu.tw/research1.htm.
The EEMD trend of phenology series at a given time t is defined as the value increment of R n from the beginning time 1981, that is: For the trends that are time-varying, their changing rates are also temporally local quantities and can be calculated as: Furthermore, we also conducted a significance test approach to test the EEMD trends. The detailed processes are as follows [34]: (1) generate 5000 white noise series with the same temporal length as the phenology series and implement EEMD to extract their secular trends; (2) divide the EEMD trends of that spatial location by the standard deviation of the corresponding phenology data; (3) find the 1.96 times standard deviation spread value of the trends of the white noise series as the 95% confidence interval; (4) judge whether the trend value is outside confidence interval.
The nonlinear trends interpreted by EEMD can be divided into five categories on the basis of their statistical significance [34].
Nonsignificant change: the trends were not significant at any year. Monotonic Decrease/Monotonic Increase (the monotonic trends): the trend was monotonically increasing/decreasing and was statistical significant for at least one year. If a place displayed Monotonic Increase/Monotonic Decrease trend, it meant that SOS/EOS experienced delayed/advanced trend or LOS experienced prolonged/shortened trend.
Increase to Decrease/Decrease to Increase (the shifted trends): the trend contained one local maximum/minimum and was statistical significant for at least one year. If a place displayed Increase to Decrease trend, it meant that SOS/EOS delayed first and then turned to advance or LOS prolonged first and then turned to shortened. And if a place displayed Decrease to Increase trend, it meant that SOS/EOS advanced first and then turned to delay or LOS shortened first and then turned to prolonged.
We calculated the average of all pixels with the same latitude or longitude or land cover class to study the phenology trends from the angle of latitude, longitude and land cover class (Section 3.3, Figures 4-6). Moreover, a moving average over 0.5 • in the latitudinal and longitudinal direction was used to eliminate possible noise. Function 'contour' in MATLAB was used to study phenology trends at different latitudes or longitudes and function 'imagesc' in MATLAB was used to study phenology trends at different land cover types.

Mann-Kendall Trend Test
In this study, Mann-Kendall trend test (M-K) was also used for trend analysis. Since this method does not require sample data to obey a certain distribution and is less affected by outliers, it can be used to effectively test the variation trend in time series [46]. The test statistic S is calculated by: The standard normal variate Z is computed by: Under the given significance level (P<0.05), there is a significant change if |Z| > Z 1−α/2 . When |Z| > 1.96, it means that the linear trend has passed the significance test of 0.05.
The linear trends interpreted by M-K were divided into three categories on the basis of their Z value.
Nonsignificant: Z value of the trend was between −1.96 and 1.96. Significant Increase/ Significant Decrease: Z value of the trend was more than 1.96/less than −1.96. If a place displayed Significant Increase/ Significant Decrease trend, it meant that SOS/EOS experienced delayed/advanced trend or LOS experienced prolonged/shortened trend.

Linear Regression
We used the linear ordinary least squares regression method to calculate the average rate of phenology trends during the whole study period: where t = 1, 2, 3, . . . , 36; y is value of phenology time series; α, β are the fitted intercept and slope, respectively; and ε refers to the error of the equation. A two-tailed t-test was used to verify the null hypothesis and examine the significance level of the trend.

Vegetation Phenology Variation at the National Scale
As shown from Figure 2, the national average timing of SOS advanced at a rate of 0.22 d/yr (P < 0.01) during the period 1981-2016. The EEMD trend showed a significant advance (P < 0.05) during the whole period, while its rate gradually raised and exceeded the linear rate in the early-2000s. Meanwhile, the national average timing of EOS delayed at a rate of 0.20 d/yr (P < 0.01). The EEMD trends in EOS showed a significant delay (P < 0.05) during the whole period but its rate gradually shrank and dropped below the level of the linear rate in the later-1980s. Since LOS was determined by SOS and EOS together, the length of growing season prolonged for about 13 days at the end of our study period. The rate also displayed a significant increase trend and exceeded the linear rate in the early 21st century.
Remote Sens. 2020, 12, x FOR PEER REVIEW 7 of 20 and higher than before, while that of EOS was just the opposite. Hence, it seemed that the increasing change rate of LOS might be mainly determined by SOS.  The linear trends in all phenological events showed the same tendency as the EEMD trends, which is an advanced SOS, a delayed EOS and a prolonged LOS. However, the rate of linear trends was stationary and that of EEMD trends was time-varying. The change rate of SOS became higher Remote Sens. 2020, 12, 572 8 of 19 and higher than before, while that of EOS was just the opposite. Hence, it seemed that the increasing change rate of LOS might be mainly determined by SOS. Figure 3 showed the spatial pattern of phenology trends, which was divided into five categories by EEMD (Nonsig, Monotonic Increase, Monotonic Decrease, Increase to Decrease and Decrease to Increase) and three by M-K (Nonsig, Significant Increase and Significant Decrease).

Spatial Pattern of Vegetation Phenology Trends
The EEMD and M-K outcome of each phenological event both suggested that the nonsignificant trends occupied the similar proportion of China (Tables 1 and 2). For EEMD method, among the four significant trend types in SOS, monotonic decrease trends prevailed and occupied the largest area, while the regions that experienced monotonic increase trends were the smallest (Figure 3a). Moreover, there were two types of trend shifts (including Increase to Decrease and Decrease to Increase) in EEMD outcome which should not be ignored. Although their distribution was relatively scattered, roughly concentrated in the Loess Plateau, North China Plain and north-eastern China, each of their area percentages was more than monotonic decrease trends and less than monotonic increase trends, but their total percentages outnumbered that of monotonic increase trends (29.37% vs. 25.81%; Table 1). For EOS and LOS, monotonic increase trends occupied the largest area, while monotonic decrease trends occupied the smallest (Tables 1 and 2). Moreover, the shifted trends in EOS and LOS were concentrated in the eastern Qinghai-Tibetan Plateau, eastern Sichuan Basin and north-eastern China, and eastern Sichuan Basin, Loess Plateau and north-eastern China, respectively (Figure 3b,c). Different from the EEMD outcome, as M-K method could not reveal the reversal of trends, there were only two significant trend types in the M-K outcome, without the two types of trend shifts. The area that experienced significant decrease trends in SOS was much larger than others (Figure 3d). The ratio of significant decrease trends was more than that of monotonic decrease trends from EEMD method. Meanwhile, the ratios of significant increase trends in EOS and LOS was much larger than others (Figure 3e,f). Considering the similar ratio of nonsignificant trends, it seems that a considerable part of shifted trends (i.e., Loess Plateau and north-eastern China) in EEMD was mistaken as significant decrease trends in M-K and thus the advance of SOS might be actually exaggerated. Similarly, there was also a considerable part of shifted trends in EOS (i.e., eastern Qinghai-Tibetan Plateau and eastern Sichuan Basin) and LOS (i.e., Loess Plateau and eastern Sichuan Basin) that was mistaken as significant increase trends, and thus the delay of EOS and the extension of LOS might be overstated.
It is clear that although considerable regions had the same phenology dynamics as the national overall trends (the advanced SOS, the delayed EOS and the prolonged LOS), the two methods agreed that the vegetation phenology trends in some regions of north-eastern China actually were the opposite ( Figure 3). Furthermore, EEMD detected that the shifted trends mainly distributed in Loess Plateau, north-eastern China, eastern Qinghai-Tibetan Plateau and eastern Sichuan Basin. However, they were misunderstood in M-K as significant decrease/increase trends and thus the change rate of phenology might be enlarged.

Evolution of Phenology Trends at Different Latitudes
As seen from Figure 4, it was fairly clear that mid-latitude and low-latitude (south of the Tropic of Cancer) (between the Tropic of Cancer and 45°N) area experienced advanced SOS, delayed EOS and prolonged LOS. Nevertheless, the phenology dynamics in the high-latitude (north of 45°N) were just the opposite.

Evolution of Phenology Trends at Different Longitudes
As can be seen from Figure 5, the amplitude of delay trends in EOS and prolonged trends in LOS increased with longitudes, while the beginning time of delay trends in EOS and prolonged trends in LOS decreased with longitudes. Meanwhile, it can be roughly divided into three regions, namely, eastern region (east of 121°E), central region (92°E-121°E) and western region (west of 92°E) based on the trends evolution of phenology varying with longitudes.
Although the eastern region experienced the same delayed EOS trends as the central and western regions, the trends in SOS and LOS were the opposite. That is that SOS delayed and LOS shortened in the eastern region, while SOS advanced and LOS prolonged in the central and western regions. Though the central and western regions experienced the same trends of phenology, their amplitudes and frequencies of trends variation with time were not the same. There were about two For SOS, the advance was starting at 40 • N in the early-1990s and moving southward and northward, and the amplitude became greater after two or three times of advance. Meanwhile, the high-latitude area began to delay since the very beginning of the study period and the delayed amplitude was getting weaker, and even changed from delay to advance. However, there was a delayed trend occurring in a narrow band centered at 20 • N and getting much stronger than before.
For EOS, considerable regions began to delay since the very beginning of the study period. But the evident delay appeared at 40 • N in the early-1990s. The delayed amplitude became much stronger after two or three times of delay. Especially since the early-2000s, there was an enhancement in delay occurring in the southern low-latitude area. Meantime, there was a continuous advance occurring in the high-latitude, although its amplitude was extremely weak.
LOS has been prolonging in the mid-latitude since the early-1990s. By and large, the higher the latitude was, the earlier the extension of LOS appeared. Moreover, since the early 2000s, the extension has further increased. Meanwhile, the length of LOS in the high-latitude changed from decrease to increase. In the most area, there were about three or four time of extension of LOS. It was clearly that the change characteristics of LOS in the mid-latitude and low-latitude were closer to the sum of SOS and EOS, indicating that it may be mainly controlled by SOS and EOS in the mid-latitude and low-latitude. Meanwhile, LOS was mainly controlled by SOS in the high-latitude.

Evolution of Phenology Trends at Different Longitudes
As can be seen from Figure    Although the eastern region experienced the same delayed EOS trends as the central and western regions, the trends in SOS and LOS were the opposite. That is that SOS delayed and LOS shortened in the eastern region, while SOS advanced and LOS prolonged in the central and western regions. Though the central and western regions experienced the same trends of phenology, their amplitudes and frequencies of trends variation with time were not the same. There were about two times of advance in SOS, four times of delay in EOS and five times of extension in LOS occurring in the western region. However, there were two times of advance in SOS, one time of delay in EOS and two times of extension in LOS occurring in the central region. The characteristics of the trends evolution of SOS and LOS in the eastern region were alike. The early characteristics of LOS in the central region was similar to those of SOS, while similar to those of EOS in the later. Meanwhile, the characteristics of LOS in the western region was similar to those of SOS and EOS. Hence, LOS was mainly induced by SOS in the east and by SOS and EOS in the west, respectively. Meanwhile, LOS in the central region was determined by SOS in the early stage and by EOS in the later stage. Figure 6 indicated that except Needleleaf forest (DNF and ENF), all other land cover types experienced the same trends as an advanced SOS, a delayed EOS and a prolonged LOS, though their evolution processes varied considerably. For SOS, most of the land cover types exhibited the advanced trends since mid-1990s and EBF was the earliest one. Moreover, except CRO, UB and ENF, other land cover types experienced several times of advance represented by different colors in Figure 6a. Meanwhile, GRA has the largest advanced amplitude. For EOS, GRA was the earliest one to exhibit obviously delayed trends, while EBF has the largest delayed days. Meanwhile, the delay of EOS began in 1990s for most land cover types, and experienced a continuously delayed trend, especially for GRA and EBF. Moreover, the ranges of other land cover types were all less than 10 days for both SOS and EOS, and mostly began to change since the late 1990s. Among them, the phenology trends of CRO was the most stable one, whose amplitude was always less than 5 days. The extension of LOS began in the later-1990s for most land cover types, except that GRA began in 1980s. For ENF, EBF, GRA and CNM, the extension of LOS was much larger than others, and experienced several times of variation. especially for GRA and EBF. Moreover, the ranges of other land cover types were all less than 10 days for both SOS and EOS, and mostly began to change since the late 1990s. Among them, the phenology trends of CRO was the most stable one, whose amplitude was always less than 5 days. The extension of LOS began in the later-1990s for most land cover types, except that GRA began in 1980s. For ENF, EBF, GRA and CNM, the extension of LOS was much larger than others, and experienced several times of variation.

Evolution of the Vegetation Phenology Trends at the National Scale
Our study showed that over the whole China, the average vegetation phenology derived from EVI2 showed that SOS has advanced by 0.22 d/yr, EOS has delayed by 0.20 d/yr and thus LOS has lengthened by 0.42 d/yr during 1981-2016 ( Figure 2). It is because that unprecedented global warming in the past several decades, especially in the past 20 years, may lead to earlier spring phenology and later autumn phenology, and thus further lengthened growing season [21]. The mean decrease rate of SOS estimated in our study was just a little smaller than Ma and Zhou (2.35 days per decade) and Ge et al. (2.75 days per decade) [20,22]. The spring phenology trend also generally agreed with that in Europe (2.5 days per decade) and in the Northern Hemisphere (2.8 days per decade) [21,47].
Ge et al. found that the advance of SOS became stronger since 1980s [22]. Our study suggested that not only the advance of SOS, but also the extension of LOS became stronger with EEMD method as time goes. Whereas the delay of EOS was getting weaker. Numerous studies suggested that warmer preseason (summer or autumn) temperatures delayed EOS and the preseason warming trends decelerated or even reversed in recent decades [15,48,49]. Hence, the rate of EOS was decreasing over the whole period. EEMD rate of SOS and LOS exceeded the linear rate in the early-2000s, while the EEMD rate of EOS dropped below the linear rate in the late-1980s. Hence, the variation of phenology trends over time seems to be largely nonlinear and more nonlinear methods like EEMD should be taken to detect the evolution of phenology trends [50].

Spatial Distribution of Different Types of Trends Evolution of Vegetation Phenology
Previous studies pointed out that the advance of SOS, delay of EOS, and extension of LOS are widely distributed in the world [5,15,22,48], which are consistent with the results of M-K method in our study. However, with EEMD method, the area that experienced Monotonic Decrease of SOS and Monotonic Increase of EOS and LOS were unexpectedly reduced. In the meanwhile, the shifted trends (including Increase to Decrease and Decrease to Increase) took up considerable area which were mostly from the main significant trends in M-K and have been ignored or misunderstood.
Our study found that North China Plain experienced Decrease to Increase of SOS and Loess Plateau experienced Increase to Decrease of SOS rather than continuous delayed trends reported by Ge et al. [22]. For EOS, we found that considerable area in north-eastern China and eastern Sichuan Basin experienced Increase to Decrease and some area in the eastern Qinghai-Tibetan Plateau experienced Monotonic Decrease and Increase to Decrease. Hence, the shifted trends mainly distributed in the north-eastern China, North China Plain, Loess Plateau, Sichuan Basin and Qinghai-Tibetan Plateau. For LOS, we found that the Increase to Decrease trends mainly occurred in eastern Sichuan Basin experienced, while the Decrease to Increase trends occurred in Loess Plateau and north-eastern China. These results show that the shifted trends did distribute widely in various regions, but were rarely mentioned in previous studies [15,22,26,50]. These discrepancies may be because that most traditional methods were linear and assumed that there were no trend shifts during the study period [34]. Without considering the potential trend shifts, the inherent trend shifts of vegetation phenology are completely ignored. Such neglect of trend shifts made the decrease trends in SOS change from 25.81% in EEMD to 49.87% in M-K the increase trends in EOS from 23.02% in EEMD to 44.89% in M-K, and the increase trends in LOS from 26.77% in EEMD to 48.24% in M-K. The nonlinear method, EEMD, can well reveal the trends evolution, and further extract the shifted trends which are usually ignored in linear methods. Therefore, we suggest that the trend shifts should be paid more attention with nonlinear methods to precisely reveal the trends of vegetation phenology including monotonic and shifted trends.

Vegetation Phenology Evolution with Different Latitudes, Longitude and Land Cover Types
Although previous studies have explored the spatial pattern of vegetation phenology [15,48,49,51], how the trends evolution of phenology along the latitudinal or longitudinal gradient was rarely reported. In our study, the amplitude and frequencies of the trends in all phenological events all evidently varied with latitude and longitude, especially with longitude.
From the angle of latitudinal gradient, the trends in the low and mid latitudes were exactly the opposite to those in the high-latitude. A delayed trend in SOS occurred in the high-latitude, while an advanced trend occurred in the low-latitude and mid-latitude weaker from 40 • N to the north and south. Moreover, the frequencies of trends evolution were also getting less from 40 • N to northward and southward. The delayed amplitude became much stronger after two or three times of delay in the mid-latitude. For EOS, an advanced trend occurred in the high-latitude, while a delayed trend occurred in the low-latitude and mid-latitude. Similarly, the frequencies of trends evolution were also getting less from 40 • N to northward and southward. These differences might be related to the effect of temperature. In the high-latitude area in China, including middle temperate zone and sub-frigid zone, vegetation may not be able to fulfill the chilling requirements in a warmer winter and thus induce a later SOS [52,53]. Liu et al. found that EOS in the high-latitude area was positively correlated with precipitation but negatively with insolation, and thus reduced precipitation and increased insolation led to an earlier EOS and a shorter LOS [48]. While in the low-latitude and mid-latitude, including warm temperate and sub-tropics zone, the higher preseason temperatures could easily thaw frozen soil and bring earlier snow first melt and lead to an earlier SOS [5]. Moreover, higher temperatures in the summer and autumn may enhance photosynthetic activities, reduce the potential frost damage and the number of freezing days and thus induce, a later EOS and a longer LOS [50].
From the angle of longitudinal gradient, it was extremely clear that there was a dividing line between the eastern region and central region, whose longitude was approximately 121 • E. The SOS outcome showed that the trends on both sides of 121 • E were exactly the opposite. An advanced trend stronger from east to west occurred in the central and western region, while a delayed trend occurred in the east. Furthermore, the EOS outcome also showed that the magnitude of delay was stronger from east to west as longitudes changes. This phenomenon might be related to the effect of water stress on vegetation growth. For arid/semi-arid regions like the western and central regions in our study, water availability is the limiting factor for vegetation growth. Limited water availability would restrain vegetation growth and photosynthetic activities [54]. Under global warmer climate condition, more snow melt and earlier snow first melt may affect the timing of spring phenology [5,55]. Moreover, increased precipitation in the summer and autumn may mitigate water stress on vegetation in the western and central regions, and then promote vegetation activities and induce a later EOS [5,55]. Nevertheless, increased winter temperature and declined radiation brought by increased precipitation in the eastern region played a negative and more dominate role in determining vegetation growth and thus induced a later SOS and an earlier EOS respectively [51]. We also noticed that in some areas, such as areas near 35 • N or 100 • E, the variation amplitude of trends was unexpectedly weaker, which was consistent with Zheng et al. [56]. They discovered that the trends of spring phenology in the Weihe Plain and the west of Henan Province barely changed and concluded that the weak variation amplitude should be attributed to the nonsignificant temperature change and even decreased temperature.
Various results have confirmed that an earlier SOS is generally associated with a longer LOS [57,58]. However, a later EOS may also prolong LOS. Zhu et al. concluded that the prolonged LOS was mainly due to the delayed EOS in North America during the period 1982-2006 [59]. Our results suggested that LOS was mainly determined by SOS in the high-latitude, while LOS was collectively affected by SOS and EOS in the low-latitude and mid-latitude. Moreover, LOS was induced by SOS in the eastern region, while LOS was affected by SOS and EOS together in the western region. However, in the central region, LOS was controlled by SOS in the beginning and after the early-2000s, LOS became a joint effort of SOS and EOS. These results revealed that the effects of SOS and EOS on LOS may vary with location and time.
Studies have found that the spring phenology for herbaceous plants showed an overall advanced trend, while the advanced trends for woody plants were much stronger [48,60]. Our study detected from the angle of different land cover types that needle-leaf forests (DNF and ENF) experienced a later SOS and DNF experienced an earlier EOS and a shorter LOS, while SOS advanced, EOS delayed and LOS prolonged in all the other land cover types since 1980s, 1990s and later-1990s, respectively. We noticed that only needle-leaf forests experienced the delayed trends in SOS. This phenomenon was mainly caused by the winter-spring temperature increasing. As needle-leaf forests were low-temperature vegetation widely concentrated in the high-latitude cold zone, they might not fulfill the chilling requirements in a warmer winter [5]. Though all the other land cover types experienced the same trends, their evolution processes varied considerably. For EBF and GRA, their amplitudes and frequencies of phenology variation were more than others, while CRO just had little small amplitude and nonsignificant trend evolution. This may be because that CRO is greatly managed by human and thus the response of CRO to climate change may be different from other vegetation. GRA mainly located in the mid-latitude temperate zone, a warmer winter may break the ecodormancy and thus an earlier SOS [5]. EBF mainly located in the mid-and low latitude subtropical zone, higher temperatures in the summer and autumn may enhance activities of photosynthetic enzymes and slow the speed of chlorophyll degradation [48,51]. These results suggested that different species showed the same trends but diverse trend evolution of phenology variation under global warming.

Limitation and Future Direction
Our results are subject to some caveats. First, the phenology dataset used in the study is at a spatial resolution of 0.05 • extracted from EVI2, which can be coarse for some areas. However, this dataset is still likely the best one we can currently use for a large-scale analysis. The dataset is composed of three widely-used satellites and hence it has covers a relative long time period, though it has a lower spatial resolution. Future development of higher-resolution satellite dataset may help us deepen understanding of the evolution of phenology trends. Second, we aimed to seek the effects of geographical factors on phenology trends and we did find that the difference exists in both latitudinal and longitudinal direction. This means that latitude and longitude are two suitable variables for analyzing the spatial pattern of vegetation phenology trends at a national scale. But considering the impact of Eastern Asia Monsoon, Westerlies and etc., the climatic factors (i.e., precipitation) are not strictly varying with geographical factors (i.e., longitude). Future studies should take both geographical factors and climatic factors and their interaction into consideration.

Conclusions
The EEMD results showed that the change rates of SOS and LOS were increasing with time, while that of EOS was decreasing. Moreover, the EEMD rates of SOS and LOS exceeded the linear rate in the early-2000s, while that of EOS exceeded the linear rate in the mid-1980s. Hence, EEMD method is able to well reveal the nonlinearity of phenology trends. Although a large regions experienced the advance of SOS, the delay of EOS and the extension of LOS, EEMD detected a large proportion of the shifted trends while M-K did not. All three phenological events showed latitudinal and longitudinal gradients, while longitudinal gradient was much stronger, especially in EOS. The area located in the mid-latitude and low-latitude displayed the same trends as the overall trends on the national scale, while the areas in the high-latitude experienced the contrary trends. The amplitudes and frequencies of phenology variation in the mid-latitude were stronger than those in the high-latitude and low-latitude. Areas in the western and central region experienced an advanced SOS trend and a prolonged LOS trend during the whole study period, which was contrary to the eastern region. The magnitude of delayed trends in EOS was stronger from east to west as longitudes change, while the beginning time decreased with longitudes. LOS was mostly determined by SOS in the high-latitude and the eastern region, while LOS was induced collectively by SOS and EOS in the mid-latitude, low-latitude and western region. Moreover, in the central region, LOS was controlled by SOS in the beginning and after the early-2000s, LOS became a joint effort of SOS and EOS. Almost all land cover types, except Needleleaf forests, experienced the advance of SOS, delay of EOS and extension of LOS. Among them, the magnitude of GRA and EBF was biggest one in SOS and EOS, respectively. Whereas, the amplitude of CRO was smallest one in all phenological events. As a time-varying analyses method, EEMD displayed the evolution of vegetation phenology dynamics and the spatial pattern of the shift trends. Thus, EEMD method is capable of better detecting temporal variation and hidden shift trends than traditional linear method.