Direct Impacts of Climate Change and Indirect Impacts of Non-Climate Change on Land Surface Phenology Variation across Northern China

Land surface phenology (LSP) is a sensitive indicator of climate change. Understanding the variation in LSP under various impacts can improve our knowledge on ecosystem dynamics and biosphere-atmosphere interactions. Over recent decades, LSP derived from remote sensing data and climate change-related variation of LSP have been widely reported at the regional and global scales. However, the smoothing methods of the vegetation index (i.e., NDVI) are diverse, and discrepancies among methods may result in different results. Additionally, LSP is affected by climate change and non-climate change simultaneously. However, few studies have focused on the isolated impacts of climate change and the impacts of non-climate change on LSP variation. In this study, four methods were applied to reconstruct the MODIS enhanced vegetation index (EVI) dataset to choose the best smoothing result to estimate LSP. Subsequently, the variation in the start of season (SOS) and end of season (EOS) under isolated impacts of climate change were analyzed. Furthermore, the indirect effects of isolated impacts of non-climate change were conducted based on the differences between the combined impact (the impacts of both climate change and non-climate change) and isolated impacts of climate change. Our results indicated that the Savitzky-Golay method is the best method of the four for smoothing EVI in Northern China. Additionally, SOS displayed an advanced trend under the impacts of both climate change and non-climate change (hereafter called the combined impact), isolated impacts of climate change, and isolated impacts of non-climate change, with mean values of −0.26, −0.07, and −0.17 days per year, respectively. Moreover, the trend of SOS continued after 2000, but the magnitudes of changes in SOS after 2000 were lower than those that were estimated over the last two decades of the twentieth century (previous studies). EOS showed a delayed trend under the combined impact and isolated impacts of non-climate change, with mean values of 0.41 and 0.43 days per year, respectively. However, EOS advanced with a mean value of −0.16 days per year under the isolated impacts of climate change. Furthermore, the absolute mean values of SOS and EOS trends under the isolated impacts of non-climate change were larger than that of the isolated impacts of climate change, indicating that the effect of non-climate change on LSP variation was larger than that of climate change. With regard to the relative contribution of climatic factors to the variation in SOS and EOS, the proportion of solar radiation was the largest for both SOS and EOS, followed by precipitation and temperature.


Introduction
Land surface phenology (LSP), which refers to the responses to inter-and intra-annual variations in climate, is a useful indicator in the study of the response of ecosystems to climate variability [1].It has an important role in regulating regional and global carbon, water, and energy cycling [2,3].Changes in the timing of LSP have strong impacts on ecological processes, including the long-term distribution of tree species [4], plant competition [5], ecosystem functions [6], and even human activities (i.e., agriculture, forestry), human health (e.g., allergies), and tourism [7,8].Moreover, the advances in spring compound shifts in the growing season length increase the survival and activity of harmful insects and pathogens [9].Therefore, the phenology products have been widely applied to climatic change, farm management, and biomass monitoring [10,11].
Currently, there are four main methods for studying vegetation phenology: ground observation, satellite-based estimation, visible digital camera imaging [12,13], and bioclimatic models [14,15].Among the four methods, the satellite-based method, which has high spatial and temporal resolution, is an unprecedented, powerful, integrative, and objective technique for monitoring and characterizing LSP and its relationship with climate factors across a variety of spatial and temporal scales [16].However, the estimation of LSP for satellite-based methods is mainly based on smoothing or filtering vegetation index (VI) (i.e., Normalized Difference Vegetation Index (NDVI)); nevertheless, there has not been a general method that has been applied universally until now [17], which increases the uncertainness of estimation of LSP.For instance, the magnitude of the start of season (SOS) in temperate China was shown to range from −0.01 to −0.19 days per year across five different methods [18].Additionally, high-magnitude biases in LSP have been found in the same region and/or over the same time period [18].Due to the difficulties in separating the impacts of climate change and impacts of non-climate change on LSP effectively [19], the variation in LSP was mostly considered to be associated with climate change (i.e., increase temperature), despite the occurrence of simultaneous impacts of climate change and impacts of non-climate change.Therefore, simply attributing LSP variation to climate change is unjustified [20].In addition, the measurement sensitivity of the LSP to climate change has mainly been based on an ordinary linear regression model [21][22][23], which ignores the multicollinearity among climatic variables.These limitations have severely hindered our understanding of LSP variation and its correlation with impact factors (i.e., climatic factors).
Northern China covers an extensive territory characterized by complex ecosystems and climatic zones.In recent decades, Northern China has experienced dramatic changes in climate, such as the strongest El Niño events ever recorded [24], frequent occurrences of severe droughts [25], and intensive human activities (i.e., reforest, urbanization) [26].These characteristics make northern China one of the most critical and sensitive regions in which to observe the variation of LSP and its correlations with impact factors.Considering the fundamental role of LSP and to improve our ability to develop predictive models, it is urgent to determine whether trends in LSP have continued since 2000 as compared with previous studies, and particularly to investigate the relationships between LSP and diverse climatic factors and non-climate change.Therefore, four smoothing methods (the Savitzky-Golay method (S-G), the Double logistic function method (D-L), the Asymmetric Gaussian function method (A-G), and the Harmonic Analysis of Time Series method (HANTS)) were applied in this study to fitting VI data, and then chose the best smoothing results to estimate LSP, in order to (1) explore the variation of LSP in Northern China during 2001-2014; (2) isolate and quantify the impacts of climate change and impacts of non-climate change on LSP variation; and, (3) characterize the relative contribution of climatic factors (i.e., mean temperature) on the variation of LSP.We hypothesize that the trends and magnitudes of LSP will change since 2000 under the impacts of climate change and impacts of non-climate change in Northern China.

Study Area
The study area is located in China's temperate zone.The climate zones include warm, middle, and cold temperate from south to north, the mean annual temperature is between −4 and 14 • C, and the total annual precipitation ranges from 200 mm in the Northwest to 1000 mm in the Southeast.The dominant vegetation types are deciduous conifer forest, deciduous broad-leaved forest, temperate steppe, and temperate desert [27].Due to the markedly seasonality and spatial heterogeneity in thermal and moisture conditions, vegetation phenology displays extensive temporal and spatial differences [27].

Data
The MODIS surface reflectance (specifically, MOD09A1, version 6) with a spatial resolution of 0.5 km and an eight-day temporal resolution from 2001-2014 were used to in this study.All image data were reprojected from Sinusoidal to Asia North Albers Equal Area Conic, then subset to the study area.After the processing procedures (i.e., reprojection and imagery clipping), EVI values were calculated and pixels with EVI value less than 0.1 were excluded to mask non-vegetation further.Finally, LSP during 2001-2014 were estimated based on the EVI value mentioned above.
Monthly meteorological data from 2001-2014, including temperature, precipitation, and solar radiation, were acquired from the China Meteorological Data Sharing Service System (downloaded from http://data.cma.cn/,accessed on 16 May 2018).All meteorological data used in this study were verified by China's Meteorological Information Center; thus, false or missing data from some of the stations were eliminated.There are 301 recordings for temperature and precipitation, and 48 recordings for solar radiation in Northern China (Figure S1) were used.The Kriging method, which was recognized as a better interpolation with lower bias than other interpolation methods [28], was then used for the spatial interpolation of climate data across the study area.The spatial resolution of climate date was 0.5 km × 0.5 km, which is consistent with MODIS EVI data.

Comparison of EVI Smoothing Approaches
Four methods namely the Savitzky-Golay method (S-G), the double logistic function method (D-L), the asymmetric Gaussian function method (A-G), and the Harmonic Analysis of Time Series method (HANTS), were applied to reconstructed EVI time-series data for each year, and the algorithms for the four methods are detailed in [29,30].
Three statistical indicators, namely the root mean square error (RMSE), Akaike's Information Criterion (AIC), and the Bayesian Information Criterion (BIC) were used to evaluate the performance of each method mentioned above.The three statistical indicators were calculated as (Equations ( 1)-( 3)); the lower the values of the RMSE, AIC, and BIC, the more preferable the method is.
where EVI(t) is the mean EVI value obtained from the four methods (assumed to be accurate) [17], EVI * (t) is the result of EVI smoothing, k is the number of free parameters, N is the number of time points, and RSS is the residual sum of squares between the mean EVI and the corresponding smoothing methods.Our results indicated that the S-G method was the best of the four (Table S1).Therefore, the EVI results from the S-G method were then applied to extract LSP.

Extraction of LSP
TIMESAT software, which is widely used on regional and global scales [18], with the S-G filtering method, and an adaptation strength of 2.0, a seasonal parameter of 0.5, an S-G window size of 2, no spike filtering, and an amplitude season start and end of 20% was applied to estimate LSP (SOS and EOS).Furthermore, the interquartile range rule was applied to mask outliers for the date of SOS and EOS.For instance, date smaller than 25th percentile minus 1.5 times interquartile or date greater than 75th percentile plus 1.5 times interquartile were excluded [31].

Change Trend Analysis
The trends of LSP (SOS and EOS) during the period 2001-2014 were calculated at pixel level using a robust, non-parametric Mann-Kendall (M-K) trend analysis [32].This method does not require the independence and normality of the time series data [17].Previous studies reported that the M-K test statistic Z was approximately normally distributed when the sample size was greater than eight [32].A positive or a negative Z value indicated an increasing or a decreasing trend, respectively, which were all monotonic [32].The trends of LSP were tested at a significance level of α = 0.05.Additionally, the rate of change of LSP was calculated based on the Theil-Sen median slope estimator, which is more appropriate for assessing the rate of change in the short or noisy time series.

Time-Lag Effect Analysis
Previous studies have observed a time-lag effect between vegetation growth and climate variability on a monthly scale [33].The effects of lagged months during the preseason and growing season (see [23] for details) with one-month steps between LSP and climatic factors were calculated based on a partial least-squares regression (PLS) model.The month with the highest determination coefficient (R 2 ) of the estimated regression model was regarded as the lagged month for each climatic factor [33].For each year of the analysis, we defined the preseason period as May to the preceding November, and the growing season was considered to be June to October, in accordance with the LSP results mentioned above.
Our results indicated that nearly 80% of the total pixels displayed a four month time-lag effect between climatic factors (temperature, precipitation, and solar radiation) and the SOS, while for the EOS, the percentage was nearly 100% (Figure S2).

Isolated Effects of Climate Change and Non-Climate Change on Trends of LSP
To identify the isolated impacts of climate change on LSP variation, the mean temperature, cumulative precipitation, and cumulative solar radiation based on a four-month time-lag effect were first calculated.Then, the first-difference method, which calculates the differences between two consecutive years and it is regarded as a common de-trending technique for reducing the impacts of long-term trends that are caused by technological improvements or other effects [34,35], was applied to calculate the difference values (∆Val = Val t+1 − Val t ) of LSP and climatic factors between two consecutive years (year t and year (t + 1)).Thirdly, a PLS model was used at the pixel level to reduce the multiple colinearity among climate factors, and the corresponding coefficients of the model were considered to represent the sensitivity of LSP to the corresponding climatic factors.The multiple regression model was established, as follows: where ∆Phe, ∆Rad, ∆Tem, and ∆Pre represent the first difference of LSP, the solar radiation, temperature, and precipitation, respectively; Sen Rad , Sen Tem , and Sen Pre represent the sensitivity of LSP to solar radiation (days/(MJ m −2 )), temperature (days/ • C), and precipitation (days/mm), respectively; and, inte is the intercept of the partial least squares regression model.
Finally, the isolated effects of climate change and the isolated effects of non-climate change on the trends of LSP were calculated directly and indirectly, respectively [21].The calculation was as follows: Tre non_cli = Tre phe − Tre cli_phe (6) where Tre cli_phe and Tre non_cli represent the trend of LSP under the isolated impacts of climate change and isolated impacts of non-climate, respectively; Tre phe represent the trend of LSP under the combined impact of climate change and non-climate change.Tre Rad , Tre Tem , and Tre Pre represent the trends of the monthly cumulative solar radiation, monthly mean temperature, and monthly cumulative precipitation, respectively.The other parameters are defined the same, as in Equation (4).

Relative Contribution of Climatic Factors to the Variation of LSP
The relative contribution of each climatic factors to the trends of LSP were calculated, as follows: where Rel i represents the relative contribution of factor i to trends of LSP; Tre i represents the trend of factor i; m is the number of climatic factors, with m = 3 in this study; and, Sen i represents the sensitivity of LSP to climate i.

Impacts of Both Climate Change and Non-Climate Change on LSP Variation
The temporal trends of LSP under the impacts of both climate change and non-climate change (hereafter called the combined impact) in Northern China in 2001-2014 are displayed (Figure 1).Under the combined impact, an advanced SOS mostly occurred in Northeastern China and the North China Plain (Figure 1a), accounting for 68.85% of the total pixels (Figure 1b).However, pixels with a positive (i.e., delayed) trend in SOS accounted for 31.15% of the total pixels, which mainly occurred in the west of Northern China (Figure 1a).The average trend for the SOS under the combined impact was −0.27 days per year, and negative (i.e., advanced) trends for SOS mainly occurred between −0.8 and 0 days per year, while the values of positive trends mainly occurred between 0 and 0.5 days per year (Figure 1b).In terms of the EOS, a delayed EOS, mainly distributed in the northwest of Northern China and the North China Plain (Figure 1c), accounted for 72.21% of the total pixels (Figure 1d).However, the negative trends of EOS were mainly scattered in the northeast of Northern China (Figure 1c).The EOS was delayed with a mean value of 0.41 days per year, and the magnitude of positive trends mainly ranged from 0 to 0.5 days per year (Figure 1d).

Isolated Impacts of Climate Change and Non-Climate Change on LSP Variation
Under the isolated impacts of climate change on SOS, approximately 53.89% of the total pixels displayed a negative (i.e., advanced) trend, and these were mostly scattered in the west and south of Northern China (Figure 2a,b).In contrast, pixels with a positive SOS trend accounted for 46.11%, and they were mainly dispersed in the northwest and northeast of Northern China.SOS advanced with a mean value of −0.07 days per year under the isolated impacts of climate change, and the absolute values of the negative trends and positive trends for SOS both mainly ranged from 0 to 0.3 days per year (Figure 2b).Over the whole study region, about 61.61% of the total pixels displayed a negative trend for the SOS under the isolated impacts of non-climate change, and they mainly occurred in the northeast of Northern China (Figure 2c,d), while positive SOS trends accounted for approximately 38% of the total pixels, and these were mainly distributed in the west of Northern China.SOS advanced with a mean value of −0.17 days per year under the isolated impacts of non-climate change, and the magnitude of negative trends for the SOS mainly ranged from −1.5 to 0 days per year, while the positive trends for the SOS were mainly between 0 and 0.5 days per year.
Under the isolated impacts of climate change on EOS, approximately 60% of the total pixels displayed a negative trend and they were mainly distributed in the east, south, and west of Northern China (Figure 3a,b).EOS advanced with a mean value of 0.16 days per year, and the negative trend values mainly ranged from −0.3 to 0 days per year (Figure 3b).During the study period, the EOS trends under the isolated impacts of non-climate change were similar to those of the combined impact, with an average delayed trend of 0.43 days per year.The positive trends in the EOS under isolated impacts of non-climate change accounted for 71.95% of the total pixels, and they mainly occurred in the west, south, and east of Northern China (Figure 3c,d).In contrast, the negative trends in the EOS were mainly distributed in the north and northeast of Northern China (Figure 3c).The magnitude of the positive trends mainly ranged from 0 to 0.5 days per year, while the magnitude of the negative trends for the EOS were mainly between −0.3 and 0 days per year (Figure 3d).

Relative Contribution of Climatic Factors to LSP Variation
The relative contribution of climatic factors to the variation in LSP are illustrated in Figure 4.Among the three climatic factors, the impact of cumulative solar radiation on SOS was larger than those of mean temperature and cumulative precipitation.Approximately 35.77% of the total pixels displayed that the relative contribution of cumulative solar radiation to the SOS was the largest, and these were mainly dispersed in the east, west, and middle of Northern China (Figure 4a).Over 32% of the total pixels showed that the relative contribution of cumulative precipitation was the largest, and these were mostly scattered around desert areas and grassland areas (Figure 4a).For the remaining areas, the relative contribution of the mean temperature was the largest, and this was mainly distributed in the northwest and south of Northern China (Figure 4a).In terms of the relative contributions of climatic factors to the EOS, the proportion and spatial pattern for each climatic factor were similar to those of the SOS.Over 36.62% of the total pixels displayed that the relative contribution of cumulative solar radiation was the largest, and these were mainly distributed in the east and middle of Northern China (Figure 4b).Approximately 30.73% of the total pixels showed that the relative contribution of cumulative precipitaion was the largest, and these were mostly scattered around desert areas and the Northeast China Plain (Figure 4b).For the remaining areas, the mean temperature was the most important, which was mostly distributed in mountainous regions, such as the Great Khingan, Lesser Khingan, Changbai Moutain, and Tianshan Mountains (Figure 4b).

Temporal Variation in LSP
The SOS displayed advanced trends under the combined impact, isolated impacts of climate change, and isolated impacts of non-climate change, with mean values of −0.26, −0.07, and −0.17 days per year, respectively.The trends for the SOS were in line with previous findings [36,37].Interestingly, the effect of non-climate change on SOS variation was larger than that of climate change.Additionally, the magnitude of changes in the SOS after 2000 (this study) was less than that estimated over the last two decades of the twentieth century (previous studies), and this is supported by previous studies [36][37][38][39], which showed that the slope of SOS variation has decreased due to the deceleration of strong warming in recent decades [27] or because of changes in winter chilling or fire regimes on a regional scale [36,39].Previous studies reported that winter chilling requirements and the photoperiod are also known to play important roles in SOS.If chilling temperatures are not required before the end of winter in the previous year, the SOS might be delayed, despite a continued warming in the current spring.Additionally, strong photoperiod control may also limit the degree to which the SOS can advance in the future [40,41].In Northern China, during the period of 2001-2014, over 51% of the total pixels showed a positive relationship between the winter temperature and the SOS (Table S2).Additionally, the winter temperature displayed a decreasing trend (Figure S3), which satisfied the winter chilling requirement for leaf-out.Furthermore, solar radiation and precipitation in spring showed increasing trends (Figure S3), meeting both the material and water requirements for vegetation growth and resulting in an advanced trend for the SOS in Northern China during the period of 2001-2014.Notably, our results showed that the trends in the SOS continued after 2000, negating our first hypothesis; however, the magnitude of changes in SOS slowed when compared with previous studies that were conducted in the last two decades, partly affirming our second hypothesis.In terms of the EOS, it showed delayed trends under the combined impact and isolated impacts of non-climate change, with mean values of 0.41 and 0.43 days per year, respectively, which is consistent with previous studies [17,23,36,42].However, the EOS advanced with a mean value of −0.16 days per year under the isolated impacts of climate change.Our results displayed negative relationships between the EOS and both the autumn temperature and autumn solar radiation, accounting for approximately 52% and 63% of the total pixels, respectively (Table S2).Therefore, the advanced trend of the EOS under the isolated impacts of climate change may be due to the increasing trends of both temperature and solar radiation in autumn (Figure S3), both of which increase evapotranspiration and limit water supply for vegetation growth and then accelerate leaf coloration changes or leaf fall.The positive relationship between autumn precipitation and EOS illustrated this conclusion (Table S2).Additionally, our results indicated that the impact of non-climate change on EOS variation is larger than that of climate change.Although an advanced trend of EOS was found under the isolated impacts of climate change, the EOS was ultimately delayed under the combined impact.
Although the variation direction of LSP in Northern China was consistent with previous studies (except for the EOS trend under the isolated impacts of climate change), the magnitudes of both the SOS and EOS were diverse when compared to previous studies [18,22,36].This may have resulted from differences in the study period and study methods.For instance, Zhu et al. [43] revealed that the SOS in North America advanced by 0.27, 0.35, and 0.13 days per year, and the EOS was delayed by 0.78, 0.42, and 0.55 days per year from 1982 to 1991, 1982 to 1999, and 1982 to 2006, respectively.Liu et al. [44] illustrated that the EOS was delayed by 0.12 days per year from 1982 to 2011 in temperate China, and, similarly to Yang et al. [36], calculated 0.13 days per year for the EOS delay in temperate China from 1982 to 2010.However, both of these are much less than that of Piao et al. [23], which illustrated a delayed EOS with a rate of 0.37 days per year in temperate China during the period 1982-1998.Moreover, the discrepancies among methods may have resulted in different magnitudes in phenophase (i.e., SOS).For instance, the magnitude of SOS in temperate China ranged from −0.01 to −0.19 days per year across five different methods [18].The SOS was delayed at a rate of 0.03 and 0.13 days per year using the TIMESAT method [45] and piecewise logistic method [43], respectively.These results indicate the presence of relatively large uncertainties in the noise filtering of the vegetation index (i.e., EVI) time-series data, and a substantial influence in threshold setting for the onset date [18,46].Therefore, it is critical to choose an appropriate method and threshold for onset date extraction in a specific region, especially for vast areas with diverse vegetation types and/or harsh environments [46].

Relative Contribution of Impact Factors to LSP
In terms of the relative contribution of climatic factors to SOS, 35.77% of the total pixels were dominated by cumulative solar radiation, and these were mainly dispersed in the east, west, and middle of Northern China.This is because an increasing trend of cumulative solar radiation (Figure S3) in winter and spring is usually accompanied by abundant sunlight intensity and sufficient materials and energy, which enhance the photosynthetic capacity and promote vegetation growth.In addition, increased solar radiation may result in the melting of seasonal snow, which might lead to higher soil moisture content to support vegetation growth [47].The negative relationship between the SOS and precipitation in both winter and spring demonstrated the conclusion mentioned above (Table S2).Precipitation is one of the critical climatic factors that regulates vegetation growth in arid and semi-arid areas.Therefore, the relative contribution of cumulative precipitation accounted for approximately 32.22% of the total pixels, and these were mostly scattered around desert areas and grassland areas, which is in line with previous studies [26,48,49].This may be related to precipitation in that area is already deficient, with root systems for grassland that are shallow and soil profiles that are superficial; thus, the ability of grassland in this area to hold soil moisture is limited [50].A decreasing trend in temperature and an increasing trend in precipitation in spring (Figure S3) decreased evapotranspiration and provided sufficient water for vegetation growth and advanced the SOS.The relative contribution of the annual mean temperature to SOS variation was lower than that of the other two climatic factors, accounting for approximately 32.01% of the total pixels, which were mainly distributed in the northwest and south of Northern China.This response pattern occurred because the thermal requirement is essential for vegetation leaf onset in the Northwest of Northern China due to the lower temperatures in spring resulting from the high elevation.Additionally, the impacts of the Siberian high pressure are coldness and dryness.In contrast, in the south of Northern China, forest is the dominant vegetation type, and chilling requirements in winter and a warmer temperature in spring are more important than other climatic factors in regulating leaf onset [8].However, the continuous decreasing trend of temperature that occurred throughout winter and spring (Figure S3) identified temperature as a critical factor in regulating forest SOS variation in the south of Northern China.
Regarding the relative contribution of climatic factors to the EOS, the proportion of cumulative solar radiation was the largest, accounting for 37% of the total pixels, which were mainly distributed in the east and middle of Northern China.This is because the increasing trend for solar radiation in autumn (Figure S3) provides sufficient materials and energy for vegetation photosynthesis.Additionally, the increasing trends for temperature and precipitation in autumn (Figure S3) mitigate the importance of low temperature and water deficiency that may restrict vegetation growth in autumn.Similar to SOS, the largest relative contribution to precipitation occurred around the desert areas.This indicates that leaf coloring or senescence may be more dependent on the water supply for desert vegetation.For the remaining areas, such as the Great Khingan, Lesser Khingan, Changbai Moutain, and Tianshan Mountains, temperature contributed the most to the regulation of vegetation senescence due to the low temperatures resulting from the high elevation, and thermal requirements were shown to be more important than other climatic factors for vegetation growth.

Limitations
The different magnitudes and trends among phenophases (i.e., SOS) at the pixel level may be related to the diverse responses and adaptations to climate factors and non-climate factors (i.e., anthropogenic activities) and spatially uneven changes in climate change and non-climate change [37,44]; however, the relatively sparse distribution of climate data, especially for solar radiation, might also result in some ambiguity and error in the variation of LSP under the isolated effects of climate change and isolated effects of non-climate change.Additionally, the diverse changes in phenophases (i.e., SOS) may also result from calculation errors (i.e., from VI smoothing) to some extent.Furthermore, this study only considered three climate variables to calculate the LSP variation under the isolated impacts of climate change; however, other factors, such as natural disturbances (i.e., wildfires, pest, plant disease, droughts and floods) and other climate-related factors (i.e., winter chilling temperature, soil temperature and moisture, sunshine duration, CO 2 concentrations, and nitrogen enrichment and deposition), may also affect LSP variation [51], and it should be considered in future studies.Moreover, cropland with two growing seasons was included in the study region, but only the first growing season was considered.Therefore, the impacts of climate change and non-climate change on the second growing season of cropland should be conducted in the future.In addition, the variation in LSP was impacted by both climate change and non-climate change factors (i.e., anthropogenic activities) simultaneously, so determination of the isolated effects of these two factors on LSP per se is difficult [20].The method that was applied in this study to isolate climate change and non-climate change was based on the hypothesis that vegetation growth is impacted by climate change and non-climate change factors; after the removal of climate change impacts, the effects of non-climate change on vegetation change could be identified [52].Obviously, this method is not effective, and it did not distinguish different types of non-climate change factors (anthropogenic activities vs. vegetation succession and competition) on vegetation growth, despite the fact that this method has been widely used to isolate the effects of climate change and non-climate change, such as anthropogenic activities, not only to measure net primary production variation, but also to measure NDVI variation and vegetation phenology variation [21,47,[53][54][55].Therefore, it is necessary to identify the effects of different types of non-climate change factors, such as species coexistence, afforestation, and deforestation on LSP variation in the future.

Conclusions
In this study, land surface phenology (LSP) in Northern China during the period 2001-2014 was estimated, and the estimation and field observations were compared.Additionally, the trends for the SOS and EOS under the combined effect, isolated effect of climate change, and isolated effect of non-climate change were analyzed, respectively.Our results showed that the Savitzky-Golay method was the best of the four for smoothing the enhanced vegetation index (EVI) in Northern China.In terms of the trends for the SOS and EOS under different impact factors, the SOS advanced under the combined impact, isolated impacts of climate change, and isolated impacts of non-climate change.The EOS was delayed under the combined impact and isolated impacts of non-climate change.However, an advanced trend for the EOS occurred under the isolated impacts of climate change.Furthermore, our results indicated that the absolute mean values of the SOS and EOS trends under the isolated impacts of non-climate change were larger than those of the isolated impacts of climate change, which indicates that the effect of non-climate change (i.e., anthropogenic activities) on LSP variation was larger than that of climate change.With regard to the relative contribution of climatic factors to the variation in the SOS and EOS, solar radiation had the largest effect on the SOS and EOS, followed by precipitation and temperature.

Figure 1 .
Figure 1.Trends of the start of season (SOS) (a) and end of season (EOS) (c) under the impacts of climate change and the impacts of non-climate change.(b,d) show the corresponding frequency distributions of SOS and EOS, respectively.

Figure 2 .
Figure 2. Isolated impacts of climatic factors (a) and non-climate change (c) on the variation in the SOS.(b,d) are the corresponding frequency distributions for climatic factors (a) and non-climate change (c), respectively.

Figure 3 .
Figure 3. Isolated impacts of climatic factors (a) and non-climate change (c) on the variation in EOS.(b,d) are the corresponding frequency distributions for climatic factors (a) and non-climate change (c), respectively.

Figure 4 .
Figure 4. Relative contributions of climatic factors to the variation in the SOS (a) and EOS (b), respectively.Tem, temperature; Pre, precipitation; Rad, solar radiation.

:
Root mean square error (RMSE), Akaike's Information Criterion (AIC), and Bayesian Information Criterion (BIC) values for the four methods of analyzing the Moderate Resolution Imaging Spectroradiometer (MODIS) enhanced vegetation index (EVI) time-series data from 2001 to 2014 in Northern China; Table S2: Percentage of partial correlation coefficients between LSP and climatic factors; Figure S1: Distribution of meteorological stations across Northern China; red points represent stations for temperature and precipitation, and blue triangles represent stations for solar radiation; Figure S2: The time lag related to climatic factors for which the highest determination coefficients (R2) were obtained between phenophases (SOS (a) and EOS (b) from top to bottom, respectively) and climatic factors in Northern China.Tem, Pre and Sad represent mean temperature, total precipitation and total solar radiation, respectively; Figure S3: Frequency distribution and the percentages of correlations (pixels at p < 0.05, in parentheses) between phenological parameters (SOS and EOS from left to right, respectively) and climate factors (temperature, precipitation and solar radiation from up to down, respectively) were shown; Figure S4: Climatic factors variation in Northern China from 2001 to 2015.Tem, mean annual temperature; Pre, annual cumulative precipitation; Rad, annual cumulative solar radiation.Author Contributions: H.Z. (Hengjun Zhang) and W.W. outlined the research topic, assisted with the manuscript writing, and coordinated the revision activities.Z.L., Q.S., and T.W. performed the data analysis, interpretation of results, manuscript writing, and coordinated the revision activities.H.Z. (Huanmu Zeng), and T.H. performed the data collection and coordinated the revision activities.

Funding:
This study is supported by the Guangdong science and technology plan (Grant No. 2017A020216003), and National key research and development plan (Grant No. 2016YFC0400703).